【原】使用Python求解微分方程（1）一阶ODE

很多仿真问题的实质都是求解微分方程，包括常微分方程（ODE）及偏微分方程（PDE）。使用Python的科学计算类型库（numpy、scipy等）可以方便地对微分方程进行求解。

对于简单的一阶ODE，这里以下落时的空气阻力问题为例进行介绍。

该问题可以由下式进行描述：

要进行求解，我们首先将其变换为下式这种形式

变换后得到：

接着我们将上式在python中表达出来（系数任给）

import numpy as npimport matplotlib.pyplot as pltimport scipy as spfrom scipy.integrate import odeintfrom scipy.integrate import solve_ivpdef dvdt(t, v):    return 2*v**2 - 3v0 = 0

求解微分方程时，scipy中有着两个常见的求解器：

• odeint：使用了 FORTRAN 的odepack库 中称为 lsoda 的解法。
  

• solve_ivp：定制化程度更高。
  

我们定义时间为1s，步长为0.01s，并使用上述两种方法分别求解。

t = np.linspace(0, 1, 100)sol_m1 = odeint(dvdt, y0=v0, t=t, tfirst=True)sol_m2 = solve_ivp(dvdt, t_span=(0,max(t)), y0=[v0], t_eval=t)

将结果出图

v_sol_m1 = sol_m1.T[0]v_sol_m2 = sol_m2.y[0]l1,=plt.plot(t, v_sol_m1,'ob')#plot返回的list对象（list of Line2D）需要解构，因此需要在line1和等号之间加一个逗号l2,=plt.plot(t, v_sol_m2, '--r')plt.ylabel('v(t)', fontsize=22)plt.xlabel('t', fontsize=22)plt.legend((l1,l2),['odeint','solve_ivp'])plt.show()

运行结果如下