解初值问题各方种法比较.docVIP

解初值问题各种方法比较 一、实验目的：掌握了解各种解初值问题的方法，体会步长对问题解的影响。 二、实验内容：给定初值问题 其精确解为 三、实验要求： 分别按 （1）欧拉法，步长； （2）改进的欧拉法，步长； （3）四阶标准龙格-库塔法，步长； 编写程序求在节点处的数值解及误差，并比较各方法的优缺点。用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与上述结果进行比较。 四、实验过程： 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; 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; 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,:), 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 [k,X,Y,fxy,wch,wucha,P]=RK4 (funfcn,fun,x0,b,C,y0,h) x=x0; y=y0;p=128; 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); x1=x+a2*h; x2=x+a3*h; x3=x+a4*h; k1=feval(funfcn,x,y); 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*