使用 solve_ivp 代替 odeint 求解初始问题值

Question

目前，我使用 odeint

求解以下 ODE 方程组

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

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

初始条件：x = 0, y = 0

但是，我想使用 solve_ivp 这似乎是解决此类问题的推荐选项，但老实说我不知道如何调整代码...

这是我在 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()

Answer 1

重要的是 solve_ivp 期望 f(t, z) 作为 ODE 的右侧。如果你不想改变你的ode函数并且还想传递你的参数u，我建议定义一个包装函数：

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)

现在可以轻松调用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()

请注意，如果不通过 t_eval=t，求解器将自动选择 tspan 内的时间点来存储解决方案。

使用 solve_ivp 代替 odeint 求解初始问题值

Using solve_ivp instead of odeint to solve initial problem value

python

scipy

ode

differential-equations