Btukfyl

Currently, I solve the following ODE system of equations using odeint

dx/dt = (-x + u)/2.0

dy/dt = (-y + x)/5.0

initial conditions: x = 0, y = 0

However, I would like to use solve_ivp which seems to be the recommended option for this type of problems, but honestly I don't know how to adapt the code...

Here is the code I'm using with odeint:

import numpy as np

from scipy.integrate import odeint, solve_ivp

import matplotlib.pyplot as plt



def model(z, t, u):

    x = z[0]

    y = z[1]

    dxdt = (-x + u)/2.0

    dydt = (-y + x)/5.0

    dzdt = [dxdt, dydt]

    return dzdt



def main():

    # initial condition

    z0 = [0, 0]



    # number of time points

    n = 401



    # time points

    t = np.linspace(0, 40, n)



    # step input

    u = np.zeros(n)

    # change to 2.0 at time = 5.0

    u[51:] = 2.0



    # store solution

    x = np.empty_like(t)

    y = np.empty_like(t)

    # record initial conditions

    x[0] = z0[0]

    y[0] = z0[1]



    # solve ODE

    for i in range(1, n):

        # span for next time step

        tspan = [t[i-1], t[i]]

        # solve for next step

        z = odeint(model, z0, tspan, args=(u[i],))

        # store solution for plotting

        x[i] = z[1][0]

        y[i] = z[1][1]

        # next initial condition

        z0 = z[1]



    # plot results

    plt.plot(t,u,'g:',label='u(t)')

    plt.plot(t,x,'b-',label='x(t)')

    plt.plot(t,y,'r--',label='y(t)')

    plt.ylabel('values')

    plt.xlabel('time')

    plt.legend(loc='best')

    plt.show()



main()

It's important that solve_ivp expects f(t, z) as right-hand side of the ODE. If you don't want to change your ode function and also want to pass your parameter u, I recommend to define a wrapper function:

def model(z, t, u):

    x = z[0]

    y = z[1]

    dxdt = (-x + u)/2.0

    dydt = (-y + x)/5.0

    dzdt = [dxdt, dydt]

    return dzdt



def odefun(t, z):

    if t < 5:

        return model(z, t, 0)

    else:

        return model(z, t, 2)

Now it's easy to call solve_ivp:

def main():

    # initial condition

    z0 = [0, 0]



    # number of time points

    n = 401



    # time points

    t = np.linspace(0, 40, n)



    # step input

    u = np.zeros(n)

    # change to 2.0 at time = 5.0

    u[51:] = 2.0



    res = solve_ivp(fun=odefun, t_span=[0, 40], y0=z0, t_eval=t)

    x = res.y[0, :]

    y = res.y[1, :]



    # plot results

    plt.plot(t,u,'g:',label='u(t)')

    plt.plot(t,x,'b-',label='x(t)')

    plt.plot(t,y,'r--',label='y(t)')

    plt.ylabel('values')

    plt.xlabel('time')

    plt.legend(loc='best')

    plt.show()



main()

Note that without passing t_eval=t, the solver will automatically choose the time points inside tspan at which the solution will be stored.