常微分方程數(shù)值解法

1、實(shí)用標(biāo)準(zhǔn)文案第八章 常微分方程的數(shù)值解法一內(nèi)容要點(diǎn)dyf ( x, y)考慮一階常微分方程初值問題 : dxy(x0 ) y0微分方程的數(shù)值解：設(shè)微分方程的解y( x) 的存在區(qū)間是 a,b ，在 a,b 內(nèi)取一系列節(jié)點(diǎn) a= x 0 x 1 x n =b , 其中 hk=xk+1-x k；( 一般采用等距節(jié)點(diǎn)， h=(b-a)/n 稱為步長) 。在每個節(jié)點(diǎn) xk 求解函數(shù) y(x) 的近似值 : yk y(x k) ，這樣 y0 , y 1 ,.,yn 稱為微分方程的數(shù)值解。用數(shù)值方法，求得 f(xk) 的近似值 yk，再用插值或擬合方法就求得y(x)的近似函數(shù)。( 一 ) 常微分方程處置問

2、題解得存在唯一性定理dy對于常微分方程初值問題 :dx f (x, y) 如果：y( x0 ) y0(1)在 x0x A yy0B 的矩形內(nèi)f ( x, y)是一個二元連續(xù)函數(shù)。,(2)f (x, y) 對于 y 滿足利普希茨條件，即dyf ( x, y1 )f (x, y2 )L y1y2 則在 x0x C 上方程 dx f ( x, y) 的解存在且唯一，這y(x0 )y0里 C=min(A- x 0), x0 +B/L),L 是利普希茨常數(shù)。定義：任何一個一步方法可以寫為yk 1ykh(xk , yk , h) ，其中( xk , yk , h) 稱為算法的增量函數(shù)。收斂性定理：若一步方

3、法滿足：(1) 是 p 解的 .(2)增量函數(shù) (xk , yk , h) 對于 y滿足利普希茨條件 .(3)初始值 y0 是精確的。則()()(hp ),kh=x- x0, 也就是有y khy xOlimyky( x)0h 0kh x x0( 一 ) 、主要算法1局部截斷誤差局部截斷誤差： 當(dāng) y( xk) 是精確解時， 由 y( xk ) 按照數(shù)值方法計算出來的y k 1 的誤差y( xk+1)- y k 1 稱為局部截斷誤差。ek+1=y( xk+1) - yk +1不同。注意： yk+1 和 yk 1 的區(qū)別。因而局部截斷誤差與誤差p+1如果局部截斷誤差是O( h), 我們就說該數(shù)值方

4、法具有p 階精度。精彩文檔實(shí)用標(biāo)準(zhǔn)文案yk 1ykhf ( xk , yk )1. 顯式歐拉格式：k=1,2, n-1 .y( x0 )y0顯式歐拉格式的特點(diǎn)：(1) 、單步方法；(2) 、顯式格式；(3) 、局部截斷誤差 O(h2 ) 因而是一階精度 。yk 1ykhf (xk 1 , yk 1 )隱式歐拉格式y(tǒng)( x0 ) y0yk 1yk 12hf (xk , yk )n-1 .2. 兩步歐拉格式：y( x0 )k=1,2,y0兩步歐拉格式的局部截斷誤差O (h3 ) ，因而是二階精度3. 梯形方法 :yk 1yk 1ykhf ( xk 1 , yk 1 ) f ( xk , yk )

5、；2y( x0 ) y03. 改進(jìn)的歐拉方法：預(yù)測值 :h f (xk , yk )yk 1yk校正值 :hyk 1ykf (xk 1, yk 1 ) f ( xk , yk ) 。2ypykhf ( xk , yk )或改寫為 :ycyk hf (xk 1, yp )1yk 1( ypyc )4、梯形方法與改進(jìn)歐拉方法的截斷誤差是O(h3 ), 具有二階精度。5、龍格 - 庫塔法的思想K 1f (xk , yk )1).二階龍格 - 庫塔法計算公式：K 2f ( xkph, ykphK 1 )yk 1ykh(1 )K 1K 2 )當(dāng):1時，得一簇龍格 - 庫塔公式，其局部截斷誤差均為3都是二

6、階精度。特別pO( h )2取1 , p 1, 就是改進(jìn)歐拉公式。2取 p1 ,1, 得二階龍格庫塔法為：2精彩文檔實(shí)用標(biāo)準(zhǔn)文案K 1f ( xk , yk )K 2f ( xkh , ykh K 1 ) ，稱為二階中點(diǎn)格式。22yk 1ykhK 22) 、經(jīng)典龍格庫塔格式 ( 也稱為標(biāo)準(zhǔn)龍格庫塔方法 ) ：K 1f ( xk , yk )K 2f ( xkh , ykh K 1 )22K 3f ( xkh , ykh K 2 )22K 4f ( xkh, ykhK 3 )yk 1ykh (K 12K2 2K3 K4)6四階龍格庫塔方法的截斷誤差為h5 ，具有四階精度。一般一階常微分方程初值問

7、題用四階龍格庫塔方法計算，其精度均滿足實(shí)際問題精度要求。3). 變步長龍格庫塔方法：從節(jié)點(diǎn) xk 出發(fā)，以步長 h 據(jù)四階龍格庫塔方法求出一個近似值 yk(h )1 , 然后以步長 h/2求出一個近似值 yk(h1/2)，得誤差事后估計式：yk 1 yk(h1/ 2) 1( yk(h 1/2)yk(h 1) )15根據(jù)yk(h1/2)yk(h )1 來選取步長 h。4). RKF格式：變步長龍格庫塔方法， 因頻繁加倍或折半步長會浪費(fèi)計算量。 Felhberg 改進(jìn)了傳統(tǒng)龍格庫塔方法，得到 RKF格式，較好解決了步長的確定，而且提高了精度與穩(wěn)定性，為 Matlab 等許多數(shù)值計算軟件采用。4/5

8、 階 RKF格式：由 4 階龍格庫塔方法與5 階龍格庫塔方法結(jié)合而成。yn 1ynh25K11408 K 32197 K41 K 5216256541015166656 K 328561 K 49 K 52 K 6yn 1ynhK 113512825564305055K1 f ( xn , yn )K 2f ( xnhh, ynK 1 )44K 3f ( xn3 h, yn3 K 19 K2)83232K 4f ( xn12h1932K 172007296, ynK2K 3 )13219721972197K 5f ( xnh, yn439 K18K 23680 K 3845 K4)216513

9、4104K 6f ( xnh , yn8 K1 2K23544 K 31859 K 411 K4)2272565410440精彩文檔實(shí)用標(biāo)準(zhǔn)文案1 / 4Felhberg 得到的最佳步長 hs，其中 h 為當(dāng)前步長， sh. 為精度要求，yn 1yn 1若 s1.5 步長加倍。6. 亞當(dāng)斯方法 (Adams)1). 顯式 Adams方法：記： fkf ( xk , yk ) ；二階顯式 Adams方法： yk 1yk三階顯式 Adams方法： yk 1yk四階顯式 Adams方法： yk 1yk2). 隱式 Adams方法二階隱式 Adams方法： yk 1yk三階隱式 Adams方法： yk

10、 1yk四階隱式 Adams方法： yk 1ykh 3 f kf k 1 ;2h 23 f k16 f k 1 5 fk 2 ;12h 9 fk 3 37 f k 2 59 f k 1 55 fk .24h 3 f kf k 1 ;2h 5 f k 18 f kfk 1 ;12h ( f k 25 f k 119 f k 9 f k 1 )243） .Adams 預(yù)報 - 校正系統(tǒng)：先用顯式格式作為預(yù)測值，再用隱式格式來校正。預(yù)測值： yk 1ykh(55 f k59 fk 137 f k29 f k 3 )24校正值：yk1ykh ( fk2 5 fk 119 f k9 f xk 1, yk

11、 1)244). 改進(jìn)的 Adams預(yù)報 - 校正系統(tǒng)：預(yù)測： pk 1ykh(55 f k59 f k 1 37 f k 29 f k 3 )24改進(jìn)： mk 1pk1251 (ckpk )270校正值：ck1ykh(9 f xk 1 ,mk 119 f k5 fk 1 f k2 )24改進(jìn)： yk 1ck119 (ck 1pk 1 )2707. 收斂與穩(wěn)定性對于固定的 xk x0 kh ，如果數(shù)值解 yk 當(dāng) h0（,同時 n）時趨向于準(zhǔn)確解 y( xk ) ，則稱該數(shù)值方法方法是收斂的。m k上產(chǎn)生的偏如果一種數(shù)值方法， 在節(jié)點(diǎn)值yk 上大小為 的擾動，于以后各節(jié)點(diǎn)值 ym( )差均不超

12、過，則稱數(shù)值方法是穩(wěn)定的.精彩文檔實(shí)用標(biāo)準(zhǔn)文案一般，隱式格式較顯式格式有較好的穩(wěn)定性。8. 剛性方程組：考慮 n 階常微分方程組： yAy ( x) ，（* ）max Re( i )若矩陣 A 的特征值 1, 2,n 的實(shí)部 Re( i )ey1=y;y=y0+h*feval(dyfun,x,y);k=k+1;if kK,error(迭代發(fā)散 );endend改進(jìn) Euler 格式functionx,y=naeulerg(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1K1=feval(dyfun,x(n),y

13、(n);y(n+1)=y(n)+h*K1;K2=feval(dyfun,x(n+1),y(n+1);y(n+1)=y(n)+h*(K1+K2)/2;endx=x;y=y;4 階經(jīng)典 Runge-Kutta格式functionx,y=nak4(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1K1=feval(dyfun,x(n),y(n);K2=feval(dyfun,x(n)+h/2,y(n)+h/2*K1);K3=feval(dyfun,x(n)+h/2,y(n)+h/2*K2);K4=feval(dyfun

14、,x(n+1),y(n)+h*K3);y(n+1)=y(n)+h*(K1+2*K2+2*K3+K4)/6;endx=x;y=y;變步長 4 階經(jīng)典 Runge-Kutta格式functionx,y=nark4v(dyfun,xspan,y0,e,h)if nargin5,h=(xspan(2)-xspan(1)/10;endn=1;x(n)=xspan(1);y(n)=y0;y1,y2=comput(dyfun,x(n),y(n),h);while x(n)ewhile abs(y2-y1)/10eh=h/2;y1,y2=comput(dyfun,x(n),y(n),h);精彩文檔實(shí)用標(biāo)準(zhǔn)文案

15、endelsewhile abs(y2-y1)/10=eh=2*h;y1,y2=comput(dyfun,x(n),y(n),h);endh=h/2;h=min(h,xspan(2)-x(n);y1,y2=comput(dyfun,x(n),y(n),h);endn=n+1;x(n)=x(n-1)+h;y(n)=y2;y1,y2=comput(dyfun,x(n),y(n),h);endx=x;y=y;function y1,y2=comput(dyfun,x,y,h);y1=rk4(dyfun,x,y,h);y21=rk4(dyfun,x,y,h/2);y2=rk4(dyfun,x+h/2,

16、y21,h/2);function y=rk4(dyfun,x,y,h);K1=feval(dyfun,x,y);K2=feval(dyfun,x+h/2,y+h/2*K1);K3=feval(dyfun,x+h/2,y+h/2*K2);K4=feval(dyfun,x+h,y+h*K3);y=y+h*(K1+2*K2+2*K3+K4)/6;五. 試驗(yàn)例題例 1不同方法的精度比較用多種方法解下述初值問題，并與其準(zhǔn)確解ye xx 進(jìn)行比較。yyx1,0 x0.6y(0)1h=, xk=khk= ,。各方法計算結(jié)果及絕對誤差見表、解：取步長,6)(1)(2)0.1(01(3) 。表（ 1）xn歐拉

17、公式改進(jìn)的歐拉公式四階標(biāo)準(zhǔn)龍格庫塔公式y(tǒng)n誤差yn誤差yn誤差0.01.000000010100.11.0000001.00500001.004837500.21.0100001.0190251.018730900.31.0290001.0412181.040818420.41.0561001.0708021.070320290.51.0904901.1070761.10653093精彩文檔實(shí)用標(biāo)準(zhǔn)文案0.61.1314411.1494041.14881193表（ 2）四階阿當(dāng)姆斯公式nxn顯示公式隱式公式 yn誤差yn誤差30.3取自準(zhǔn)確解1.0408180140.41.070322921.

18、0703196650.51.106535481.1065304160.61.148814811.14881101備注y1y2,y3 取自準(zhǔn)確解y1,y2 取自準(zhǔn)確解, 本題關(guān)于 y 為線性，代入隱式公式， 可解出 yn+1，因此可直接求隱式公式之解。表（ 3）四階阿當(dāng)姆斯預(yù)測校正公式nxn顯示公式隱式公式 Yn誤差yn誤差40.41.070319921.0703201450.51.106530271.1065307760.61.148811031.14881175備注y1,y2 y3 由四階標(biāo)準(zhǔn)龍格庫塔公式提供,對比以上各表數(shù)據(jù)知，在相同步長下求解同一問題時，方法階數(shù)越高，解的精度也越高，一階

19、的歐拉公式精度最低， 四階標(biāo)準(zhǔn)龍格庫塔公式的精度大大高于改進(jìn)的歐拉公式。四階阿當(dāng)姆斯方法、 顯式的精度略低于隱式。 同是阿當(dāng)姆斯預(yù)測校正系統(tǒng)，帶誤差修正的精度又高于不帶誤差修正的。從計算量看，四階龍格庫塔法計算量最大，每前進(jìn)一步要計算4 次函數(shù)值f ，而帶誤差修正的阿當(dāng)姆斯預(yù)測校正法與它精度差不多的，每前進(jìn)一步只計算兩次函數(shù)值 f ，所以后者是可取的。但它是四步法，必須用其它方法提供出始值，程序略復(fù)雜些。例 2 ( 步長的計算結(jié)果的影響 ) 用歐拉公式求下述初值問題在 x = 1 處的近似解，并與準(zhǔn)確解 y (1)=1 比較。解分別取步長h=10-1 、h=10-2 、h=10-3 、 h=1

20、0-4 計算，結(jié)果見下表。由表中數(shù)據(jù)可見，當(dāng)時結(jié)果完全失真，而取比計算量增加了十倍，但解的精度卻基本一樣，可見取太浪費(fèi)計算量了。精彩文檔實(shí)用標(biāo)準(zhǔn)文案步長 hy (1)的計算值上述結(jié)果差異很大的原因在于歐拉公式的絕對穩(wěn)定區(qū)間為（2, 0 ）步長 h 應(yīng)滿足f h( 2,0) ，對本題，f ( x, y)1000( yx2 )2xf1000 ，故應(yīng)取h 滿yy足即。可見取時歐拉公式是數(shù)值不穩(wěn)定的，導(dǎo)致結(jié)果失真，而取和都滿足穩(wěn)定性要求，可用于求解。此例說明，求解微分方程數(shù)值解一定選取注意步長，過大則導(dǎo)致解的失真，過小又會使計算量大增。究竟取多大步長才合適，不僅取決于所采用的數(shù)值方法， 還決定于待解微

21、分方程本身的特性。例 3： (Matlab不同命令之差異) 取 h=0.2 分別用不同的Matlab 命令解方程y y - 2t/y其準(zhǔn)確解為 : y1 2ty(1)(0tt,y=ode45(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n); n,eans =45.000000000000000.00020257808813t,y=ode23(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n); n,eans =13.000000000000000.1904836011

22、3013t,y=ode113(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n); n,eans =18.000000000000000.00973277413962t,y=ode23t(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n); n,eans =18.000000000000000.03919668151270t,y=ode23s(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n); n,eans

23、=81.000000000000002.54374270001824精彩文檔實(shí)用標(biāo)準(zhǔn)文案t,y=ode23tb(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =15.000000000000000.24308870757312t,y=ode15s(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =22.000000000000000.45509043060378可見 ode45 精度較高，計算量較大， ode23 計算量小，但誤差大，

24、 ode113 適中，用剛性方程組揭發(fā)解非剛性問題不適合，特別 ode23s 計算量大且但誤差大。例 4：解剛性方程組：y0.01y - 99.99z,y(1)2z -100z,z(0)1其準(zhǔn)確解為 : y e 0.01xe 100 x , ze 100 x 就 h=0.01、 h=0.02使用改進(jìn)歐拉格式及方程組，并在一張圖上繪出其圖形。functionx,y=naeuler2s(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y=zeros(length(y0),length(x);y(:,1)=y0(:);for n=1:length(x)-1K1=feva

25、l(dyfun,x(n),y(:,n);y(:,n+1)=y(:,n)+h*K1;K2=feval(dyfun,x(n+1),y(:,n+1);y(:,n+1)=y(:,n)+h*(K1+K2)/2;endx=x;y=y; %用歐拉格式解常微分方程組的m文件function f=dyfun(x,y)f(1)=-0.01*y(1)-99.99*y(2);f(2)=-100*y(2);f=f(:); %函數(shù) m文件clear;x,y=naeuler2s(dyfun,0 500,2,1,0.01);plot(x,y),axis(-50 500 -0.5 2)gtext(y(x) h=0.01), g

26、text(z(x) h=0.01) %步長為 0.01時解方程并畫圖clear;x,y=naeuler2s(dyfun,0 500,2,1,0.02);hold onplot(x,y),gtext(y(x) h=0.02), gtext(z(x) h=0.02) %步長為 0.02 時解方程并畫圖。clear;x=0:0.01:500;y=exp(-0.01*x)+exp(-100*x);z=exp(-100*x);hold onplot(x,y,r,x,z,r),gtext(y(x)真根 ), gtext(z(x)真根 ) % 畫出真根圖精彩文檔實(shí)用標(biāo)準(zhǔn)文案形運(yùn)行結(jié)果真根與 h=0,01 時

27、結(jié)果重合， h=0.02 時結(jié)果是錯誤的。故剛性方程組應(yīng)采用穩(wěn)定性較好方法。例 6：數(shù)值實(shí)驗(yàn)（摘自蔡）觀察歐拉顯示方法的收斂性。內(nèi)容：取 h=0.1, 0.01,用歐拉顯式方法求解： 用歐拉法解方程yxy1/3y(1)(0x1).1計算到 y(2)，并與精確解 y( x)( x22)/33/ 2比較。MATLAB程序t=1; y=1;h=0.1for k=1:10;y=y+h*t*y(1/3)t=t+hend運(yùn)行結(jié)果hy0.10.010.0010.0001t111111.11.11.11751.10681.10681.21.21361.23931.22771.22791.31.34161.37

28、621.36391.36411.41.48491.52941.51621.51651.51.64471.69981.68571.68611.61.82171.88831.87341.87391.72.0172.09622.08042.08101.82.23192.32432.30762.30831.92.46712.57392.55632.55712.02.72392.81772.82742.8283歐拉顯示方法是收斂的。例 7： 實(shí)驗(yàn)?zāi)康?：比較不同算法用于“非光滑”解時的結(jié)果。y=(t2+2)/3) 3/2(準(zhǔn)確解 )11.106816606308381.227879593566201.3

29、64135990288361.516564538686041.686170601183731.873981856902572.081044689573002.308421169608422.557186539930162.82842712474619內(nèi)容：用歐拉顯式方法和 2 階龍格 庫塔算法求解： y sinkx , y(0) 1, 計算 y(2 ) 并在一張圖上畫出解曲線。精彩文檔實(shí)用標(biāo)準(zhǔn)文案k=6 時的程序如下：cleart=0;y=1;h=pi/60;for k=0.1:0.1:2;y=y+h*abs(sin(6*t)t=t+hplot(t,y,o), hold on,endt0=0;

30、tf=1/3*pi;x0=1;k=pi/60t,x=ode23(ex545f,t0:k:tf,1)plot(t,x), hold on,for i=0:pi/60:pi/3if i T,Y=ode45(rigid,0 12,0 1 1);程序見數(shù)學(xué)建模高教出版社數(shù)學(xué)建模與plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+)數(shù)學(xué)實(shí)驗(yàn) To Matlab （ ff5 ）3、結(jié)果如圖 .【例11】采用最簡潔格式的ODE文件和解算指令，研究圍繞地球旋轉(zhuǎn)的衛(wèi)星軌道。（ 1）問題的形成軌道上運(yùn)動的衛(wèi)星，在Newton 第二定律Fmad 2rm2，和萬有引力定律dtFG mM Er

31、作用下，有r 3d 2 rGM EGM Exyx2y2adt2r3 r 。即 axr3 ， ayGM E3 ，而 r。這里rG6.672 10 11 (N m2/ kg 2 ) 是引力常數(shù)， M E5.97 1024 ( kg) 是地球的質(zhì)量。又假定衛(wèi)星以初速度 vy (0)4000( m / s) 在 x(0)4.2107(m) 處進(jìn)入軌道。（ 2）構(gòu)成一階微分方程組令 Yy1y2 y3y4Tx y vxvyTx y x y T ，則y1y3y4y2y1Y (t)GM E( x2(5.14.2.1-5)y3y 2 ) 3/ 2y4GM E( x2y2y 2 ) 3/ 2初始條件為Y(0)x(

32、0)00 v y ( 0) T(5.14.2.1-6)（ 1）根據(jù)式 (5.14.2.1-5) 編寫最簡潔格式的 ODE文件DYDt2.mfunctionYd=DYDt2(t,Y,flag,G,ME)% flag按 ODE文件格式規(guī)定，必須是第三輸入宗量。對它的賦值由 ode45指令自動產(chǎn)生。第 4、 5宗量是被傳遞的參數(shù)switchflag精彩文檔實(shí)用標(biāo)準(zhǔn)文案case %按規(guī)定：這里必須是空串。在此為 真 時，完成以下導(dǎo)數(shù)計算。X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.2);Yd=V;-G*ME*X/r3; otherwiseerror(Unknown flag flag.);end（ 2）對微分方程進(jìn)行解算G=6.672e-11;ME=5.97e24;vy0=40