数值分析作业8.docx

《数值分析作业8.docx》由会员分享，可在线阅读，更多相关《数值分析作业8.docx（7页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、9-1 使用四阶Adams预测校正算法和经典RK算法。 RK 算法function yrk4=rk4(fxy,a,b,y0,n)format long;h=(b-a)/n;x=a:h:b;yrk4=zeros(n,1);yrk4(1)=y0;for i=1:n-1 k1=fxy(x(i),yrk4(i); k2=fxy(x(i)+h/2,yrk4(i)+h*k1/2); k3=fxy(x(i)+h/2,yrk4(i)+h*k2/2); k4=fxy(x(i)+h,yrk4(i)+h*k3); yrk4(i+1)=yrk4(i)+h*(k1+2*k2+2*k3+k4)/6;endfunction

2、 z=dm(x,y)z=-y+2*cos(x); yrk4=rk4(dm,0,pi,1,10)yrk4 = 1.000000000000000 1.260037091060021 1.396736267543184 1.396715573348376 1.259976385352878 0.999903214083919 0.641953487423769 0.221165565757292 -0.221271082024303 -0.642047809096457精确值x=0:pi/10:pi;y=cos(x)+sin(x) deey = Columns 1 through 4 1.0000

3、00000000000 1.260073510670101 1.396802246667421 1.396802246667421 Columns 5 through 8 1.260073510670101 1.000000000000000 0.642039521920206 0.221231742082474 Columns 9 through 11 -0.221231742082474 -0.642039521920206 -1.0000000000000009-2. function re(a,b,c,d,max)f1=(t,x,y,z)-y-z;f2=(t,x,y,z)x+a*y;f

4、3=(t,x,y,z)(x-c)*z+b;x(1)=d;y(1)=d;z(1)=d;h=0.1;i=1;for t=0:h:max-h k1(1)=f1(t,x(i),y(i),z(i); k1(2)=f2(t,x(i),y(i),z(i); k1(3)=f3(t,x(i),y(i),z(i); k2(1)=f1(t+h,x(i)+h*k1(1),y(i)+h*k1(1),z(i)+h*k1(1); k2(2)=f2(t+h,x(i)+h*k1(2),y(i)+h*k1(2),z(i)+h*k1(2); k2(3)=f3(t+h,x(i)+h*k1(3),y(i)+h*k1(3),z(i)+h

5、*k1(3); x(i+1)=x(i)+h/2*(k1(1)+k2(1); y(i+1)=y(i)+h/2*(k1(2)+k2(2); z(i+1)=z(i)+h/2*(k1(3)+k2(3); i=i+1;endfigure (1)plot3(0:h:max,x,y,-);figure (2)plot(x,y);endre (0.2,0.8,4.5,1,1000)9-3 function de(h,max)f1=(t,x,y,z,m,n,p)m;f2=(t,x,y,z,m,n,p)n;f3=(t,x,y,z,m,n,p)p;f4=(t,x,y,z,m,n,p)2*x*z*m+3*x2*y*t

6、2;f5=(t,x,y,z,m,n,p)exp(y)*n+4*x*z*t2;f6=(t,x,y,z,m,n,p)2*t*p+2*t*exp(x*y);x(1)=3;y(1)=3;z(1)=3;m(1)=2;n(1)=2;p(1)=2;i=1;for t=1:h:max-h k1(1)=f1(t,x(i),y(i),z(i),m(i),n(i),p(i); k1(2)=f2(t,x(i),y(i),z(i),m(i),n(i),p(i); k1(3)=f3(t,x(i),y(i),z(i),m(i),n(i),p(i); k1(4)=f4(t,x(i),y(i),z(i),m(i),n(i),p

7、(i); k1(5)=f5(t,x(i),y(i),z(i),m(i),n(i),p(i); k1(6)=f6(t,x(i),y(i),z(i),m(i),n(i),p(i); k2(1)=f1(t+h,x(i)+h*k1(1),y(i)+h*k1(1),z(i)+h*k1(1),m(i)+h*k1(1),n(i)+h*k1(1),p(i)+h*k1(1); k2(2)=f2(t+h,x(i)+h*k1(2),y(i)+h*k1(2),z(i)+h*k1(2),m(i)+h*k1(2),n(i)+h*k1(2),p(i)+h*k1(2); k2(3)=f3(t+h,x(i)+h*k1(3),y

8、(i)+h*k1(3),z(i)+h*k1(3),m(i)+h*k1(3),n(i)+h*k1(3),p(i)+h*k1(3); k2(4)=f4(t+h,x(i)+h*k1(4),y(i)+h*k1(4),z(i)+h*k1(4),m(i)+h*k1(4),n(i)+h*k1(4),p(i)+h*k1(4); k2(5)=f5(t+h,x(i)+h*k1(5),y(i)+h*k1(5),z(i)+h*k1(5),m(i)+h*k1(5),n(i)+h*k1(5),p(i)+h*k1(5); k2(6)=f6(t+h,x(i)+h*k1(6),y(i)+h*k1(6),z(i)+h*k1(6)

9、,m(i)+h*k1(6),n(i)+h*k1(6),p(i)+h*k1(6); x(i+1)=x(i)+h/2*(k1(1)+k2(1); y(i+1)=y(i)+h/2*(k1(2)+k2(2); z(i+1)=z(i)+h/2*(k1(3)+k2(3); m(i+1)=m(i)+h/2*(k1(4)+k2(4); n(i+1)=n(i)+h/2*(k1(5)+k2(5); p(i+1)=p(i)+h/2*(k1(6)+k2(6); i=i+1;endfigure (1)plot3(x,y,z,-);ii=1;for t=1:h:max fprintf(t=%.6f,x=%.6f,y=%.

10、6f,z=%.6fn,t,x(ii),y(ii),z(ii) ii=ii+1;endendde(1e-5,1.03325)t=1.033080,x=3.280523,y=3.151594,z=45.355426t=1.033090,x=3.280953,y=3.151698,z=45.523123t=1.033100,x=3.281384,y=3.151802,z=45.695588t=1.033110,x=3.281816,y=3.151906,z=45.873292t=1.033120,x=3.282249,y=3.152011,z=46.056803t=1.033130,x=3.2826

11、84,y=3.152115,z=46.246814t=1.033140,x=3.283120,y=3.152220,z=46.444196t=1.033150,x=3.283557,y=3.152324,z=46.650070t=1.033160,x=3.283996,y=3.152429,z=46.865926t=1.033170,x=3.284437,y=3.152534,z=47.093846t=1.033180,x=3.284878,y=3.152639,z=47.336910t=1.033190,x=3.285321,y=3.152744,z=47.600105t=1.033200,x=3.285766,y=3.152849,z=47.892665t=1.033210,x=3.286212,y=3.152954,z=48.236208t=1.033220,x=3.286659,y=3.153059,z=48.713819t=1.033230,x=3.287108,y=3.153165,z=51.049529t=1.033240,x=3.287558,y=3.153270,z=77662871239455918000000.000000t=1.033250,x=3.288010,y=3.153376,z=Inf