MATLAB改進(jìn)歐拉法與四階龍格-庫塔求解一階常微分方程

姓名：樊元君 學(xué)號：2012200902 日期：2012.11.06一、實驗?zāi)康恼莆誐ATLAB語言、C/C++語言編寫計算程序的方法、掌握改進(jìn)歐拉法與四階龍格-庫塔求解一階常微分方程的初值問題。掌握使用MATLAB程序求解常微分方程問題的方法。、實驗內(nèi)容1、分別寫出改進(jìn)歐拉法與四階龍格-庫塔求解的算法，編寫程序上機(jī)調(diào)試出結(jié)果,要求所編程序適用于任何一階常微分方程的數(shù)值解問題，即能解決這一類問題，而不是某一個冋題。實驗中以下列數(shù)據(jù)驗證程序的正確性。&&-*~ /，步長h=0.252、實驗注意事項yr=一科亠応—二 -必-匯的精確解為- ，通過調(diào)整步長，觀察結(jié)果的精度的變化三、程序流程圖:?改進(jìn)歐拉格式流程圖:?四階龍格庫塔流程圖:開始初!|點CxO,yO);(K.Y)xl=xO+h;kl=f(xO,yO)；k2=f(x0+h/2ty0+h/2*kl)；k3=f(xO+h/2.y0+h/2*k2)；k4=f(xl,yO+h*k3)；yl=yO+h/6*(k!+2*k2+2*k3+k4).i=i+l；xO=xl；yO=y1四、源程序:?改進(jìn)后歐拉格式程序源代碼:function[]=GJOL(h,x0,y0,X,Y)formatlongh=input('h=');x0=input('x0=');y0=input('y0=');disp('輸入的范圍是：');X=input('X=');Y=input('Y=');n=round((Y-X)/h);i=1;x1=0;yp=0;yc=0;fori=1:1:nx1=x0+h;yp=yO+h*(-xO*(yOF2);%yp=yO+h*(yO-2*xO/yO);%yc=y0+h*(-x1*(ypf2);%yc=y0+h*(yp-2*x1/yp);%y1=(yp+yc)/2;x0=x1;y0=y1;y=2/(1+x0A2);%y=sqrt(1+2*x0);%fprintf('結(jié)果=%.3f,%.8f,%.8f\n',x1,y1,y);endend?四階龍格庫塔程序源代碼：function[]=LGKT(h,x0,y0,X,Y)formatlongh=input('h=');x0=input('x0=');y0=input('y0=');disp('輸入的范圍是：');X=input('X=');Y=input('Y=');n=round((Y-X)/h);i=1;x1=0;k1=0;k2=0;k3=0;k4=0;fori=1:1:nx1=x0+h;k1=-xO*yOA2;%k1=yO-2*xO/yO;%k2=(-(x0+h/2)*(y0+h/2*k1)A2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);%k3=(-(x0+h/2)*(y0+h/2*k2)A2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);%k4=(-(x1)*(y0+h*k3f2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);%y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);%x0=x1;y0=y1;y=2/(1+x0A2);%y=sqrt(1+2*x0);%fprintf('結(jié)果=%.3f,%.7f,%.7f\n',x1,y1,y);endend五、運行結(jié)果：改進(jìn)歐拉格式結(jié)果：?GJOLh=0?25xO=Oy0=2輸入的范國是：X=0Y=5結(jié)果刊?250,k87600000,1.88235294結(jié)果二0.500,1.59389108,1.60000000結(jié)果二0.750,L2823900S,1-28000000結(jié)果叫.000,1,00962125,L00000000結(jié)果=1.250,0.79318809,0.73048780結(jié)果“?500,0.62915123,0.61538462結(jié)果=1.750,0-50372S54,0.4923076S結(jié)果二N000,0.40966655,0.40000000結(jié)果二2.250,0*33736499,0.3298&691結(jié)果二近500,0.2823574S,0.275S6207結(jié)果W?750,0.23885673,0.23357564結(jié)果=王000,0.20429990,0.20000000結(jié)果250,0.17648S88,0.17297297結(jié)果二500,0.15383629,0.15094340結(jié)果二3.「50,61351747S,0.13278008結(jié)果二4?000,0.11964242,0.11764706結(jié)果二4?250,0.10659158,0.10-191803結(jié)果二4?500,0.09553028.0*09411765結(jié)果匕.為山（L03608040,0.O848SO64結(jié)果二弓.000,0.077948Q7J0.07692308四階龍格庫塔結(jié)果:?LGKTh=0,25i0=0y0=2輸入的萄Hl是：X=0Y=5結(jié)^0,250,1.8823030,1.8823529結(jié)果=0.500,1.5998962,1.6000000結(jié)果=0.750,1.2799478,1.2800000結(jié)果=1.000,L0000271,1-0000000結(jié)果=1.250,0-7805556,0.7804378結(jié)果二1.500j0.6154594j0.6153846Slltl.750,0.4923742：,0.4923077結(jié)果=2.000,0.4000543,0.4000000結(jié)果-2.250,0.3299396,0.3298969結(jié)果=2.500,0.2758952,0.2758621結(jié)果=2.750,0.2336023,0.2335766^^=3,000,0,2000200,0,2000000結(jié)果=3?25Q衛(wèi)17298S6.0,1728730結(jié)果=1500,0.1509558,0-1509434結(jié)果=3-750,0.1327899,0.1327801結(jié)果=4.000,0.1176550,0.1176471結(jié)果=4.250,0.1049245,0.1049180^^4.500,0.0941229,0,0941176結(jié)果=4.750,0.0848850,0.0848806結(jié)果=5.000,0.0769267,0.0769231

?LGKTX1xO=Oy€=2輸入的范園是；1=0Y=5結(jié)果100,1.98019?8,1,9801980結(jié)果=0.200,1.9230763,1.923076&結(jié)果刃*300,1.8348612」.8348624結(jié)果刃.400,L7241364,L724137&結(jié)果二Ck500#1.5999^84,1.6000000結(jié)果=0,60^1.4705870,1*4705802結(jié)果=0.700,1.3422312,L342281&結(jié)果-0.800fL2195122f1.2195122結(jié)果刃*900,1.1049731,L1049724結(jié)果000,1.0000012,1,0000000結(jié)果=1-100,0.9049790,0.9049774結(jié)果=1,200,0,8196739,0.8196721結(jié)果=】*300,0.7434963,0,7434944結(jié)果打*400,0.6756776,0.6756757結(jié)果=1*500,0.6153865,山6153346結(jié)果600,0.5617996,認(rèn)5617978結(jié)果=1.F00?0.5141405,0.5141388結(jié)果=1.800,0,4716997,0,4716981結(jié)果=】*900,0.4338409,0,4338395結(jié)賽2?000,Q.400001%1.4000000結(jié)果=2.100,0.3696870,0-3696858結(jié)果200,0.3424668,仇3424658

結(jié)果啞.300,0.3179660,0.3179650結(jié)果二工400,0.2958589,0.2958580結(jié)果=2.500,0,2*58629,0.2758621結(jié)果弍*600,0.2577327,0.2577320結(jié)果=監(jiān)700,0.2412552,0.^12545結(jié)果=2.800,0.2262449,0.2262443結(jié)果=2.900,0.2125404,0.2125399結(jié)果=J.：：\LI.2000uu57'：.：：'U|jijij0結(jié)果-3.100^0.1835018,0.1885014結(jié)果=3.200,0.1779363,0.1779359結(jié)果=3?300,0.16S2OS9,0.1682086結(jié)果=3.400,0.1592360,0,1592357結(jié)果二3.500」.1609437,0-1509434結(jié)果=3,600,0.143266^0.1432665結(jié)果書*700,0,1361473,0,136U70結(jié)果=3.800,0.1295330,0.1295337結(jié)果=1900,0.1233808,0,1233806結(jié)果=4.000,0.1176472,0.1176471拮UM?M0.1122966,0,1122965結(jié)果X*200,0.1072963,0.1072961結(jié)果=4.300,0.1026169,0.1026167結(jié)果=4.400,0.0982320,(L0982313結(jié)果=4.500,0.09411F8?O.0941176結(jié)果二4*600,0..0902628,0,0902527結(jié)果二4.7?,0.0866177,0,0866176絃肚4?800,0.0831948,0.0831947^$=4.900,0.0799681,0.0799680結(jié)果翡.000,0.0769232,0,0769231步長分別為:0.25和0.1時，不冋結(jié)果顯示驗證了步長減少，對于精度的提高起到很大作用，有效數(shù)字位數(shù)明顯增加。六、 實驗小結(jié)：通過這次實驗學(xué)習(xí)，首先第一點對改進(jìn)歐拉格式和四階龍格庫塔的原理推導(dǎo)有了深入的理解，改進(jìn)歐拉格式采用（預(yù)報+校正）模式得到較精確的原函數(shù)數(shù)值解；而四階龍格庫塔則采用多預(yù)報幾個點的斜率值，采用加權(quán)平均作為平均斜率的近似值的思想達(dá)到更高精度的數(shù)值解，二階龍格庫塔的特例就是改進(jìn)后的歐拉格式。七、 思考題