为了账号安全,请及时绑定邮箱和手机立即绑定

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

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

斯蒂芬大帝 2021-08-24 16:55:00
目前,我使用 odeint 求解以下 ODE 方程组dx/dt = (-x + u)/2.0dy/dt = (-y + x)/5.0初始条件:x = 0,y = 0但是,我想使用 solve_ivp 这似乎是此类问题的推荐选项,但老实说我不知道如何调整代码......这是我与 odeint 一起使用的代码:import numpy as npfrom scipy.integrate import odeint, solve_ivpimport matplotlib.pyplot as pltdef model(z, t, u):    x = z[0]    y = z[1]    dxdt = (-x + u)/2.0    dydt = (-y + x)/5.0    dzdt = [dxdt, dydt]    return dzdtdef 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()
查看完整描述

1 回答

  • 1 回答
  • 0 关注
  • 409 浏览
慕课专栏
更多

添加回答

举报

0/150
提交
取消
意见反馈 帮助中心 APP下载
官方微信