常微分方程數(shù)值解法第1頁，共72頁，2023年，2月20日，星期一§1引言1.0基本概念1.常微分方程的初值問題：稱為具有初值(1.2)的常微分方程.①若f(x,y)在{axb,|y|<+}上連續(xù)，且關(guān)于y滿足Lip條件：常數(shù)L使|f(x,y1)–f(x,y2)|

L|y1–y2|則初值問題(1.1)(1.2)存在唯一連續(xù)可微解y(x).注：以下總假設(shè)f滿足Lip條件.第2頁，共72頁，2023年，2月20日，星期一§1引言1.0基本概念1.常微分方程的初值問題：稱為具有初值(1.2)的常微分方程.②(1.1)(1.2)等價于微分方程：

(1.3)注：一般無初等解(解析解)，即使有形式也復(fù)雜.第3頁，共72頁，2023年，2月20日，星期一§1引言1.0基本概念2.初值問題的數(shù)值解

設(shè)(1.1)(1.2)的解y(x)在節(jié)點xi處的近似解值為

yi

y(xi),a<x1<x2<…<xn=b則稱yi(i=1,2,…,n)為(1.1)(1.2)的數(shù)值解，又稱y(xi)的計算值.第4頁，共72頁，2023年，2月20日，星期一§1引言1.0基本概念3.數(shù)值方法①兩種轉(zhuǎn)化：由微分出發(fā)的數(shù)值方法.

由積分出發(fā)的數(shù)值方法.②計算方法步進法：從初始條件出發(fā)，逐步求y1,y2,…,yn.

又有兩種：單步法，多步法.注：采用等距節(jié)點：第5頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.

(1.6)第6頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.1.前進歐拉公式

(1.6)的前半部分為：令yi+1=yi

+

hf(xi,yi)

(1.7)其中yi=y(xi),則yi+1

y(xi+1)第7頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.1.前進歐拉公式

令yi+1=

yi

+

hf(xi,yi)

(1.7)其中yi=y(xi),則yi+1

y(xi+1)記(1.8)則稱(1.7)為前進歐拉求解公式.簡稱為歐拉公式或歐拉法.(1.8)稱為歐拉公式的余項：ei+1(h)=y(xi+1)–yi+1

第8頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.2.后退歐拉公式

(1.6)的后半部分令yi+1=yi

+hf(xi+1,yi+1)

(1.9)其中yi=y(xi),則yi+1

y(xi+1)

第9頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.2.后退歐拉公式令yi+1=yi

+hf(xi+1,yi+1)

(1.9)其中yi=y(xi),則yi+1

y(xi+1)

注：①(1.9)中f(xi+1,yi+1)

f(xi+1,y(xi+1))

∴余項

(1.10)第10頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.2.后退歐拉公式令yi+1=yi

+hf(xi+1,yi+1)

(1.9)其中yi=y(xi),則yi+1

y(xi+1)

注：②稱(1.9)為后退歐拉公式(后退歐拉法).

稱(1.10)為后退歐拉法的誤差近似值.③歐拉法與后退歐拉公式的區(qū)別：(1.7)為直接計算公式稱顯式公式.(1.9)為關(guān)于函數(shù)方程稱隱式公式.第11頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.【例1】取h=0.1求解初值問題：

(1.11).解：，xi=ih=0.1i,(i=0,1,2,…,10)

①歐拉法：第12頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.【例1】取h=0.1求解初值問題：

(1.11).解：，xi=ih=0.1i,(i=0,1,2,…,10)

②后退歐拉法：

第13頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.注：為避免求解函數(shù)方程，采用顯式與隱式結(jié)合的方法：

此方法稱為預(yù)測——校正系統(tǒng).求解過程為：第14頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.預(yù)測——校正系統(tǒng)：【例2】利用預(yù)測——校正系統(tǒng)求解例1.第15頁，共72頁，2023年，2月20日，星期一§1引言1.1基于數(shù)值微分的求解公式.預(yù)測——校正系統(tǒng)：注：顯式比隱式方便，但有時隱式效果比顯式好.(§4介紹).第16頁，共72頁，2023年，2月20日，星期一§1引言1.2截斷誤差定義1.1稱ek(h)=y(xk)–yk為計算yk的公式第k步的局部截斷誤差.注：①“局部”是指在計算第k步時，假定前面yi=y(xi)

(i<k).而yk

y(xk)

②

——歐拉法.——后退歐拉法.③一般根據(jù)y"(xk)對y"(k),y"(k)做估計.第17頁，共72頁，2023年，2月20日，星期一§1引言1.2截斷誤差定義1.2設(shè)ei(h)(i=1,2,…,k)為求解公式第i步的局部截斷誤差.稱為該求解公式在點上的整體截斷誤差.注：①局部截斷誤差ek(h)與yk有關(guān).

整體截斷誤差Ek(h)與y1,y2,…,yk有關(guān).②所有ek(h)都與h有關(guān).第18頁，共72頁，2023年，2月20日，星期一§1引言1.2截斷誤差定義1.3若局部截斷誤差e(h)=O(hp+1)，則稱該求解公式具有p階精度.注：歐拉法具有一階精度.(精度越高越好)第19頁，共72頁，2023年，2月20日，星期一§1引言作業(yè)P2081，2，3.第20頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式

(1.13)若已知y(xk)=yk,則計算積分可求出y(xk+1)

.

如用矩形公式求積分則有y(xk+1)=y(xk)+hf(xk,yk)令yk+1=y(xk)+hf(xk,yk)即為歐拉公式.故歐拉公式又稱矩形法.第21頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式

(1.13)考慮1.梯形公式記(1.14)第22頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式1.梯形公式記(1.14)稱(1.14)為梯形(求解)公式.簡稱梯形法.第23頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式1.梯形公式梯形(求解)公式,簡稱梯形法:

(1.14)注：①梯形公式的余項：

故是二階精度.第24頁，共72頁，2023年，2月20日，星期一1.3基于數(shù)值積分的求解公式1.梯形公式

(1.14)②梯形公式為隱式公式.<顯化>預(yù)測——校正系統(tǒng)

(1.15)稱(1.15)為改進的歐拉公式，也可記為§1引言第25頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式1.梯形公式

(1.14)③可以證明，改進歐拉公式也具有二階精度.第26頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式【例3】用歐拉法，梯形法以及改進歐拉法求解取h=0.1.計算到x=0.5.解：f(x,y)=x–y+1,a=x0=0,b=0.5,y0=1,n=5(Euler法)求解公式：yk=yk–1+h(xk–1–yk–1+1)=hxk–1+(1–h)yk–1+h

=0.1xk–1+0.9yk–1+0.1

第27頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式【例3】用歐拉法，梯形法以及改進歐拉法求解解：f(x,y)=x–y+1,a=x0=0,b=0.5,y0=1,n=5(梯形法)求解公式：yk=yk–1+h[(xk–1–yk–1+1)+(xk–yk+1)]/2解出yk，得方程第28頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式【例3】用歐拉法，梯形法以及改進歐拉法求解解：f(x,y)=x–y+1,a=x0=0,b=0.5,y0=1,n=5(改進Euler法)求解公式：yk=yk–1+h[(xk–1–yk–1+1)+xk–(yk+h(xk–yk+1))+1]/2得=0.905yk–1+0.045xk–1+0.05xk+0.095方程第29頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式2.辛卜生公式

記

(1.17)第30頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式2.辛卜生公式記

(1.17)其余項第31頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式2.辛卜生公式記

(1.17)將[xk–1,xk]

對分：調(diào)整下標(biāo)為[xi–2,xi]：xi–2=xk–1,xi–1=xk–1+h1,xi=xk–1+2h1=xk則(1.17)化為

(1.19)稱(1.19)為辛卜生求解公式，其中fk–2=f(xk–2,y(xk–2))，fk–1=f(xk–1,y(xk–1))，fk=f(xk,y(xk))第32頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式2.辛卜生公式記

(1.17)

(1.19)稱(1.19)為辛卜生求解公式，其中fi–2=f(xi–2,y(xi–2))，fi–1=f(xi–1,y(xi–1))，fi=f(xi,y(xi))注：①(1.19)的誤差：第33頁，共72頁，2023年，2月20日，星期一§1引言1.3基于數(shù)值積分的求解公式2.辛卜生公式記

(1.17)

(1.19)稱(1.19)為辛卜生求解公式，其中fi–2=f(xi–2,y(xi–2))，fi–1=f(xi–1,y(xi–1))，fi=f(xi,y(xi))注：②隱式(需顯化)多步——將在§3中討論.第34頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.0原理

其中K=f(,y())=y'()稱為y在[xi–1,xi]上的平均斜率.歐拉法：改進歐拉法：(2.1)第35頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.0原理

其中K=f(,y())=y'()稱為y在[xi–1,xi]上的平均斜率.對(1.17)顯化：辛卜生：

(2.4)第36頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.0原理其中K=f(,y())=y'()稱為y在[xi–1,xi]上的平均斜率.設(shè)想：在中多計算(預(yù)測)幾個點上的值然后可加權(quán)取平均值作為的近似值可能構(gòu)成更高階的公式.一階二階三階第37頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.1Runge-Kutta公式

(*)其中0

j

1，yi–1+jh是y(xi–1+jh)的預(yù)測值.

稱(*)為R-K公式<R-K法>注：①(2.1)(2.4)分別稱為二階，三階R-K公式.②j，j，j為待定系數(shù).使(*)的階數(shù)盡量高.第38頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.1Runge-Kutta公式參數(shù)的確定，以m=2為例.

欲求1，2，2

.原則:使ei(h)=y(xi)–yi的階數(shù)盡可能高第39頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.1Runge-Kutta公式展開展開

原則:使ei(h)=y(xi)–yi的階數(shù)盡可能高第40頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.1Runge-Kutta公式

原則:使ei(h)=y(xi)–yi的階數(shù)盡可能高第41頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.1Runge-Kutta公式

欲求截斷誤差ei(h)=y(xi)–yi關(guān)于h的階數(shù)盡可能高，應(yīng)使無窮多解，從而有許多2階R-K公式第42頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.1Runge-Kutta公式應(yīng)使注：①取1=2=1/2，2=1，即為改進歐拉公式.第43頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.1Runge-Kutta公式應(yīng)使注：②取1=0，2=1，2=1/2，即為中點公式第44頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.1Runge-Kutta公式應(yīng)使注：③二階R-K公式的截斷誤差為故為二階方法.相仿可得更高階的R-K公式.第45頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.2經(jīng)典R-K公式在4解R-K公式中最重要的是經(jīng)典R-K公式<標(biāo)準(zhǔn)R-K公式>.(2.6)注：①(2.6)為4階方法.第46頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.2經(jīng)典R-K公式在4解R-K公式中最重要的是經(jīng)典R-K公式<標(biāo)準(zhǔn)R-K公式>.(2.6)注：②R-K法對4階以上不一定能提高整數(shù)階.第47頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法2.2經(jīng)典R-K公式【例4】使用三階，四階R-K法求解初值問題：

的部分計算值y1，y2，y3，其中h=0.1.解①使用三階R-K法第48頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法【例4】使用三階，四階R-K法求解初值問題：

的部分計算值y1，y2，y3，其中h=0.1.解①使用三階R-K法第49頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法【例4】使用三階，四階R-K法求解初值問題：

的部分計算值y1，y2，y3，其中h=0.1.解②使用四階R-K法第50頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法【例4】使用三階，四階R-K法求解初值問題：

的部分計算值y1，y2，y3，其中h=0.1.解②使用四階R-K法第51頁，共72頁，2023年，2月20日，星期一§2Runge-Kutta法注使用R-K法要求具備較好的光滑性，否則效果不如低階的.<P189>作業(yè)P2098②③9，10.第52頁，共72頁，2023年，2月20日，星期一§3線性多步法單步法的優(yōu)點：簡單，計算yk+1只用yk.缺點:沒有充分利用前面的信息且計算y(xk+h)較困難回顧Simpson：(1.19)考慮：(3.1)兩種插值求積：①將[xk–1,xk]增加內(nèi)部節(jié)點，改為[xk–2,xk]導(dǎo)出的公式稱為閉型求解公式.線性多步第53頁，共72頁，2023年，2月20日，星期一§3線性多步法考慮：(3.1)兩種插值求積：①將[xk–1,xk]增加內(nèi)部節(jié)點，改為[xk–2,xk]導(dǎo)出的公式稱為閉型求解公式.②在[xk–1,xk]外增加插值節(jié)點，導(dǎo)出的公式稱為開型求解公式.開型有顯和隱，閉型也有顯和隱.第54頁，共72頁，2023年，2月20日，星期一§3線性多步法3.1開型求解公式1.亞當(dāng)斯顯式求解公式取節(jié)點xk–3,xk–2,xk–1，在[xk–3,xk]上作F(x)=f(x,y(x))的插值多項式.第55頁，共72頁，2023年，2月20日，星期一§3線性多步法3.1開型求解公式1.亞當(dāng)斯顯式求解公式取節(jié)點xk–3,xk–2,xk–1，在[xk–3,xk]上記xk–i=xk–ih,x=xk+th，則

第56頁，共72頁，2023年，2月20日，星期一§3線性多步法3.1開型求解公式1.亞當(dāng)斯顯式求解公式取節(jié)點xk–3,xk–2,xk–1，記xk–i=xk–ih,x=xk+th，則

代入(3.1)得第57頁，共72頁，2023年，2月20日，星期一§3線性多步法3.1開型求解公式1.亞當(dāng)斯顯式求解公式取節(jié)點xk–3,xk–2,xk–1，記xk–i=xk–ih,x=xk+th，則

令

(3.4)稱(3.4)為亞當(dāng)斯顯式求解公式(線性多步).第58頁，共72頁，2023年，2月20日，星期一§3線性多步法3.1開型求解公式1.亞當(dāng)斯顯式求解公式取節(jié)點xk–3,xk–2,xk–1，記xk–i=xk–ih,x=xk+th，則余項：∵∴從而(3.4)具有3階精度.稱為3階亞當(dāng)斯求解公式.第59頁，共72頁，2023年，2月20日，星期一§3線性多步法3.1開型求解公式1.亞當(dāng)斯顯式求解公式類似地取xk–4,xk–3,xk–2,xk–1

在[xk–4,xk]上作F(x)=f(x,y(x))的插值多項式，可導(dǎo)出4階亞當(dāng)斯顯式求解公式：

(3.6)

(3.7)——4階精度第60頁，共72頁，2023年，2月20日，星期一§3線性多步法3.1開型求解公式2.亞當(dāng)斯隱式求解公式取xk–3,xk–2,xk–1,xk，在[xk–3,xk]上作F(x)=f(x,y(x))的插值多項式用上述方法可導(dǎo)出：

(3.8)

(3.9)稱為亞當(dāng)斯隱式求解公式.第61頁，共72頁，2023年，2月20日，星期一§3線性多步法3.1開型求解公式2.亞當(dāng)斯隱式求解公式

(3.8)

(3.9)稱為亞當(dāng)斯隱式求解公式.注：利用4階公式(3.6)顯化之：

(3.10)稱(3.10)為亞當(dāng)斯預(yù)測——校正系統(tǒng).第62頁，共72頁，2023年，2月20日，星期一§3線性多步法3.2閉型求解系統(tǒng)將[xk–1,xk]擴充為[xk–4,xk]，取xk–4，xk–3，xk–2，xk–1為節(jié)點，作F(x)=f(x,y(x))的牛頓前插多項式.

則第63頁，共72頁，2023年，2月20日，星期一§3線性多步法3.2閉型求解系統(tǒng)將[xk–1,xk]擴充為[xk–4,xk]，取xk–4，xk–3，xk–2，xk–1為節(jié)點，作F(x)=f(x,y(x))的牛頓前插多項式.則令x=xk+(t–4)h

則

第64頁，共72頁，2023年，2月20日，星期一§3線性多步法3.2閉型求解系統(tǒng)令x=xk+(t–4)h

則由第65頁，共72頁，2023年，2月20日，星期一