《解初值问题各种方法比较》由会员分享，可在线阅读，更多相关《解初值问题各种方法比较（5页珍藏版）》请在金锄头文库上搜索。
1、本文格式为Word版，下载可任意编辑解初值问题各种方法比较 解初值问题各种方法对比 一、测验目的：掌管了解各种解初值问题的方法，体会步长对问题解的影 响。 二、测验内容：给定初值问题 ?dy22x?y?xe,1?x?2 ?dxx?y(0)?1其精确解为y?x2(ex?e) 三、测验要求： 分别按 （1）欧拉法，步长h?0.025,h?0.1； （2）提升的欧拉法，步长h?0.05,h?0.1； （3）四阶标准龙格-库塔法，步长h?0.1； 编写程序求在节点xk?1?0.1k(k?1,2,?10)处的数值解及误差，并对比各方法的优缺点。用MATLAB中的内部函数dsolve求此常微分方程初值问题
2、的解并与上述结果举行对比。 四、测验过程： 1、（1）编写主函数。开启Editor编辑器，输入欧拉法法主程序语句： function h,k,X,Y,P=Qeuler1(funfcn,x0,y0,b,n,tol) x=x0; h=(b-x)/n; X=zeros(n,1); y=y0; Y=zeros(n,1); k=1; X(k)=x; Y(k)=y; for k=2:n+1 fxy=feval(funfcn,x,y); delta=norm(h*fxy,inf); wucha=tol*max(norm(y,inf),1.0); if delta=wucha x=x+h; y=y+h*fxy
3、; X(k)=x;Y(k)=y; end plot(X,Y,rp) grid xlabel(自变量 X), ylabel(因变量 Y) title(用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解) end P=X,Y; 以文件名Qeuler1.m保存。 （2）编写主函数。开启Editor编辑器，输入提升的欧拉法法主程序语句： function X,Y,n,P= odtixing1(funfcn,x0,b,y0,h,tol) n=fix(b-x0)/h); X=zeros(n+1,1); Y=zeros(n+1,1); k=1; X(k)=x0;
4、Y(k,:)=y0; Y1(k,:)=y0; % 绘图. clc,x0,h,y0 % 产生初值. for i=2:n+1 X(i)=x0+h; fx0y0=feval(funfcn,x0,y0); Y(i,:)=y0+h*fx0y0; fxiyi=feval(funfcn,X(i),Y(i,:); Y1(i,:)=y0+h*(fxiyi+fx0y0)/2; % 主循环. Wu=abs(Y1(i,:)-Y(i,:); while Wutol p=Y1(i,:), fxip=feval(funfcn,X(i),p); Y1(i,:)=y0+h*(fx0y0+fxip)/2,P1=Y1(i,:),
5、Y(i,:)=p1; end x0=x0+h; y0=Y1(i,:); Y(i,:)=y0; plot(X,Y,ro) grid on xlabel(自变量 X), ylabel(因变量 Y) title(用梯形公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解) end X=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=n,X,Y 以文件名odtixing1.m保存。 （3）编写主函数。开启Editor编辑器，输入四阶标准龙格-库塔法主程序语句： function (funfcn,fun,x0,b,C,y0,h) x=x0; y=y0;p=128;
6、 n=fix(b-x0)/h); fxy=zeros(p,1); wucha=zeros(p,1); wch=zeros(p,1); X=zeros(p,1); Y=zeros(p,length(y); k=1; X(k)=x; Y(k,:)=y; % 绘图. clc,x,h,y %计算 %fxy=fxy(:); for k=2:n+1 x=x+h;a2=C(5); a3=C(6); a4=C(7);b21=C(8); b31=C(9); b32=C(10); b41=C(11); b42=C(12); b43=C(13);c1=C(1); c2=C(2); c3=C(3); c4=C(4);
7、 x1=x+a2*h; x2=x+a3*h; x3=x+a4*h; k1=feval(funfcn,x,y); k,X,Y,fxy,wch,wucha,P=RK4 y1=y+b21*h*k1; k2=feval(funfcn,x1,y1); y2=y+b31*h*k1+b32*h*k2; k3=feval(funfcn,x2,y2); y3=y+b41*h*k1+b42*h*k2+b43*h*k3; k4=feval(funfcn,x3,y3); fxy(k)=feval(fun,x); y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4); X(k)=x; Y(k,:)=y; k=
8、k+1; plot(X,Y,rp,X,fxy,bo), grid xlabel(自变量 X), ylabel(因变量 Y) legend(用四阶龙格-库塔方法计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解,y/dx=f(x,y)，y(x0)=y0的精确解y=f(x) end %计算误差. for k=2:n+1 wucha(k)=norm(Y(k-1)-Y(k); wch(k)=norm(fxy(k)-Y(k); end X=X(1:k); Y=Y(1:k,:);fxy=fxy(1:k,:); n=1:k;wucha=wucha(1:k,:); wch=wch(1:k,:)
9、; P=n,X,Y,fxy,wch,wucha; 以文件名RK4 .m保存。 2、运行程序。 （1）在MATLAB命令窗口输入： x=1:0.05:2; Subplot(2,1,1) x0=1;y0=0;b=2-1.e-4;n=40;tol=1.e-4;h,k,x1,y1P=Qeuler1(funfcn,x0,y0,b,n,tol) hold on s1=x1.2.*(exp(x1)-exp(1);plot(x1,s1,b-) legend(h=0.025时,dy/dx=2/x*y+x2*exp(x),y(1)=0在1,2上的数值解,dy/dx=2/x*y+x2*exp(x),y(1)=0在1
10、,2上的精确解) title(用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解） x=1:0.05:2; Subplot(2,1,1) x0=1;y0=0;b=2-1.e-4;n=10;tol=1.e-4;h,k,x1,y1,P=Qeuler1(funfcn,x0,y0,b,n,tol) hold on s1=x1.2.*(exp(x1)-exp(1);plot(x1,s1,b-) legend(h=0.1时,dy/dx=2/x*y+x2*exp(x),y(1)=0在1,2上的数值解,dy/dx=2/x*y+x2*exp(x),y(1)=0在1,2
11、上的精确解) title(用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解） hold off juwY1=S1-Yt; xiwYt=juwYt./Yt; 回车得到的结果经过整理后如下表1： 表1 取h=0.025时，初值问题在X1处的精确解y、数值解Y1、十足解jwY1和相对误差xwY1 k 自变量X1 数值解Y1 精确解y 十足误差jwY1 相对误差xwY1 NaN 0.0600 0.0597 0.0594 0.0591 0.0588 0.0585 0 1 2 3 4 5 6 1.0000 1.0250 1.0500 1.0750 1.1000 1.1250 1.1500 0 0.0680 0.1445 0.2301 0.3254 0.4311 0.5477 0 0.0723 0.1536 0.2446 0.3459 0.4580 0.5817 0 0.0043 0.0092 0.0145 0.0204 0.0269 0.0340 5
