资源描述:
《解初值问题各种方法比较.doc》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、解初值问题各种方法比较一、实验目的:掌握了解各种解初值问题的方法,体会步长对问题解的影响.二、实验内容:给定初值问题其精确解为三、实验要求:分别按(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=ze
2、ros(n,1);y=y0;Y=zeros(n,1);k=1;X(k)=x;Y(k)=y';fork=2:n+1fxy=feval(funfcn,x,y);delta=norm(h*fxy,’inf’);wucha=tol*max(norm(y,’inf'),1。0);ifdelta>=wuchax=x+h;y=y+h*fxy;X(k)=x;Y(k)=y’;endplot(X,Y,'rp')gridxlabel('自变量X’),ylabel(’因变量Y')title('用向前欧拉(Euler)公式计算dy/dx=f(x,y),y(x0)=y0在[x0,b]上的数值解')endP=[
3、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%产生初值.fori=2:n+1X(i)=x0+h;fx0y0=feval(funfcn,x0,y0);Y(i,:)=y0+h*fx0y0;fxiyi=feval(funfcn,X(i),
4、Y(i,:));Y1(i,:)=y0+h*(fxiyi+fx0y0)/2;%主循环。Wu=abs(Y1(i,:)-Y(i,:));whileWu〉tolp=Y1(i,:),fxip=feval(funfcn,X(i),p);Y1(i,:)=y0+h*(fx0y0+fxip)/2,P1=Y1(i,:),Y(i,:)=p1;endx0=x0+h;y0=Y1(i,:);Y(i,:)=y0;plot(X,Y,'ro’)gridonxlabel('自变量X'),ylabel(’因变量Y')title(’用梯形公式计算dy/dx=f(x,y),y(x0)=y0在[x0,b]上的数值解’)end
5、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’;%绘
6、图.clc,x,h,y%计算%fxy=fxy(:);fork=2:n+1x=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
7、+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=k+1;plot(X,Y,'rp',X,fxy,'bo'),gridxlabel(’自变量X’),ylabel('因变量Y’)legend('用四阶龙格-库塔方法计算dy/dx=f(x,y),y(x0)=y0在[x0,b]上的数值