《数值分析matlab版第五章.doc》由会员分享,可在线阅读,更多相关《数值分析matlab版第五章.doc(11页珍藏版)》请在金锄头文库上搜索。
1、 实验报告五题目:常微分方程数值解法摘要:熟悉常微分方程的数值解法的基本原理。掌握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+y
2、c).四阶经典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
3、(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.3
4、000 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预测-校正法解常微分方程);le
5、gend(准确值,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
6、(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,-
7、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*
8、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
9、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(
10、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
11、(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(