Python⼩⽩的数学建模课-09微分⽅程模型
⼩⽩往往听到微分⽅程就觉得害怕，其实数学建模中的微分⽅程模型不仅没那么复杂，⽽且很容易写出⾼⽔平的数模论⽂。
本⽂介绍微分⽅程模型的建模与求解，通过常微分⽅程、常微分⽅程组、⾼阶常微分⽅程 3个案例⼿把⼿教你搞定微分⽅程。
通过⼆阶 RLC 电路问题，学习微分⽅程模型的建模、求解和讨论。
欢迎关注系列，每周持续更新
1. 微分⽅程
1.1 基本概念
微分⽅程是描述系统的状态随时间和空间演化的数学⼯具。物理中许多涉及变⼒的运动学、动⼒学问题，如空⽓的阻⼒为速度函数的落体运动等问题，很多可以⽤微分⽅程求
解。微分⽅程在化学、⼯程学、经济学和⼈⼝统计等领域也有⼴泛应⽤。
具体来说，微分⽅程是指含有未知函数及其导数的关系式。
微分⽅程按⾃变量个数分为：只有⼀个⾃变量的常微分⽅程（Ordinary Differential Equations）和包含两个或两个以上独⽴变量的偏微分⽅程（Partial Differential
Equations）。
微分⽅程按阶数分为：⼀阶、⼆阶、⾼阶，微分⽅程的阶数取决于⽅程中最⾼次导数的阶数。
微分⽅程还可以分为：（⾮）齐次，常（变）系数，（⾮）线性，初值问题/边界问题...
以上内容看看就算了，看多了就吓跑了。
1.2 微分⽅程的数学建模
微分⽅程的数学建模其实并不复杂，基本过程就是分析题⽬属于哪⼀类问题、可以选择什么微分⽅程模型，然后如何使⽤现有的微分⽅程模型建模。
在数学、⼒学、物理、化学等各个学科领域的课程中，针对该学科的各种问题都会建⽴适当的数学模型。在中学课程中，各学科的数学模型主要是线性或⾮线性⽅程，⽽在⼤学
物理和各专业的课程中，越来越多地出现⽤微分⽅程描述的数学模型。
数学建模中的微分⽅程问题，通常还是这些专业课程中相对简单的模型，专业课程的教材在介绍⼀个模型时，往往都做了⾮常详细的讲解。只要搞清楚问题的类型、选择好数学
模型，建模和求解并不是很难，⽽且在撰写论⽂时对问题背景、使⽤范围、假设条件、求解过程有⼤量现成的内容可以复制参考。
⼩⽩之所以害怕，⼀是看到微分⽅程就⼼⾥发怵，⼆是缺乏专业背景，不知道从哪⾥查资料、不能判断问题的类型、不知道选择什么模型、不善于从题⽬内容得出模型参数，也
不知道如何编程求解。所以，⽼师说，⼀看这就是××问题，显然就可以⽤××模型。⼩⽩说，我们还是换 B题吧。
本系列将会从简单的微分⽅程模型⼊⼿，重点介绍微分⽅程数值解法的编程实现，并通过分析问题、建⽴模型的案例帮助⼩⽩树⽴信⼼和动⼒。
希望你在学习本系列之后，会发现微分⽅程模型是数学建模中最容易的题型：模型教材，建模例题，求解有例程，讨论有套路，论⽂够档次。
1.3 微分⽅程的数值解法
在学习专业课程时，经常会推导和求解微分⽅程的解析解，⼩⽩对微分⽅程模型的恐惧就是从⾼等数学“微分⽅程”开始，经过专业课的不断强化⽽形成的。实际上，只有很少的
微分⽅程可以解析求解，⼤多数的微分⽅程只能采⽤数值⽅法进⾏求解。
微分⽅程的数值求解是先把时间和空间离散化，然后将微分化为差分，建⽴递推关系，然后反复进⾏迭代计算，得到任意时间和空间的值。
如果你还是觉得头晕⽬眩，我们可以说的更简单⼀些。建模就是把专业课教材上的公式抄下来，求解就是把公式的参数输⼊到 Python 函数中。
我们先说求解。求解常微分⽅程的基本⽅法，有欧拉法、龙格库塔法等，可以详见各种教材，撰写数模竞赛论⽂时还是可以抄⼏段的。本⽂沿⽤“编程⽅案”的概念，不涉及这些
算法的具体内容，只探讨如何使⽤ Python 的⼯具包、库函数，零基础求解微分⽅程模型。
我们的选择是 Python 常⽤⼯具包三剑客：Scipy、Numpy 和 Matplotlib：
Scipy 是 Python 算法库和数学⼯具包，包括最优化、线性代数、积分、插值、特殊函数、傅⾥叶变换、信号和图像处理、常微分⽅程求解等模块。有⼈介绍 Scipy 就是
Python 语⾔的 Matlab，所以⼤部分数学建模问题都可以⽤它搞定。
Numpy 提供了⾼维数组的实现与计算的功能，如线性代数运算、傅⾥叶变换及随机数⽣成，另外还提供了与 C/C++ 等语⾔的集成⼯具。
Matplotlib 是可视化⼯具包，可以⽅便地绘制各种数据可视化图表，如折线图、散点图、直⽅图、条形图、箱形图、饼图、三维图，等等。
顺便说⼀句，还有⼀个 Python 符号运算⼯具包 SymPy，以解析⽅式求解积分、微分⽅程，也就是说给出的结果是微分⽅程的解析解表达式。很⽜，但只能求解有解析解的微分
⽅程，所以，你知道就可以了。
2. SciPy 求解常微分⽅程（组）
2.1 ⼀阶常微分⽅程（组）模型
给定初始条件的⼀阶常微分⽅程（组）的标准形式是：
\[\begin{cases} \begin{aligned} &\frac{dy}{dt} = f(y,t)\\ &y(t_0) = y_0 \end{aligned} \end{cases} \]
式中的 y 在常微分⽅程中是标量，在常微分⽅程组中是数组向量。
2.2 scipy.integrate.odeint() 函数
SciPy 提供了两种⽅式求解常微分⽅程：基于odeint函数的 API ⽐较简单易学，基于ode类的⾯向对象的 API 更加灵活。
**scipy.integrate.odeint() **是求解微分⽅程的具体⽅法，通过数值积分来求解常微分⽅程组。在odeint函数内部使⽤ FORTRAN 库 odepack 中的 lsoda，可以求解⼀阶刚性系统
和⾮刚性系统的初值问题。官⽹介绍详见：。
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=Fa odeint 的主要参数：
求解标准形式的微分⽅程（组）主要使⽤前三个参数：
func: callable(y, t, …) 导数函数 \(f(y,t)\) ，即 y 在 t 处的导数，以函数的形式表⽰
y0: array： 初始条件 \(y_0\)，对于常微分⽅程组 \(y_0\) 则为数组向量
小白学python买什么书t: array： 求解函数值对应的时间点的序列。序列的第⼀个元素是与初始条件 \(y_0\) 对应的初始时间 \(t_0\)；时间序列必须是单调递增或单调递减的，允许重复值。
其它参数简介如下：
args: 向导数函数 func 传递参数。当导数函数 \(f(y,t,p1,p2,..)\) 包括可变参数 p1,p2.. 时，通过 args =(p1,p2,..) 可以将参数p1,p2.. 传递给导数函数 func。argus 的⽤法参见
2.4 中的实例2。
Dfun: func 的雅可⽐矩阵，⾏优先。如果 Dfun 未给出，则算法⾃动推导。
col_deriv: ⾃动推导 Dfun的⽅式。
printmessg: 布尔值。控制是否打印收敛信息。
其它参数⽤于控制求解算法的参数，⼀般情况可以忽略。
odeint 的主要返回值：
y: array 数组，形状为 (len(t),len(y0)，给出时间序列 t 中每个时刻的 y 值。
3. 实例1：Scipy 求解⼀阶常微分⽅程
3.1 例题 1：求微分⽅程的数值解
\[\begin{cases} \begin{aligned} &\frac{dy}{dt} = sin(t^2)\\ &y(-10) = 1 \end{aligned} \end{cases} \]
3.2 常微分⽅程的编程步骤
以该题为例讲解 scipy.integrate.odeint() 求解常微分⽅程初值问题的步骤：
1. 导⼊ scipy、numpy、matplotlib 包；
2. 定义导数函数 \(f(y,t)=sin(t^2)\) ；
3. 定义初值 \(y_0\) 和 \(y\) 的定义区间 \([t_0,\ t]\)；
4. 调⽤ odeint() 求 \(y\) 在定义区间 \([t_0,\ t]\) 的数值解。
3.3 Python 例程
eof函数vfp# 1. 求解微分⽅程初值问题(scipy.integrate.odeint)
from scipy.integrate import odeint # 导⼊ scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
def dy_dt(y, t): # 定义函数 f(y,t)
return np.sin(t**2)
y0 = [1] # y0 = 1 也可以
t = np.arange(-10,10,0.01) # (start,stop,step)
y = odeint(dy_dt, y0, t) # 求解微分⽅程初值问题
# 绘图
plt.plot(t, y)
plt.title("scipy.integrate.odeint")
plt.show()
3.4 Python 例程运⾏结果
4. 实例2：Scipy 求解⼀阶常微分⽅程组
4.1 例题 2：求洛伦兹（Lorenz）⽅程的数值解
洛伦兹（Lorenz）混沌吸引⼦的轨迹可以由如下的 3个微分⽅程描述：
\[\begin{cases} \begin{aligned} &\frac{dx}{dt} = \sigma (y-x)\\ &\frac{dy}{dt} = x (\rho-z) - y\\ &\frac{dz}{dt} = xy - \beta z\\ \end{aligned} \end{cases} \]
洛伦兹⽅程将⼤⽓流体运动的强度 x 与⽔平和垂直⽅向的温度变化 y 和 z 联系起来，进⾏⼤⽓对流系统的模拟，现已⼴泛应⽤于天⽓预报、空⽓污染和全球⽓候变化的研究。参数 \(\sigma\) 称为普兰特数，\(\rho\) 是规范化的瑞利数，\(\beta\) 和⼏何形状相关。洛伦兹⽅程是⾮线性微分⽅程组，⽆法求出解析解，只能使⽤数值⽅法求解。
4.2 洛伦兹（Lorenz）⽅程问题的编程步骤
以该题为例讲解 scipy.integrate.odeint() 求解常微分⽅程初值问题的步骤：
新生活cms订货系统下载1. 导⼊ scipy、numpy、matplotlib 包；
2. 定义导数函数 lorenz(W, t, p, r, b)
注意 odeint() 函数中定义导数函数的标准形式是 \(f(y,t)\) ，对于微分⽅程组 y 表⽰向量。
为避免混淆，我们记为 \(W=[x,y,z]\)，函数 lorenz(W,t) 定义导数函数 \(f(W,t)\) 。
⽤ p,r,b 分别表⽰⽅程中的参数 \(\sigma、\rho、\beta\)，则对导数定义函数编程如下：
# 导数函数，求 W=[x,y,z] 点的导数 dW/dt
def lorenz(W,t,p,r,b):
x, y, z = W # W=[x,y,z]
dx_dt = p*(y-x) # dx/dt = p*(y-x), p: sigma
dy_dt = x*(r-z) - y # dy/dt = x*(r-z)-y, r:rho
dz_dt = x*y - b*z # dz/dt = x*y - b*z, b;beta
return np.array([dx_dt,dy_dt,dz_dt])
3. 定义初值 \(W_0\) 和 \(W\) 的定义区间 \([t_0,\ t]\)；
4. 调⽤ odeint() 求 \(W\) 在定义区间 \([t_0,\ t]\) 的数值解。
注意例程中通过 args=paras 或 args = (10.0,28.0,3.0) 将参数 (p,r,b) 传递给导数函数 lorenz(W,t,p,r,b)。参数 (p,r,b) 当然也可以不作为函数参数传递，⽽是在导数函数lorenz() 中直接设置。但例程的参数传递⽅法，使导数函数结构清晰、更为通⽤。另外，对于可变参数问题，使⽤这种参数传递⽅式就⾮常⽅便。
4.3 洛伦兹（Lorenz）⽅程问题 Python 例程
# 2. 求解微分⽅程组初值问题(scipy.integrate.odeint)
from scipy.integrate import odeint # 导⼊ scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 导数函数, 求 W=[x,y,z] 点的导数 dW/dt
def lorenz(W,t,p,r,b): # by youcans
x, y, z = W # W=[x,y,z]
dx_dt = p*(y-x) # dx/dt = p*(y-x), p: sigma
dy_dt = x*(r-z) - y # dy/dt = x*(r-z)-y, r:rho
dz_dt = x*y - b*z # dz/dt = x*y - b*z, b;beta
return np.array([dx_dt,dy_dt,dz_dt])
t = np.arange(0, 30, 0.01) # 创建时间点 (start,stop,step)
paras = (10.0, 28.0, 3.0) # 设置 Lorenz ⽅程中的参数 (p,r,b)
# 调⽤ode对lorenz进⾏求解, ⽤两个不同的初始值 W1、W2 分别求解
W1 = (0.0, 1.00, 0.0) # 定义初值为 W1
track1 = odeint(lorenz, W1, t, args=(10.0, 28.0, 3.0)) # args 设置导数函数的参数
W2 = (0.0, 1.01, 0.0) # 定义初值为 W2
track2 = odeint(lorenz, W2, t, args=paras) # 通过 paras 传递导数函数的参数
# 绘图
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2], color='magenta') # 绘制轨迹 1
ax.plot(track2[:,0], track2[:,1], track2[:,2], color='deepskyblue') # 绘制轨迹 2
ax.set_title("Lorenz attractor by scipy.integrate.odeint")
plt.show()
4.4 洛伦兹（Lorenz）⽅程问题 Python 例程运⾏结果
js能实现的简单效果5. 实例3：Scipy 求解⾼阶常微分⽅程
⾼阶常微分⽅程，必须做变量替换，化为⼀阶微分⽅程组，再⽤ odeint 求数值解。
5.1 例题 3：求⼆阶 RLC 振荡电路的数值解
零输⼊响应的 RLC 振荡电路可以由如下的⼆阶微分⽅程描述：
\[\begin{cases} \begin{aligned} &\frac{d^2 u}{dt^2} + \frac{R}{L} * \frac{du}{dt} + \frac{1}{LC}*u = 0\\ &u(0) = U_0\\ &u'(0)= 0 \end{aligned} \end{cases} \]
令 $ \alpha = R/2L\(、\)\omega_0^2=1/LC$，在零输⼊响应 \(u_s=0\) 时上式可以写成：
\[\begin{cases} \begin{aligned} &\frac{d^2 u}{dt^2} + 2 \alpha \frac{du}{dt} + \omega_0^2 u = 0\\ &u(0) = U_0\\ &u'(0)= 0 \end{aligned} \end{cases} \]
对⼆阶微分⽅程问题，引⼊变量 \(v = {du}/{dt}\)，通过变量替换就把原⽅程化为如下的微分⽅程组：
\[\begin{cases} \begin{aligned} &\frac{du}{dt} = v \\ &\frac{dv}{dt} = -2\alpha v - \omega_0^2 u\\ &u(0)=U_0\\ &v(0)=0 \end{aligned} \end{cases} \]
这样就可以⽤上节求解微分⽅程组的⽅法来求解⾼阶微分⽅程问题。
5.2 ⼆阶微分⽅程问题的编程步骤
以RLC 振荡电路为例讲解 scipy.integrate.odeint() 求解⾼阶常微分⽅程初值问题的步骤：
1. 导⼊ scipy、numpy、matplotlib 包；
2. 定义导数函数 deriv(Y, t, a, w)
注意 odeint() 函数中定义导数函数的标准形式是 \(f(y,t)\) ，本问题中 y 表⽰向量，记为 \(Y=[u,v]\)
导数定义函数 deriv(Y, t, a, w) 编程如下，其中 a, w 分别表⽰⽅程中的参数 \(\alpha、\omega\)：
# 导数函数，求 Y=[u,v] 点的导数 dY/dt
def deriv(Y, t, a, w):
u, v = Y # Y=[u,v]
dY_dt = [v, -2*a*v-w*w*u]
return dY_dt
round函数怎么按出来3. 定义初值 \(Y_0=[u_0,v_0]\) 和 \(Y\) 的定义区间 \([t_0,\ t]\)；
4. 调⽤ odeint() 求 \(Y=[u,v]\) 在定义区间 \([t_0,\ t]\) 的数值解。
例程中通过 args=paras 将参数 (a,w) 传递给导数函数 deriv(Y, t, a, w) 。本例要考察不同参数对结果的影响，这种参数传递⽅法使⽤⾮常⽅便。
5.3 ⼆阶微分⽅程问题 Python 例程
# 3. 求解⼆阶微分⽅程初值问题(scipy.integrate.odeint)
# Second ODE by scipy.integrate.odeint
from scipy.integrate import odeint # 导⼊ scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
# 导数函数，求 Y=[u,v] 点的导数 dY/dt
def deriv(Y, t, a, w):
u, v = Y # Y=[u,v]
dY_dt = [v, -2*a*v-w*w*u]
return dY_dt
t = np.arange(0, 20, 0.01) # 创建时间点 (start,stop,step)
# 设置导数函数中的参数 (a, w)
paras1 = (1, 0.6) # 过阻尼：a^2 - w^2 > 0
paras2 = (1, 1) # 临界阻尼：a^2 - w^2 = 0
paras3 = (0.3, 1) # ⽋阻尼：a^2 - w^2 < 0
# 调⽤ode对进⾏求解, ⽤两个不同的初始值 W1、W2 分别求解
Y0 = (1.0, 0.0) # 定义初值为 Y0=[u0,v0]
Y1 = odeint(deriv, Y0, t, args=paras1) # args 设置导数函数的参数
Y2 = odeint(deriv, Y0, t, args=paras2) # args 设置导数函数的参数
Y3 = odeint(deriv, Y0, t, args=paras3) # args 设置导数函数的参数
# W2 = (0.0, 1.01, 0.0) # 定义初值为 W2
# track2 = odeint(lorenz, W2, t, args=paras) # 通过 paras 传递导数函数的参数
# 绘图
plt.plot(t, Y1[:, 0], 'r-', label='u1(t)')
plt.plot(t, Y2[:, 0], 'b-', label='u2(t)')
plt.plot(t, Y3[:, 0], 'g-', label='u3(t)')
plt.plot(t, Y1[:, 1], 'r:', label='v1(t)')
plt.plot(t, Y2[:, 1], 'b:', label='v2(t)')
plt.plot(t, Y3[:, 1], 'g:', label='v3(t)')
plt.axis([0, 20, -0.8, 1.2])
plt.legend(loc='best')
plt.title("Second ODE by scipy.integrate.odeint")
plt.show()
5.4 ⼆阶⽅程问题 Python 例程运⾏结果
结果讨论：
RLC串联电路是典型的⼆阶系统，在零输⼊条件下根据 \(\alpha\) 与 \(\omega\) 的关系，电路的输出响应存在四种情况：
1. 过阻尼： \(\alpha^2 - \omega^2>0\) ，有 2 个不相等的负实数根；
2. 临界阻尼： \(\alpha^2 - \omega^2 = 0\)，有 2 个相等的负实数根；
3. ⽋阻尼： \(\alpha^2 - \omega^2 <0\)，有⼀对共轭复数根；
4. ⽆阻尼：\(R=0\)，有⼀对纯虚根。
例程中所选择的 3 组参数分别对应过阻尼、临界阻尼和⽋阻尼的条件，微分⽅程的数值结果很好地体现了不同情况的相应曲线。
6. ⼩结
1. ⼩⽩⾸先要有信⼼，⽤ Scipy ⼯具包求解标准形式的微分⽅程（组），编程实现是⾮常简单、容易掌握的。
2. 其次要认识到，由于微分⽅程的解的特性是与模型参数密切相关的，不同参数的解可能具有完全不同的形态。这就涉及模型稳定性、灵敏度的分析，很容易使论⽂写的⾮常
丰富和精彩。
3. 不会从问题建⽴微分⽅程模型怎么办，不会展开参数对稳定性、灵敏度的影响进⾏讨论怎么办？谁让你⾃⼰做呢，当然是先去相关专业的教材、论⽂，从中选择⽐较接
近、⽐较简单的理论和模型，然后通过各种假设强⾏将题⽬简化为模型中的条件，这就可以照猫画虎了。
【本节完】
版权声明：
欢迎关注原创作品
原创作品，转载必须标注原⽂链接：
编程软件名称Copyright 2021 Youcans, XUPT
Crated：2021-06-08
欢迎关注，每周更新数模笔记
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论