第三篇 第十章 常微分方程（组）求解
options=odeset('RelTol', 1e-1,'AbsTol',[1e-1,1e-1]); [t,y1]= ode15s(@funfcn,[0 10],[2 3],options) yf=2*exp(-t)+sin(t); zf=2*exp(-t)+cos(t);
plot(t,y1(:,1),'mo',t,y1(:,2),'rp'); hold on
plot(t,yf,'b-',t,zf,'g-.') hold off grid
xlabel('自变量 X'), ylabel('因变量 Y')
legend( '用ode15s函数解刚性方程的y=f(x)数值解','用ode15s函数解
刚性方程的z=g(x)数值解','精确解y=f(x)','精确解z=g(x)')
运行后屏幕显示计算的结果（略），并画出精确解和用ode15s求的数值解的图形.
下面计算时间，编写并保存tode.m 的程序
function tode tic; p1=flops;
options=odeset('RelTol', 1e-1,'AbsTol',[1e-1,1e-1]); [t,y1]= ode45(@funfcn,[0 10],[2 3],options); p2=flops; tode45=toc, pode45=p2-p1, tic;,p3=flops;
options=odeset('RelTol', 1e-1,'AbsTol',[1e-1,1e-1]); [t,y2]= ode15s(@funfcn,[0 10],[2 3],options); p4=flops; tode15s=toc, pode15s=p4-p3 输入命令
>> tode
运行后屏幕显示计算结果
tode45 = tode15s =
4.5100 0.0500
由此可以看出，对于刚性的方程组，用ode15s求数值解比用ode45速度快许多. 习题五：
1，分别用二阶数值方法求下列初值问题
2x?dy?dy?y?,0?x?2,???x?y,0?x?1,y (1)?dx(2)?dx?y(0)?1;??y(0)?1?的数值解，并计算与精确解的误差，画出精确解和数值解的图形。 2，先判别方程组
?dy?dx??3y?2z?sinx,??dz?998y?999z?999(cosx?sinx), ?dx??y(0)?2,?z(0)?3?是否刚性的。再分别用函数和函数求在上的数值解，精确到10，将计算结果与精确值比较，并画出精确解和数值解的图形。 3，用表10-12中的库函数求下列初值问题
?1 205．
第三篇 第十章 常微分方程（组）求解
?dy3y?dy23?,0?x?1,???x?x,0?x?2, (1)?dx1?x(2)?dx???y(0)?1;?y(1)?1的数值解，并计算与精确解的误差，画出精确解和数值解的图形。
4，选择表10-12中合适的库函数解初值问题
?dy?dx??y?z?2sinx,??dz?998y?999z?999(cosx?sinx), ?dx??y(0)?2,?z(0)?3?5，选择表10-12中合适的库函数解计算积分I(x)??0etdt在x?1附近的近似值，取
h?0.2，计算结果保留四位小数。
6，选择表10-12中合适的库函数解下列初值问题
x2??dy2dy223?(1?x)?2xy?4x,0?x?2,??x?x,0?x?2, (1)?(2)?dxdx???y(0)?1;?y(1)?1?6的数值解，要求精确到10，并计算与精确解的误差，画出精确解和数值解的图形。
7，选择表10-12中合适的库函数解初值问题
2ty?dy,0?t?2,??1? 1?t2?dx??y(0)?0.取h?0.5，小数点后至少保留六位。
4x,0?x?1,y8，选择表10-12中合适的库函数求初值问题的数值解。
y(0)?1;h?0.25.?xy?9，选择表10-12中合适的库函数解下列初值问题的数值解。
1??y???y?x?1,0?x?1,?y??(y2?y),1?x?3, (1)?(2)?x?y(0)?1;h?0.1.??y(0)??2;h?0.5.
10.6 线性多步法及其MATLAB程序
10.6.2 亚当斯（Adams）显式公式及其MATLAB程序 取k?4,?1??2??3???1?0，由线性多步法的方程可得：
5559379?0?1,?0?,?1??,?2?,?3??,
24242424代入线性多步法公式，得四阶亚当斯显式公式，即四阶亚当斯外插公式为
yn?1?yn?h?55fn?59fn?1?37fn?2?9fn?3?, (10.50) 24局部截断误差为
206．
第三篇 第十章 常微分方程（组）求解
35h5??(5)5Rn?1??1??(?i)?i?5?(?i)4?i?yn?O(h6)
5!?i?1i??1?2515(5)hyn?O(h).6 (10.51) ?720用四阶亚当斯显式公式求解常微分方程初值问题（10.5）的数值解的MATLAB
主程序
输入量：funfcn是函数f(x,y)，fun是初值问题的精确解函数，x0和y0是初值问
题y(x0)?y0,b是自变量x的最大值，h是步长。
输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用(10.50)式求出常微分方程初值问题(10.5)在X的元素处的数值解，n是自变量x取值的序号，wucha(k+1)=norm(Y(k+1)-Y(k)),并画出数值解和精确解的图形。
function [k,X,Y,wucha,P]= Adams4x(funfcn,x0,b,y0,h) x=x0; y=y0;p=128; n=fix((b-x0)/h); if n<5,return,end; X=zeros(p,1);
Y=zeros(p,length(y)); f=zeros(p,1); k=1; X(k)=x; Y(k,:)=y'; for k=2:4
c1=1/6;c2=2/6;c3=2/6; c4=1/6;a2=1/2; a3=1/2; a4=1;b21=1/2;b31=0;
b32=1/2; b41=0;b42=0;b43=1; x1=x+a2*h; x2=x+a3*h;
x3=x+a4*h; k1=feval(funfcn,x,y); y1=y+b21*h*k1; x=x+h; k2=feval(funfcn,x1,y1); y2=y+b31*h*k1+b32*h*k2; k3=feval(funfcn,x2,y2);
y3=y+b41*h*k1+b42*h*k2+b43*h*k3; k4=feval(funfcn,x3,y3);
y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4); X(k)=x; Y(k,:)=y; end
X;Y;f(1:4)=feval(funfcn,X(1:4),Y(1:4)); for k=4:n
f(k)=feval(funfcn,X(k),Y(k)); X(k+1)=X(1)+h*k;
Y(k+1)=Y(k)+(h/24)*((f(k-3:k))'*[-9 37 -59 55]');
f(k+1)= feval(funfcn,X(k+1),Y(k+1));f(k)=f(k+1); k=k+1; end
for k=2:n+1
wucha(k)=norm(Y(k)-Y(k-1)); k=k+1; end
X=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1,
wucha=wucha(1:n,:); P=[n',X,Y,wucha'];
例10.6.1 先利用常用的四阶龙格-库塔公式求初值问题
2xy?dy,0?x?2,??1? ?dx1?x2??yx?0?0,的几个点的数值解，再利用四阶亚当斯显式公式求解常微分方程初值问题，h=1/15,并计
算它与精确解的误差，在同一图形窗口画出精确解和数值解的图形.
解 （1）编写并保存名为funfcn.m和fun.m的文件如下 function f=funfcn(x,y)
207．
第三篇 第十章 常微分方程（组）求解
f=1-(2.*x.*y)./(1+x.^2); function y=fun(x)
y=(x+1/3.*x.^3)/(1+x.^2);
(2) 在MATLAB工作窗口输入下面的程序
>> y=dsolve('Dy=1-(2*x*y)/(1+x^2)','x') >> x0=0;b=2;y0=0;h=1/15;
[k,X,Y,wucha,P]=Adams4x(@funfcn,x0,b,y0,h),
y=(X+1/3*X.^3)./(1+X.^2); b31=0; b41=0;b42=0;b43=1; c1=1/6;c2=2/6;c3=2/6;c4=1/6;a2=1/2; a3=1/2; a4=1;b21=1/2; b32=1/2;
C=[c1,c2,c3, c4,a2, a3, a4,b21,b31,b32,b41,b42,b43]; [k,X,Y1,fxy,wch,wucha,P]=RK4(@funfcn,@fun,x0,b,C,y0,h) plot(X,Y,'gh',X,Y1,'mp',X,y,'bo'),
grid,xlabel('自变量 X'), ylabel('因变量 Y')
legend('用四阶亚当斯显式公式计算dy/dx=1-(2xy)/(1+x^2)，y(0)=0
在[0,2]上的数值解','用常用的四阶龙格-库塔公式计算dy/dx=1-(2xy)/(1+x^2)，y(0)=0在[0,2]上的数值解
','dy/dx=1-(2xy)/(1+x^2)的精确解y=(x+1/3x^3)/(1+x^2)')
wchY=abs(y-Y), wchY1=abs(y-Y1), m=zeros(1,k),
for n=1:k,m(1,n)=n-1,end, [m',X,y,Y,Y1,wchY,wchY1],
运行后屏幕显示图和计算结果（略）. 由此可见，四阶亚当斯显式公式在前4次迭代的数值解与精确解的绝对误差较小，但是从第5次迭代开始到第30次迭代为止，随着自变量的增大，四阶亚当斯显式公式计算的数值解与精确解的绝对误差逐渐增大，大约在x?1左右又逐渐减少，但是误差始终比Runge?Kutta公式的误差大，尤其是在xn?0.9333时，误差最大。
10.6.3 亚当斯隐式公式及其MATLAB程序 （一）亚当斯隐式公式及其MATLAB程序
取?1??2??3??3?0，由线性多步法的方程可得：
91951?0?1,??1?,?0?,?1??,?2?,
24242424代入线性多步法公式，得四阶亚当斯隐式公式，即四阶亚当斯内插公式为
yn?1?yn?h?9fn?1?19fn?5fn?1?fn?2?, (10.52) 24局部截断误差为
Rn?1??195(5)hyn?O(h6). (10.53) 720
例10.6.2 先利用常用的四阶龙格-库塔公式求初值问题
?dy?x?y,0?x?1,? ?dx??yx?0?0,的几个点的数值解，再利用四阶亚当斯隐式公式求解常微分方程初值问题，h=1/10,并计算它与精确解y?x?1?e的误差，在同一图形窗口画出精确解和数值解的图形.
解 （1） 编写并保存名为Adams4y.m的MATLAB程序
function [k,X,Y,wucha,P]=Adams4y(x0,b,y0,h) x=x0; y=y0;p=128; n=fix((b-x0)/h); if n<5,
208．
?x第三篇 第十章 常微分方程（组）求解
return, end;
X=zeros(p,1);
Y=zeros(p,length(y)); f=zeros(p,1); k=1; X(k)=x; Y(k,:)=y'; for k=2:3
x1=x+h/2; x2=x+h/2; x3=x+h;k1=x-y; y1=y+h*k1/2;
x=x+h;k2=x1-y1;
y2=y+h*k2/2;k3=x2-y2; y3=y+h*k3;k4=x3-y3;
y=y+h*(k1+2*k2+2*k3+k4)/6; X(k)=x;Y(k,:)=y; k=k+1; end X,Y,
for k=3:n
X(k+1)=X(1)+h*k;
Y(k+1)=(1/24.9)*(0.24*k+0.12+(Y(k-2:k))'*[-0.1 0.5
22.1]'),
k=k+1, end
for k=2:n+1
wucha(k)=norm(Y(k)-Y(k-1)); end
X=X(1:n+1);Y=Y(1:n+1,:);n=1:n+1,
wucha=wucha(1:n,:);P=[n',X,Y,wucha']; （2）编写并保存名为funfcn.m和fun.m的M文件如下 function f=funfcn(x,y) f=x-y;
function y=fun(x) y=x-1+exp(-x);
（3）在MATLAB工作窗口输入下面的程序
>> x0=0;b=1;y0=0;h=1/10;
[k,X,Y,wucha,P]= Adams4y (x0,b,y0,h) y=X-1+exp(-X); b31=0; b41=0; b42=0;b43=1; c1=1/6;c2=2/6; c3=2/6;c4=1/6; a2=1/2; a3=1/2; a4=1;b21=1/2; b32=1/2;
C=[c1,c2,c3, c4,a2, a3, a4,b21,b31,b32,b41,b42,b43]; [k,X,Y1,fxy,wch,wucha,P]=RK4(@funfcn,@fun,x0,b,C,y0,h) plot(X,Y,'gh',X,Y1,'mp',X,y,'bo'), grid
xlabel('自变量 X'), ylabel('因变量 Y')
legend('用四阶亚当斯隐式公式计算dy/dx=x-y，y(0)=0在[0,1]上的数
值解','用常用的四阶龙格-库塔公式计算dy/dx=x-y，y(0)=0在[0,1]上的数值解',' dy/dx=x-y，y(0)=0的精确解y=x-1+exp(-x)')
wchY=abs(y-Y),
wchY1=abs(y-Y1), m=zeros(1,k), for n=1:k m(1,n)=n-1 end
[m',X,y,Y,Y1,wchY,wchY1],
运行后屏幕显示图和计算结果（略）.
209．