做自媒体网站需要注册什么公司,用糖做的网站,wordpress能注册,网站开发学习课程小白往往听到微分方程就觉得害怕#xff0c;其实数学建模中的微分方程模型不仅没那么复杂#xff0c;而且很容易写出高水平的数模论文。
本文介绍微分方程模型的建模与求解#xff0c;通过常微分方程、常微分方程组、高阶常微分方程 3个案例手把手教你搞定微分方程。
通过…
小白往往听到微分方程就觉得害怕其实数学建模中的微分方程模型不仅没那么复杂而且很容易写出高水平的数模论文。
本文介绍微分方程模型的建模与求解通过常微分方程、常微分方程组、高阶常微分方程 3个案例手把手教你搞定微分方程。
通过二阶 RLC 电路问题学习微分方程模型的建模、求解和讨论。
更多微分方程数学模型案例参见 【新冠疫情 模型系列】
欢迎关注『Python小白的数学建模课 Youcans』系列每周持续更新 1. 微分方程
1.1 基本概念
微分方程是描述系统的状态随时间和空间演化的数学工具。物理中许多涉及变力的运动学、动力学问题如空气的阻力为速度函数的落体运动等问题很多可以用微分方程求解。微分方程在化学、工程学、经济学和人口统计等领域也有广泛应用。
具体来说微分方程是指含有未知函数及其导数的关系式。
微分方程按自变量个数分为只有一个自变量的常微分方程Ordinary Differential Equations和包含两个或两个以上独立变量的偏微分方程Partial Differential Equations。微分方程按阶数分为一阶、二阶、高阶微分方程的阶数取决于方程中最高次导数的阶数。微分方程还可以分为非齐次常变系数非线性初值问题/边界问题…
以上内容看看就算了看多了就吓跑了。 欢迎关注『Python小白的数学建模课 Youcans』系列每周持续更新 Python小白的数学建模课-01.新手必读 Python小白的数学建模课-02.数据导入 Python小白的数学建模课-03.线性规划 Python小白的数学建模课-04.整数规划 Python小白的数学建模课-05.0-1规划 Python小白的数学建模课-06.固定费用问题 Python小白的数学建模课-07.选址问题 Python小白的数学建模课-09.微分方程模型 Python小白的数学建模课-10.微分方程边值问题 Python小白的数学建模课-B2. 新冠疫情 SI模型 Python小白的数学建模课-B3. 新冠疫情 SIS模型 Python小白的数学建模课-B4. 新冠疫情 SIR模型 Python小白的数学建模课-B5. 新冠疫情 SEIR模型 Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型 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 一阶常微分方程组模型
给定初始条件的一阶常微分方程组的标准形式是 {dydtf(y,t)y(t0)y0\begin{cases} \begin{aligned} \frac{dy}{dt} f(y,t)\\ y(t_0) y_0 \end{aligned} \end{cases} ⎩⎨⎧dtdyf(y,t)y(t0)y0
式中的 y 在常微分方程中是标量在常微分方程组中是数组向量。 2.2 scipy.integrate.odeint() 函数
SciPy 提供了两种方式求解常微分方程基于 odeint 函数的 API 比较简单易学基于 ode 类的面向对象的 API 更加灵活。
**scipy.integrate.odeint() **是求解微分方程的具体方法通过数值积分来求解常微分方程组。在 odeint 函数内部使用 FORTRAN 库 odepack 中的 lsoda可以求解一阶刚性系统和非刚性系统的初值问题。官网介绍详见 scipy.integrate.odeint — SciPy v1.6.3 Reference Guide 。
scipy.integrate.odeint(func, y0, t, args(), DfunNone, col_deriv0, full_output0, mlNone, muNone, rtolNone, atolNone, tcritNone, h00.0, hmax0.0, hmin0.0, ixpr0, mxstep0, mxhnil0, mxordn12, mxords5, printmessg0, tfirstFalse)odeint 的主要参数
求解标准形式的微分方程组主要使用前三个参数
func: callable(y, t, …) 导数函数 f(y,t)f(y,t)f(y,t) 即 y 在 t 处的导数以函数的形式表示y0: array 初始条件 y0y_0y0对于常微分方程组 y0y_0y0 则为数组向量t: array 求解函数值对应的时间点的序列。序列的第一个元素是与初始条件 y0y_0y0 对应的初始时间 t0t_0t0时间序列必须是单调递增或单调递减的允许重复值。
其它参数简介如下 args: 向导数函数 func 传递参数。当导数函数 f(y,t,p1,p2,..)f(y,t,p1,p2,..)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. 实例1Scipy 求解一阶常微分方程
3.1 例题 1求微分方程的数值解
{dydtsin(t2)y(−10)1\begin{cases} \begin{aligned} \frac{dy}{dt} sin(t^2)\\ y(-10) 1 \end{aligned} \end{cases} ⎩⎨⎧dtdysin(t2)y(−10)1 3.2 常微分方程的编程步骤
以该题为例讲解 scipy.integrate.odeint() 求解常微分方程初值问题的步骤 导入 scipy、numpy、matplotlib 包 定义导数函数 f(y,t)sin(t2)f(y,t)sin(t^2)f(y,t)sin(t2) 定义初值 y0y_0y0 和 yyy 的定义区间 [t0,t][t_0,\ t][t0, t] 调用 odeint() 求 yyy 在定义区间 [t0,t][t_0,\ t][t0, t] 的数值解。 3.3 Python 例程
# 1. 求解微分方程初值问题(scipy.integrate.odeint)
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as pltdef 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. 实例2Scipy 求解一阶常微分方程组
4.1 例题 2求洛伦兹Lorenz方程的数值解
洛伦兹Lorenz混沌吸引子的轨迹可以由如下的 3个微分方程描述
{dxdtσ(y−x)dydtx(ρ−z)−ydzdtxy−βz\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} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧dtdxσ(y−x)dtdyx(ρ−z)−ydtdzxy−βz
洛伦兹方程将大气流体运动的强度 x 与水平和垂直方向的温度变化 y 和 z 联系起来进行大气对流系统的模拟现已广泛应用于天气预报、空气污染和全球气候变化的研究。参数 σ\sigmaσ 称为普兰特数ρ\rhoρ 是规范化的瑞利数β\betaβ 和几何形状相关。洛伦兹方程是非线性微分方程组无法求出解析解只能使用数值方法求解。 4.2 洛伦兹Lorenz方程问题的编程步骤
以该题为例讲解 scipy.integrate.odeint() 求解常微分方程初值问题的步骤 导入 scipy、numpy、matplotlib 包 定义导数函数 lorenz(W, t, p, r, b) 注意 odeint() 函数中定义导数函数的标准形式是 f(y,t)f(y,t)f(y,t) 对于微分方程组 y 表示向量。 为避免混淆我们记为 W[x,y,z]W[x,y,z]W[x,y,z]函数 lorenz(W,t) 定义导数函数 f(W,t)f(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: sigmady_dt x*(r-z) - y # dy/dt x*(r-z)-y, r:rhodz_dt x*y - b*z # dz/dt x*y - b*z, b;betareturn np.array([dx_dt,dy_dt,dz_dt])定义初值 W0W_0W0 和 WWW 的定义区间 [t0,t][t_0,\ t][t0, t] 调用 odeint() 求 WWW 在定义区间 [t0,t][t_0,\ t][t0, t] 的数值解。 注意例程中通过 argsparas 或 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 youcansx, y, z W # W[x,y,z]dx_dt p*(y-x) # dx/dt p*(y-x), p: sigmady_dt x*(r-z) - y # dy/dt x*(r-z)-y, r:rhodz_dt x*y - b*z # dz/dt x*y - b*z, b;betareturn 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, argsparas) # 通过 paras 传递导数函数的参数# 绘图
fig plt.figure()
ax Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2], colormagenta) # 绘制轨迹 1
ax.plot(track2[:,0], track2[:,1], track2[:,2], colordeepskyblue) # 绘制轨迹 2
ax.set_title(Lorenz attractor by scipy.integrate.odeint)
plt.show()4.4 洛伦兹Lorenz方程问题 Python 例程运行结果 5. 实例3Scipy 求解高阶常微分方程
高阶常微分方程必须做变量替换化为一阶微分方程组再用 odeint 求数值解。
5.1 例题 3求二阶 RLC 振荡电路的数值解
零输入响应的 RLC 振荡电路可以由如下的二阶微分方程描述
{d2udt2RL∗dudt1LC∗u0u(0)U0u′(0)0\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} ⎩⎪⎪⎪⎨⎪⎪⎪⎧dt2d2uLR∗dtduLC1∗u0u(0)U0u′(0)0
令 αR/2L\alpha R/2LαR/2L、ω021/LC\omega_0^21/LCω021/LC在零输入响应 us0u_s0us0 时上式可以写成
{d2udt22αdudtω02u0u(0)U0u′(0)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} ⎩⎪⎪⎪⎨⎪⎪⎪⎧dt2d2u2αdtduω02u0u(0)U0u′(0)0 对二阶微分方程问题引入变量 vdu/dtv {du}/{dt}vdu/dt通过变量替换就把原方程化为如下的微分方程组
{dudtvdvdt−2αv−ω02uu(0)U0v(0)0\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} ⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧dtduvdtdv−2αv−ω02uu(0)U0v(0)0
这样就可以用上节求解微分方程组的方法来求解高阶微分方程问题。 5.2 二阶微分方程问题的编程步骤
以RLC 振荡电路为例讲解 scipy.integrate.odeint() 求解高阶常微分方程初值问题的步骤 导入 scipy、numpy、matplotlib 包 定义导数函数 deriv(Y, t, a, w) 注意 odeint() 函数中定义导数函数的标准形式是 f(y,t)f(y,t)f(y,t) 本问题中 y 表示向量记为 Y[u,v]Y[u,v]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定义初值 Y0[u0,v0]Y_0[u_0,v_0]Y0[u0,v0] 和 YYY 的定义区间 [t0,t][t_0,\ t][t0, t] 调用 odeint() 求 Y[u,v]Y[u,v]Y[u,v] 在定义区间 [t0,t][t_0,\ t][t0, t] 的数值解。 例程中通过 argsparas 将参数 (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_dtt 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, argsparas1) # args 设置导数函数的参数
Y2 odeint(deriv, Y0, t, argsparas2) # args 设置导数函数的参数
Y3 odeint(deriv, Y0, t, argsparas3) # args 设置导数函数的参数
# W2 (0.0, 1.01, 0.0) # 定义初值为 W2
# track2 odeint(lorenz, W2, t, argsparas) # 通过 paras 传递导数函数的参数# 绘图
plt.plot(t, Y1[:, 0], r-, labelu1(t))
plt.plot(t, Y2[:, 0], b-, labelu2(t))
plt.plot(t, Y3[:, 0], g-, labelu3(t))
plt.plot(t, Y1[:, 1], r:, labelv1(t))
plt.plot(t, Y2[:, 1], b:, labelv2(t))
plt.plot(t, Y3[:, 1], g:, labelv3(t))
plt.axis([0, 20, -0.8, 1.2])
plt.legend(locbest)
plt.title(Second ODE by scipy.integrate.odeint)
plt.show()5.4 二阶方程问题 Python 例程运行结果 结果讨论
RLC串联电路是典型的二阶系统在零输入条件下根据 α\alphaα 与 ω\omegaω 的关系电路的输出响应存在四种情况
过阻尼 α2−ω20\alpha^2 - \omega^20α2−ω20 有 2 个不相等的负实数根临界阻尼 α2−ω20\alpha^2 - \omega^2 0α2−ω20有 2 个相等的负实数根欠阻尼 α2−ω20\alpha^2 - \omega^2 0α2−ω20有一对共轭复数根无阻尼R0R0R0有一对纯虚根。
例程中所选择的 3 组参数分别对应过阻尼、临界阻尼和欠阻尼的条件微分方程的数值结果很好地体现了不同情况的相应曲线。 6. 小结
小白首先要有信心用 Scipy 工具包求解标准形式的微分方程组编程实现是非常简单、容易掌握的。其次要认识到由于微分方程的解的特性是与模型参数密切相关的不同参数的解可能具有完全不同的形态。这就涉及模型稳定性、灵敏度的分析很容易使论文写的非常丰富和精彩。不会从问题建立微分方程模型怎么办不会展开参数对稳定性、灵敏度的影响进行讨论怎么办谁让你自己做呢当然是先去找相关专业的教材、论文从中选择比较接近、比较简单的理论和模型然后通过各种假设强行将题目简化为模型中的条件这就可以照猫画虎了。
更多微分方程数学模型案例参见 新冠疫情 模型系列 Python小白的数学建模课-B2. 新冠疫情 SI模型 Python小白的数学建模课-B3. 新冠疫情 SIS模型 Python小白的数学建模课-B4. 新冠疫情 SIR模型 Python小白的数学建模课-B5. 新冠疫情 SEIR模型 Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型 【本节完】 版权声明
欢迎关注『Python小白的数学建模课 Youcans』 原创作品
原创作品转载必须标注原文链接(https://blog.csdn.net/youcans/article/details/117702996)。
Copyright 2021 Youcans, XUPT
Crated2021-06-08 欢迎关注 『Python小白的数学建模课 Youcans』 系列持续更新 Python小白的数学建模课-01.新手必读 Python小白的数学建模课-02.数据导入 Python小白的数学建模课-03.线性规划 Python小白的数学建模课-04.整数规划 Python小白的数学建模课-05.0-1规划 Python小白的数学建模课-06.固定费用问题 Python小白的数学建模课-07.选址问题 Python小白的数学建模课-09.微分方程模型 Python小白的数学建模课-10.微分方程边值问题 Python小白的数学建模课-12.非线性规划 Python小白的数学建模课-15.图论的基本概念 Python小白的数学建模课-16.最短路径算法 Python小白的数学建模课-17.条件最短路径算法 Python小白的数学建模课-18.最小生成树问题 Python小白的数学建模课-19.网络流优化问题 Python小白的数学建模课-20.网络流优化案例 Python小白的数学建模课-A1.国赛赛题类型分析 Python小白的数学建模课-A2.2021年数维杯C题探讨 Python小白的数学建模课-A3.12个新冠疫情数模竞赛赛题及短评 Python小白的数学建模课-B2. 新冠疫情 SI模型 Python小白的数学建模课-B3. 新冠疫情 SIS模型 Python小白的数学建模课-B4. 新冠疫情 SIR模型 Python小白的数学建模课-B5. 新冠疫情 SEIR模型 Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型 Python数模笔记-PuLP库 Python数模笔记-StatsModels统计回归 Python数模笔记-Sklearn Python数模笔记-NetworkX Python数模笔记-模拟退火算法
