數值分析(8)NumericalAnalysisWenjianYu2第八章常微分方程初值問題

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

WenjianYu4

常微分方程

WenjianYu5

初值問題:

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

度較大時擺錘的軌跡呈現混沌現象WenjianYu7求解微分方程初值問題,得到

擺錘的運動規律Matlab演示swinger

常微分方程

WenjianYu8

線性齊次常系數微分方程

實際的問題基本上都是穩定的!

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

WenjianYu9

局部穩定簡單方法與有關概念WenjianYu10簡單的初值問題數值解法

初值問題的數值解法

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

步長h=0.1,和0.05步長小的更準

數值解法的穩定性

WenjianYu13

-10

歐拉法解模型問題的穩定區域

數值解法的穩定性

WenjianYu14-10

歐拉法穩定

數值解法的穩定性

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這里設的h太大!計算結果如下表:數值解法的局部截斷誤差

WenjianYu16整體誤差

穩定的問題，整體誤差小于局部誤差之和不穩定的問題呢？一般僅能控制局部誤差整體誤差？~局部誤差數值解法的局部截斷誤差

WenjianYu17

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

梯形

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

WenjianYu19

準確解

01穩定區域

無條件穩定(unconditionallystable)!向后歐拉法

WenjianYu20

具有1階準確度!向后歐拉法與梯形法

WenjianYu21

穩定的條件是:

思考無條件穩定!

具有2階準確度簡單方法與有關概念WenjianYu22Runge-Kutta方法在歐拉法基礎上改進再增加一次函數求值:數值積分的中矩形或梯形公式利用歐拉法算半個步長的結果,估算中點處被積函數值先用歐拉法估計區間終點處斜率,再用它與

起始點斜率的平均值算一整步Runge-Kutta方法

中矩形梯形公式(中點公式)(Heun方法/改進的歐拉法)都比歐拉法準確均為2級R-K公式WenjianYu23

Runge-Kutta方法只能估算

,令

WenjianYu24

Runge-Kutta方法其他,用所有前面點的信息

WenjianYu25積分節點………幾種顯式R-K公式參數的值不按具體數值

積分公式設置,而根據

準確度階數要求設置2級公式：改進歐拉法、

中點公式經典4級、4階Runge-Kutta法

(1905)Runge-Kutta方法WenjianYu26

Runge-Kutta方法此時局部截斷誤差只要對非模型問題也有相同結論！WenjianYu27如

例8.7:用2階改進歐拉、3階Ralston、

4階經典R-K解問題,h=0.1,算到y(2)精確解為Runge-Kutta方法r4對應的r級R-K公式有r階準確度高于4階的公式很少單獨使用WenjianYu28r>4對應的r級R-K公式達不到r階準確度二階Heun三階Ralston四階Runge-Kutta準確值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方法(此時)一般顯式單步法

恰好與歐拉法一樣顯格式,都不是無條件穩定的局部截斷誤差判斷單步法收斂性的簡便方法WenjianYu29

簡單方法與有關概念WenjianYu30多步法

多步法WenjianYu31

(線性m步法)

固定步長hTaylor展開法求線性m步法的系數

這些系數應該等于0例8.8:求兩步法公式中參數值多步法WenjianYu32

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

,

有二階準確度

多步法WenjianYu33

同例8.8的結果!

多步法公式中包含p個待定參數，至少可達到p-1階準確度

多步法WenjianYu34

Vandermonde陣T,非奇異

插值節點函數值

常用的多步法公式Adams公式的推導：用插值多項式近似被積函數例8.10:推導m=4對應的顯式Adams公式多步法WenjianYu35

可證明滿足最高準確度階數

類似地算其他系數，得

顯式四階Adams-

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

階數穩定閾值誤差常數11

-21/223/2-1/2

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

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

階數穩定閾值誤差常數11

--1/221/21/2

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

-6-1/2449/2419/24-5/241/24-3-19/720向后歐拉法梯形法并非無條件穩定!

多步法WenjianYu37

用Matlab求解初值問題WenjianYu38用Matlab解ODE-IVP

Matlab中的ODE-IVP求解器

WenjianYu39[T,

Y,

TE,

YE,

IE]=solver(odefun,

tspan,

y0,

options)

求解單個ODE火焰燃燒問題當點燃一根火柴時,火焰迅速增大直到一個臨界體積,然后維持這一體積不變,此時火焰內部燃燒耗費的氧氣和其表面現存的氧氣達到了一種平衡.火焰(近似為球)半徑y滿足ODE>>f=@(t,y)y^2-y^3;>>ode23(f,[0,2.0e4],1e-4)WenjianYu40設初始半徑=0.0001例8.15嘗試用ode23s求解它求解ODE方程組

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

雙聯擺問題的求解WenjianYu42隱格式非線性常微分方程初值

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

odeset可以設置質量矩陣質量矩陣swinger_solve.m很大時混沌現象:確定的但不可預測(非線性科學)

小結WenjianYu43函數名內部算法說明ode23顯格式，單步法，采用BS23算法的自動變步長R-K方法對于精度要求不高的情況，效率好于ode45ode45顯格式，單步法，包含一個4階和5階公式的自動變步長R-K方法一般情況下，首先嘗試使用它來求解ode113顯格式，多步法，采用變階數的Adam