常微分方程的數(shù)值解法(歐拉法、改進(jìn)歐拉法、泰勒方法和龍格庫塔法)

1、頁眉內(nèi)容例1用歐拉方法與改進(jìn)的歐拉方法求初值問題在區(qū)間0,1上取h 0.1的數(shù)值解。解歐拉方法的計(jì)算公式為使用excel表格進(jìn)行運(yùn)算，相應(yīng)結(jié)果如下例一:歐拉法nxy精確解001110.111.00332220.21.0066671.01315930.31.0198241.02914240.41.0390541.05071850.51.0637541.07721760.61.0932111.10793270.71.1266811.14216580.81.1634431.17927490.91.2028451.2186891011.2443141.259921現(xiàn)用matlab編程，程序如下x0=0

2、;y0=1;x(1)=0.1;y(1)=y0+0.1*2*x0/(3*y0A2);for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+0.1*2*x(n)/(3*y(n42);end;xy結(jié)果為x =Columns 1 through 8 0.10000.20000.30000.40000.50000.60000.70000.8000Columns 9 through 10 0.90001.0000y =Columns 1 through 81.00001.00671.01981.03911.06381.09321.12671.1634Columns 9 through

3、 10 1.20281.2443改進(jìn)的歐拉方法其計(jì)算公式為本題的精確解為y(x) 31 x2 ,可用來檢驗(yàn)數(shù)值解的精確度，列出計(jì)算結(jié)果。使用excel表格進(jìn)行運(yùn)算，相應(yīng)如下例一:改進(jìn)的歐拉法nx預(yù)測y校正y精確解0011110.111.0033331.00332220.21.0099561.013181.01315930.31.0261691.0291711.02914240.41.0480541.0507511.05071850.51.0749041.0772521.07721760.61.1059761.1079651.10793270.71.1405491.1421941.1421658

4、0.81.1779651.1792971.17927490.91.2176461.2187061.2186891011.2591031.259931.259921現(xiàn)用matlab編程，程序如下x0=0;y0=1;x(1)=0.1;ya(1)=y0+0.1*2*x0/(3*y0A2);y(1)=y0+0.05*(2*x0/(3*y0A2)+2*x0/(3*yaA2); for n=1:9x(n+1)=0.1*(n+1);ya(n+1)=ya(n)+0.1*2*x(n)/(3*ya(n)A2);y(n+1)=y(n)+0.05*(2*x(n)/(3*y(n)A2)+2*x(n+1)/(3*ya(n

5、+142); end; x y 結(jié)果為 x =Columns 1 through 8 0.10000.20000.30000.40000.50000.60000.70000.8000Columns 9 through 10 0.90001.0000y =Columns 1 through 81.00001.00991.02611.04791.07481.10591.14071.1783Columns 9 through 10 1.21831.2600例2用泰勒方法解分別用二階、四階泰勒方法計(jì)算點(diǎn)xn =0.1, 0.2，1.0處的數(shù)值解，并與精確解進(jìn)行比較。解：二階泰勒方法對于本題故使用exc

6、el表格進(jìn)行運(yùn)算，相應(yīng)結(jié)果如下nxy精確解001110.11.0033331.00332220.21.0132231.01315930.31.0292911.02914240.41.0509691.05071850.51.0775751.07721760.61.1083881.10793270.71.1427041.14216580.81.1798781.17927490.91.219341.2186891011.2606011.259921現(xiàn)用matlab編程，程序如下x0=0;y0=1;x(1)=0.1;y(1)=y0+0.1/(3*y0A2)*(2*x0+0.1*(1-4*x0A2/(3

7、*y0A3);for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+0.1/(3*y(n)A2)*(2*x(n)+0.1*(1-4*x(n)A2/(3*y(n)A3);end; x y 結(jié)果為 x =Columns 1 through 90.10000.20000.30000.40000.50000.60000.70000.80000.9000Column 101.0000y =Columns 1 through 91.00331.01321.02931.05101.07761.10841.14271.17991.2193Column 101.2606四階泰勒方法使用ex

8、cel表格進(jìn)行運(yùn)算，相應(yīng)結(jié)果如下例二：四階泰勒方法nxyy的一階導(dǎo)數(shù)y的二階導(dǎo)數(shù)y的三階導(dǎo)數(shù)y的四階導(dǎo)數(shù)精確解00100.6666670-1.33333110.11.0033280.0662250.653509-0.13402-1.183061.00332220.21.0131910.1298840.616121-0.2711-0.790411.01315930.31.0292110.1888080.560087-0.40991-0.29471.02914240.41.0508230.2414960.49274-0.543790.1597871.05071850.51.0773460.2871

9、890.421266-0.663420.4828091.07721760.61.1080630.3257850.351405-0.760550.6515451.10793270.71.1422740.3576560.286967-0.830570.6901671.14216580.81.1793390.3834610.229962-0.872960.6414121.17927490.91.2186920.4039830.181038-0.89030.546371.2186891011.259850.4200210.13996-0.886940.4355851.259921現(xiàn)用matlab編程，

10、程序如下x0=0;y0=1;ya0=2*x0/(3*y0A2); %一階導(dǎo)數(shù)yb0=2/(3*y0A2)-8*x0A2/(9*y0A5); %二階導(dǎo)數(shù) yc0=-4*x0/(3*y0A5)-80*x0A3/(27*y0A8); %三階導(dǎo)數(shù) yd0=-4/(3*y0A5)+40*x0A2/(3*y0A8)-1280*x0A4/(81*y0A11)；% 四階導(dǎo)數(shù)x(1)=0.1;y(1)=y0+0.1*ya0+0.01/2*yb0+0.001/6*yc0+0.0001/24*yd0;ya(1)=2*x(1)/(3*y(1)A2); %一階導(dǎo)數(shù) yb(1)=2/(3*y(1)A2)-8*x(1)A2

11、/(9*y(1)A5); %二階導(dǎo)數(shù) yc(1)=-4*x(1)/(3*y(1)A5)-80*x(1)A3/(27*y(1)A8); %三階導(dǎo)數(shù) yd(1)=-4/(3*y(1)A5)+40*x(1)A2/(3*y(1)A8)-1280*x(1)A4/(81*yA11)；% 四階導(dǎo)數(shù)for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+0.1*ya(n)+0.01/2*yb(n)+0.001/6*yc(n)+0.0001/24*yd(n);ya(n+1)=2*x(n+1)/(3*y(n+1)A2); %一階導(dǎo)數(shù) yb(n+1)=2/(3*y(n+1)A2)-8*x(n+1

12、)A2/(9*y(n+1)A5);% 二階導(dǎo)數(shù)yc(n+1)=-4*x(n+1)/(3*y(n+1)A5)-80*x(n+1)A3/(27*y(n+1)A8)；%三階導(dǎo)數(shù)yd(n+1)=-4/(3*y(n+1)A5)+40*x(n+1)A2/(3*y(n+1)A8)-1280*x(n+1)A4/(81*y(n+1)A11)；% 四階導(dǎo)數(shù)end; x Y 結(jié)果為 x =Columns 1 through 8 0.10000.20000.30000.40000.50000.60000.70000.8000Columns 9 through 10 0.90001.0000y =Columns 1 t

13、hrough 8 1.00331.01321.02921.05081.07731.10811.14231.1793Columns 9 through 10 1.21871.2598例3用標(biāo)準(zhǔn)四階R- K方法求在區(qū)間0,1上，h0.1的數(shù)值解以及在區(qū)間1,10止，h 1的數(shù)值解，并與精確解進(jìn)行比較。解：對于本題使用excel表格進(jìn)行運(yùn)算，相應(yīng)結(jié)果如下例三:標(biāo)準(zhǔn)四階R K方法nxyk1k2k3k4精確解00100.0033330.0033220.006623110.11.0033220.0066230.0098690.0098370.0129891.00332220.21.0131590.0129

14、890.016030.0159830.0188831.01315930.31.0291430.0188830.0216320.0215750.0241541.02914240.41.0507180.0241540.026560.02650.0287261.05071850.51.0772170.0287260.0307720.0307150.0325861.07721760.61.1079320.0325860.0342860.0342340.0357721.10793270.71.1421650.0357720.0371550.0371110.038351.14216580.81.17927

15、40.038350.0394540.0394170.0403981.17927490.91.2186890.0403990.0412640.0412350.0419971.2186891011.2599210.0419970.0426630.0426410.0432221.259921現(xiàn)用matlab編程，程序如下x0=0;y0=1;k10=2*0.1*x0/(3*y0A2);k20=2*0.1*(x0+0.05)/(3*(y0+k10/2)A2);k30=2*0.1*(x0+0.05)/(3*(y0+k20/2)A2);k40=2*0.1*(x0+0.1)/(3*(y0+k30)A2);x(

16、1)=0.1;y(1)=y0+(k10+2*k20+2*k30+k40)/6;k1(1)=2*0.1*x(1)/(3*y(1)A2);k2(1)=2*0.1*(x(1)+0.05)/(3*(y(1)+k1(1)/2)A2);k3(1)=2*0.1*(x(1)+0.05)/(3*(y(1)+k2(1)/2)A2);k4(1)=2*0.1*(x(1)+0.1)/(3*(y+k3(1)F2);for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+(k1(n)+2*k2(n)+2*k3(n)+k4(n)/6;k1(n+1)=2*0.1*x(n+1)/(3*y(n+1)A2);k2(n+1)=2*0.1*(x(n+1)+0.05)/(3*(y(n+1)+k1(n+1)/2)A2);k3(n+1)=2*0.1*(x(n+1)+0.05)/(3*(y(n+1)+k2(n+1)/2)A2);k4(n+1)=