ode15s_old.m

国外经典书籍MULTIVARIABLE FEEDBACK CONTROL-多变量反馈控制 的源码

function [tout,yout,stats] = ode15s(ydot,tspan,y0,options)
%ODE15S	Solve stiff differential equations, variable order method.
%	[T,Y] = ODE15S('ydot',TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates
%	the system of first order differential equations y' = ydot(t,y) from
%	time T0 to TFINAL with initial conditions Y0.  Function ydot(t,y)
%	must return a column vector.  Each row in solution matrix Y
%	corresponds to a time returned in column vector T.  To obtain
%	solutions at the specific times T0, T1, ..., TFINAL (all increasing
%	or all decreasing), use TSPAN = [T0 T1 ... TFINAL].
%	
%	[T,Y] = ODE15S('ydot',TSPAN,Y0,OPTIONS) solves as above with default
%	integration parameters replaced by values in OPTIONS, an argument
%	created with the ODESET function.  See ODESET for details.  Commonly
%	used options are scalar relative error tolerance 'rtol' (1e-3 by
%	default) and vector of absolute error tolerances 'atol' (all
%	components 1e-6 by default).
%	
%	It is possible to specify tspan, y0 and options in ydot.  If TSPAN
%	or YO is empty, or if ODE15S is invoked as ODE15S('ydot'), ODE15S
%	calls [tspan,y0,options] = ydot([],[]) to obtain any values not
%	supplied at the command line.  TYPE CHM6EX to see how this is coded.
%	
%	The Jacobian matrix J(t,y) is critical to the reliability and
%	efficiency of the integration.  If J(t,y) is constant and/or sparse,
%	use the 'constantJ' and/or 'sparseJ' options (see B5EX, BRUSSEX).
%	If ydot is coded so that ydot([t1 t2 ...],[y1 y2 ...]) returns
%	[ydot(t1,y1) ydot(t2,y2) ...], setting 'vectorized' true may speed
%	up the computation of J (see VDPEX).  If an M-file function that
%	evaluates analytically J(t,y) is available, use the 'analyticJ'
%	option (see VDPJAC, BRUSSJAC).
%	
%	As an example, the command
%	
%	    ode15s('vdpex',[0 3000],[2 0]);
%	
%	solves the system y' = vdpex(t,y) with the default relative error
%	tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
%	component.  When called with no output arguments, as in this
%	example, ODE15S calls the default output function ODEPLOT to plot
%	the solution as it is computed.
%	
%	ODE15S also solves problems M(t)*y' = ydot(t,y) with a mass matrix
%	M(t) that is nonsingular and (usually) sparse.  Use the 'mass'
%	option to supply the name of a function of t that returns M(t) (see
%	FEM1EX).  If M(t) is a constant matrix M, use the 'constantM'
%	option, or supply M as the value of 'mass' (see FEM2EX).
%	
%	See also
%	    other ODE solvers:   ODE23S, ODE45, ODE23, ODE113
%	    options handling:    ODESET, ODEGET
%	    output functions:    ODEPLOT, ODEPHAS2, ODEPHAS3
%	    ydot examples:       VDPEX, BRUSSEX, B5EX, CHM6EX
%	    ydot Jacobians:      VDPJAC, BRUSSJAC
%	    Jacobian functions:  NUMJAC, COLGROUP
%	    mass matrices:	 FEM1EX, FEM1MASS, FEM2EX, FEM2MASS

%	ODE15S is a quasi-constant step size implementation in terms of
%	backward differences of the Klopfenstein-Shampine family of
%	Numerical Differentiation Formulas of orders 1-5.  The natural
%	"free" interpolants are used.  Local extrapolation is not done.  By
%	default, Jacobians are generated numerically.  Details are to be
%	found in The MATLAB ODE Suite, L. F. Shampine and M. W. Reichelt,
%	Rept. 94-6, Math. Dept., SMU, Dallas, TX, 1994.

%	Mark W. Reichelt and Lawrence F. Shampine, 6-13-94
%	Copyright (c) 1984-95 by The MathWorks, Inc.

false = 0;
true = ~false;

nsteps = 0; 				% stats
nfailed = 0; 				% stats
nfevals = 0; 				% stats
npds = 0; 				% stats
ndecomps = 0; 				% stats
nsolves = 0; 				% stats

if nargin == 1
  tspan = []; y0 = []; options = [];
elseif nargin == 2
  y0 = []; options = [];
elseif nargin == 3
  options = [];
end

% Get default tspan and y0 from ydot if none are specified.
if (length(tspan) == 0) | (length(y0) == 0)
  [ydot_tspan,ydot_y0,ydot_options] = feval(ydot,[],[]);
  if length(tspan) == 0
    tspan = ydot_tspan;
  end
  if length(y0) == 0
    y0 = ydot_y0;
  end
  if length(options) == 0
    options = ydot_options;
  else
    options = odeset(ydot_options,options);
  end
end

% Test that tspan is internally consistent.
tspan = tspan(:);
ntspan = length(tspan);
if ntspan == 1
  t = 0;
  next = 1;
else
  t = tspan(1);
  next = 2;
end
tfinal = tspan(ntspan);
if t == tfinal
  error('The last entry in tspan must be different from the first entry.');
end
tdir = sign(tfinal - t);
if any(tdir * (tspan(2:ntspan) - tspan(1:ntspan-1)) <= 0)
  error('The entries in tspan must strictly increase or decrease.');
end

y = y0(:);
neq = length(y);
zerosneq = zeros(neq,1);

% Get options, and set defaults.
rtol = odeget(options,'rtol',1e-3);
atol = odeget(options,'atol',1e-6);
atol = atol(:);
if (rtol <= 0) | any(atol <= 0)  
  error('Tolerance rtol and all entries of atol must be positive.');
end
if rtol < 100 * eps 
  rtol = 100 * eps;
  fprintf('Warning: rtol has been increased to %e\n',rtol);
end
if length(atol) == 1
  atol = atol + zerosneq;
elseif length(atol) ~= neq
  error(['Vector atol must be same length as y0 vector (' num2str(neq) ').']);
end
threshold = atol ./ rtol;

userhmax = abs(odeget(options,'hmax',0.1*abs(tfinal-t)));
hmax = min(userhmax, abs(tfinal-t));

hmin = 4 * eps * abs(t);
userhmin = abs(odeget(options,'hmin',hmin));
hmin = min(max(hmin, userhmin), hmax);

stopfun = odeget(options,'stopfun');
if nargout ~= 0
  outfun = odeget(options,'outfun');
else
  outfun = odeget(options,'outfun','odeplot');
end

refine = odeget(options,'refine',1);

printstats = odeget(options,'printstats',false);

vectorized = odeget(options,'vectorized',false);
constantJ = odeget(options,'constantJ',false);
Js = odeget(options,'sparseJ');
analyticJ = odeget(options,'analyticJ');
if length(analyticJ) == 0
  notanalyticJ = true;
else
  notanalyticJ = false;
end

mass = odeget(options,'mass');
if length(mass) ~= 0
  if isstr(mass)
    Mt = feval(mass,t);
    constantM = odeget(options,'constantM',false);
  else
    Mt = mass;
    constantM = odeget(options,'constantM',true);
    if ~constantM
      error('Mass matrix has constant value, and yet constantM is set false.');
    end
  end
  [L,U] = lu(Mt);
else
  Mt = sparse((1:neq)',(1:neq)',1,neq,neq);
  L = Mt;
  U = Mt;
  constantM = true;
end
if constantM
  Mtnew = Mt;
end

maxk = odeget(options,'maxorder',5);
bdf = odeget(options,'bdf',false);

% Initialize the output function.
if length(outfun) == 0
  haveoutfun = false;
else
  haveoutfun = true;
  feval(outfun,[t tfinal],y,'init');
end

% Allocate memory if we're generating output.
if nargout ~= 0
  if (ntspan == 1) & (refine == 0) 	% only 1 output at tfinal
    nout = 0;
  else
    if ntspan > 2 			% output only at tspan points
      tout = zeros(ntspan,1);
      yout = zeros(ntspan,neq);
    else 				% alloc in chunks
      chunk = max(ceil(128 / neq),refine);
      tout = zeros(chunk,1);
      yout = zeros(chunk,neq);
    end
    nout = 1;
    tout(nout) = t;
    yout(nout,:) = y.';
  end
end

% Initialize method parameters.
G = [1; 3/2; 11/6; 25/12; 137/60];
if bdf
  alpha = [0; 0; 0; 0; 0];
else
  alpha = [-37/200; -1/9; -0.0823; -0.0415; 0];
end
invGa = 1 ./ (G .* (1 - alpha));
erconst = alpha .* G + (1 ./ (2:6)');
difU = [ -1, -2, -3, -4,  -5; 		% difU is its own inverse!
          0,  1,  3,  6,  10;
          0,  0, -1, -4, -10;
          0,  0,  0,  1,   5;
          0,  0,  0,  0,  -1 ];
maxK = 1:maxk;
[kJ,kI] = meshgrid(maxK,maxK);
difU = difU(maxK,maxK);
maxit = 4;

% Set the output flag.
if (ntspan > 2) | (refine == 0)
  outflag = 1; 				% output only at tspan points
elseif refine == 1
  outflag = 2; 				% computed points, no refinement
else
  outflag = 3; 				% computed points, with refinement
  ones1r = ones(1,refine-1);
  S = 1 / refine;
  S = cumsum(S(ones1r)) - 1;
  p = cumprod((S(ones(maxk,1),:)+kI(:,1:refine-1)-1) ./ kI(:,1:refine-1));
end

% Compute the partial derivatives dfdt and dfdy.  The error of
% BDF1 is 0.5*h^2*y''(t).  Use this to determine the optimal h.
f0 = feval(ydot,t,y);
[m,n] = size(f0);
if n > 1
  error(['Function ' ydot '(t,y) must return a column vector.'])
elseif m ~= neq
  error(['Vector ' ydot '(t,y) must be same length as initial conditions.']);
end
F0 = U \ (L \ f0);
wt = abs(y) + threshold;
absh = min(hmax, abs(tspan(next) - t));
rh = norm(F0 ./ wt,inf) / (0.8 * sqrt(rtol));
if absh * rh > 1
  absh = 1 / rh;
end
absh = max(absh, hmin);
h = tdir * absh;

tdel = (t + tdir*min(sqrt(eps)*max(abs(t),abs(t+h)),absh)) - t;
f1 = feval(ydot,t+tdel,y);
nfevals = nfevals + 2; 			% stats
dfdt = (f1 - f0) ./ tdel;

if notanalyticJ
  [dfdy,fac,g,nF] = numjac(ydot,t,y,f0,atol,[],vectorized,Js,[]);
  nfevals = nfevals + nF; 		% stats
else
  dfdy = feval(analyticJ,t,y);
end
npds = npds + 1; 			% stats
JMcurrent = true;

absh = min(hmax, abs(tspan(next) - t));
rh = 1.25 * sqrt(0.5*norm((U \ (L \ (dfdt + dfdy*F0))) ./ wt,inf) / rtol);
if absh * rh > 1
  absh = 1 / rh;
end
absh = max(absh, hmin);
h = tdir * absh;

% Initialize.
hnew = absh;
k = 1; 					% start at order 1 with BDF1
K = 1; 					% K = 1:k
newhk = false;

dif = zeros(neq,maxk+2);
dif(:,1) = h * F0;

hinvGak = h * invGa(k);
nconhk = 0; 				% steps taken with current h and k
[L,U] = lu(Mt - hinvGak * dfdy);
ndecomps = ndecomps + 1; 		% stats
havrate = false;

% THE MAIN LOOP

last = false; 				% is it the last step?
while ~last
  
  hmin = max(userhmin, 4 * eps * abs(t));
  if 1.1*hnew >= abs(tfinal - t)
    hnew = abs(tfinal - t);
    newhk = true;
    last = true;
  end

  if newhk
    newhk = false;
    
    difRU = cumprod((kI - 1 - kJ*(hnew/absh)) ./ kI) * difU;
    dif(:,K) = dif(:,K) * difRU(K,K);
    absh = hnew;
    h = tdir * absh;