% % %
switch flag case '' %
varargout{1} = f(t,Y,G,ME); % case 'init' %
[varargout{1:3}] = fi(tspan,Y0); % case 'events' %
[varargout{1:3}] = fev(t,Y,Y0); otherwise
error(['Unknown flag ''' flag '''.']); end
% ------------------------------------------------------------------ function Yd = f(t,Y,G,ME) %
X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V; -G*ME*X/r^3];
% ------------------------------------------------------------------ function [ts,y0,options] = fi(tspan,Y0) %
ts=tspan;y0 = Y0;
options = odeset('Events','on','Reltol',1e-5,'Abstol',1e-4); %
% ------------------------------------------------------------------ function [value,isterminal,direction] = fev(t,Y,Y0) %
dDSQdt = 2 * ((Y(1:2)-Y0(1:2))' * Y(3:4)); YSQdt
value = [dDSQdt; dDSQdt]; % direction = [1; -1]; % isterminal = [1; 0]; %
（2）
G=6.672e-11;ME=5.97e24;vy0=4000; x0=-4.2e7;t0=0;tf=60*60*24*9; tspan=[t0,tf];Y0=[x0;0;0;vy0];
[t,YY,Te,Ye,Ie]=ode45('DYDt3',[],[],[],G,ME,tspan,Y0); % <3> X=YY(:,1);Y=YY(:,2);
plot(X,Y,'b','Linewidth',2);hold on text(0,6e7,'轨道','Color','b') axis('image'); %
plot(Ye(1,1),0.4e7+Ye(1,2),'r^','MarkerSize',10) plot(Ye(2,1),0.4e7+Ye(2,2),'bv','MarkerSize',10) plot(Ye(3,1),-0.4e7+Ye(3,2),'b^','MarkerSize',10) %
text(0.8*Ye(3,1),-2e7+Ye(3,2),['t3=' sprintf('%6.0f',Te(3))]) text(0.8*Ye(2,1),1.5e7+Ye(2,2),['t2=' sprintf('%6.0f',Te(2))]) %
[XE,YE,ZE] = sphere(10);RE=0.64e7;XE=RE*XE;YE=RE*YE;ZE=0*ZE; mesh(XE,YE,ZE)
text(1e7,1e7,'地球','Color','r'), hold off
36
图 4.14-3
4.14.3 关于ODE文件的说明
4.14.4 关于解算指令选项options的属性设置 4.14.4.1 options的属性域名
4.14.4.2 options属性处理和输出函数使用演示
【例4.14.4.2-1】仍以卫星轨道问题为例（原题见例4.14.2.1-1, 4.14.2.1-2 , 4.14.2.2-1）。本例演示如何通过对options域的直接设置，借助微分方程解算输出指令，表现解算的中间结果。具体目标是：画出解向量Y??y1y2y3y4?T?xyvxvyT中由x,vx构成
??的相平面。相平面的绘制是在微分方程解算中间逐步完成的。 （1）
[DYDt4.m]
function varargout=DYDt4(t,Y,flag,G,ME,tspan,Y0) switch flag case ''
varargout{1} = f(t,Y,G,ME); case 'init'
[varargout{1:3}] = fi(tspan,Y0); otherwise
error(['Unknown flag ''' flag '''.']); end
% ------------------------------------------------------------------ function Yd = f(t,Y,G,ME)
X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V; -G*ME*X/r^3];
% ------------------------------------------------------------------ function [ts,y0,options] = fi(tspan,Y0)
37
ts=tspan;y0 = Y0; %
options.RelTol=1e-5;options.AbsTol=1e-4; options.OutputFcn='odephas2'; options.OutputSel=[1 3];
（2）
[odeexp4.m] %odeexp4.m
G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9; tspan=[t0,tf];Y0=[x0;0;0;vy0]; %
clf,set(gca,'xlim',[-5 25]*1e7,'ylim',[-3 3]*1e3);% <5> box on hold on;
ode45('DYDt4',[],[],[],G,ME,tspan,Y0);hold off
（3）
shg,odeexp4
图 4.14-4
4.14.5 常微分方程的边值问题解
4.14.5.1 bvp4c求解边值问题的基本思路 4.14.5.2 求解边值问题的基本配套指令
【例4.14.5.2-1】求二阶方程z???c?z?0满足边界条件z(0)?0,z(4)??2的解。在此，
取c=1。本例的目的：（A）完整演示解题步骤；（B）微分方程和边界条件的标准写法；（C）导数函数文件和边界残差函数文件的编写；（D）已知参数的传递；（E）初始猜测网、近似解、插值解的形态。 （1）
（2）
38
[twoode.m]
function dydx=twoode(x,y,c) dydx=[ y(2) , -c*abs(y(1))]';
[twobc.m]
function res = twobc(ya,yb,c) % %
res=[ ya(1), yb(1)+2]';
（3）
sinit=bvpinit(linspace(0,4,5),[1;0]);
（4）
c=1;
sol=bvp4c(@twoode,@twobc,sinit,[ ],c);
%<3>
（5）
x=linspace(0,4,100); y = bvpval(sol,x);
plot(x,y(1,:),'b-',sol.x,sol.y(1,:),'ro',sinit.x,sinit.y(1,:),'ks') legend('\\fontname{隶书}\\fontsize{20}插值后的解曲线','解点','猜测解点',0)
2.521.510.50-0.5-1-1.5-200.51??ó???ú????????1.522.533.54图 4.14-5
4.14.5.3 改善边值问题的求解性能
【例4.14.5.3-1】求解微分方程z???(l?2qcos2x)z?0，在此令q?15。边界条件为
z(0)?1,z?(0)?0,z?(?)?0。本例的目的：（A）演示含未知参数微分方程边值问题的求
解；（B）边界条件的数学表达和M函数文件的编写；（C）函数型猜测解的构成方式；（D）使用bvpset改变选项属性；（E）采用构架直接赋值法改变选项属性。（F）演示“接续”求解法。 （1）
（2）
[mat4ode.m]
function dydx=mat4ode(x,y,lmb)
39
q=5;
dydx=[y(2); -(lmb – 2*q*cos(2*x))*y(1)];
[mat4bc.m]
function res=mat4bc(ya,yb,lmb) res=[ ya(1)-1 ; ya(2) ; yb(2) ];
（3）
clear
x00=linspace(0,pi,4); y00=inline('[cos(4*x);-4*sin(4*x)]'); lmb=15; s0=bvpinit(x00,y00,lmb);
%
% <3> % %
（4）
opts=bvpset('AbsTol',0.5,'RelTol',0.38,'Stats','on');% <6> s1= bvp4c(@mat4ode,@mat4bc,s0,opts); % Lambda1=s1.parameters
The solution was obtained on a mesh of 7 points. The maximum residual is 1.263e-001.
There were 356 calls to the ODE function. There were 67 calls to the BC function. Lambda1 = 18.0803
（5）
opts.AbsTol=1e-6;opts.RelTol=1e-3; s2= bvp4c(@mat4ode,@mat4bc,s1,opts); Lambda2=s2.parameters
The solution was obtained on a mesh of 39 points. The maximum residual is 4.915e-004.
There were 1443 calls to the ODE function. There were 52 calls to the BC function. Lambda2 = 17.0973
（6）
plot(s0.x,s0.y(1,:),'ks--',s1.x,s1.y(1,:),'bo:',s2.x,s2.y(1,:),'r*-') legend('\\fontname{隶书}\\fontsize{16}猜测解','第一近似解','第二近似解',0) axis([0,pi,-1,1.2]),xlabel('x'),ylabel('solution y') 10.80.60.4??????ü?????ü??solution y0.20-0.2-0.4-0.6-0.8-100.511.5x22.53图 4.14-6 40