常微分方程初值問題數(shù)值解法的比較

1、常微分方程初值問題數(shù)值解法的比較數(shù)值計(jì)算實(shí)踐一課程設(shè)計(jì)報(bào)告課題名稱常微分方程初值問題數(shù)值解法的比較完成時(shí)間2013-1-17姓名班級(jí)學(xué)號(hào)成績一.實(shí)驗(yàn)?zāi)康募皟?nèi)容1實(shí)驗(yàn)?zāi)康模?1)了解常微分方程初值問題的理論背景以及初值問題穩(wěn)定性、收斂性的研究；熟練掌握歐拉法、改進(jìn)歐拉法、龍格-庫塔法以及截?cái)嗾`差分析；比較歐拉法、改進(jìn)歐拉法及龍格-庫塔法，能夠選擇合適的方法進(jìn)行問題的研究計(jì)算； ,2實(shí)驗(yàn)內(nèi)容：求微分方程yy (0 X 1)(歐拉法求解)y(0) = 1_ 2X求微分方程yy y (0 x 0,使1得| f (x, y1) | -f (x, y2) L | y1 - y2 |, yyy y2 w R

2、,則稱f關(guān)于y滿足利普希茨條件，L稱為利普希茨常數(shù)。對(duì)于常微分方程初值問題y = f (t, y), t _ 10 , 一、,00 ,考慮初值y0的擾動(dòng)是問題的解y(t)發(fā)生偏差的情y(t0) = v。形。若t t 笛時(shí)y(t)的偏差被控制在有界范圍內(nèi)，則稱該初值問題是穩(wěn)定的，否則該初值問題 不穩(wěn)定的。特別地，若t tm時(shí)y(t )的偏差收斂于零，則稱該初值問題是漸進(jìn)穩(wěn)定的。對(duì)于初值問題Jy - y ,t _ 10 _ ，一(t 3),、,0穩(wěn)定性的研究，易知其準(zhǔn)確解為 y(t) = y留*0),假定初值經(jīng)過y(t0) = v。擾動(dòng)后變?yōu)閥0十Ay0 ,對(duì)于擾動(dòng)后的解為y(t) = (y 0

3、+ Ay0)e*,0)因此帶來的擾動(dòng)誤差為y(t) = Ay0e其上0),因此考慮t t 餡時(shí)Ay(t )的值，它取決于 九。易知，若 九w 0 ,則原問題是穩(wěn)定的；若 九a 0,則原問題是不穩(wěn)定的；若九 0,則原問題是漸進(jìn)穩(wěn)定的。實(shí)際遇到的大多數(shù)常微分方程初值問題都是穩(wěn)定的，因此在后面的討論數(shù)值解法時(shí)這常常是默認(rèn)條件。依據(jù)：積分曲線上一點(diǎn)(x, y)的切線斜率等于函數(shù)值。方法：推進(jìn)法，初始點(diǎn) p0(x0, y0)出發(fā)，依照方向場在改點(diǎn)的方向推進(jìn)到xi,x2，xn向前歐拉法的得到：y(x)y(xn i) = y(xn h) = y(xn)y (xn)hy,(xn)h2y (xn)h2!十取h的

4、線性部分，得Yn 1 二 Yn 八(。*)2 ) 將初值問題中得導(dǎo)數(shù)用差商來代替y(xn 1) - y(xn)y(Xn 1)- y(xn)xn 1 - xn=f(xn,yn)二ynhf(xn,yn)公dy將丁 =f(x,y)兩邊同 dx時(shí)對(duì) x 的區(qū)間xn,xn+i上 積xn 1 y dxxnxn 1-x n _1=p(xn, yn)dx,即y(xn書)-y(xn) = Jf (xn, y(x)dx 對(duì)右端x nxnxn 1ff (xn, y(x)dx 用左xn矩形公式得yn書=yn +hf(xn,yn),此方法稱向前歐拉法，也叫 顯示歐拉法。Xn 1(4)對(duì)右端 jf (xn, y(x)dx

5、用右矩形公式得yn卡=yn + hf (xn#, y5)，也 xn叫隱式歐拉法。誤差分析：1.稱y(xnQ _yn4為計(jì)算yn邛時(shí)的局部截?cái)嗾`差； TOC o 1-5 h z 2.如果數(shù)值方法的局部截?cái)嗾`差為 qhp*),那么稱這種數(shù)值方法的階數(shù)是p,其實(shí)p為非負(fù)整數(shù)。通常情況下，步長 h越小，p越高，則局部截?cái)嗾`差越小，計(jì)算精,，, y (xn)hy (xn)h HYPERLINK l bookmark34 o Current Document .、生叱由療y(xn1) =y(xnh)二y(xn)y(xn)h-nn .y(x)在xn初泰勒展開有八 八n /八 y ”2!32.=ynhf(xn

6、, yn) qh )則有y(xn卡) -yn4i = Qh2)可見歐拉方法是一階方法，精度不是很高。2.改進(jìn)歐拉方法：xn .1梯形公式：對(duì)右端Jf (xn, y(x)dx用梯形公式得+顯然梯形公式是隱式公式。 xn改進(jìn)歐拉公式： 先用歐拉公式求的一個(gè)初步的近似值yn書，成為預(yù)測值，預(yù)測值 yn卡的精度可能達(dá)不到要求，在用梯形公式將他校正一次，記為yn41,這個(gè)結(jié)果成為校正值。預(yù)測：Yn .1 = Yn , hf(xn,yn)h校正：yn l = yn 21f(xn, yn) f(xn i, yn l)誤差分析：記Tn中為改進(jìn)歐拉公式在xn 4i處的截?cái)嗾`差，Tn+1 = y( xn+)- y

7、(xn )丁一，、 h記 Tn 書=y(xn/ - y(xn) - - f ( x n , y( x n ) + f (x n 41 , y( x n+1 )因此有h 一Tn4=y(xn書)-y(xn)=，書+萬f ( x n如y ( x n由)f ( x n由，y( x n+),八中表不在*門書出的局到相應(yīng)的計(jì)算公式。歐拉公式可以寫為Ki = f (Xn，yn)h 、y n 4i =yn ( Ki + K2 )改進(jìn)歐拉公式可以寫成 Ki = f ( Xn, y n)因此推出一般的推出廣式K2 = f(Xn+i,yn + hKi)yn 平=yn +h(CiKi +C2K2 + . +CPKP)

8、Ki = f (Xn，yn)jK2 =f(xn +a2h,yn + hb21Kl)稱為p階龍格-庫塔方法，簡稱p階R-K方法。.pKP = f(Xn +aph, yn + h bp K )1i金r(Xn，yn，h) = Ci Ki， i i因?yàn)?Ki =f(Xn,yn),這里的6,兒，也均為常數(shù)。i JKi = f(Xn -h,ynhij Kj )j i因?yàn)榻o定的系數(shù)不唯一，因此這里的常數(shù)有無窮多個(gè)解，下面是特殊情況下和一般情況下得結(jié)果(一階龍格-庫塔)當(dāng)r=i時(shí)，這就是歐拉法。,i,(二階龍格-庫塔)當(dāng)r=2時(shí)，Ci = , % = = i ,這就是改進(jìn)歐拉法。2j三階和四階龍格-庫塔也只是

9、在一般情況下得結(jié)果。hyn 書=yn + -(Ki +K2 + K3) 6Ki =f(Xn，yn)三階K2 =f(Xn + , +2Ki)K3 =f(Xn +h,y0- hKi +2h9)yn書=yn十h(K +26 +2K3 + K4) 6Ki =f(Xn，yn)rhh四階TK2 =f(Xn + 2，yn + 2 Ki)K2 =f(Xn + -2，yn + hKz) 2K3 =f(Xn + h,yn 十 h()Tn由=y(Xn中)-y(Xn) -h(4-iki + %kz)其局部截?cái)嗾`差是:三.程序代碼歐拉法代碼如下：(1)向前歐拉法：function y=Euler1(fun,x0,y0,

10、xN,N)%fun為一階微分方程，x0,y0為初始條件，xN為取值范圍的一個(gè)端點(diǎn)，N為區(qū)間個(gè)數(shù)x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N; % h為區(qū)間步長for (n=1:N)x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n);%艮據(jù)向前歐拉公式計(jì)算y 值endT=x,y(2)向后歐拉法：function y=Euler2(fun,x0,y0,xN,N)% fun為一階微分方程，x0 , y0為初始條件，xN為取值范圍的一個(gè)端點(diǎn)，N為區(qū)間個(gè)數(shù)x=zeros(1,N+1);x=z

11、eros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;% h為區(qū)間步長for (n=1:N)x(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n);for (k=1:3)z1=y(n)+h*feval(fun,x(n+1),z0);if (abs(z1-z0)1e-3) break ;endz0=z1;endy(n+1)=z1;%艮據(jù)向后歐拉公式計(jì)算y值endT=x,y改進(jìn)歐拉代碼如下：function x,y=Gaijineuler(f,x0,y0,xZ,h)%f為一階微分方程的函數(shù)；x0,y0為初始條件；xZ為取值范圍的一個(gè)端點(diǎn)，h

12、為區(qū)間步長n=fix(xZ-x0)/h);%計(jì)算分點(diǎn)數(shù)y(1)=y0;x(1)=x0;for i=1:nx(i+1)=x0+i*h;yp=y(i)+h*feval(f,x(i),y(i);yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=(yp+yc)/2;endx=x;y=y;1.3.龍格-庫塔代碼如下： (三階龍格-庫塔)function R=Longgekuta3(f,a,b,aZ,h)%a, b為端點(diǎn)，h為步長，aZ為初值 n=(b-a)/h;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=aZ;for j=1:nk1=fev

13、al(f,T(j),Y(j);k2=feval(f,TQ)+h/2,Y(j)+k1*h/2);k3=feval(f,T(j)+h,Y(j)-h*k1+k2*2*h);Y(j+1)=Y(j)+(k1+4*k2+k3)*h/6;endR=T Y;(四階龍格-庫塔)function R=Longgekuta4(f,a,b,aZ,h)% a, b為端點(diǎn)，h為步長，aZ為初值 n=(b-a)/h;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=aZ;for j=1:nk1=feval(f,T(j),Y(j);k2=feval(f,T(j)+h/2,Y(j)+k1*h

14、/2);k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);k4=feval(f,T(j)+h,Y(j)+k3*h);Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)*h/6;endR=T Y;四.數(shù)值結(jié)果：% 艮據(jù)改進(jìn)歐拉公式計(jì)算結(jié)果%定義向量%計(jì)算各個(gè)分點(diǎn)%初值賦予%根據(jù)公式計(jì)算%定義向量%計(jì)算各個(gè)分點(diǎn)%初值賦予%根據(jù)公式計(jì)算輸入：定義M文件并保存ffx.mEuler1(ffx,0,1,1,10)結(jié)果：T =01.00000.10001.10000.20001.21000.30001.33100.40001.46410.50001.61050.60001.77160

15、.70001.94870.80002.14360.90002.35791.00002.5937ans =Columns 1 through 71.77161.00001.10001.21001.33101.46411.6105Columns 8 through 111.94872.14362.35792.5937向前歐拉公式結(jié)果：輸入：Euler2(ffx,0,1,1,10)結(jié)果T =01.00000.10001.11100.20001.23440.30001.37160.40001.52400.50001.69330.60001.88140.70002.09040.80002.32270.9

16、0002.58071.00002.8674ans =Columns 1 through 71.00001.11101.23441.37161.52401.69331.8814Columns 8 through 112.09042.32272.58072.8674改進(jìn)歐拉公式結(jié)果：輸入：x,y=Gaijineuler(f,0,1,1,0.1)結(jié)果：x =00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000y =1.00001.09591.18411.26621.34341.41641.48601.55251.61651.678

17、21.7379龍格-庫塔計(jì)算結(jié)果：(三階)輸入：Longgekuta3(ff,0,1,1,0.1)結(jié)果：ans =Longgekuta3(ff,0,1,1,0.1)8ans =00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000(四階)輸入：結(jié)果:00.10000.20000.30000.40000.50000.60001.00001.10481.21881.34111.47111.60821.75201.90212.05802.21952.3863Longgekuta4(ff,0,1,1,0.1) ans =1.00001

18、.10481.21881.34111.47111.60821.75200.70001.90210.80002.05800.90002.21951.00002.3863五.計(jì)算結(jié)果分析方法顯示歐拉簡單；精度低隱式歐拉穩(wěn)定性最好精度低，計(jì)算量大梯形公式精度提高計(jì)算量大中點(diǎn)公式精度提高，顯式:多一個(gè)初值，可能影響精度1.收斂性：若某算法對(duì)于任意固定的 x = Xi = x0 + i h,當(dāng)20 (同時(shí)i t8)時(shí)有yiT y( xi),則稱該算法是收斂的。以下討論的都是單步法(指在計(jì)算 yn 卡時(shí)只用到它前一步的信息yn)歐拉法，龍格-庫塔法都 是單步法的例子； V =,t 之 t C例：就初值問題

19、 y 0考察歐拉顯式格式的收斂性。了(t ) = y解：該問題的精確解為y(t) = y0e *工)歐拉公式為yi中=yi+h，y =(1十h九)yixi1對(duì)于固定的 X = Xi = X0+ i h,有 v =y0(1+Kh)h = y0( 1 + Ih)沖NNoe ，= y( Xi )設(shè)單步法具有中階精度，滿足利普希茨條件，其整體的截?cái)嗾`差：y(Xn) - yn = qhp)判斷單步法的收斂性，歸結(jié)為驗(yàn)證增量函數(shù)邛能否滿足利普希茨條件。對(duì)于歐拉法，由于其增量函數(shù) 中就是f (x, y),因此當(dāng)f (x, y)關(guān)于y滿足利普希茨條件 時(shí)候它就是收斂的。對(duì)于改進(jìn)歐拉法和龍格-庫塔方法，找到增量

20、函數(shù)，滿足利普希茨條件的時(shí)候，找到利普希 茨常數(shù)，這些方法都是收斂的。2.穩(wěn)定性:r 例：考察初值問題 ”=一 y,在區(qū)間0,0.5上的解，分別用歐拉顯，隱式格式，改進(jìn)歐拉格式計(jì) 7(0) = 1算數(shù)值解：10節(jié)點(diǎn)歐拉顯式改進(jìn)歐拉法精 確 解30 xy = e-01.00001.00001.00001.00000.1-2.00002.5*10A(-1)2.50004.9787*10A(-2)0.24.00006.2500*10A(-2)6.25002.4788*10A(-3)0.3-8.00001.5625*10A(-2)1.5625*10A(1)1.2341*10A(-4)0.4.6000*

21、10A(1)3.9063*10A(-3)3.9063*10A(1)16.1442*10A(-6)0.513.2000*10A(1)9.7656*10A(-1)9.7656*10A(1)3.0590*10A(-7)則稱該算法是絕對(duì)穩(wěn)若某算法在計(jì)算過程中任一步產(chǎn)生的誤差在以后的計(jì)算中都逐步衰減, 定的.i 1yi 1 = V h yi =(1h)y, ； = yyi 1=(1+ h)i卡斯,由此可見,要保證初始誤差a以后逐步衰減h =Kh必須滿足:| 1 h | ： 1考察后退歐拉法：yi卅=yih yi 1玩由此可見,1梯形公式：模型y=y, yn .1yn，h,、+ -(y n + yn +)

22、,穩(wěn)定區(qū)域?yàn)閨 2 + Ah |42 九h |相對(duì)穩(wěn)定區(qū)間為-二：h改進(jìn)歐拉公式:yn相對(duì)穩(wěn)定區(qū)域?yàn)閨 1四階龍格-庫塔:相對(duì)穩(wěn)定區(qū)域:1相對(duì)穩(wěn)定區(qū)間:yn 1yn 1ynhf(Xn,yn) = (1h)yn,hf(Xn，yn) f (Xn 1, yn 1 2,1, 2)=1 h 2( h)2yn1 ,. 、2 ,-, 八+ -(油)2 |M 1穩(wěn)定區(qū)間為2 Kh 021c 1a 1“H1 h 2( h)26( h)324( -h)4yn121314h 2(h)6(h)24(h)L- 2. 78 :二 h ：二 01 - h J據(jù)對(duì)絕對(duì)穩(wěn)定絕對(duì)穩(wěn)定 絕對(duì)穩(wěn)定的區(qū)域?yàn)閨 1 - h | . 1節(jié)

23、點(diǎn)改進(jìn)歐拉法四節(jié)龍格-庫塔法準(zhǔn)確解y(0) = 1-12x,、，“在區(qū)間0,1上的解考察初值問題y1101110.21.1866671.1832291.1832160.41.3483121.3416671.3416410.61.4937041.4832811.4832400.81.6278611.6125141.61245211.7542051.7321421.732051變步長的龍格-庫塔法：實(shí)際計(jì)算中，步長h的選擇是一個(gè)十分重要的問題，若步長h取的太小，將增加許多不必要的的計(jì)算量，若步長取的太大，其計(jì)算競速很難保證，因此，與數(shù)值積分一樣，微分方程的數(shù)值解也有一個(gè)選擇步長的問題。由四階龍格-

24、庫塔法的性質(zhì)，其局部截?cái)嗾`差為y(xn+1) y：?定ch5這h()一個(gè)近似值yn.,n里的c與y5(x)在xn,xn中內(nèi)的值有關(guān)系。然后將步長折半，即取為步長，從xn跨兩步到xn中，在求h 5、一22)5 ,比較之后h 5(；)每跨一步的誤差為ci-J ,因此y(xn中)_yn1一誤差的大約減少16116h (夕) y(xn .1)- yn2.1y( x n .1) - y n 六.計(jì)算時(shí)出現(xiàn)的問題，解決方法及體會(huì)：對(duì)于初值問題，在用不同方法計(jì)算離散點(diǎn)的函數(shù)值的時(shí)候，有精確值給定的情況下，不同方法跟精確值 的差別很大，也就是不是一種好方法，這與實(shí)際問題的穩(wěn)定性有很大的關(guān)系，只要所求的初值問題

25、，能夠落 在使用方法的穩(wěn)定區(qū)間內(nèi)，所得到的結(jié)果和精確值的離散點(diǎn)的函數(shù)值相差不會(huì)很大，并且在穩(wěn)定區(qū)間內(nèi)，改 進(jìn)歐拉法優(yōu)于歐拉法，三階之后的龍格-庫塔方法優(yōu)于改進(jìn)歐拉法，離散點(diǎn)的函數(shù)值更加接近精確值。 所謂的歐拉法就是一階的龍格-庫塔方法，改進(jìn)歐拉法也是在顯示歐拉法預(yù)測和梯形公式校正的基礎(chǔ)上作出的改進(jìn)， 也是二階龍格-庫塔的特殊情況。所謂的三階，四階甚至n階龍格-庫塔方法，也是在待定系數(shù)有無數(shù)解得情況下的一般情況。一階常微分初值問題的收斂性，只要其增量函數(shù)關(guān)于y能滿足利普希茨條件，這就收斂。而起在計(jì)算過程過有避免不了的誤差， 歐拉法精度相對(duì)來說比較低，龍格-庫塔的階數(shù)越高，精度就越高，所得的結(jié)果和

26、精確值越接近，截?cái)嗾`差越小。同一個(gè)問題，想要減少截?cái)嗾`差，我們可以采用變步長的的方法，將步長變?yōu)閔h,這樣得出的結(jié)果，根據(jù)不同的方法，截?cái)嗾`差會(huì)有相應(yīng)的改變。2而在每種方法的程序設(shè)計(jì)方面，主要思想還是根據(jù)每種方法的公式，找到步長，分點(diǎn)數(shù)，用循環(huán)語句， 得出相應(yīng)的結(jié)果，只是在龍格 -庫塔方法里面的三階，四階的公式里面的待定系數(shù)有無數(shù)個(gè)解，我們給出的也 是一般情況的三階，四階龍格 -庫塔方法，精度，誤差相對(duì)于歐拉法，改進(jìn)歐拉法，梯形公式都優(yōu)于他們。數(shù)值計(jì)算實(shí)踐一課程設(shè)計(jì)原始數(shù)據(jù)記錄實(shí)驗(yàn)名稱：常微分方程初值問題數(shù)值解法研究 實(shí)驗(yàn)時(shí)間： 2013年 1月 14 日12輸入：定義M文件并保存ffx.mE

27、uler1(ffx,0,1,1,10)結(jié)果：T =01.00000.10001.10000.20001.21000.30001.33100.40001.46410.50001.61050.60001.77160.70001.94870.80002.14360.90002.35791.00002.5937 ans =Columns 1 through 71.00001.10001.21001.33101.46411.61051.7716Columns 8 through 111.94872.14362.35792.5937向前歐拉公式結(jié)果：輸入：Euler2(ffx,0,1,1,10)結(jié)果T =

28、01.00000.10001.11100.20001.23440.30001.37160.40001.52400.50001.69330.60001.88140.70002.09040.80002.32270.90002.58071.00002.8674 ans =13Columns 1 through 71.00001.11101.23441.37161.52401.69331.8814Columns 8 through 112.09042.32272.58072.8674改進(jìn)歐拉公式結(jié)果：輸入：x,y=Gaijineuler(f,0,1,1,0.1)結(jié)果：x =00.10000.20000.30000.40000.50000.600