第5講：常微分方程數(shù)值解

常微分方程數(shù)值解《數(shù)值分析》第五講第五章：常微分方程數(shù)值解§5.1引言1、常微分方程與解為n階常微分方程。如果函數(shù)在區(qū)間[a,b]內(nèi)n階可導，稱方程滿足方程的函數(shù)稱為微分方程的解。則如為任意常數(shù)）為方程的解,一般稱為方程的通解。如果則有為方程滿足定解條件的解。第五章：常微分方程數(shù)值解方程的通解滿足定解條件的解微分關系（方程）解的圖示第五章：常微分方程數(shù)值解本教材重點討論定解問題(初值問題）定解條件（初始條件）是否能夠找到定解問題的解取決于僅有極少數(shù)的方程可以通過“常數(shù)變易法”、“可分離變量法”等特殊方法求得初等函數(shù)形式的解，絕大部分方程至今無法理論求解。如等等2、數(shù)值解的思想第五章：常微分方程數(shù)值解（1）將連續(xù)變量離散為（2）用代數(shù)的方法求出解函數(shù)在點的近似值*數(shù)學界關注工程師關注如果找不到解函數(shù)數(shù)學界還關注：解的存在性解的唯一性解的光滑性解的振動性解的周期性解的穩(wěn)定性解的混沌性解的分岔性……§5.2Euler方法第五章：常微分方程數(shù)值解第一步：連續(xù)變量離散化第二步：用直線步進·····1、Euler格式Euler格式例P106第五章：常微分方程數(shù)值解初值問題Bernoulli型方程常數(shù)變易法第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解令將代入Euler格式計算得x精確值Euler方法Euler方法誤差01.00000001.00000000.00000000.11.09544511.10000000.00455490.21.18321601.19181820.00860220.31.26491111.27743780.01252680.41.34164081.35821260.01657180.51.41421361.43513290.02091940.61.48323971.50896630.02572660.71.54919331.58033820.03114490.81.61245151.64978340.03733190.91.67332011.71777930.044459311.73205081.78477080.0527200第五章：常微分方程數(shù)值解Euler值2、Euler格式的誤差分析第五章：常微分方程數(shù)值解事實上Euler格式的每一步都存在誤差，為了方便討論算法的好壞，假定第n步準確的前提下分析第n+1步的誤差，稱為局部截斷誤差。定義1定義2第五章：常微分方程數(shù)值解Euler格式的誤差即Euler格式具有一階精度將在點Taylor展開的計算格式第五章：常微分方程數(shù)值解3、Euler格式與后退Euler格式Euler格式的值后退Euler格式的值隱式格式第五章：常微分方程數(shù)值解4、后退格式的精度局部誤差分析的要求第五章：常微分方程數(shù)值解對照即后退Euler格式具有1階精度第五章：常微分方程數(shù)值解3、梯形格式得梯形格式Euler格式的值后退Euler格式的值第五章：常微分方程數(shù)值解梯形格式幾何解釋第五章：常微分方程數(shù)值解4、梯形格式的精度局部誤差分析的要求第五章：常微分方程數(shù)值解梯形格式5、改進的Euler格式預測校正第五章：常微分方程數(shù)值解隱式格式為方便計算，一般用以下改進格式計算用改進格式計算例5.1的結(jié)果為x精確值Euler方法改進Euler方法Euler方法誤差改進Euler誤差01.00000001.00000001.00000000.00000000.00000000.11.09544511.10000001.09590910.00455490.00046400.21.18321601.19181821.18409660.00860220.00088060.31.26491111.27743781.26620140.01252680.00129030.41.34164081.35821261.34336020.01657180.00171940.51.41421361.43513291.41640190.02091940.00218840.61.48323971.50896631.48595560.02572660.00271590.71.54919331.58033821.55251410.03114490.00332080.81.61245151.64978341.61647480.03733190.00402320.91.67332011.71777931.67816640.04445930.004846311.73205081.78477081.73786740.05272000.0058166第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解6、預測精度的改進—兩步Euler格式第五章：常微分方程數(shù)值解如果令第五章：常微分方程數(shù)值解整理得兩步Euler格式具有2階精度。因此可得具有2階精度的Euler改進格式預測校正第五章：常微分方程數(shù)值解問題：精度還能提高嗎？§5.3Lunge-Kutta方法1、二階Lunge-Kutta方法（P113-P115)第五章：常微分方程數(shù)值解依據(jù)精度要求的待定系數(shù)法令加權(quán)平均斜率

在點的斜率

在點的斜率用Euler格式預測第五章：常微分方程數(shù)值解對照第五章：常微分方程數(shù)值解可解得得改進的Euler格式3、三階Lunge-Kutta方法第五章：常微分方程數(shù)值解補充第五章：常微分方程數(shù)值解分別將第五章：常微分方程數(shù)值解由得一個三階精度的Runge-Kutta格式4、四階Lunge-Kutta方法第五章：常微分方程數(shù)值解補充點第五章：常微分方程數(shù)值解同理可得一個經(jīng)典四階精度的Runge-Kutta格式§5.4幾種方法的數(shù)值計算例5.1P106改進的Euler格式Euler格式第五章：常微分方程數(shù)值解四階經(jīng)典Lunge-Kutta方法例5.1P106第五章：常微分方程數(shù)值解x精確值Euler方法改進Euler方法四階Lunge-KuttaEuler方法誤差改進Euler誤差四階Lunge-Kutta誤差01.00000001.00000001.00000001.00000000.00000000.00000000.00000000.11.09544511.10000001.09590911.09544550.00455490.00046400.00000040.21.18321601.19181821.18409661.18321670.00860220.00088060.00000080.31.26491111.27743781.26620141.26491220.01252680.00129030.00000120.41.34164081.35821261.34336021.34164240.01657180.00171940.00000160.51.41421361.43513291.41640191.41421560.02091940.00218840.00000200.61.48323971.50896631.48595561.48324220.02572660.00271590.00000250.71.54919331.58033821.55251411.54919650.03114490.00332080.00000310.81.61245151.64978341.61647481.61245530.03733190.00402320.00000380.91.67332011.71777931.67816641.67332470.04445930.00484630.000004611.73205081.78477081.73786741.73205640.05272000.00581660.0000056幾種方法的結(jié)果與誤差第五章：常微分方程數(shù)值解參考程序-Eulerx=0:0.01:1;y=sqrt(1+2.*x);a=0.0;b=1.0;n=10;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;fork=1:10y0(k+1)=y0(k)+h*(y0(k)-2*x0(k)/y0(k));endfori=1:10y1(1)=1.0;y1(i+1)=y1(i)+h*(y1(i)-2*x0(i)/y1(i));y1(i+1)=y1(i)+h*(y1(i)-2*x0(i)/y1(i)+y1(i+1)-2*x0(i+1)/y1(i+1))/2;endplot(x,y,'b');holdon;plot(x0,y0,'or');holdon;plot(x0,y1,'*');第五章：常微分方程數(shù)值解參考程序-Lunge_Kuttax=0:0.01:1;y=sqrt(1+2.*x);a=0.0;b=1.0;n=10;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;fork=1:10k1=y0(k)-2*x0(k)/y0(k);k2=y0(k)+h*k1/2-(2*x0(k)+h)/(y0(k)+h*k1/2);k3=y0(k)+h*k2/2-(2*x0(k)+h)/(y0(k)+h*k2/2);k4=y0(k)+h*k3-2*(x0(k)+h)/(y0(k)+h*k3);y0(k+1)=y0(k)+h*(k1+2*k2+2*k3+k4)/6;endholdon;plot(x,y,'b');holdon;plot(x0,y0,'or');第五章：常微分方程數(shù)值解§5.5線性多步方法1、Adams顯式格式第五章：常微分方程數(shù)值解Euler公式梯形公式····參看P27(2.5.12)第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解r=3時第五章：常微分方程數(shù)值解當r=1時，（5.5.6）式的誤差分析第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解例5.1P106第五章：常微分方程數(shù)值解Adams程序x=0:0.01:1;y=sqrt(1+2.*x);a=0.0;b=1.0;n=10;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;y0(2)=1.0954;y0(3)=1.1832;y0(4)=1.2649;fork=4:10y0(k+1)=y0(k)+h*(55*(y0(k)-2*x0(k)/y0(k))-59*(y0(k-1)-2*x0(k-1)/y0(k-1))+37*(y0(k-2)-2*x0(k-2)/y0(k-2))-9*(y0(k-3)-2*x0(k-3)/y0(k-3)))/24;endholdon;plot(x,y,'b');holdon;plot(x0,y0,'or');第五章：常微分方程數(shù)值解x精確值四階Lunge-KuttaAdams四階Lunge-Kutta誤差Adams誤差01.00000001.00000001.00000000.00000000.00000000.11.09544511.09544551.09540000.00000040.00004510.21.18321601.18321671.18320000.00000080.00001600.31.26491111.26491221.26490000.00000120.00001110.41.34164081.34164241.34153280.00000160.00010800.51.41421361.41421561.41402400.00000200.0001896