odefile.m

这本书是matlab经典书籍

function varargout = odefile(t,y,flag,p1,p2)
% odefile.m	是MATLAB提供的模板文件。用户可根据需要，使用任何其他文件名。
% varargout	是MATLAB专门设计的"变长度"输出宗量。它由元胞数组构成，因此
%					可适应任意多个、任意形式的输出宗量。请参看第7.5.2节。
% t					自变量。不管在下面指令中是否出现，自变量t必须作为第一输入宗量。
% y					一阶微分方程组的列向量形式函数。它必须作为第二输入宗量。
% flag			切换变量。它必须处在第三输入宗量位置。用户无须也不要对它直接赋
%					值；它的赋值由微分方程解算指令（solver），自动产生。
% p1, p2		是被传递的参数。这里，作为示意，仅列出两个。
%					根据需要，传递参数的数目不受限制。

% 以下是switch-case构成的多向选择控制。应注意：
%	（A）各情况的情况表达式是MATLAB规定的，不得任意改变。
% 	（B）各情况所执行的任务也是规定的，不得任意改变。
%	（C）各情况内，变长输出宗量元胞数是规定的，不得任意改变。
%	（D）各情况内，（等式右边的）函数名称可以改变，但相应子函数名称要一致。
%	（E）各情况内，（等式右边的）函数中所包含的t , y 是必须的，
%	     不管它们是否以显式出现相应子函数体中。
%	（F）各情况内，（等式右边的）函数中所包含的p1, p2 是传递参数的示意性表示。
%	     具体参数视需要而定，只要该参数已经传入odefile.m 函数内存。

 switch flag
 case ''         	% 规定空字符串 '' 情况：专管一阶导数 dy/dt = f(t,y)的计算
   varargout{1} = f(t,y,p1,p2);
 case 'init'  	% 规定 'init' 情况：专管三个宗量 [tspan , y0 , options]的设置。
   [varargout{1:3}] = init(p1,p2);
 case 'jacobian' 	% 规定 'jacobian' 情况：专管计算解析的 Jacobian 矩阵 df/dy。
   varargout{1} = jacobian(t,y,p1,p2);
 case 'jpattern' 	% 规定 'jpattern' 情况：专管计算稀疏的数值 Jacobian 矩阵 df/dy。  
   varargout{1} = jpattern(t,y,p1,p2);
 case 'mass    ' 	% 规定 'mass' 情况：专管计算质量矩阵（mass matrix）。
   varargout{1} = mass(t,y,p1,p2);
 case 'events'  	% 规定'events'情况：专管事件定义和判断
   [varargout{1:3}] = events(t,y,p1,p2);
 otherwise
   error(['Unknown flag ''' flag '''.']);
 end
 
% 以下是odefile.m 的子函数，它们分别与"父"函数中各情况对应。
%---------------------------------------------------------------------------
  function dydt = f(t,y,p1,p2)
% 空字符串 '' 情况内调用的子函数：专管一阶导数 dy/dt的计算
%	dydt		一阶导数。该变量名可由用户按需要取名。
% 下面函数内容，由用户自己编写，最后产生输出dydt 。
dydt =〈由用户编写〉
% ---------------------------------------------------------------------------
  function [tspan,y0,options] = init(p1,p2)
% 'init' 情况内调用的子函数：专管三个宗量 [tspan , y0 , options]的设置。
% 三个输出宗量的名称可由用户自起，但各位置上宗量的性质不能改变。
% 下面函数内容，由用户自己编写，最后产生输出三个输出 。
tspan =〈积分的时间区段或时间点集〉
 y0 =〈初值〉
 options =〈或用 odeset 设置算法参数，或用空阵符 [ ] 表示缺省设置〉
% ---------------------------------------------------------------------------
  function dfdy = jacobian(t,y,p1,p2)
% 'jacobian' 情况调用的子函数：专管计算解析的 Jacobian 矩阵 df/dy。
% 输出宗量的名称可由用户自起。
% 下面函数内容，由用户自己编写，最后产生输出输出dfdy 。
dfdy =〈 Jacobian 矩阵〉
% ---------------------------------------------------------------------------
  function S = jpattern(t,y,p1,p2)
% 'jpattern' 情况调用的子函数：专管计算解析的 Jacobian 矩阵 df/dy。
% 输出宗量的名称可由用户自起。
% 下面函数内容，由用户自己编写，最后产生输出S 。
 S =〈稀疏形式的数值 Jacobian 矩阵〉
% ---------------------------------------------------------------------------
  function M = mass(t,y,p1,p2)
% 'mass' 情况调用的子函数：专管计算质量矩阵（mass matrix）。解刚性方程时用。
% 输出宗量的名称可由用户自起。
% 下面函数内容，由用户自己编写，最后产生输出M 。
 M =〈质量矩阵〉
 % --------------------------------------------------------------------------
  function [value,isterminal,direction] = events(t,y,p1,p2)
% 'events' 情况内调用的子函数：专管三个宗量 [value,isterminal,direction]的设置。
% 三个输出宗量的名称可由用户自起，但各位置上宗量的性质不能改变。
% 下面函数内容，由用户自己编写，最后产生输出三个输出 。
value =〈各元素是标量事件函数，即函数穿越0点为事件〉
direction =〈各元素定义事件函数穿越方式：1为向正穿越0；-1为向负穿越0；0不管方向〉
isterminal =〈各元素定义事件发生后计算是否继续：1为终止计算；0为继续计算〉
 % --------------------------------------------------------------------------