如何用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 ) = 0 b = 0.25 ; c = 5 θ ( 0 ) = π − 0.1 ; θ ′ ( 0 ) = 0 heta''(t)+b heta'(t)+csin heta(t)=0\ b=0.25;quad c=5\ heta(0)=pi-0.1;quad heta'(0)=0 θ′′(t)+bθ′(t)+csinθ(t)=0b=0.25;c=5θ(0)=π−0.1;θ′(0)=0
将其中的二阶导数项用一个新变量替代, ω ( t ) = θ ′ ( t ) omega(t)= heta'(t) ω(t)=θ′(t),则常微分方程可拆分成微分方程组
θ ′ ( t ) = ω ( t ) ω ′ ( t ) = − b ω ( t ) − c sin θ ( t ) egin{aligned} heta'(t)&=omega(t)\ omega'(t)&=-bomega(t)-csin heta(t) end{aligned} θ′(t)ω′(t)=ω(t)=−bω(t)−csinθ(t)
令
y
=
[
θ
,
ω
]
y=[ heta, omega]
y=[θ,ω],则
y
′
=
[
θ
′
,
ω
′
]
y'=[ heta', 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使用tkinter组件Label显示简单数学公式
- 内网渗透之DCOM横向移动
- 以目标为导向的语义交流的共同语言——一个课程学习框架
- python爬虫前奏【成信笔记】
- HTML 5 File API:文件拖放上传功能
- 教你快速创建 Python 虚拟环境
- pyenv 实现Python多版本自由切换
- 用 Python 对 Excel文件进行批量操作
- Python - 接入钉钉机器人
- Python - 抓取 iphone13 pro 线下店供货信息并发送到钉钉机器人,最后设置为定时任务
- crontab - 解决 mac 下通过 crontab 设置了 Python 脚本的定时任务却无法运行
- [源码解析] PyTorch分布式(5) ------ DistributedDataParallel 总述&如何使用
- Python科普系列——类与方法(上篇)
- SAP对STO的交货单执行PGI,报错 -Fld selectn for mvmt type 643 acct 400020 differs
- Spring Boot 实现通用 Auth 认证的 4 种方式
- 盘点4种使用Python批量合并同一文件夹内所有子文件夹下的Excel文件内所有Sheet数据
- OushuDB 学习经验分享(三):技术特点
- Java和Python思维方式的不同之处
- Python中日志记录新技能
- 奥比中光Gemini OpenCV—Python使用