dydt3.m

精通matlab

function varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0)
% DYDt3.m 	供主程序调用的ODE函数文件
%					本文件自带三个子函数：f,fi,fev。
%		t, Y		分别是自变量和一阶函数向量，是最基本的输入宗量。
%		flag		第三输入宗量，它专供解算指令（如ode45）作调用通知。
%					在运行中，解算指令会根据需要向flag发不同的字符串。
% varargout 是"变"输出宗量。它由变维的元胞数组构成。每个元胞中可以存放指令
%					所产生的任意形式的数据。
switch lower flag
case ''      	% 必须用空串符。情况为"真"时，计算导数 dY/dt = f(t,Y)。
  varargout{1} = f(t,Y,G,ME);		%输出为一个元胞，容纳f子函数的一个输出Yd
case 'init'  	% 必须用'init'，情况为"真"时，传送计算区间、初值、设置参数。
  [varargout{1:3}] = fi(tspan,Y0);	%输出为三个元胞，容纳fi子函数的三个输出量
case 'events'	% 必须用'events'，情况为"真"时，设置事件性质。
  [varargout{1:3}] = fev(t,Y,Y0);	%输出三个元胞，容纳fev子函数的三个输出量
otherwise
  error(['Unknown flag ''' flag '''.']);
end case
% ------------------------------------------------------------------
function Yd = f(t,Y,G,ME)
%计算导数子函数，被"父"函数DYDt3调用。
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)
%设置时间区间、初值、算法参数子函数，被"父"函数DYDt3调用。
ts=tspan;y0 = Y0;
options = odeset('Events','on','Reltol',1e-5,'Abstol',1e-4);
						%开动ode45的"事件判断"功能，设置相对误差和绝对误差。
% ------------------------------------------------------------------
function [value,isterminal,direction] = fev(t,Y,Y0)
%事件判断子函数，被"父"函数DYDt3调用。
dDSQdt = 2 * ((Y(1:2)-Y0(1:2))' * Y(3:4));
%dDSQdt是离初始点的二乘距离u的时间导数du/dt，而u=(x-x0)^2+(y-y0)^2 。
value = [dDSQdt; dDSQdt];	%定义两个穿越0的事件
direction = [1; -1];	%第一事件：以渐增方式穿越0。第二事件：以渐减方式穿越0
isterminal = [1; 0];	%第一事件发生后，终止计算；而第二事件发生后，继续计算。