《數(shù)值分析與算法》第八講-常微分方程初值問題第一頁(yè)，共44頁(yè)。WenjianYu2第八章常微分方程初值問題

(主要是前3節(jié))第二頁(yè)，共44頁(yè)。常微分方程初值問題常微分方程基本概念3第三頁(yè)，共44頁(yè)。常微分方程

4第四頁(yè)，共44頁(yè)。常微分方程

初值問題:

5第五頁(yè)，共44頁(yè)。常微分方程–例子1含電容元件的電路問題通過節(jié)點(diǎn)分析法得到微分方程組問題的解反映了電容充/放電過程C1R2R1R3C2電流/電壓關(guān)系節(jié)點(diǎn)電流方程t6第六頁(yè)，共44頁(yè)。常微分方程–例子2雙聯(lián)擺的運(yùn)動(dòng)兩個(gè)擺錘(重物),剛性桿的重量可忽略不考慮摩擦力,運(yùn)動(dòng)不會(huì)停止,且在初始角

度較大時(shí)擺錘的軌跡呈現(xiàn)混沌現(xiàn)象求解微分方程初值問題,得到

擺錘的運(yùn)動(dòng)規(guī)律Matlab演示swinger

7第七頁(yè)，共44頁(yè)。常微分方程

線性齊次常系數(shù)微分方程

實(shí)際的問題基本上都是穩(wěn)定的!

(由于歷史原因)8第八頁(yè)，共44頁(yè)。常微分方程

局部穩(wěn)定9第九頁(yè)，共44頁(yè)。簡(jiǎn)單方法與有關(guān)概念簡(jiǎn)單的初值問題數(shù)值解法

10第十頁(yè)，共44頁(yè)。初值問題的數(shù)值解法

否則為多步法否則為隱格式方法11第十一頁(yè)，共44頁(yè)。歐拉法

“左矩形”求積公式

h=0.1h=0.050.11.0000001.0048370.051.0000000.31.0350920.21.0100001.0187310.11.0025000.351.0483370.31.0290001.0408180.151.0073750.41.0634200.41.0561001.0703200.21.0145060.451.0802490.51.0904901.1065310.251.0237810.51.098737

步長(zhǎng)h=0.1,和0.05步長(zhǎng)小的更準(zhǔn)

12第十二頁(yè)，共44頁(yè)。數(shù)值解法的穩(wěn)定性

-10

歐拉法解模型問題的穩(wěn)定區(qū)域

13第十三頁(yè)，共44頁(yè)。數(shù)值解法的穩(wěn)定性

-10

歐拉法穩(wěn)定

14第十四頁(yè)，共44頁(yè)。數(shù)值解法的穩(wěn)定性

00.0250.050.0750.10.1250.151-1.52.25-3.3755.0625-7.5937511.390610.0820850.0067380.0005534.5410-53.7310-63.0610-7這里設(shè)的h太大!計(jì)算結(jié)果如下表:15第十五頁(yè)，共44頁(yè)。數(shù)值解法的局部截?cái)嗾`差

整體誤差

穩(wěn)定的問題，整體誤差小于局部誤差之和不穩(wěn)定的問題呢？一般僅能控制局部誤差整體誤差？~局部誤差16第十六頁(yè)，共44頁(yè)。數(shù)值解法的局部截?cái)嗾`差

歐拉法是一階方法我們討論的所有方法都至少有1階準(zhǔn)確度數(shù)值解法的收斂性：隨著h0，誤差017第十七頁(yè)，共44頁(yè)。向后歐拉法與梯形法從數(shù)值積分的角度推導(dǎo)向后歐拉法:梯形法:兩者均為單步、隱格式方法,每步計(jì)算要求解(非線性)方程例8.6:用向后歐拉法求解右矩形

梯形

00.0250.050.0750.10.1250.1510.0066630.0019040.00054410.0820850.0067380.0005534.5410-53.7310-63.0610-718第十八頁(yè)，共44頁(yè)。向后歐拉法

準(zhǔn)確解

01穩(wěn)定區(qū)域

無條件穩(wěn)定(unconditionallystable)!19第十九頁(yè)，共44頁(yè)。向后歐拉法

具有1階準(zhǔn)確度!20第二十頁(yè)，共44頁(yè)。向后歐拉法與梯形法

穩(wěn)定的條件是:

思考無條件穩(wěn)定!

具有2階準(zhǔn)確度21第二十一頁(yè)，共44頁(yè)。簡(jiǎn)單方法與有關(guān)概念Runge-Kutta方法22第二十二頁(yè)，共44頁(yè)。在歐拉法基礎(chǔ)上改進(jìn)再增加一次函數(shù)求值:數(shù)值積分的中矩形或梯形公式利用歐拉法算半個(gè)步長(zhǎng)的結(jié)果,估算中點(diǎn)處被積函數(shù)值先用歐拉法估計(jì)區(qū)間終點(diǎn)處斜率,再用它與

起始點(diǎn)斜率的平均值算一整步Runge-Kutta方法

中矩形梯形公式(中點(diǎn)公式)(Heun方法/改進(jìn)的歐拉法)都比歐拉法準(zhǔn)確均為2級(jí)R-K公式

23第二十三頁(yè)，共44頁(yè)。

Runge-Kutta方法只能估算

,令

24第二十四頁(yè)，共44頁(yè)。

Runge-Kutta方法其他,用所有前面點(diǎn)的信息

積分節(jié)點(diǎn)………25第二十五頁(yè)，共44頁(yè)。幾種顯式R-K公式參數(shù)的值不按具體數(shù)值

積分公式設(shè)置,而根據(jù)

準(zhǔn)確度階數(shù)要求設(shè)置2級(jí)公式：改進(jìn)歐拉法、

中點(diǎn)公式經(jīng)典4級(jí)、4階Runge-Kutta法

(1905)Runge-Kutta方法

26第二十六頁(yè)，共44頁(yè)。

Runge-Kutta方法此時(shí)局部截?cái)嗾`差只要對(duì)非模型問題也有相同結(jié)論！如

27第二十七頁(yè)，共44頁(yè)。例8.7:用2階改進(jìn)歐拉、3階Ralston、

4階經(jīng)典R-K解問題,h=0.1,算到y(tǒng)(2)精確解為Runge-Kutta方法r4對(duì)應(yīng)的r級(jí)R-K公式有r階準(zhǔn)確度高于4階的公式很少單獨(dú)使用r>4對(duì)應(yīng)的r級(jí)R-K公式達(dá)不到r階準(zhǔn)確度二階Heun三階Ralston四階Runge-Kutta準(zhǔn)確值10.40.40.40.41.10.4756410.4746260.47463830.47463821.20.5834080.5813640.58138680.58138671.30.7281350.7250340.72506630.72506621.40.9153290.9111370.91117730.91117711.51.1511101.1457851.14583361.14583331.61.4421691.4356641.43572031.43572001.71.7957381.7880041.78806741.78806711.82.2195782.2105612.21063152.21063111.92.7219612.7116062.71168362.71168322.03.3116653.2999163.30000043.300000028第二十八頁(yè)，共44頁(yè)。

Runge-Kutta方法(此時(shí))一般顯式單步法

恰好與歐拉法一樣顯格式,都不是無條件穩(wěn)定的局部截?cái)嗾`差判斷單步法收斂性的簡(jiǎn)便方法

29第二十九頁(yè)，共44頁(yè)。簡(jiǎn)單方法與有關(guān)概念多步法

30第三十頁(yè)，共44頁(yè)。

多步法

(線性m步法)

固定步長(zhǎng)h31第三十一頁(yè)，共44頁(yè)。Taylor展開法求線性m步法的系數(shù)

這些系數(shù)應(yīng)該等于0例8.8:求兩步法公式中參數(shù)值多步法

滿足它們才可能收斂(相容性)

,

有二階準(zhǔn)確度

32第三十二頁(yè)，共44頁(yè)。

多步法

同例8.8的結(jié)果!

多步法公式中包含p個(gè)待定參數(shù)，至少可達(dá)到p-1階準(zhǔn)確度

33第三十三頁(yè)，共44頁(yè)。

多步法

Vandermonde陣T,非奇異

插值節(jié)點(diǎn)函數(shù)值

34第三十四頁(yè)，共44頁(yè)。常用的多步法公式Adams公式的推導(dǎo)：用插值多項(xiàng)式近似被積函數(shù)例8.10:推導(dǎo)m=4對(duì)應(yīng)的顯式Adams公式多步法

可證明滿足最高準(zhǔn)確度階數(shù)

類似地算其他系數(shù)，得

顯式四階Adams-

Bashforth公式(單項(xiàng)式函數(shù)代入法)35第三十五頁(yè)，共44頁(yè)。Adams公式幾種顯式公式幾種隱式公式多步法

階數(shù)穩(wěn)定閾值誤差常數(shù)11

-21/223/2-1/2

-15/12323/12-16/125/12

-6/113/8455/24-59/2437/24-9/24-3/10251/720歐拉法

階數(shù)穩(wěn)定閾值誤差常數(shù)11

--1/221/21/2

--1/1235/128/12-1/12

-6-1/2449/2419/24-5/241/24-3-19/720向后歐拉法梯形法并非無條件穩(wěn)定!36第三十六頁(yè)，共44頁(yè)。

多步法

37第三十七頁(yè)，共44頁(yè)。用Matlab求解初值問題用Matlab解ODE-IVP

38第三十八頁(yè)，共44頁(yè)。Matlab中的ODE-IVP求解器

[T,

Y,

TE,

YE,

IE]=solver(odefun,

tspan,

y0,

options)

39第三十九頁(yè)，共44頁(yè)。求解單個(gè)ODE火焰燃燒問題當(dāng)點(diǎn)燃一根火柴時(shí),火焰迅速增大直到一個(gè)臨界體積,然后維持這一體積不變,此時(shí)火焰內(nèi)部燃燒耗費(fèi)的氧氣和其表面現(xiàn)存的氧氣達(dá)到了一種平衡.火焰(近似為球)半徑y(tǒng)滿足ODE>>f=@(t,y)y^2-y^3;>>ode23(f,[0,2.0e4],1e-4)設(shè)初始半徑=0.0001例8.15嘗試用ode23s求解它40第四十頁(yè)，共44頁(yè)。求解ODE方程組

例8.16functionydot=myode2(t,y);ydot=[y(2);6*t];%列向量

41第四十一頁(yè)，共44頁(yè)。雙聯(lián)擺問題的求解隱格式非線性常微分方程初值

問題,初值為[,,0,0]T

odeset可以設(shè)置質(zhì)量矩陣質(zhì)量矩陣swinger_solve.m很大時(shí)混沌現(xiàn)象:確定的但不可預(yù)測(cè)(非線性科學(xué))

42第四十二頁(yè)，共44頁(yè)。小結(jié)函數(shù)名內(nèi)部算法說明ode23顯格式，單步法，采用BS23算法的自動(dòng)變步長(zhǎng)R-K方法對(duì)于精度要求不高的情況，效率好于ode45ode45顯格式，單步法，包含一個(gè)4階和5階公式的自動(dòng)變步長(zhǎng)R-K方法一般情況下，首先嘗試使用它來求解ode113顯格式，多步法，采用變階數(shù)的Adams-Bashforth-Moulton預(yù)測(cè)校正公式ode15s