《Numerical Methods Using MATLAB》（第四版）第九章：常微分方程求解

前言

这章主要介绍解决常微分方程，微分方程组合和边值问题的方法。

学习过程

<1>初值问题 initial value problem $I.V.P.$

常微分方程的一般形式： $\frac{dy}{dx}=f(x,y)$
方程会因为初值不同有变化。

<2>Lipschitz条件

给出矩形区域 $R=\{(t,y):a \leq t \leq b,c \leq y \leq \}$ ，假设 $f(t,y)$ 在 $R$ 上连续，且存在一个常量 $L>0$ 满足性质 $|f(t,y_1)-f(t,y_2)| \leq L|y_1-y_2|$ ，其中任意 $(t,y_1),(t,y_2)∈R$ ，那么 $f$ 可以说是在 $R$ 上满足 $Lipschitz$ 条件。而 $L$ 也可以叫作 $Lipschitz$ 常量。

为了更好的判断，我们一般使用下面的方法。

$f$ 在 $R$ 上定义，如果这里存在常量 $L>0$ 使得所有 $(t,y)∈R,|f_y(t,y)| \leq L$ ，则 $f$ 在 $R$ 上满足 $Lipschitz$ 条件。

如果 $f$ 满足 $Lipschitz$ 条件且有初值，那么 $y'=f(t,y)$ 在子区间 $t_0 \leq t \leq t_0+\delta$ 上有唯一解 $y=y(t)$ 。

<3>Euler法

微分方程不总是那么容易求解的，为此，我们会在接下来提出很多方法来求解，其中之一就是Euler法。但是Euler法实际上使用是有限的因为它有很大误差，但它具有教育意义。
我们提出一个函数满足初值问题 $I.V.P.$ ，即在区间 $[a,b]$ 下， $y'=f(t,y),y(a)=y_0$ 。在我们无法找出函数 $y=y(t)$ 的情况下，我们可以找出一系列点来逼近函数。
设 $h=\frac{b-a}{M}$ ，我们将区间分为 $M$ 个子区间，并假设 $y(t)$ ， $y'(t)$ ， $y''(t)$ 在区域连续，然后使用泰勒定理，我们可以得到以下公式：
$y(t)=y(t_0)+y'(t_0)(t-t_0)+\frac{y''(c_1)(t-t_0)^2}{2}.$
把 $t_1$ 代入，且 $h$ 十分小时，二次项可以忽略不计，公式就会变为：
$y_1=y_0+hf(t_0,y_0).$
递推下去就是：
$y_{k+1}=y_k+hf(t_k,y_k),k=0,1,...,M-1.$
这就是Euler法。

我们把微分方程求解的方法叫做微分方法或者是离散变量方法。这些方法都是找出一系列逼近函数的离散点集。其中一步法有这种形式 $y_{k+1}=y_k+h\Phi (t_k,y_k)$ ， $\Phi$ 叫做增量函数。
当使用这种离散变量方法来求解初值问题时，误差会有两种来源：离散化，四舍五入。

假设点集 $\{(t_k,y_k)\}_{k=0}^M$ ，初值问题有唯一解 $y=y(t)$ 。
全局离散误差： $e_k=y(t_k)-y_k,k=0,1,...,M.$
局部离散误差： $\epsilon_{k+1}=y(t_{k+1})-y_k-h\Phi(t_k,y_k),k=0,1,...,M-1.$

Euler法中，在满足以上初值问题条件，且 $y(t)∈C^2[t_0,b]$ 时，有：
$|e_k|=|y(t_k)-y_k|=O(h),$

$|\epsilon_{k+1}|=|y(t_{k+1})-y_k-hf(t_k,y_k)|=O(h^2).$

最终全局误差 $F.G.E. =E(y(b),h)=|y(b)-y_M|=O(h).$

<4>Heun法

我们使用积分的方法来解决初值问题。一开始有 $\int_{t_0}^{t_1}f(t,y(t))dt=y(t_1)-y(t_0)$ 。于是我们利用梯形规则假设 $y(t_1)\approx y(t_0)+\frac{h}{2}(f(t_0,y(t_0))+f(t_1,y(t_1)))$ 。再结合Euler法。
即Heun法为：
$p_{k+1}=y_k+hf(t_k,y_k),t_{k+1}=t_k+h,$

$y_{k+1}=y_k+\frac{h}{2}(f(t_k,y_k)+f(t_{k+1},p_{k+1})).$

Heun法中，在满足以上初值问题条件，且 $y(t)∈C^3[t_0,b]$ 时，有：
$|e_k|=|y(t_k)-y_k|=O(h^2),$

$|\epsilon_{k+1}|=|y(t_{k+1})-y_k-h\Phi(t_k,y_k)|=O(h^3).$

最终全局误差 $F.G.E. =E(y(b),h)=|y(b)-y_M|=O(h^2).$

<5>泰勒序列方法

假设 $y(t)∈C^{N+1}[t_0,b]$ ，并且 $y(t)$ 在 $t=t_k∈[t_0,b]$ 上有 $N$ 阶展开：
$y(t_k+h)=y(t_k)+hT_N(t_k,y(t_k))+O(h^{N+1}),$

$T_N(t_k,y(t_k))=\sum_{j=1}^N\frac{y^{(j)}(t_k)}{j!}h^{j-1}$

泰勒序列法中，在满足以上初值问题条件，且 $y(t)∈C^{N+1}[t_0,b]$ 时，有：
$|e_k|=|y(t_k)-y_k|=O(h^N),$

$|\epsilon_{k+1}|=|y(t_{k+1})-y_k-hT_N(t_k,y_k)|=O(h^{N+1}).$

最终全局误差 $F.G.E. =E(y(b),h)=|y(b)-y_M|=O(h^N).$

<6>Runge-Kutta法

在该方法中，我们假设一个区间内再取出多个点，计算它们的加权平均斜率。并且使得它们精度和泰勒序列前几项一致。这个方法的好处就是不需要计算更高阶的导数。
$N$ 阶通式：
$y_{i+1}=y_i+\Phi h,$

$\Phi=w_1k_1+w_2k_2+...+w_nk_n,$

$k_1=f(t_i,y_i),$

$k_2=f(t_i+p_1h,y_i+q_{11}k_1h),$

$k_3=f(t_i+p_2h,y_i+q_{21}k_1h+q_{22}k_2h),$

$k_n=f(t_i+p_{n-1}h,y_i+q_{n-1,1}k_1h+...),$

$w$ 之和为1， $p$ 等于 $q$ 之和。
Runge-Kutta法中，在满足以上初值问题条件，且 $y(t)∈C^{5}[t_0,b]$ 时，即4阶，有：
$|e_k|=|y(t_k)-y_k|=O(h^4),$

$|\epsilon_{k+1}|=|y(t_{k+1})-y_k-hT_N(t_k,y_k)|=O(h^{5}).$

最终全局误差 $F.G.E. =E(y(b),h)=|y(b)-y_M|=O(h^4).$

例子：
RK4：模仿了四阶泰勒
$y_{k+1}=y_k+\frac{h(f_1+2f_2+2f_3+f_4)}{6},$
$f_1=f(t_k,y_k),$
$f_2=f(t_k+\frac{h}{2},y_k+\frac{h}{2}f_1),$
$f_3=f(t_k+\frac{h}{2},y_k+\frac{h}{2}f_2),$
$f_4=f(t_k+h,y_k+hf_3).$
RKF45：可以知道步长越小，误差会越小，但是这会导致大量的计算。因此我们提出了RKF45，它是找到一个恰当步长的程序。在每步中，对解的两种不同逼近会用来比较，如果两个解答案很近，则步长是可以接受的。如果两个解答案准确率不太一致，则继续缩小步长。当需要更多有效位数的解时，则增加步长。
上述的方法都是一步法。

<7>Adams-Bashforth-Moulton 法

该方法是基于积分公式来推的：
$y(t_{k+1})=y(t_k)+\int_{t_k}^{t_{k+1}}f(t,y(t))dt$
在 $[t_k,t_{k+1}]$ 内积分。
先使用Adams-Bashforth预测，在 $[t_k,t_{k+1}]$ 内积分：
$p_{k+1}=y_k+\frac{h}{24}(-9f_{k-3}+37f_{k-2}-59f_{k-1}+55f_k)$
再用Adams-Moulton矫正，在 $[t_k,t_{k+1}]$ 内积分：
$y_{k+1}=y_k+\frac{h}{24}(f_{k-2}-5f_{k-1}+19f_k+9f_{k+1})$
其中预测和矫正的误差都是 $O(h^5)$ 。

当需要更精确时，就缩小步长，而缩小步长，就需要新的初始值。

<8>Milne-Simpson法

该方法是基于积分公式来推的：
$y(t_{k+1})=y(t_{k-3})+\int_{t_{k-3}}^{t_{k+1}}f(t,y(t))dt$
在 $[t_{k-3},t_{k+1}]$ 内积分。
先使用Milne预测，在 $[t_{k-3},t_{k+1}]$ 内积分：
$p_{k+1}=y_{k-3}+\frac{4h}{3}(2f_{k-2}-f_{k-1}+2f_k)$
再用Simpson矫正，在 $[t_{k-1},t_{k+1}]$ 内积分：
$y_{k+1}=y_{k-1}+\frac{h}{3}(f_{k-1}+4f_{k}+f_{k+1})$
其中预测和矫正的误差都是 $O(h^5)$ 。

<9>微分方程组

在本小节，我们考虑到以下的初值问题：
$\begin{cases}\frac{dx}{dt}=f(t,x,y) \\ \frac{dy}{dt}=g(t,x,y)\end{cases} with \begin{cases}x(t_0)=x_0 \\ y(t_0)=y_0\end{cases}$
如果使用Euler法来解决这个问题的话，有：
$dx=f(t,x,y)dt,dy=g(t,x,y)dt$
其中 $dt=t_{k+1}-t_k,dx=x_{k+1}-x_k,dy=y_{k+1}-y_k$
于是，
$x_{k+1}-x_k \approx f(t_k,x_k,y_k)(t_{k+1}-t_k)$
$y_{k+1}-y_k \approx g(t_k,x_k,y_k)(t_{k+1}-t_k)$
分为步长为 $h$ 的 $M$ 个子区间后，
$t_{k+1}=t_k+h$
$x_{k+1}=x_k+hf(t_k,x_k,y_k)$
$y_{k+1}=y_k+hg(t_k,x_k,y_k),k=0,1,...,M-1$
遇到更高阶的微分方程组时，我们可以设置多个未知量对应不同阶的导数，相当于n维求欧拉公式。在这里，龙格库塔公式也是可以的，一般4阶的效果比较好。

<10>边值问题

<11>初值问题

词汇学习

modification 修改
trajectory 弹道
conjecture 推测
interpretation 解释
mesh 网眼
discrete 离散的
increment 增量
predominate 占主导地位
trapezoidal 梯形的
trade-off 交易
elaborate 详尽的
omit 忽略