我一直在使用 sympy 和 scipy,但无法找到或弄清楚如何求解耦合微分方程组(非线性,一阶)。
那么有什么方法可以求解耦合微分方程吗?
方程的形式为:
V11'(s) = -12*v12(s)**2
v22'(s) = 12*v12(s)**2
v12'(s) = 6*v11(s)*v12(s) - 6*v12(s)*v22(s) - 36*v12(s)
初始条件为 v11(s)、v22(s)、v12(s)。
有关使用 scipy 的 ODE 的数值解,请scipy.integrate.solve_ivp
参阅 scipy.integrate.odeint
或scipy.integrate.ode。
SciPy Cookbook中给出了一些示例(向下滚动到“常微分方程”部分)。
除了已经提到的 SciPy 方法之外,它现在还具有odeint
更新且通常更方便的方法。一个完整的例子,编码为一个数组:ode
solve_ivp
[v11, v22, v12]
v
from scipy.integrate import solve_ivp
def rhs(s, v):
return [-12*v[2]**2, 12*v[2]**2, 6*v[0]*v[2] - 6*v[2]*v[1] - 36*v[2]]
res = solve_ivp(rhs, (0, 0.1), [2, 3, 4])
(0, 0.1)
这解决了具有初始值的区间上的系统[2, 3, 4]
。结果具有自变量(您的符号中的 s)为res.t
:
array([ 0. , 0.01410735, 0.03114023, 0.04650042, 0.06204205,
0.07758368, 0.0931253 , 0.1 ])
这些值是自动选择的。可以提供t_eval
在所需点评估解决方案:例如,t_eval=np.linspace(0, 0.1)
。
因变量(我们正在寻找的函数)位于res.y
:
array([[ 2. , 0.54560138, 0.2400736 , 0.20555144, 0.2006393 ,
0.19995753, 0.1998629 , 0.1998538 ],
[ 3. , 4.45439862, 4.7599264 , 4.79444856, 4.7993607 ,
4.80004247, 4.8001371 , 4.8001462 ],
[ 4. , 1.89500744, 0.65818761, 0.24868116, 0.09268216,
0.0345318 , 0.01286543, 0.00830872]])
使用 Matplotlib,此解决方案被绘制为(如果我按上述plt.plot(res.t, res.y.T)
提供,绘图会更平滑)。t_eval
最后,如果系统涉及高于 1 的阶方程,则需要使用归约到一阶系统。