Numerical Methods Using MATLAB(4版)-09-ODE

前言

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

学习过程

<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 忽略

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 199,902评论 5赞 468
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 84,037评论 2赞 377
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 146,978评论 0赞 332
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 53,867评论 1赞 272
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 62,763评论 5赞 360
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,104评论 1赞 277
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,565评论 3赞 390
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,236评论 0赞 254
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,379评论 1赞 294
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,313评论 2赞 317
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,363评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,034评论 3赞 315
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,637评论 3赞 303
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,719评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,952评论 1赞 255
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,371评论 2赞 346
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 41,948评论 2赞 341