如何用Python求解微分方程组

Posted 微小冷

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了如何用Python求解微分方程组相关的知识,希望对你有一定的参考价值。

文章目录

odeint简介

scipy文档中将odeint函数和ode, comples_ode这两个类称为旧API,是scipy早期使用的微分方程求解器,但由于是Fortran实现的,尽管使用起来并不方便,但速度没得说,所以有的时候还挺推荐使用的。

其中,odeint的参数如下

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)

其中func为待求解函数;y0为初值;t为自变量列表,其他参数都有默认选项,可以不填,而且这些参数非常多,其中常用的有

  • args func中除了t之外的其他变量
  • Dfun func的梯度函数,当此参数不为None时,若将col_deriv设为True,则可提升效率。
  • full_output 如果为True,则额外返回一个参数字典
  • ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5,
  • printmessgTrue时打印信息。
  • tfirst 当为False时,func的格式为func(y,t...),否则格式为func(t, y...)

示例

对于常微分方程

θ ′ ′ ( t ) + b θ ′ ( t ) + c sin ⁡ θ ( t ) = 0 b = 0.25 ; c = 5 θ ( 0 ) = π − 0.1 ; θ ′ ( 0 ) = 0 \\theta''(t)+b\\theta'(t)+c\\sin\\theta(t)=0\\\\ b=0.25;\\quad c=5\\\\ \\theta(0)=\\pi-0.1;\\quad \\theta'(0)=0 θ′′(t)+bθ(t)+csinθ(t)=0b=0.25;c=5θ(0)=π0.1;θ(0)=0

将其中的二阶导数项用一个新变量替代, ω ( t ) = θ ′ ( t ) \\omega(t)=\\theta'(t) ω(t)=θ(t),则常微分方程可拆分成微分方程组

θ ′ ( t ) = ω ( t ) ω ′ ( t ) = − b ω ( t ) − c sin ⁡ θ ( t ) \\beginaligned \\theta'(t)&=\\omega(t)\\\\ \\omega'(t)&=-b\\omega(t)-c\\sin\\theta(t) \\endaligned θ(t)ω(t)=ω(t)=(t)csinθ(t)

y = [ θ , ω ] y=[\\theta, \\omega] y=[θ,ω],则 y ′ = [ θ ′ , ω ′ ] y'=[\\theta', \\omega'] y=[θ,ω],据此可设计函数func

import numpy as np
def pend(y, t, b, c):
    th, om = y
    dydt = [om, -b*om - c*np.sin(th)]
    return dydt

然后调用并求解

from scipy.integrate import odeint
y0 = [np.pi-0.1, 0]
t = np.linspace(0, 10, 101)
sol = odeint(pend, y0, t, args=(0.25, 5))

然后绘制一下结果

import matplotlib.pyplot as plt
plt.plot(t, sol[:,0], label="theta")
plt.plot(t, sol[:,1], label="omega")
plt.legend()
plt.show()

这个形状还是比较离奇的。

以上是关于如何用Python求解微分方程组的主要内容,如果未能解决你的问题,请参考以下文章

如何用matlab求解微分方程并画图

如何用matlab求解微分方程组的数值解

如何用python进行符号求解?

如何用matlab求解常微分方程?matlab解常微分方程之符号解法介绍

如何用matlab求解二阶微分方程,以及程序实例

如何用ode45求解matlab中的耦合微分方程