资源预览内容
第1页 / 共11页
第2页 / 共11页
第3页 / 共11页
第4页 / 共11页
第5页 / 共11页
第6页 / 共11页
第7页 / 共11页
第8页 / 共11页
第9页 / 共11页
第10页 / 共11页
亲,该文档总共11页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述
实验报告五题目:常微分方程数值解法摘要:熟悉常微分方程的数值解法的基本原理。掌握Euler法,改进Euler法,后退Euler法,梯形法,四阶Runge-Kutta法,四阶显式Adams法和四阶隐式Adams法等基本算法。原理:Euler法:预测:y*n+1=yn+hf(xn,yn),校正:yn+1=yn+h/2*f(xn,yn)+f(xn+1,y*n+1),这一计算格式亦可以表示为: Yn+1=yn+h/2*f(xn,yn)+f(xn+h,yn+hf(xn,yn),或表示为下列平均化形式: Yp=yn+hf(xn,yn), Yc=yn+hf(xn+1,yp), Yn+1=1/2*(yp+yc).四阶经典Runge-Kutta方法四阶显式Adams四阶隐式Adams习题3。(改进的Euler法,Euler法)F=x2+x-y;a=0;b=1;h=0.1;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;for i=2:n+1 x=X(i-1); y=Y(i-1); Y(i)=Y(i-1)+eval(F)*h;endY1=zeros(1,n+1);Y1(1)=1;for i=2:n+1 x=X(i-1); y=Y1(i-1); ty=Y1(i-1)+eval(F)*h;Y1(i)=Y1(i-1)+h/2*eval(F);x=X(i); y=ty; Y1(i)=Y1(i)+h/2*eval(F);endtemp=;f=dsolve(Dy=x2+x-y,y(0)=0,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 Euler法 Euler预测-校正公式 准确值); disp(X,Y,Y1,df); Untitled5 步长 Euler法 Euler预测-校正公式 准确值 0 1.0000 1.0000 0 0.1000 0.9000 0.9105 0.0052 0.2000 0.8210 0.8410 0.0213 0.3000 0.7629 0.7914 0.0492 0.4000 0.7256 0.7617 0.0897 0.5000 0.7090 0.7521 0.1435 0.6000 0.7131 0.7624 0.2112 0.7000 0.7378 0.7926 0.2934 0.8000 0.7830 0.8429 0.3907 0.9000 0.8487 0.9131 0.5034 1.0000 0.9349 1.0033 0.6321 figure;plot(X,df,k-,X,Y,-r,X,Y1,.-b);grid on;title(Euler法和Euler预测-校正法解常微分方程);legend(准确值,Euler法,Euler预测-校正法);习题6(1)(四阶经典RungeKutta)F=x+y;a=0;b=1;h=0.2;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;for i=1:n x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x; y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;endtemp=;f=dsolve(Dy=x+y,y(0)=1,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 四阶经典R-K法 准确值);disp(X,Y,df);Untitled3 步长 四阶经典R-K法 准确值 0 1.0000 1.0000 0.2000 1.2428 1.2428 0.4000 1.5836 1.5836 0.6000 2.0442 2.0442 0.8000 2.6510 2.6511 1.0000 3.4365 3.4366 figure;plot(X,df,k*,X,Y,-r);grid on;title(四阶经典R-K法解常微分方程);legend(准确值,四阶经典R-K法)(2)四阶经典RungeKutta)F=3*y/(1+x);a=0;b=1;h=0.2;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;for i=1:n x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x; y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;endtemp=;f=dsolve(Dy=3*y/(1+x),y(0)=1,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 四阶经典R-K法 准确值);disp(X,Y,df); Untitled4 步长 四阶经典R-K法 准确值 0 1.0000 1.0000 0.2000 1.7275 1.7280 0.4000 2.7430 2.7440 0.6000 4.0942 4.0960 0.8000 5.8292 5.83201.0000 7.9960 8.0000 figure;plot(X,df,k*,X,Y,-r);grid on;title(四阶经典R-K法解常微分方程);legend(准确值,四阶经典R-K法);习题9(1)(四阶显式Adams法)F=1-y;a=0;b=1;h=0.2;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=0;for i=1:3 x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x; y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;endfor i=4:n x=X(i-3); y=Y(i-3); f1=eval(F); x=X(i-2); y=Y(i-2); f2=eval(F); x=X(i-1); y=Y(i-1); f3=eval(F); x=X(i); y=Y(i); f4=eval(F); Y(i+1)=Y(i)+h*(55*f4-59*f3+37*f2-9*f1)/24;endtemp=;f=dsolve(Dy=1-y,y(0)=0,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 Adams四步四阶显式法 准确值);disp(X,Y,df);figure;plot(X,df,k*,X,Y,-r);grid on;title(Adams四步四阶显式法解常微分方程);legend(准确值,Adams四步四阶显式法); diwu9 步长 Adams四步四阶显式法 准确值 0 0 0 0.2000 0.1813 0.1813 0.4000 0.3297 0.3297 0.6000 0.4512 0.4512 0.8000 0.5506 0.5507 1.0000 0.6320 0.6321(2) (四阶隐式Adams法)F=1-y;a=0;b=1;h=0.2;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=0;for i=1:3 x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x; y=Y(i)+K2/2; K3=h*eval(
收藏 下载该资源
网站客服QQ:2055934822
金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号