Scipy 微分方程type
status
date
slug
summary
tags
category
icon
password
Property
基本概念
微分方程是描述系统的状态随时间和空间演化的数学工具。物理中许多涉及变力的运动学、动力学问题，如空气的阻力为速度函数的落体运动等问题，很多可以用微分方程求解。微分方程在化学、工程学、经济学和人口统计等领域也有广泛应用。
具体来说，微分方程是指含有未知函数及其导数的关系式。
- 微分方程按自变量个数分为：只有一个自变量的常微分方程（Ordinary Differential Equations）和包含两个或两个以上独立变量的偏微分方程（Partial Differential Equations）。
- 微分方程按阶数分为：一阶、二阶、高阶，微分方程的阶数取决于方程中最高次导数的阶数。
- 微分方程还可以分为：（非）齐次，常（变）系数，（非）线性，初值问题/边界问题…
微分方程的数值解法
只有很少的微分方程可以解析求解，大多数的微分方程只能采用数值方法进行求解。Python有个符号运算工具包SymPy，以解析方式求解积分、微分方程，也就是说给出的结果是微分方程的解析解表达式，很牛，但只能求解有解析解的微分方程
微分方程的数值求解是先把时间和空间离散化，然后将微分化为差分，建立递推关系，然后反复进行迭代计算，得到任意时间和空间的值。
求解常微分方程的基本方法，有欧拉法、龙格库塔法等。
微分方程初值
一阶常微分方程（组）模型
给定初始条件的一阶常微分方程（组）的标准形式是：
式中的在常微分方程中是标量，在常微分方程组中是数组向量。
scipy.integrate.odeint() 函数
SciPy 提供了两种方式求解常微分方程：基于
odeint 函数的 API 比较简单易学，基于 ode 类的面向对象的 API 更加灵活。scipy.integrate.odeint() 是求解微分方程的具体方法，通过数值积分来求解常微分方程组。在 odeint 函数内部使用 FORTRAN 库 odepack 中的 lsoda，可以求解一阶刚性系统和非刚性系统的初值问题。官网介绍详见： odeint 的主要参数：
求解标准形式的微分方程（组）主要使用前三个参数：
- func: callable(y, t, …) 导数函数 ，即 y 在 t 处的导数，以函数的形式表示
- y0: array： 初始条件，对于常微分方程组 则为数组向量
- t: array： 求解函数值对应的时间点的序列。序列的第一个元素是与初始条件对应的初始时间 ；时间序列必须是单调递增或单调递减的，允许重复值。
其它参数简介如下：
- args: 向导数函数 func 传递参数。当导数函数 包括可变参数 p1,p2… 时，通过 args =(p1,p2,…) 可以将参数p1,p2… 传递给导数函数 func。
- Dfun: func 的雅可比矩阵，行优先。如果 Dfun 未给出，则算法自动推导。
- col_deriv: 自动推导 Dfun的方式。
- printmessg: 布尔值。控制是否打印收敛信息。
- 其它参数用于控制求解算法的参数，一般情况可以忽略。
odeint 的主要返回值：
- y: array 数组，形状为 (len(t),len(y0)，给出时间序列t中每个时刻的y值。
例题 1：求微分方程的数值解
例题 2：求洛伦兹（Lorenz）方程的数值解
洛伦兹（Lorenz）混沌吸引子的轨迹可以由如下的 3个微分方程描述：
洛伦兹方程将大气流体运动的强度 x 与水平和垂直方向的温度变化 y 和 z 联系起来，进行大气对流系统的模拟，现已广泛应用于天气预报、空气污染和全球气候变化的研究。参数 σ 称为普兰特数，ρ 是规范化的瑞利数，β 和几何形状相关。洛伦兹方程是非线性微分方程组，无法求出解析解，只能使用数值方法求解。
注意 odeint() 函数中定义导数函数的标准形式是 ，对于微分方程组 y 表示向量。
为避免混淆，我们记为 ，函数 lorenz(W,t) 定义导数函数 。
用 p,r,b 分别表示方程中的参数 ，
通过 args=paras 或 args = (10.0,28.0,3.0) 将参数 (p,r,b) 传递给导数函数 lorenz(W,t,p,r,b)。参数 (p,r,b) 当然也可以不作为函数参数传递，而是在导数函数 lorenz() 中直接设置。但参数传递方法，使导数函数结构清晰、更为通用。另外，对于可变参数问题，使用这种参数传递方式就非常方便。
Scipy 求解高阶常微分方程
高阶常微分方程，必须做变量替换，化为一阶微分方程组，再用 odeint 求数值解。
例题 3：求二阶 RLC 振荡电路的数值解
零输入响应的 RLC 振荡电路可以由如下的二阶微分方程描述：
令 、 ，在零输入响应时上式可以写成：
对二阶微分方程问题，引入变量 ，通过变量替换就把原方程化为如下的微分方程组：
这样就可以用之前求解微分方程组的方法来求解高阶微分方程问题。
RLC串联电路是典型的二阶系统，在零输入条件下根据与的关系，电路的输出响应存在四种情况：
- 过阻尼： 有 2 个不相等的负实数根；
- 临界阻尼： ，有 2 个相等的负实数根；
- 欠阻尼： ，有一对共轭复数根；
- 无阻尼： ，有一对纯虚根。
所选择的 3 组参数分别对应过阻尼、临界阻尼和欠阻尼的条件，微分方程的数值结果很好地体现了不同情况的相应曲线。
微分方程边值
常微分方程边值问题的数学模型
只含边界条件作为定解条件的常微分方程求解问题，称为常微分方程的边值问题（boundary value problem）。
一般形式的二阶常微分方程边值问题：
有三种情况的边界条件：
（1）第一类边界条件（两点边值问题）：
（2）第二类边界条件：
（3）第三类边界条件：
其中：
常微分方程边值问题的数值解法
简单介绍求解常微分方程边值问题的数值解法，常用方法有：打靶算法、有限差分法和有限元法。打靶算法把边值问题转化为初值问题求解，是根据边界条件反复迭代调整初始点的斜率，使初值问题的数值解在边界上“命中”问题的边值条件。有限差分法把空间离散为网格节点，用差商代替微商，将微分方程离散化为线性或非线性方程组来求解。 有限元法将微分方程离散化，有限元就是指近似连续域的离散单元，对每一单元假定一个近似解，然后推导求解域满足条件，从而得到问题的解。
SciPy 求解常微分方程边值问题
BVP 问题的标准形式
Scipy 用
solve_bvp() 函数求解常微分方程的边值问题，定义微分方程的标准形式为：因此要将第一类边界条件 改写为：
scipy.integrate.solve_bvp() 函数
scipy.integrate.solve_bvp() 是求解微分方程边值问题的具体方法，可以求解一阶微分方程（组）的两点边值问题（第一类边界条件）。在 odeint 函数内部使用 FORTRAN 库 odepack 中的 lsoda，可以求解一阶刚性系统和非刚性系统的初值问题。官网介绍详见： solve_bvp 的主要参数：
求解标准形式的微分方程（组）主要使用前 4 个参数：
- func: callable fun(x, y, …) 导数函数 ， y 在 x 处的导数，以函数的形式表示。可以带有参数 p。
- bc: callable bc(ya, yb, …) 边界条件，y 在两点边界的函数，以函数的形式表示。可以带有参数 p。
- x: array： 初始网格的序列，shape (m,)。必须是单调递增的实数序列，起止于两点边界值 xa，xb。
- y: array： 网格节点处函数值的初值，shape (n,m)，第 i 列对应于 x[i]。
- p: array： 可选项，向导数函数 func、边界条件函数 bc 传递参数。
其它参数用于控制求解算法的参数，一般情况可以忽略。
solve_bvp 的主要返回值：
- sol: PPoly 通过 PPoly （如三次连续样条函数）插值求出网格节点处的 y 值。
- x: array 数组，形状为 (m,)，最终输出的网格节点。
- y: array 二维数组，形状为 (n,m)，输出的网格节点处的 y 值。
- yp: array 二维数组，形状为 (n,m)，输出的网格节点处的 y’ 值。
例题 1：求常微分方程边值问题的数值解
引入变量 ，通过变量替换就把原方程化为如下的标准形式的微分方程组：
这样就可以用
solve_bvp()求解该常微分方程的边值问题。注意本问题中 y 表示向量，记为 ，导数定义函数 dydx(x,y) 如下
定义边界条件函数 boundCond(ya,yb)
例题 2：水滴横截面形状问题
水平面上的水滴横截面形状，可以用如下的微分方程描述：
引入变量 ，通过变量替换就把原方程化为如下的标准形式的微分方程组：
这样就可以用
solve_bvp()求解该常微分方程的边值问题。
例题 3：Mathieu 方程的特征函数
Mathieu 在研究椭圆形膜的边界值问题时，导出了一个二阶常微分方程，其形式为：
用这种形式的数学方程可以描述自然中的物理现象，包括振动椭圆鼓、四极质谱仪和四极离子阱、周期介质中的波动、强制振荡器参数共振现象、广义相对论中的平面波解决方案、量子摆哈密顿函数的本征函数、旋转电偶极子的斯塔克效应。
式中 是两个实参数，方程的系数是以 π 或 2π 为周期的，但只有在 满足一定关系时 Mathieu 方程才有周期解。
引入变量 ，通过变量替换就把原方程化为如下的标准形式的微分方程组：
这样就可以用
solve_bvp()求解该常微分方程的边值问题。需要注意的是：（1）本案例涉及一个待定参数 λ 需要通过 solve_bvp(fun, bc, x, y, p=None) 中的可选项 p 传递到导数函数和边界条件函数，（2）本案例涉及 3 个边界条件，要注意边界条件函数的定义。
初值 A从 0~3.0 变化时，y-x 曲线（图中虚线）几乎不变，但 y’-x 的振幅增大；当 A 再稍微增大，系统就进入不稳定区， y-x 曲线振荡发散