第三篇 第十章 常微分方程（组）求解
hold off jdwuc1=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)
hold on
S1 = 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]上的精确解')
hold off
jwY2=S1-Y2;xwY2=jwY2./S1;k1=1:n1;k=[0,k1]; P2=[k',x2,Y2,S1,jwY2,xwY2]
运行后屏幕显示分别取n?10,100时，所给的初值问题在[0,1]上的自变量x的值构成的数组Xi (i=1,2)，利用向前欧拉公式（10.8）求出的与Xi (i=1,2)对应的数值解Yi (i=1, 2)，Yi的绝对误差jwYi (i=1,2)和相对误差xwYi (i=1,2)（略），每个子区间[xk,xk?1]上的局部截断误差公式Ren2 =1/200*dy2， 近似解Y与精确解的图形.
例10.3.3 分别用向前欧拉公式(10.8)、一、四阶泰勒逼近及用Matlab函数double求解
?dy??1?ycosx,0?x?2,初值问题?dx分别取n?40,80时，并将计算结果与精确解做比
??y(0)?0.较，写出在每个子区间?xn,xn?1?上的局部截断误差公式，画出数值解与精确解在区间[0，
2]上的图形。
解：(1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=1-y*cos(x); (2)输入程序
>> S1=dsolve('Dy=1-y*cos(x)','y(0)=0','x')
取一阶泰勒逼近
y?e?sinx(1?x?cosx).
同理，取四阶泰勒逼近，输入程序 >> syms X
y=exp(-sin(X))*int(1+sin(u)+((sin(u))^2)/2+((sin(u))^3)/6+((sin(u))^4)/24,u,0,X)
根据运行后，将屏幕显示结果整理得四阶泰勒逼近为
?108111?1017??y?e?sinx??x?cosx??sinx?sin2x?sin3x??.
1896?964???964（3）取n?40时，输入程序
>> syms u
h=2/40,X=0:h:2;
S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S); subplot(4,1,1)
plot(X,S1,'bo'),grid
legend('用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解')
subplot(4,1,2)
y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,'g-'),grid
legend('y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形') subplot(4,1,3)
y1=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X).^2.*cos(X)-1/96*sin(X).^3.*cos(X)+10/9);plot(X,y1,'m-'),grid
legend('四阶泰勒逼近的解在[0,2]上的图形') subplot(4,1,4)
180．
第三篇 第十章 常微分方程（组）求解
x0=0;y0=0;b=2-1.e-4;n=40;tol=1.e-4; [h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol) legend('n=40时，用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解')
syms u
X=0:0.05:2;S=exp(-sin(X))*int(exp(sin(u)),u,0,X); S1=double(S);jwS1y1=y1'-S1',
xwS1y1=jwS1y1./S1',jwYy1=y1'-Y,xwYy1=jwYy1./Y, k1=1:n;k=[0,k1];p=[k',X',Y,y1',jwYy1,xwYy1], pS=[k',X',S1',jwS1y1],xwS1y1
运行后屏幕显示取n=40时，分别用向前欧拉公式(10.8)和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y，y，y1，S1及其图形，Y与y1，S1与y1的绝对误差向量jwYy1，xwS1y1和相对误差向量xwYy1，xwS1y1（略）.向前欧拉公式在每个子区间
?xk,xk?1?上的局部截断误差公式为
Ren?0.0012y??(?k),?k??xk,xk?1?,k?1,2,...,40。
取n?40时，输入程序 >> syms u
X=0:0.05:2;
S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S); y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,'g-'),grid x0=0;y0=0;b=2-1.e-4;n=40;tol=1.e-4;
[h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol)
y4=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X).^2.*cos(X)-1/96*sin(X).^3.*cos(X)+10/9);plot(X,S1,'mo',X,y,'g-.',X,Y,'rp',X,y4,'b-'),grid
legend('用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解','y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形','n=40时，用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解
','y4=exp(-sin(X))(10/9+81/64X-cos(X)(10/9+17/64sin(X)+1/18sin(X)^2+1/96sin(X)^3))在[0,2]上的图形')
运行后屏幕显示取n=40时，分别用向前欧拉公式(10.8)和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y，y，y1，S1及其图形。 取n=40时，向前欧拉公式和四阶泰勒多项式逼近法所求初值问题的数值解的曲线几乎重合，二者的绝对误差和相对误差最小。一阶泰勒多项式逼近法所求初值问题的数值解的曲线与前两条曲线较近。而用函数double所求初值问题的数值解的曲线与前三条曲线很远，且绝对误差和相对误差很大，尤其是相对误差，当x=0时，Y，y1，和y与真值0相等，而S1=0.6912e-093不等于真值0，由此推测，用向前欧拉公式和一、四阶泰勒多项式逼近法所求初值问题的值都接近与真值，是有实际意义的，而用函数double计算的值严重失真，没有实际意义。 下面再取n=80验证。
（4）取n=80时，输入程序 >> syms u
X=0:0.0250:2;
S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S); subplot(4,1,1)
plot(X,S1,'bo'),grid
legend('用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解')
subplot(4,1,2)
y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,'g-'),grid
legend('y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形') subplot(4,1,3)
y1=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X).^2.*cos(X)-1/96*sin(X).^3.*cos(X)+10/9);plot(X,y1,'m-'),grid
legend('四阶泰勒逼近的解在[0,2]上的图形') subplot(4,1,4)
x0=0;y0=0;b=2-1.e-4;n=80;tol=1.e-4;
[h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol)
181．
第三篇 第十章 常微分方程（组）求解
legend('n=80时，用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解')
syms u
X=0:0.0250:2;S=exp(-sin(X))*int(exp(sin(u)),u,0,X); S1=double(S);jwS1y=S1'-y',
xwS1y=jwS1y./S1',jwYy=Y-y',xwYy=jwYy./y', k1=1:n;k=[0,k1];p=[k',X',Y,y',jwYy,xwYy], pS=[k',X',S1',jwS1y],xwS1y
运行后屏幕显示取n=80时，分别用向前欧拉公式(10.8)和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y，y，y1，S1及其图形，Y与y，S1与y的绝对误差向量jwYy，xwS1y和相对误差向量xwYy，xwS1y（略）.向前欧拉公式在每个子区间
?xk,xk?1?上的局部截断误差公式为
Ren?0.0003y??(?k),?k??xk,xk?1?,k?1,2,...,80。
?xk,xk?1?上的局部截断误差公式为
由运行的结果可见，当x=0时，Y，y1，和y与真值0相等，但是S1却随着n的改变而发生变化，但是不等于真值0，由此推测，用向前欧拉公式和一、四阶泰勒多项式逼近法所求初值问题的值都接近与真值，是有实际意义的，而用函数double计算的值严重失真，没有实际意义。下面再取n=100时，向前欧拉公式在每个子区间
Ren?1.9998e?004y??(?k),?k??xk,xk?1?,k?1,2,...,100。
这说明，节点数n越大，向前欧拉公式在每个子区间?xk,xk?1?上的局部截断误差越小，所以函数double计算的值严重失真。
(二) 向前欧拉公式的MATLAB主程序2
用向前欧拉公式（10.8）求解常微分方程初值问题（10.5）的数值解及其截断误差公式的MATLAB主程序2
输入量：funfcn是函数f(x,y)，x0和y0是初值问题y(x0)?y0,b是自变量x的最
大值，h是步长。
输出量：k是自变量x取值的个数，向量X的元素是自变量x的取值，数组Y的元素是利用向前欧拉公式(10.8)求出常微分方程初值问题(10.5)在X的元素处的数值解，REn是在每个子区间[xx,xk?1]上的局部截断误差公式，h是步长，其中dy2是
y(x)的2阶导数y??(x)。
function [k,X,Y,P,REn]=Qeuler2(funfcn,x0,y0,b,h) x=x0; n=fix((b-x)/h); X=zeros(n+1,1); y=y0; Y=zeros(n+1,1); k=1; X(k)=x; Y(k)=y'; for k=2:n+1
X(k)=x+(k-1)*h,
fxy=feval(funfcn,x,y), Y(k)=y+h*fxy
y=Y(k), k=k+1,plot(X,Y,'rp'),
grid,xlabel('自变量 X'), ylabel('因变量 Y') title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]
上的数值解')
end
k1=1:n;k=[0,k1]; P=[k',X,Y]; syms dy2,
REn=0.5*dy2*h^2;
例10.3.4 用向前欧拉公式（10.8）求解初值问题
dyx?y?,y(0)?1，取dx3y 182．
第三篇 第十章 常微分方程（组）求解
n?10,100，并将计算结果与精确解作比较，写出在每个子区间[xk,xk?1]上的局部截断误差公式，画出数值解与精确解在区间[0,2]上的图形.
解 (1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=y-x/(3*y);
(2)建立并保存名为Qeuler2.m的M文件函数。 (3)输入程序
>> S1=dsolve('Dy=y-x/(3*y)','y(0)=1','x') (4)输入程序
>> subplot(2,1,1)
x0=0;y0=1;b=2;n=10;h=2/10;
[k,X,Y,P,REn]=Qeuler2(@funfcn,x0,y0,b,h) hold on
S1=1/6*(6+12*X+30*exp(2*X)).^(1/2); plot(X,S1,'b-')
title('用向前欧拉公式计算dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值
解')
legend('n=10时，dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值解','
dy/dx =y-x/(3y)，y(0)=1在[0,2]上的精确解')
hold off
jdwucY=S1-Y; jwY=S1-Y;
xwY=jwY./Y;k1=1:n;k=[0,k1]; P1=[k',X,Y,S1,jwY,xwY] subplot(2,1,2) n1=100;h1=2/100;
[k,X1,Y1,P1,Ren1]=Qeuler2(@funfcn,x0,y0,b,h1) hold on
S2=1/6*(6+12*X1+30*exp(2*X1)).^(1/2); plot(X1,S2,'b-') legend('n=100时，dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值解','
dy/dx =y-x/(3y)，y(0)=1在[0,2]上的精确解')
hold off
jwY1=S2-Y1;xwY1=jwY1./Y1;k1=1:n1;k=[0,k1]; P2=[k',X1,Y1,S2,jwY1,xwY1]
运行后屏幕显示取n?10,100时，此初值问题在[0,2]上的自变量x的向量X，X1，利用向前欧拉公式（10.8）求出的与X，X1对应的数值解Y，Y1及其绝对误差向量jwY，jwY1和相对误差向量xwY，xwY1（略），每个子区间[xk,xk?1]上的局部截断误差公式Ren2 =1/200*dy2， 数值解Y与精确解的图形.
(三) 自适应向前欧拉公式的MATLAB主程序
用自适应向前欧拉公式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序2
输入量：funfcn是函数f(x,y)，x0和y0是初值问题y(x0)?y0,b是自变量x的最大值，tol是精度。
输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用向前欧拉公式(10.8)求出常微分方程初值问题(10.5)在X的元素处的数值解，向量H是步长h序列，n是自变量x取值的序号。并画出数值解向量Y的图形。
function [H,X,Y,k,h,P]=QEuler(funfcn,x0,b,y0,tol) %初始化. pow=1/3;
if nargin < 5 | isempty(tol), tol = 1.e-6; end;
if nargin < 6 | isempty(trace),
183．
第三篇 第十章 常微分方程（组）求解
trace = 0; end;
x=x0; h=0.0078125*(b-x); y=y0(:);p=128; H=zeros(p,1); X=zeros(p,1); Y=zeros(p,length(y)); k=1; X(k)=x;
Y(k,:)=y';
% 绘图.
clc,x,h,y end
% 主循环.
while (x
h=b-x; end
%计算斜率.
fxy=feval(funfcn,x,y); fxy=fxy(:); %计算误差，设定可接受误差.
delta=norm(h*fxy,'inf');
wucha=tol*max(norm(y,'inf'),1.0);
% 当误差可接受时重写解. if delta<=wucha
x=x+h; y=y+h*fxy; k=k+1; if k>length(X)
X=[X;zeros(p,1)]; Y=[Y;zeros(p,length(y))]; H=[H;zeros(p,1)];
end
H(k)=h;X(k)=x;Y(k,:)=y'; plot(X,Y,'rp'), grid xlabel('自变量 X'), ylabel('因变量 Y') title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0
在[x0,b]上的数值解')
end
%更新步长.
if delta~=0.0
h=min(h*8,0.9*h*(wucha/delta)^pow); end end
if (x disp('Singularity likely.'), x end H=H(1:k); X=X(1:k); Y=Y(1:k,:); n=1:k; P=[n',H,X,Y] 例 10.3.5 用向前欧拉公式（10.8）求解在区间[0,2]上的初值问题 dy??3y?8x?7,y(0)?1，取精度10?1并与精确解作比较，并在同一个坐标系中作出dx图形。 解：(1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=8*x-3*y-7; (2)建立并保存名为Qeuler2.m的M文件函数。 (3)输入程序 >> S1=dsolve('Dy=8*x-3*y-7','y(0)=1','x') (4)输入程序 >> subplot(2,1,1) >> x0=0;y0=1;b=2;tol=1.e-1; [H,X,Y,k,h,P]=Qeuler2(@funfcn,x0,b,y0,tol) 184．