十常微分方程组求解Word下载.docx

1、的数值解，分别取，并计算误差，画出精确解和数值解的图形.解 编写并保存名为Eulerli1.m的MATLAB计算和画图的主程序如下function P=Eulerli1（x0,y0,b,hn=（b-x0/h。 X=zeros（n,1Y=zeros（n,1 k=1。X（k=x0。 Y（k=y0。for k=1:nX（k+1=X（k+h。Y（k+1=Y（k+h*（X（k-Y（k k=k+1。endy=X-1+2*exp（-X plot（X,Y,mp,X,y,b-gridxlabel（自变量 X, ylabel（因变量 Ytitle（用向前欧拉公式求dy/dx=x-y，y（0=1在0,1上的数值解和

2、精确解y=x-1+2 exp（-xlegend（h=0.075时，dy/dx=x-y，y（0=1在0,1上的数值解精确解y=x-1+2 exp（-xjwY=y-Y。xwY=jwY./y。k1=1:n。k=0,k1。P=k,X,Y,y,jwY,xwY。在MATLAB工作窗口输入下面的程序x0=0。y0=1。b=1。h=0.0750。P=Eulerli1（x0,y0,b,hh1=0.0075。 P1=Eulerli1（x0,y0,b,h1h1=0.0075时，dy/dx=x-y，y（010.3.2 向前欧拉方法的三种MATLAB程序（一 向前欧拉公式的MATLAB主程序1用向前欧拉公式10.8）求

3、解常微分方程初值问题x=x0。 h=（b-x/n。 y=y0。 X（k=x。=yfor k=2:n+1 fxy=feval（funfcn,x,ydelta=norm（h*fxy,inf wucha=tol*max（norm（y,1.0if delta=wuchax=x+h。 y=y+h*fxy。Y（kplot（X,Y,rplabel（用向前欧拉，y（x0=y0在x0,b上的数值解P=X,Y。syms dy2,REn=0.5*dy2*h2。例10.3.2 用向前欧拉公式 x0=0。b=1-1.e-4。n=100。tol=1.e-4。h1,k1,x1,Y1,P1,Ren1=QEuler1（funf

4、cn,x0,y0,b,n,tolhold onS1= 8/3*x1-29/9+38/9*exp（-3*x1, plot（x1,S1,用向前欧拉公式计算dy/dx=8x-3y-7，y（0n=100时，dy/dx=8x-3y-7，y（0 dy/dx=8x-3y-7，y（0=1在0,1上的精确解hold offjdwuc1=S1-Y1。 jwY1=S1-Y1。xwY1=jwY1./S1。P1=k,x1,Y1,S1,jwY1,xwY1subplot（2,1,2n1=10。 h2,k2,x2,Y2,P2,Ren2=QEuler1（funfcn,x0,y0,b,n1,tolS1 = 8/3*x2-29/9

5、+38/9*exp（-3*x2, plot（x2,S1,n=10时，dy/dx=8x-3y-7，y（0jwY2=S1-Y2。xwY2=jwY2./S1。n1。 P2=k,x2,Y2,S1,jwY2,xwY2运行后屏幕显示分别取时，所给的初值问题在上的自变量的值构成的数组Xi （i=1,2，利用向前欧拉公式对应的数值解Yi（i=1, 2，Yi的绝对误差jwYi （i=1,2和相对误差xwYi （i=1,2 向前欧拉公式的MATLAB主程序210.5）的数值解及其截断误差公式的MATLAB主程序2function k,X,Y,P,REn=Qeuler2（funfcn,x0,y0,b,h n=fix

6、（b-x/hX=zeros（n+1,1Y=zeros（n+1,1=x+（k-1*h,fxy=feval（funfcn,x,y, =y+h*fxyy=Y（k, k=k+1,plot（X,Y,grid,xlabel（ P=k,X,Y。例10.3.4 用向前欧拉公式S1=1/6*（6+12*X+30*exp（2*X.（1/2 plot（X,S1,用向前欧拉公式计算dy/dx=y-x/（3y，y（0=1在0,2上的数值解n=10时，dy/dx=y-x/（3y dy/dx=y-x/（3y=1在0,2上的精确解jdwucY=S1-Y。 jwY=S1-Y。xwY=jwY./Y。,X,Y,S1,jwY,xwY

7、n1=100。h1=2/100。k,X1,Y1,P1,Ren1=Qeuler2（funfcn,x0,y0,b,h1S2=1/6*（6+12*X1+30*exp（2*X1 plot（X1,S2,n=100时，dy/dx=y-x/（3yjwY1=S2-Y1。xwY1=jwY1./Y1。P2=k,X1,Y1,S2,jwY1,xwY1运行后屏幕显示取时，此初值问题在的向量X，X1，利用向前欧拉公式10.8）求出的与X，X1对应的数值解Y，Y1及其绝对误差向量jwY，jwY1和相对误差向量xwY，xwY1 自适应向前欧拉公式的MATLAB主程序用自适应向前欧拉公式求解常微分方程初值问题%初始化.pow=

8、1/3。if nargin tol = 1.e-6。end。 6 | isempty（tracetrace = 0。 h=0.0078125*（b-x y=y0（:p=128。 H=zeros（p,1X=zeros（p,1 Y=zeros（p,length（y Y（k,:% 绘图.clc,x,h,y % 主循环.while （x&（x+hx if x+hbh=b-x。 end %计算斜率. fxy=fxy（:%计算误差，设定可接受误差.% 当误差可接受时重写解.if deltalength（XX=X。zeros（p,1。 Y=Y。zeros（p,length（yH=H。H（k=h。Y（k,:,

9、 grid%更新步长.if delta=0.0h=min（h*8,0.9*h*（wucha/deltapowif （xX=X（1:Y=Y（1:k,:n=1:k。 P=n,H,X,Y10.3.4 向后欧拉方法的MATLAB程序用向后欧拉方法求解常微分方程初值问题n=fix（b-x0 X=zeros（n+1,1 Y=zeros（n+1,1k=1。 Y1（k,:clc,x0,h,y0 % 产生初值.for i=2:X（i=x0+h。 Y（i,:=y0+h*feval（funfcn,x0,y0Y1（i,:=y0+h*feval（funfcn,X（i,Y（i,:Wu=abs（Y1（i,:-Y（i,:wh

10、ile Wutol p=Y1（i,:,pY（i,:=p。x0=x0+h。 y0=Y1（i,:rogrid on用向后欧拉公式计算dy/dx=f（x,yn+1 Y=Y（1:n+1,: n=1:n+1。,X,Y例10.3.6 用向后欧拉公式求解区间上的初值问题的数值解，取步长，并与精确解作比较，在同一个坐标系中作出图形.然后再取，观察数值解与精确解误差的变化，说明与误差的关系.解输入程序S1=dsolve（Dy=8*x-3*y-7y（0=1xtol=1.e-1。h1=0.01。X1,Y1,n,P1=Heuler1（funfcn,x0,b,y0,h1,tolS2= 8/3*X1-29/9+38/9*

11、exp（-3*X1, plot（X1,S2,h=0.01用向后欧拉公式计算dy/dx=8x-3y-7，y（0dy/dx=8x-3y-7，y（0juwY1=S2-Y1。xiwY1=juwY1./Y1。L=P1,S2,juwY1,xiwY1h=0.05。X,Y,n,P=Heuler1（funfcn,x0,b,y0,h,tolS1 = 8/3*X-29/9+38/9*exp（-3*Xplot（X,S1,h=0.05用向后欧拉公式计算dy/dx=8x-3y-7，y（0juwY=S1-Y。xiwY=juwY./Y。L=P,S1,juwY,xiwY运行后屏幕显示用向后欧拉公式计算此初值问题在0,2上的自变

12、量X处数值解Y和精确解S1及其图形，步长H, Y的相对误差xiwY和绝对误差juwY（略 .10.4 改进的欧拉方法的MATLAB程序10.4.2 梯形公式的两种MATLAB程序（一） 梯形公式的MATLAB程序用梯形公式求解常微分方程初值问题n+1 fx0y0=feval（funfcn,x0,y0=y0+h*fx0y0。fxiyi=feval（funfcn,X（i=y0+h*（fxiyi+fx0y0/2。fxip=feval（funfcn,X（i=y0+h*（fx0y0+fxip/2,P1=Y1（i,:, Y（i,:=p1。用梯形公式计算dy/dx=f（x,y例10.4.1 用梯形公式求解区

13、间，取步长，精度为，并与精确解作比较，在同一个坐标系中画出图形. b=2。 tol=0.1。 h=0.05。X,Yt,n,Pt=odtixing1（funfcn,x0,b,y0,h,tolS1=8/3*X-29/9+38/9*exp（-3*X, hold offh=0.05,用梯形公式计算dy/dx=8x-3y-7，y（0juwYt=S1-Yt。 xiwYt=juwYt./Yt。 Lt=Pt,S1,juwYt,xiwYt运行后屏幕显示取精度为，分别用梯形公式和向前欧拉公式求解此初值问题在区间上的自变量X处数值解Yi（i=t,q和精确解S1，步长H, Yi的相对误差xiwYi和绝对误差juwYi

14、（略 及其数值解和精确解的图形. 自适应梯形公式的MATLAB程序用自适应梯形公式求解常微分方程初值问题 n=fix（b-x0if x+h%计算斜率. y1=y+h*fxy。 fxy1=feval（funfcn,x,y1y2=y+h*fxy1。 y=（y1+y2go用自适应梯形公式计算dy/dx=f（x,y例10.4.2 分别用自适应梯形公式和向前欧拉公式分别求解区间，取精度为，并与精确解作比较，在同一个坐标系中作出图形. tol=1.e-1。Ht,X,Yt,k,h,Pt= odtixing2（funfcn,x0,b,y0,tol,hold on,Hq,X,Yq,k,h,Pq=QEuler（funfcn,x0,b,y0,tol用自适应梯形公式计算dy/dx=8x-3y-7，y（0dy/dx= 8x-3y-7， y（0用自适应梯形公式计算dy/dx=8x-3y-7,y（0和精确解S1及其图形，步长H等 （略.10.4.4 改进的欧拉公式的MATLAB程序用自适应改进的欧拉公式求解常微分方程初值问题 trace = 0。y=y0（:H=zeros（p,1 X=zeros（p,1Y=