數(shù)值分析作業(yè) 歐拉 龍格庫塔

1、9.2對初值問題試用歐拉方法取步長h=0.1和h=0.2計算其近似解，并與精確解比較。解: 歐拉方法：h=0.1;a=0;b=1;n=(b-a)/h;1y=zeros(1,n+1);xx=0;yy=0;y(1)=yy+h*(1.0/(1+xx.2)-2*(yy.2);xx=xx+h;yy=y(1);for i=1:ny(i+1)=y(i)+h*(1.0/(1+xx.2)-2*(yy.2);xx=xx+h;yy=y(i+1);end; yy = 0.1000 0.1970 0.2854 0.3609 0.4210 0.4656 0.4957 0.5137 0.5219 0.5227 0.5180

2、當h=0.2時y = 0.2000 0.3763 0.4921 0.5423 0.5466 0.5271精確計算的結(jié)果：xx=0;h=0.1;a=0;b=1;n=(b-a)/h;y=zeros(1,n);for i=1:nxx=xx+h;y(i)=xx/(1+xx.2);end; yy =0.0990 0.1923 0.2752 0.3448 0.4000 0.4412 0.4698 0.4878 0.4972 0.5000誤差的發(fā)展情況: h=0.1;jsy= 0.1970 ;0.2854 ; 0.3609; 0.4210 ; 0.4656 ; 0.4957 ;0.5137 ; 0.5219

3、; 0.5227;0.5180;jqy=0.0990;0.1923;0.2752 ; 0.3448 ; 0.4000 ; 0.4412 ; 0.4698; 0.4878;0.4972 ; 0.5000; x=zeros(10,1); xx=0; x(1)=xx+h;for i=2:10x(i)=x(i-1)+h; end;yc=zeros(10,1);for j=1:10yc(j)=jsy(j)-jqy(j);end; plot(x,yc,b)誤差yc= 0.0980 0.0931 0.0857 0.0762 0.0656 0.0545 0.0439 0.0341 0.0255 0.0180誤

4、差的走勢圖為：9.3利用歐拉預(yù)測-校正法和經(jīng)典龍格庫塔方法，取步長為h=0.1，求解方程并與精確解進行比較。解: 歐拉預(yù)測-矯正方法h=0.1;a=0;b=1;n=(b-a)/h;yy=1;x=zeros(1,n+1);y=zeros(1,n);for i=2:11x(i)=x(i-1)+h;end;yb=yy+h*(x(1)+yy);y(1)=yy+(h/2)*(x(1)+yy+x(2)+yb);for j=2:nyb=y(j-1)+h*(x(j)+y(j-1);y(j)=y(j-1)+(h/2)*(x(j)+y(j-1)+x(j+1)+yb);end; yy =1.1100 1.2421

5、1.3985 1.5818 1.7949 2.0409 2.3231 2.6456 3.0124 3.4282四階經(jīng)典龍格庫塔方法：h=0.1;y0=1;a=0;b=1;n=(b-a)/h;y=zeros(1,n);x=zeros(1:10);k1=h*(0+1); k2=h*(0+h/2+1+k1/2); k3=h*(0+h/2+1+k2/2); k4=h*(0+h+1+k3);y(1)=y0+(k1+2*k2+2*k3+k4)/6; x(1)=0.1;for i=2:10x(i)=x(i-1)+h;end;for j=2:10k1=h*(x(j-1)+y(j-1); k2=h*(x(j-1

6、)+h/2+y(j-1)+k1/2);k3=h*(x(j-1)+h/2+y(j-1)+k2/2); k4=h*(x(j-1)+h+y(j-1)+k3);y(j)=y(j-1)+(k1+2*k2+2*k3+k4)/6;end; yy =1.1103 1.2428 1.3997 1.5836 1.7974 2.0442 2.3275 2.6511 3.0192 3.4366精確計算： x=0;h=0.1;a=0;b=1;n=(b-a)/h;y=zeros(1,n);for i=1:nx=x+h;y(i)=-x-1+2*exp(x);end; yy =1.1103 1.2428 1.3997 1.5

7、836 1.7974 2.0442 2.3275 2.6511 3.0192 3.4366誤差發(fā)展：h=0.1;jsy= 1.1100 1.2421 1.3985 1.5818 1.7949 2.0409 2.3231 2.6456 3.0124 3.4282;jqy=1.1103 1.2428 1.3997 1.5836 1.7974 2.0442 2.3275 2.6511 3.0192 3.4366; x=zeros(10,1); xx=0; x(1)=xx+h;for i=2:10x(i)=x(i-1)+h; end;yc=zeros(10,1);for j=1:10yc(j)=abs

8、(jsy(j)-jqy(j);end; ycplot(x,yc,b)ans =0.0003 0.0007 0.0012 0.0018 0.0025 0.0033 0.0044 0.0055 0.0068 0.0084誤差的走勢圖：9.5用Adams預(yù)測校正法求解并與精確解進行比較。解：當取步長為0.01時:h=0.01;y0=1;a=0;b=1;n=(b-a)/h;y=zeros(1,n);x=zeros(1,n);k1=h*(-1); k2=h*(-(1+k1/2).2); k3=h*(-(1+k2/2).2); k4=h*(-(1+k3).2);y(1)=y0+(k1+2*k2+2*k3+

9、k4)/6; x(1)=0.01;for i=2:nx(i)=x(i-1)+h;end;for j=2:nk1=h*(-y(j-1)*y(i-1); k2=h*(-(y(j-1)+k1/2).2);k3=h*(-(y(j-1)+k2/2).2); k4=h*(-(y(j-1)+k3).2);y(j)=y(j-1)+(k1+2*k2+2*k3+k4)/6;end;yb=zeros(1,100);yb(4)=y(3)+(h/24)*(55*(-y(3)*y(3)-59*(-y(3-1).2)+37*(-y(3-2).2)-9*(-(y0).2);y(4)=y(3)+(h/24)*(9*(-yb(4

10、).2)+19*(-y(3).2)-5*(-y(3-1).2)+(-y(3-2).2);for k=4:99yb(k+1)=y(k)+(h/24)*(55*(-y(k)*y(k)-59*(-y(k-1).2)+37*(-y(k-2).2)-9*(-y(k-3).2);y(k+1)=y(k)+(h/24)*(9*(-yb(k+1).2)+19*(-y(k).2)-5*(-y(k-1).2)+(-y(k-2).2);end;y(:,1:10) %100個校正值的前十個輸出;y(:,91:100) %100個校正值的最后是個輸出;ans = 0.9901 0.9820 0.9740 0.9646 0

11、.9554 0.9464 0.9375 0.9288 0.9202 0.9118ans =0.5245 0.5217 0.5190 0.5163 0.5137 0.5111 0.5085 0.5059 0.5034 0.5008精確解：h=0.01;y0=1;a=0;b=1;n=(b-a)/h;y=zeros(1,n);x=zeros(1,n);x(1)=0.01;for i=2:nx(i)=x(i-1)+h;end; h=0.01;y0=1;a=0;b=1;n=(b-a)/h;y=zeros(1,n);x=zeros(1,n);x(1)=0.01;for i=2:nx(i)=x(i-1)+h;end;for j=1:ny(j)=1/(x(j)+1);end;y(:,1:10) %100個校正值的前十個輸出;y(:,91:1