《數(shù)值分析與算法》第八講-常微分方程初值問(wèn)題

數(shù)值分析(8)NumericalAnalysisWenjianYu2第八章常微分方程初值問(wèn)題

(主要是前3節(jié))常微分方程初值問(wèn)題WenjianYu3常微分方程基本概念常微分方程

WenjianYu4

常微分方程

WenjianYu5

初值問(wèn)題:

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

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

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

常微分方程

WenjianYu8

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

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

(由于歷史原因)常微分方程

WenjianYu9

局部穩(wěn)定簡(jiǎn)單方法與有關(guān)概念WenjianYu10簡(jiǎn)單的初值問(wèn)題數(shù)值解法

初值問(wèn)題的數(shù)值解法

WenjianYu11

否則為多步法否則為隱格式方法歐拉法

WenjianYu12“左矩形”求積公式

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)

數(shù)值解法的穩(wěn)定性

WenjianYu13

-10

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

數(shù)值解法的穩(wěn)定性

WenjianYu14-10

歐拉法穩(wěn)定

數(shù)值解法的穩(wěn)定性

WenjianYu15

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é)果如下表:數(shù)值解法的局部截?cái)嗾`差

WenjianYu16整體誤差

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

WenjianYu17

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

梯形

00.0250.050.0750.10.1250.1510.0066630.0019040.00054410.0820850.0067380.0005534.5410-53.7310-63.0610-7向后歐拉法

WenjianYu19

準(zhǔn)確解

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

無(wú)條件穩(wěn)定(unconditionallystable)!向后歐拉法

WenjianYu20

具有1階準(zhǔn)確度!向后歐拉法與梯形法

WenjianYu21

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

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

具有2階準(zhǔn)確度簡(jiǎn)單方法與有關(guān)概念WenjianYu22Runge-Kutta方法在歐拉法基礎(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公式WenjianYu23

Runge-Kutta方法只能估算

,令

WenjianYu24

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

WenjianYu25積分節(jié)點(diǎn)………幾種顯式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方法WenjianYu26

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

例8.7:用2階改進(jìn)歐拉、3階Ralston、

4階經(jīng)典R-K解問(wèn)題,h=0.1,算到y(tǒng)(2)精確解為Runge-Kutta方法r4對(duì)應(yīng)的r級(jí)R-K公式有r階準(zhǔn)確度高于4階的公式很少單獨(dú)使用WenjianYu28r>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.3000000

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

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

簡(jiǎn)單方法與有關(guān)概念WenjianYu30多步法

多步法WenjianYu31

(線性m步法)

固定步長(zhǎng)hTaylor展開(kāi)法求線性m步法的系數(shù)

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

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

,

有二階準(zhǔn)確度

多步法WenjianYu33

同例8.8的結(jié)果!

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

多步法WenjianYu34

Vandermonde陣T,非奇異

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

常用的多步法公式Adams公式的推導(dǎo)：用插值多項(xiàng)式近似被積函數(shù)例8.10:推導(dǎo)m=4對(duì)應(yīng)的顯式Adams公式多步法WenjianYu35

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

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

顯式四階Adams-

Bashforth公式(單項(xiàng)式函數(shù)代入法)Adams公式幾種顯式公式幾種隱式公式多步法WenjianYu36

階數(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ú)條件穩(wěn)定!

多步法WenjianYu37

用Matlab求解初值問(wèn)題WenjianYu38用Matlab解ODE-IVP

Matlab中的ODE-IVP求解器

WenjianYu39[T,

Y,

TE,

YE,

IE]=solver(odefun,

tspan,

y0,

options)

求解單個(gè)ODE火焰燃燒問(wèn)題當(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)WenjianYu40設(shè)初始半徑=0.0001例8.15嘗試用ode23s求解它求解ODE方程組

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

雙聯(lián)擺問(wèn)題的求解WenjianYu42隱格式非線性常微分方程初值

問(wèn)題,初值為[,,0,0]T

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

小結(jié)WenjianYu43函數(shù)名內(nèi)部算法說(shuō)明ode23顯格式，單步法，采用BS23算法的自動(dòng)變步長(zhǎng)R-K方法對(duì)于精度要求不高的情況，效率好于ode45ode45顯格式，單步法，包含一個(gè)4階和5階公式的自動(dòng)變步長(zhǎng)R-K方法一般情況下，首先嘗試使用它來(lái)求解ode113顯格式，多步法，采用變階數(shù)的Adam