微分方程数值

微分方程数值解指导教师：李晓峰学院：应用科学学院班级：信息与计算科学131801姓名：王慧兰学号：201318030120微分方程数值解法实验报告目录一、实验名称......................................错误!未定义书签。二、实验目的......................................错误!未定义书签。三、实验内容......................................错误!未定义书签。四、实验要求......................................错误!未定义书签。五、实验步骤、源程序及运行截图....................错误!未定义书签。1.Euler法（向前）................................错误!未定义书签。1.1源程序..........................................错误!未定义书签。1.2实验步骤........................................错误!未定义书签。1.3运行截图........................................错误!未定义书签。2.Euler法（向后）................................错误!未定义书签。2.1源程序..........................................错误!未定义书签。2.2实验步骤........................................错误!未定义书签。2.3运行截图错误!未定义书签。3.梯形迭代法公式..................................错误!未定义书签。3.1源程序..........................................错误!未定义书签。3.2实验步骤........................................错误!未定义书签。3.3运行截图........................................错误!未定义书签。4.改进的Euler公式................................错误!未定义书签。4.1源程序.........................................................104.2实验步骤........................................错误!未定义书签。4.3运行截图.......................................................11微分方程数值解法实验报告一、实验名称Euler法、梯形法求解初值问题。二、实验目的掌握求解常微分方程初值问题的单步法；比较不同算法所得的数值结果，体会步长缩小对问题解得影响。三、实验内容求解初值问题1)0(,783yxydxdy的数值解。四、实验要求编写程序，步长取h=0.1，分别用Euler法、梯形法迭代格式、4阶RK方法求解该问题。整理比较各节点处的数值解、误差、相对误差、体会算法对数值解梯度的改进。（1）改变步长h=0.2和h=0.5，通过比较各节点处数值解，在不同的步长下，算法优劣结论一致。（2）按同一算法不同步长整理数据，比较在同一节点按不同步长计算所得的数值解的误差变化趋势，体会缩小步长对数值解的影响。（3）改进Euler法程序使其适用任意初值问题。五、源程序、实验核心步骤运行截图1、Euler法（向前）10(,),(0,1,2...)nnnnnyyhfxynxxnh1.1源程序QEuler.mfunction[h,k,X,Y,P,REn]=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';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')gridlabel('自变量X'),ylabel('因变量Y')title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')endP=[X,Y];symsdy2,REn=0.5*dy2*h^2;1.2实验步骤（1）建立并保存以funfcn.m文件命名的M文件函数functionf=funfcn(x,y)f=8*x-3*y-7;（2）建立并保存以QEuler.m文件命名的M文件函数。（3）输入程序：subplot(2,1,1)x0=0;y0=1;b=1-1.e-4;n=100;tol=1.e-4;[h1,k1,x1,Y1,P1,Ren1]=QEuler1(@funfcn,x0,y0,b,n,tol)holdonS1=8/3*x1-29/9+38/9*exp(-3*x1),plot(x1,S1,'b-')title('用向前欧拉公式计算dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解')legend('n=100时，dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解','dy/dx=8x-3y-7，y(0)=1在[0,1]上的精确解')holdoffjdwuc1=S1-Y1;jwY1=S1-Y1;xwY1=jwY1./S1;k1=1:n;k=[0,k1];P1=[k',x1,Y1,S1,jwY1,xwY1]subplot(2,1,2)n1=10;[h2,k2,x2,Y2,P2,Ren2]=QEuler1(@funfcn,x0,y0,b,n1,tol)holdonS1=8/3*x2-29/9+38/9*exp(-3*x2),plot(x2,S1,'b-')legend('n=10时，dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解','dy/dx=8x-3y-7，y(0)=1在[0,1]上的精确解')holdoffjwY2=S1-Y2;xwY2=jwY2./S1;k1=1:n1;k=[0,k1];P2=[k',x2,Y2,S1,jwY2,xwY2]1.3运行截图微分方程数值解法实验报告图1-3n=10,100时，1)0(,783yxydxdy在[0,1]上的数值解和精确解2.向后Euler公式：10(,),(0,1,2...)nnnnnyyhfxynxxnh2.2源程序：Heuler1.mfunction[X,Y,n,P]=Heuler1(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;Y(i,:)=y0+h*feval(funfcn,x0,y0);Y1(i,:)=y0+h*feval(funfcn,X(i),Y(i,:));%主循环.Wu=abs(Y1(i,:)-Y(i,:));whileWutolp=Y1(i,:);Y1(i,:)=y0+h*feval(funfcn,X(i),p);微分方程数值解法实验报告Y(i,:)=p;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]上的数值解')endX=X(1:n+1);Y=Y(1:n+1,:);n=1:n+1;P=[n',X,Y]2.1实验步骤：（1）建立并保存以funfcn.m文件命名的M文件函数functionf=funfcn(x,y)f=8*x-3*y-7;（2）建立并保存以Heuler1.m文件命名的M文件函数。（3）输入程序S1=dsolve('Dy=8*x-3*y-7','y(0)=1','x')x0=0;y0=1;b=2;tol=1.e-1;subplot(2,1,1)h1=0.01;[X1,Y1,n,P1]=Heuler1(@funfcn,x0,b,y0,h1,tol)holdonS2=8/3*X1-29/9+38/9*exp(-3*X1),plot(X1,S2,'b-')legend('h=0.01用向后欧拉公式计算dy/dx=8x-3y-7，y(0)=1在[0,2]上的数值解','dy/dx=8x-3y-7，y(0)=1在[0,2]上的精确解')holdoffjuwY1=S2-Y1;xiwY1=juwY1./Y1;L=[P1,S2,juwY1,xiwY1]subplot(2,1,2)h=0.05;[X,Y,n,P]=Heuler1(@funfcn,x0,b,y0,h,tol)holdonS1=8/3*X-29/9+38/9*exp(-3*X),plot(X,S1,'b-')legend('h=0.05用向后欧拉公式计算dy/dx=8x-3y-7，y(0)=1在[0,2]上的数值解','dy/dx=8x-3y-7，y(0)=1在[0,2]上的精确解')holdoffjuwY=S1-Y;xiwY=juwY./Y;微分方程数值解法实验报告111[(,)(,)],0,1,2...2nnnnnnhyyfxyfxynL=[P,S1,juwY,xiwY]2.3实验结果图：图10–10取h=0.05时，用向后欧拉公式求初值问题的数值解和精确解的图形3、梯形迭代法公式3.1源程序odtixing1.mfunction[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),Y(i,:));Y1(i,:)=y0+h*(fxiyi+fx0y0)/2;%主循环.Wu=abs(Y1(i,:)-Y(i,:));whileWutolp=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]上的数值解')endX=X(1:n+1);Y=Y(1:n+1,:);n=1:n+1;P=[n',X,Y]3.2实验步骤（1）建立并保存以funfcn.m文件命名的M文件函数functionf=funfcn(x,y)f=8*x-3*y-7;（2）建立并保存以QEuler.m文件命名的M文件函数。（3）输入程序：x0=0;y0=1;b=2;tol=0.1;h=0.05;[X,Yt,n,Pt]=odtixing1(@funfcn,x0,