python 求微分（如何用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, printmessg 为True时打印信息


 


 


 


 


 


 

。 tfirst 当为False时 
 
 
 
 
 
 
 
 
 
 
 
 
，func的格式为func(y,t...)


















，否则格式为func(t, y...)

示例

对于常微分方程

θ

′

′

(

t

)

+

b

θ

′

(

t

)

+

c

sin

⁡

θ

(

t

)

=

b

=

0.25

;

c

=

5

θ

(

)

=

π

−

0.1

;

θ

′

(

)

=

\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

)

\begin{aligned} \theta(t)&=\omega(t)\\ \omega(t)&=-b\omega(t)-c\sin\theta(t) \end{aligned}

θ′(t)ω′(t)​=ω(t)=−bω(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()

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