1、姓名：樊元君 學(xué)號(hào)： 日期：2012.11.06 一、實(shí)驗(yàn)?zāi)康?掌握MATLAB語(yǔ)言、C/C+語(yǔ)言編寫(xiě)計(jì)算程序的方法、掌握改進(jìn)歐拉法與四階龍格-庫(kù)塔求解一階常微分方程的初值問(wèn)題。掌握使用MATLAB程序求解常微分方程問(wèn)題的方法。二、實(shí)驗(yàn)內(nèi)容1、分別寫(xiě)出改進(jìn)歐拉法與四階龍格-庫(kù)塔求解的算法，編寫(xiě)程序上機(jī)調(diào)試出結(jié)果，要求所編程序適用于任何一階常微分方程的數(shù)值解問(wèn)題，即能解決這一類問(wèn)題，而不是某一個(gè)問(wèn)題。 實(shí)驗(yàn)中以下列數(shù)據(jù)驗(yàn)證程序的正確性。 求，步長(zhǎng)h=0.25。2、實(shí)驗(yàn)注意事項(xiàng)的精確解為，通過(guò)調(diào)整步長(zhǎng)，觀察結(jié)果的精度的變化三、程序流程圖： 改進(jìn)歐拉格式流程圖： 四階龍格庫(kù)塔流程圖：四、源程序： 改

2、進(jìn)后歐拉格式程序源代碼：function = GJOL(h,x0,y0,X,Y)format longh=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;for i=1:1:n x1=x0+h;yp=y0+h*(-x0*(y0)2);%yp=y0+h*(y0-2*x0/y0);%yc=y0+h*(-x1*(yp)2);%yc=y0+h*(yp-2*x1/yp);%y1=(yp+yc)/2;x0=x1;y0=y1;y=2/

3、(1+x02);%y=sqrt(1+2*x0);%fprintf(結(jié)果=%.3f,%.8f,%.8fn,x1,y1,y);end end 四階龍格庫(kù)塔程序源代碼：function = LGKT(h,x0,y0,X,Y)format longh=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;for i=1:1:n x1=x0+h;k1=-x0*y02;%k1=y0-2*x0/y0;%k2=(-(x0+

4、h/2)*(y0+h/2*k1)2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);%k3=(-(x0+h/2)*(y0+h/2*k2)2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);%k4=(-(x1)*(y0+h*k3)2);%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+x02);%y=sqrt(1+2*x0);%fprintf(結(jié)果=%.3f,%.7f

5、,%.7fn,x1,y1,y);end end五、運(yùn)行結(jié)果：改進(jìn)歐拉格式結(jié)果：四階龍格庫(kù)塔結(jié)果：步長(zhǎng)分別為：0.25和0.1時(shí)，不同結(jié)果顯示驗(yàn)證了步長(zhǎng)減少，對(duì)于精度的提高起到很大作用，有效數(shù)字位數(shù)明顯增加。六、實(shí)驗(yàn)小結(jié)： 通過(guò)這次實(shí)驗(yàn)學(xué)習(xí)，首先第一點(diǎn)對(duì)改進(jìn)歐拉格式和四階龍格庫(kù)塔的原理推導(dǎo)有了深入的理解，改進(jìn)歐拉格式采用（預(yù)報(bào)+校正）模式得到較精確的原函數(shù)數(shù)值解；而四階龍格庫(kù)塔則采用多預(yù)報(bào)幾個(gè)點(diǎn)的斜率值，采用加權(quán)平均作為平均斜率的近似值的思想達(dá)到更高精度的數(shù)值解，二階龍格庫(kù)塔的特例就是改進(jìn)后的歐拉格式。七、思考題： 如何對(duì)四階龍格-庫(kù)塔法進(jìn)行改進(jìn)，以保證結(jié)果的精度。答：可以通過(guò)計(jì)算結(jié)果的精度處理步長(zhǎng)來(lái)保證結(jié)果的精度。（1）步長(zhǎng)折半。對(duì)于給定