微分方程数值.docx

1、微分方程数值微分方程数值解 指导教师：李晓峰学 院： 应用科学学院 班 级： 信息与计算科学131801 姓 名： 王慧兰 学 号： 201318030120 目录一、实验名称 错误!未定义书签。二、实验目的 错误!未定义书签。三、实验内容 错误!未定义书签。四、实验要求 错误!未定义书签。五、实验步骤、源程序及运行截图 错误!未定义书签。1. Euler法（向前） 错误!未定义书签。1.1源程序 错误!未定义书签。1.2实验步骤 错误!未定义书签。1.3运行截图 错误!未定义书签。2. Euler法（向后） 错误!未定义书签。2.1源程序 错误!未定义书签。2.2实验步骤 错误!未定义书签。

2、2.3运行截图 错误!未定义书签。3.梯形迭代法公式 错误!未定义书签。3.1源程序 错误!未定义书签。3.2实验步骤 错误!未定义书签。3.3运行截图 错误!未定义书签。4.改进的Euler公式 错误!未定义书签。4.1源程序 104.2实验步骤 错误!未定义书签。4.3运行截图 11一、实验名称Euler法、梯形法求解初值问题。二、实验目的掌握求解常微分方程初值问题的单步法；比较不同算法所得的数值结果，体会步长缩小对问题解得影响。三、实验内容求解初值问题的数值解。四、实验要求编写程序，步长取h=0.1，分别用Euler法、梯形法迭代格式、4阶RK方法求解该问题。整理比较各节点处的数值解、误

3、差、相对误差、体会算法对数值解梯度的改进。（1）改变步长h=0.2和h=0.5，通过比较各节点处数值解，在不同的步长下，算法优劣结论一致。（2）按同一算法不同步长整理数据，比较在同一节点按不同步长计算所得的数值解的误差变化趋势，体会缩小步长对数值解的影响。（3）改进Euler法程序使其适用任意初值问题。五、源程序、实验核心步骤运行截图1、Euler法（向前）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

4、; 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=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; syms dy2,REn=0.5*dy2*h2;1.2实验步骤（1）建立

5、并保存以funfcn.m文件命名的M文件函数function f = 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) hold onS1= 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上的数值解)lege

6、nd(n=100时，dy/dx=8x-3y-7，y(0)=1在0,1上的数值解, dy/dx=8x-3y-7，y(0)=1在0,1上的精确解)hold offjdwuc1=S1-Y1; jwY1=S1-Y1;xwY1=jwY1./S1;k1=1:n;k=0,k1;P1=k,x1,Y1,S1,jwY1,xwY1subplot(2,1,2)n1=10; h2,k2,x2,Y2,P2,Ren2=QEuler1(funfcn,x0,y0,b,n1,tol) hold onS1 = 8/3*x2-29/9+38/9*exp(-3*x2), plot(x2,S1,b-)legend(n=10时，dy/dx

7、=8x-3y-7，y(0)=1在0,1上的数值解, dy/dx=8x-3y-7，y(0)=1在0,1上的精确解)hold offjwY2=S1-Y2;xwY2=jwY2./S1;k1=1:n1;k=0,k1; P2=k,x2,Y2,S1,jwY2,xwY2 1.3运行截图图1-3 =10,100时，在0,1上的数值解和精确解2.向后Euler公式： 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

8、(k,:)=y0; Y1(k,:)=y0; % 绘图.clc,x0,h,y0 % 产生初值.for i=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,:);while Wutol p=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)grid onxlabel(自变量

9、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,Y2.1实验步骤：（1）建立并保存以funfcn.m文件命名的M文件函数function f = 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;

10、X1,Y1,n,P1=Heuler1(funfcn,x0,b,y0,h1,tol) hold onS2= 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上的精确解)hold offjuwY1=S2-Y1; xiwY1=juwY1./Y1; L=P1,S2,juwY1,xiwY1subplot(2,1,2)h=0.05; X,Y,n,P=Heuler1(funfcn,x0,b,y0,h,tol) hold

11、onS1 = 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上的精确解)hold offjuwY=S1-Y; xiwY=juwY./Y; L=P,S1,juwY,xiwY2.3实验结果图：图1010 取h=0.05时，用向后欧拉公式求初值问题的数值解和精确解的图形3、梯形迭代法公式3.1源程序odtixing1.mfunction X,Y,n,P= odtixing1(funfcn,x0,b,y0,h,tol)

12、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(

13、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)grid onxlabel(自变量 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,Y3.2实验步骤（1）建立并保存以funfcn.m文件命名的M文件函数function f = funfcn(x,y)f =

14、 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,b,y0,h,tol)hold onS1=8/3*X-29/9+38/9*exp(-3*X); plot(X,S1,b-), hold offlegend(h=0.05,用梯形公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解,dy/dx=8x-3y-7，y(0)=1在0,2上的精确解)juwYt=S1-Yt; xiwYt=juwYt./Yt; Lt=Pt,S1,

15、juwYt,xiwYt3.3运行截图4、改进的Euler公式4.1源程序： gaiEuler.mfunction H,X,Y,k,h,P=gaiEuler(funfcn,x0,b,y0,tol)%初始化.pow=1/3; if nargin 5 | isempty(tol), tol = 1.e-6; end;if nargin 6 | isempty(trace), trace = 0; end;x=x0; h=0.0078125*(b-x); y=y0(:);p=128; n=fix(b-x0)/h);H=zeros(p,1); X=zeros(p,1); Y=zeros(p,length

16、(y); k=1; X(k)=x; Y(k,:)=y;% 绘图.clc,x,h,y end% 主循环.while (xx) if x+hbh=b-x; end %计算斜率. fxy=feval(funfcn,x,y); fxy=fxy(:);%计算误差，设定可接受误差.delta=norm(h*fxy,inf); wucha=tol*max(norm(y,inf),1.0);% 当误差可接受时重写解.if deltalength(X)X=X;zeros(p,1); Y=Y;zeros(p,length(y);H=H;zeros(p,1); endH(k)=h;X(k)=x;Y(k,:)=y;

17、plot(X,Y,mh), gridxlabel(自变量 X), ylabel(因变量 Y)title(用改进的欧拉公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解)end%更新步长.if delta=0.0h=min(h*8,0.9*h*(wucha/delta)pow);endendif (xb)disp(Singularity likely.), x endH=H(1:k); X=X(1:k); Y=Y(1:k,:); n=1:k; P=n,H,X,Y4.2实验步骤：（1）建立并保存以funfcn.m文件命名的M文件函数function f = funfcn(x,y

18、)f = 8*x-3*y-7;（2）建立并保存以gaiEuler.m文件命名的M文件函数。（3）输入程序 x0=0;y0=1; b=2;tol=1.e-1; Ht,X,Yt,k,h,Pt=odtixing2(funfcn,x0,b,y0,tol)hold onS1 = 8/3*X-29/9+38/9*exp(-3*X),plot(X,S1,b-)hold off hold onH,X,Y,k,h,P=gaiEuler(funfcn,x0,b,y0,tol)hold off legend(用梯形公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解,dy/dx= 8x-3y-7，y(0)=1在0,2上的精确解, 用改进的欧拉公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解),juwY=S1-Y; xiwY=juwY./Y; L=P,S1,juwY,xiwY3.3运行截图