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

前言

這章主要介紹解決常微分方程障般，微分方程組合和邊值問題的方法备埃。

學(xué)習過程

<1>初值問題 initial value problem $I.V.P.$

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

<2>Lipschitz條件

給出矩形區(qū)域 $R=\{(t,y):a \leq t \leq b,c \leq y \leq \}$ 嘴拢，假設(shè) $f(t,y)$ 在 $R$ 上連續(xù)，且存在一個常量 $L>0$ 滿足性質(zhì) $|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)$ 在子區(qū)間 $t_0 \leq t \leq t_0+\delta$ 上有唯一解 $y=y(t)$ 吞加。

<3>Euler法

微分方程不總是那么容易求解的，為此，我們會在接下來提出很多方法來求解衔憨，其中之一就是Euler法叶圃。但是Euler法實際上使用是有限的因為它有很大誤差，但它具有教育意義巫财。
我們提出一個函數(shù)滿足初值問題 $I.V.P.$ 盗似，即在區(qū)間 $[a,b]$ 下， $y'=f(t,y),y(a)=y_0$ 平项。在我們無法找出函數(shù) $y=y(t)$ 的情況下赫舒，我們可以找出一系列點來逼近函數(shù)。
設(shè) $h=\frac{b-a}{M}$ 闽瓢，我們將區(qū)間分為 $M$ 個子區(qū)間接癌，并假設(shè) $y(t)$ ， $y'(t)$ 扣讼， $y''(t)$ 在區(qū)域連續(xù)缺猛，然后使用泰勒定理，我們可以得到以下公式：
$y(t)=y(t_0)+y'(t_0)(t-t_0)+\frac{y''(c_1)(t-t_0)^2}{2}.$
把 $t_1$ 代入椭符，且 $h$ 十分小時荔燎，二次項可以忽略不計，公式就會變?yōu)椋?br> $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法销钝。

我們把微分方程求解的方法叫做微分方法或者是離散變量方法有咨。這些方法都是找出一系列逼近函數(shù)的離散點集蒸健。其中一步法有這種形式 $y_{k+1}=y_k+h\Phi (t_k,y_k)$ 座享， $\Phi$ 叫做增量函數(shù)。
當使用這種離散變量方法來求解初值問題時似忧，誤差會有兩種來源：離散化渣叛，四舍五入。

假設(shè)點集 $\{(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)$ 瓶籽。于是我們利用梯形規(guī)則假設(shè) $y(t_1)\approx y(t_0)+\frac{h}{2}(f(t_0,y(t_0))+f(t_1,y(t_1)))$ 匠童。再結(jié)合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>泰勒序列方法

假設(shè) $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法

在該方法中挤牛，我們假設(shè)一個區(qū)間內(nèi)再取出多個點莹痢，計算它們的加權(quán)平均斜率。并且使得它們精度和泰勒序列前幾項一致墓赴。這個方法的好處就是不需要計算更高階的導(dǎo)數(shù)竞膳。
$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：可以知道步長越小，誤差會越小藕届，但是這會導(dǎo)致大量的計算挪蹭。因此我們提出了RKF45，它是找到一個恰當步長的程序休偶。在每步中梁厉，對解的兩種不同逼近會用來比較，如果兩個解答案很近椅贱，則步長是可以接受的懂算。如果兩個解答案準確率不太一致只冻，則繼續(xù)縮小步長庇麦。當需要更多有效位數(shù)的解時，則增加步長喜德。
上述的方法都是一步法山橄。

<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}]$ 內(nèi)積分。
先使用Adams-Bashforth預(yù)測舍悯，在 $[t_k,t_{k+1}]$ 內(nèi)積分：
$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}]$ 內(nèi)積分：
$y_{k+1}=y_k+\frac{h}{24}(f_{k-2}-5f_{k-1}+19f_k+9f_{k+1})$
其中預(yù)測和矯正的誤差都是 $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}]$ 內(nèi)積分朴艰。
先使用Milne預(yù)測观蓄，在 $[t_{k-3},t_{k+1}]$ 內(nèi)積分：
$p_{k+1}=y_{k-3}+\frac{4h}{3}(2f_{k-2}-f_{k-1}+2f_k)$
再用Simpson矯正，在 $[t_{k-1},t_{k+1}]$ 內(nèi)積分：
$y_{k+1}=y_{k-1}+\frac{h}{3}(f_{k-1}+4f_{k}+f_{k+1})$
其中預(yù)測和矯正的誤差都是 $O(h^5)$ 祠墅。

<9>微分方程組

在本小節(jié)侮穿，我們考慮到以下的初值問題：
$\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$ 個子區(qū)間后亲茅，
$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$
遇到更高階的微分方程組時，我們可以設(shè)置多個未知量對應(yīng)不同階的導(dǎo)數(shù)狗准，相當于n維求歐拉公式克锣。在這里，龍格庫塔公式也是可以的驶俊，一般4階的效果比較好娶耍。

<10>邊值問題

<11>初值問題

詞匯學(xué)習

modification 修改
trajectory 彈道
conjecture 推測
interpretation 解釋
mesh 網(wǎng)眼
discrete 離散的
increment 增量
predominate 占主導(dǎo)地位
trapezoidal 梯形的
trade-off 交易
elaborate 詳盡的
omit 忽略

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末，一起剝皮案震驚了整個濱河市饼酿，隨后出現(xiàn)的幾起案子榕酒，更是在濱河造成了極大的恐慌，老刑警劉巖故俐，帶你破解...
沈念sama閱讀 211,194評論 6贊 490
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件想鹰，死亡現(xiàn)場離奇詭異，居然都是意外死亡药版，警方通過查閱死者的電腦和手機辑舷，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 90,058評論 2贊 385
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來槽片，“玉大人何缓，你說我怎么就攤上這事』顾ǎ” “怎么了碌廓？”我有些...
開封第一講書人閱讀 156,780評論 0贊 346
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長剩盒。經(jīng)常有香客問我谷婆，道長，這世上最難降的妖魔是什么辽聊？我笑而不...
開封第一講書人閱讀 56,388評論 1贊 283
?港島之戀（遺憾婚禮）
正文為了忘掉前任纪挎，我火速辦了婚禮，結(jié)果婚禮上跟匆，老公的妹妹穿的比我還像新娘异袄。我一直安慰自己，他們只是感情好玛臂，可當我...
茶點故事閱讀 65,430評論 5贊 384
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布烤蜕。她就那樣靜靜地躺著埠帕，像睡著了一般。火紅的嫁衣襯著肌膚如雪玖绿。梳的紋絲不亂的頭發(fā)上敛瓷，一...
開封第一講書人閱讀 49,764評論 1贊 290
城市分裂傳說
那天，我揣著相機與錄音斑匪，去河邊找鬼呐籽。笑死，一個胖子當著我的面吹牛蚀瘸，可吹牛的內(nèi)容都是我干的狡蝶。我是一名探鬼主播，決...
沈念sama閱讀 38,907評論 3贊 406
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼贮勃，長吁一口氣：“原來是場噩夢啊……” “哼贪惹！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起寂嘉，我...
開封第一講書人閱讀 37,679評論 0贊 266
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤奏瞬，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后泉孩，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體硼端，經(jīng)...
沈念sama閱讀 44,122評論 1贊 303
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 36,459評論 2贊 325
?白月光啟示錄
正文我和宋清朗相戀三年寓搬，在試婚紗的時候發(fā)現(xiàn)自己被綠了珍昨。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
茶點故事閱讀 38,605評論 1贊 340
活死人
序言：一個原本活蹦亂跳的男人離奇死亡句喷，死狀恐怖镣典，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情唾琼，我是刑警寧澤兄春，帶...
沈念sama閱讀 34,270評論 4贊 329
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站父叙，受9級特大地震影響神郊，放射性物質(zhì)發(fā)生泄漏肴裙。R本人自食惡果不足惜趾唱，卻給世界環(huán)境...
茶點故事閱讀 39,867評論 3贊 312
男人毒藥：我在死后第九天來索命
文/蒙蒙一、第九天我趴在偏房一處隱蔽的房頂上張望蜻懦。院中可真熱鬧甜癞，春花似錦、人聲如沸宛乃。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,734評論 0贊 21
一樁弒父案蒸辆，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽。三九已至析既，卻和暖如春躬贡，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背眼坏。一陣腳步聲響...
開封第一講書人閱讀 31,961評論 1贊 265
情欲美人皮
我被黑心中介騙來泰國打工拂玻，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留，地道東北人宰译。一個月前我還...
沈念sama閱讀 46,297評論 2贊 360
代替公主和親
正文我出身青樓檐蚜，卻偏偏與公主長得像，于是被迫代替她去往敵國和親沿侈。傳聞我的和親對象是個殘疾皇子闯第，可洞房花燭夜當晚...
茶點故事閱讀 43,472評論 2贊 348