數(shù)值分析試驗報告之常微分方程數(shù)值解匯總

1、沙理工大數(shù)學與計算科學學院實驗報告實驗項目名稱常微分方程數(shù)值解所屬課程名稱數(shù)值方法B實驗類型驗證實驗日期2013.11.11班級學號姓名一、實驗概述：【實驗目的】1 .掌握求解常微分方程的歐拉法；2 .掌握求解常微分方程的預估校正法；3 .掌握求解常微分方程的經(jīng)典的四階龍格庫塔法；4 .能用C語言或MATLAB將上述三種算法用程序運行出來；5 .將算法實例化，并得出三種算法的相關關系，如收斂性、精度等;6 .附帶書中例題的源程序見附錄1?！緦嶒炘怼?.歐拉格式(1)顯式歐拉格式：ynjH=yn+hf(xn,yn)h2.h2c局部截斷反差：y(Xn1)-yn1y(八一y(Xn)=o(h)22(

2、2)隱式歐拉格式：yn書=yn+hf(Xn如yn書)h2.c局部截斷反差：y(Xn1)-yn1：-y(Xn)=o(h)22.預估校正法預估：yn1=yn.hf(Xn,yn)h校正：yn-yn2f(Xn,Yn)f(Xn1,Yn1)統(tǒng)一格式：yn+=yn+llf(Xn,yn)+f(Xn+h,yn+hf(Xn,yn)yp=yn+hf區(qū)*),平均化格式：-yc=yn+hf%,yp),1/、yn中=2(yp+yc).*LJ3.四階龍格庫塔方法的格式(經(jīng)典格式)h%4=丫0+72+2&+2七十(),K1=f(xn,yn),一hhK2=f(Xno，YnoK1),22K3nfn+gK4=f(Xnh,y

3、nhK3).【實驗環(huán)境】1 .硬件環(huán)境：HPMicrosoft76481-640-8834005-23929HPCorporationIntel(R)Core(TM)I5-2400CPU3.10GHz3.09GHz,3.16GB的內(nèi)存2 .軟件環(huán)境：MicrosoftWindowsXPProfessional版本2002ServicePack3二、實驗內(nèi)容:【實驗方案】方案一：用歐拉法，預估校正法，經(jīng)典的四階龍格庫塔方法求解下列ODE問題：例題：在區(qū)間【0,1上以h=0.1用歐拉法，預估校正法，經(jīng)典的四階龍格庫塔法求解微分方程dy/dx=-y+x+1,初值y(0)=1；其精確解為y=x+exp

4、(-x),且將計算結(jié)果與精確解進行比較，對三個算法的收斂性的進行分析比較。用歐拉法，預估校正法，經(jīng)典的四階龍格庫塔方法求解初值問題dy/dx=xe'-y,初值y(0)=1;將計算結(jié)果與精確解為-(x2+2)e"比較在區(qū)2間0,1上分別取步長h=0.1;0.05時進行計算。對三個算法的收斂性進行分析比較，【實驗過程】(實驗步驟、記錄、數(shù)據(jù)、分析)注：以下圖形是通過Excel表格處理數(shù)據(jù)得出，并未通過MATLAB編程序所得!1、年一y(0)=1由題可知精確解為：y=x+e",當x=0時，y(x)=0。h=0.1表1h=0.1時三個方法與精確值的真值表JkEuler法預估

5、校正法經(jīng)典四階庫精確值0.11.0100001.0050001.0048381.2490800.21.0290001.0190251.0187311.0554550.31.0561001.0412181.0408181.0912170.41.0904901.0708021.0703201.1318030.51.1314411.1070761.1065311.1768510.61.1782971.1494041.1488121.2260250.71.2304671.1972111.1965861.2790160.81.2874211.2499751.2493291.3355360.91.3486

6、781.3072281.3065701.3953221.01.4138111.3685411.3678801.4581270.8預估校正法門匚T一經(jīng)典四階庫U.Df-精確值0.40.2000.10.20.30.40.50.60.70.80.91圖1h=0.1時三個方法走勢圖h=0.05(此時將源程序中i的范圍進行擴大，即for(i=0;i<20;i+)表2h=0.05時三個方法與精確值的真值表JkEuler法預估校正法經(jīng)典四階庫精確值0.051.0025001.0012501.0012291.0117210.101.0073751.0048771.0048371.0249080.151.

7、0145061.0107641.0107081.0395040.201.0237811.0188021.0187311.0554550.251.0350921.0288851.0288011.0727100.301.0483371.0409151.0408181.0912170.351.0634211.0547951.0546881.1109310.401.0802501.0704361.0703201.1318010.451.0987371.0877521.0876281.1537910.501.1188001.1066621.1065311.1768510.551.1403601.1270

8、871.1269501.2009420.601.1633421.1489541.1488121.2260250.651.1876751.1721931.1720461.2520620.701.2132911.1967361.1965851.2790160.751.2401271.2225201.2223671.3068520.801.2681211.2494851.2493291.3355360.851.2972151.2775721.2774151.3650370.901.3273541.3067281.3065701.3953220.951.3584861.3369001.3367411.

9、4263621.001.3905621.3680391.3678801.4581271,60.8預估校正法一經(jīng)典四階庫T-精確值0.40.20npIIIIII|IfIIiii1iiii00A0203040.566070£0.91圖2h=0.05時三個方法走勢圖dy斗xe-y2、ddxy(0)=12x由題可知精確解為:(x+2)e,當x=0時，y(x)=0。2h=0.1表3h=0.1時三個方法與精確值的真值表JkEuler法預估校正法經(jīng)典四階庫精確值:0.10.9000000.9096250.9094280.9295330.20.8192490.8359270.8355930.8725

10、6410.30.7544330.7760810.7756550.826822F0.40.7027260.7276710.7271890.7903480.50.6617260.6886360.6881270.7614570.60.6293960.6572250.6567110.738709r0.7P0.6040180.6319570.6314530.720874r0.80.5841470.6115820.6111000.7069080.90.5685750.5950500.5945990.6959271.00.5562970.5814870.5810720.687191T-Eulei法預估校工法

11、經(jīng)典四階庫一精確值圖3h=0.1時三個方法走勢圖h=0.05(此時將源程序中i的范圍進行擴大，即for(i=0;i<20;i+)表4h=0.05時三個方法與精確值的真值表JkEuler法預估校正法經(jīng)典四階庫精確值0.050.9500000.9524520.9524270.9629240.100.9049040.9094740.9094280.92953310.150.86428410.8706700.8706060.8995110.200.8277410.8356710.8355920.872564m.250.794908I0.8041370.8040470.8484190.300.76

12、54470.7757550.7756550.8268220.350.7390430.7502320.7501250.807538|0.400.715407(0.7273020.7271890.7903480.450.69427210.7067150.7065990.7750500.500.6753940.6882450.6881260.7614570.550.6585460.6716820.6715610.74939710.600.64351910.6568300.6567100.73870910.650.63012410.6435140.6433950.7292470.700.6181850

13、.6315700.6314530.72087470.750.6075410.6208480.6207330.7134660.800.5980460.6112110.6111000.7069080.850.5895650.6025350.6024260.7010940.900.5819760.5947030.5945990.6959270.950.5751670.5876120.5875120.691320：1.000.5690350.5811670.5810710.687191EulerS-預估校正法T-經(jīng)典四階庫精蠲值圖4h=0.05時三個方法走勢圖【實驗結(jié)論】(結(jié)果)1 .預估校正法的精確

14、度比經(jīng)典的四階庫法的精確度高，庫塔法最低；2 .從表中數(shù)據(jù)可知三個算法所得數(shù)據(jù)與精確值相比，可得出以下結(jié)論(針對方案二)：(1)Euler法所得值偏離精確值最大，因此可知其精度相對來說最差；(2)經(jīng)典的四階庫塔法所得值與精確值距離較近，因此可知對于Euler來說,該法更加有效；(3)預估校正法的數(shù)據(jù)時距離精確值最近的，具驟減幅度較小，因此對精度上的考慮而言，預估校正法應屬于最佳解法；(4)由數(shù)據(jù)可知，上述兩個方程的解的光滑性都比較差，從而導致四階龍格庫塔法的精度低十預估校正法的精度。3.由圖形可知，三個算法所得數(shù)據(jù)均呈遞減趨勢，對于他們的收斂性有以下結(jié)論：(1)用上述三法得到的結(jié)果大致趨近于0

15、.581,相對于精確值來說，還是存在較大的誤差；(2)就誤差最小，應首選預估校正法對問題進行求解，它與經(jīng)典的四階庫塔法所得值比較接近，因此誤差不會相差太大?！緦嶒炐〗Y(jié)】(收獲體會)由此次試驗，心方面強化了自己的編程能力，另一方面也對(1)庫塔法，(2)預估校正法，(3)經(jīng)典的四階龍格庫塔法有了全新的認識，并能巧妙的將他們運用到數(shù)學建模中，解決一些追求高精度的問題。其次在使用上述三種方法時要充分考慮方程的解的光滑性質(zhì)，并對照光滑性的好壞選擇相應的求解方法，以達到想要的精度的目的。三、指導教師評語及成績：評語評語等級優(yōu)良中及格/、及格1.實驗報告按時完成，字跡清楚，文字敘述流暢，邏輯性強2.實驗方

16、案設計合理3.實驗過程（實驗步驟詳細，記錄完整，數(shù)據(jù)合理，分析透徹）4實驗結(jié)論正確.成績：指導教師簽名：批閱日期：附錄1:源程序1.書中例題:精確化#include<stdio.h>#include<math.h>main()(inti,p;floaty0,x0,yi,xi,yp,xp,h;xi=0.0;printf("請輸入步長h:");scanf("%f",&h);for(p=0;p<=10;p+)(printf("p=%d",p);10if(p=0)(xp=0.0;yp=1.0;)yp=sq

17、rt(1+2*xp);printf("x%d=%f,y%d=%fn",p,xp,p,yp)xp+=h;)Euler格式：#include<stdio.h>main()(inti,p;floaty0,x0,yi,xi,yp,xp,h;xi=0.0;printf("請輸入步長h:");scanf("%f',&h);for(i=1;i<=10;i+)(p=i-1;printf("i=%d",i);if(p=0)(xp=0.0;yp=1.0;)yi=yp+h*(yp-2*xp/yp);printf(

18、"t=%f",xp/yp);yp=yi;xi+=h;xp=xi;printf("x%d=%f,y%d=%fn",i,xi,i,yi);)預估校正格式：#include<stdio.h>main()(inti,t;floatyi,xi,m,xp,n,yt,xt,h;11xi=0.0;printf("請輸入步長h:");scanf("%f',&h);printf("t=0x0=0.000000,y0=1.000000n")for(i=0;i<10;i+)t=i+1;print

19、f("t=%d",t);if(i=0)xi=0.0;yi=1.0;xt=xi+h;m=yi+h*(yi-2*xi/yi);n=yi+h*(m-2*xt/m);yt=(m+n)/2.0;yi=yt;xi+=h;printf("x%d=%f,y%d=%fn",t,xt,t,yt);經(jīng)典的四階龍格庫塔方法的格式：#include<stdio.h>main()inti,t;floatyi,xi,K1,K2,K3,K4,xp,yt,xt,h;xi=0.0;printf("請輸入步長h:");scanf("%f",

20、&h);printf("t=0x0=0.000000,y0=1.000000n")for(i=0;i<10;i+)t=i+1;printf("t=%d",t);if(i=0)xi=0.0;yi=1.0;K1=yi-2*xi/yi;K2=yi+h*K1/2.0-(2*xi+h)/(yi+h*K1/2.0);12K3=yi+h*K2/2.0-(2*xi+h)/(yi+h*K2/2.0);K4=yi+h*K3-2*(xi+h)/(yi+h*K3);yt=yi+h*(K1+2*K2+2*K3+K4)/6.0;xt=xi+h;xi+=h;yi=yt;

21、printf("x%d=%f,y%d=%fn",t,xt,t,yt);)2.精確解：y=x+e，(方案一)#include<stdio.h>#include<math.h>#definee2.1828main()inti,p;floaty0,x0,yi,xi,yp,xp,h;xi=0.0;printf("請輸入步長h:");scanf("%f',&h);for(p=0;p<=10;p+)printf("p=%d",p);if(p=0)xp=0.0;yp=1.0;)yp=xp+po

22、w(e,-xp);printf("x%d=%f,y%d=%fn",p,xp,p,yp);xp+=h;)Euler格式：#include<stdio.h>main()inti,p;floaty0,x0,yi,xi,yp,xp,h;13xi=0.0;printf("請輸入步長h:");scanf("%f',&h);for(i=1;i<=11;i+)(P=i-1;printf("p=%d",p);if(p=0)(xp=0.0;yp=1.0;yi=yp+h*(-yp+xp+1);yp=yi;prin

23、tf("x%d=%f,y%d=%fn",p,xi,p,yi);xi+=h;xp=xi;預估校正格式：#include<stdio.h>main()(inti,t;floatyi,xi,m,xp,n,yt,xt,h;xi=0.0;printf("請輸入步長h:");scanf("%f",&h);printf("t=0x0=0.000000,y0=1.000000n");for(i=0;i<20;i+)(t=i+1;printf("t=%d",t);if(i=0)(xi=0

24、.0;yi=1.0;xt=xi+h;m=yi+h*(-yi+xi+1);n=yi+h*(-m+xt+1);yt=(m+n)/2.0;yi=yt;14xi+=h;printf("x%d=%f,y%d=%fn",t,xt,t,yt);)經(jīng)典的四階龍格庫塔方法的格式：#include<stdio.h>main()inti,t;floatyi,xi,K1,K2,K3K4,yt,xt,h;xi=0.0;printf("請輸入步長h:");scanf("%f',&h);printf("t=0x0=0.000000,y0

25、=1.000000n");for(i=0;i<10;i+)t=i+1;printf("t=%d",t);if(i=0)xi=0.0;yi=1.0;)K1=-yi+xi+1;K2=-(yi+h*K1/2.0)+xi+h/2.0+1;K3=-(yi+h*K2/2.0)+xi+h/2.0+1;K4=-(yi+h*K3)+xi+h+1;yt=yi+h*(K1+2*K2+2*K3+K4)/6.0;xt=xi+h;xi+=h;yi=yt;printf("x%d=%f,y%d=%fn",t,xt,t,yt);)3.精確解：1(x2+2)e"(

26、方案二)2#include<stdio.h>#include<math.h>#definee2.1828main()inti,p;15floaty0,x0,yi,xi,yp,xp,h;xi=0.0;printf("請輸入步長h:");scanf("%f',&h);for(p=0;p<=10;p+)printf("p=%d",p);if(p=0)xp=0.0;yp=1.0;yp=(xp*xp+2)*pow(e,-xp)/2.0;printf("x%d=%f,y%d=%fn",p,x

27、p,p,yp);xp+=h;Euler格式：#include<stdio.h>#include<math.h>#definee2.1828main()inti,p;floaty0,x0,yi,xi,yp,xp,h;xi=0.0;printf("請輸入步長h:");scanf("%f",&h);printf("i=0x0=0.000000,y0=1.000000n");for(i=1;i<=10;i+)p=i-1;printf("i=%d",i);if(p=0)xp=0.0;yp

28、=1.0;yi=yp+h*(xp*pow(e,-xp)-yp);yp=yi;xi+=h;xp=xi;printf("x%d=%f,y%d=%fn",i,xi,i,yi);16)預估校正格式：#include<stdio.h>#include<math.h>#definee2.1828main()(inti,t;floatyi,xi,m,xp,n,yt,xt,h;xi=0.0;printf("請輸入步長h:");scanf("%f",&h);printf("i=0x0=0.000000,y0=1.000000n");for(i=0;i<20;i+)(t=i+1;printf("t=%d",t);if(i=0)(xi=0.0;yi=1.0;)xt=xi+h;m=yi+h*(xi*pow(e,-xi)-yi);n=yi+h*(xt*pow(e,-xt)-m);yt=(m+n)/2.0;yi=yt;xi+=h;printf("x%d