www.gusucode.com > funfun工具箱matlab源码程序 > funfun/ode15i.m
function varargout = ode15i(ode,tspan,y0,yp0,options,varargin)
%ODE15I Solve fully implicit differential equations, variable order method.
% [TOUT,YOUT] = ODE15I(ODEFUN,TSPAN,Y0,YP0) with TSPAN = [T0 TFINAL]
% integrates the system of differential equations f(t,y,y') = 0 from time
% T0 to TFINAL, with initial conditions Y0,YP0. Function ODE15I solves ODEs
% and DAEs of index 1. The initial conditions must be "consistent", meaning
% that f(T0,Y0,YP0)=0. Use the function DECIC to compute consistent initial
% conditions close to guessed values. ODEFUN is a function handle. For a
% scalar T and vectors Y and YP, ODEFUN(T,Y,YP) must return a column vector
% corresponding to f(t,y,y'). Each row in the solution array YOUT
% corresponds to a time returned in the column vector TOUT. To obtain
% solutions at specific times T0,T1,...,TFINAL (all increasing or all
% decreasing), use TSPAN = [T0 T1 ... TFINAL].
%
% [TOUT,YOUT] = ODE15I(ODEFUN,TSPAN,Y0,YP0,OPTIONS) solves as above with
% default integration properties replaced by values in OPTIONS, an argument
% created with the ODESET function. See ODESET for details. Commonly used
% options are scalar relative error tolerance 'RelTol' (1e-3 by default)
% and vector of absolute error tolerances 'AbsTol' (all components 1e-6
% by default).
%
% The Jacobian matrices df/dy and df/dy' are critical to reliability and
% efficiency. Use ODESET to set 'Jacobian' to a function handle FJAC if
% [DFDY,DFDYP] = FJAC(T,Y,YP) returns the matrices df/dy and df/dy'.
% If DFDY = [], df/dy is approximated by finite differences and similarly
% for DFDYP. If the 'Jacobian' option is not set (the default), both
% matrices are approximated by finite differences.
% Set 'Vectorized' {'on','off'} if the ODE function is coded so that
% ODEFUN(T,[Y1 Y2 ...],YP) returns [ODEFUN(T,Y1,YP) ODEFUN(T,Y2,YP) ...].
% Set 'Vectorized' {'off','on'} if the ODE function is coded so that
% ODEFUN(T,Y,[YP1 YP2 ...]) returns [ODEFUN(T,Y,YP1) ODEFUN(T,Y,YP2) ...].
% If df/dy or df/dy' is a sparse matrix, set 'JPattern' to the sparsity
% patterns, {SPDY,SPDYP}. A sparsity pattern of df/dy is a sparse
% matrix SPDY with SPDY(i,j) = 1 if component i of f(t,y,yp)
% depends on component j of y, and 0 otherwise. Use SPDY = [] to
% indicate that df/dy is a full matrix. Similarly for df/dy' and
% SPDYP. The default value of 'JPattern' is {[],[]}.
%
% [TOUT,YOUT,TE,YE,IE] = ODE15I(ODEFUN,TSPAN,Y0,YP0,OPTIONS) with the 'Events'
% property in OPTIONS set to a function handle EVENTS, solves as above while
% also finding where functions of (T,Y,YP), called event functions, are zero.
% For each function you specify whether the integration is to terminate at
% a zero and whether the direction of the zero crossing matters. These are
% the three column vectors returned by EVENTS: [VALUE,ISTERMINAL,DIRECTION] =
% EVENTS(T,Y,YP). For the I-th event function: VALUE(I) is the value of the
% function, ISTERMINAL(I)=1 if the integration is to terminate at a zero of
% this event function and 0 otherwise. DIRECTION(I)=0 if all zeros are to
% be computed (the default), +1 if only zeros where the event function is
% increasing, and -1 if only zeros where the event function is
% decreasing. Output TE is a column vector of times at which events occur.
% Rows of YE are the corresponding solutions, and indices in vector IE
% specify which event occurred.
%
% SOL = ODE15I(ODEFUN,[T0 TFINAL],Y0,YP0,...) returns a structure that
% can be used with DEVAL to evaluate the solution or its first derivative
% at any point between T0 and TFINAL. The steps chosen by ODE15I are returned
% in a row vector SOL.x. For each I, the column SOL.y(:,I) contains
% the solution at SOL.x(I). If events were detected, SOL.xe is a row vector
% of points at which events occurred. Columns of SOL.ye are the corresponding
% solutions, and indices in vector SOL.ie specify which event occurred.
%
% Example
% t0 = 1;
% y0 = sqrt(3/2);
% yp0 = 0;
% [y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0);
% uses a helper function DECIC to hold fixed the initial value for y(t0)
% and compute a consistent initial value for y'(t0) for the Weissinger
% implicit ODE. The ODE is solved using ODE15I and the numerical solution
% is plotted against the analytical solution
% [t,y] = ode15i(@weissinger,[1 10],y0,yp0);
% ytrue = sqrt(t.^2 + 0.5);
% plot(t,y,t,ytrue,'o');
%
% See also ODE15S, ODE23S, ODE23T, ODE23TB, ODE45, ODE23, ODE113, DECIC,
% ODESET, ODEGET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
% ODEEXAMPLES, IHB1DAE, IBURGERSODE, FUNCTION_HANDLE.
% Jacek Kierzenka and Lawrence F. Shampine
% Copyright 1984-2011 The MathWorks, Inc.
solver_name = 'ode15i';
% Check inputs
if nargin < 5
options = [];
if nargin < 4
error(message('MATLAB:ode15i:NotEnoughInputs'));
end
end
% Stats
nsteps = 0;
nfailed = 0;
nfevals = 0;
npds = 0;
ndecomps = 0;
nsolves = 0;
% Output
FcnHandlesUsed = true; % No MATLAB v. 5 legacy.
output_sol = (FcnHandlesUsed && (nargout==1)); % sol = odeXX(...)
output_ty = (~output_sol && (nargout > 0)); % [t,y,...] = odeXX(...)
% There might be no output requested...
sol = []; kvec = [];
if output_sol
sol.solver = solver_name;
sol.extdata.odefun = ode;
sol.extdata.options = options;
sol.extdata.varargin = varargin;
end
odeArgs = varargin;
% Handle solver arguments -- pass yp0 as first extra parameter.
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, ~, odeFcn, ...
options, threshold, rtol, ~, ~, hmax, htry, htspan] = ...
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, ...
options, [{yp0(:)} varargin]);
nfevals = nfevals + 1;
% Non-negative solution components
nonNegative = ~isempty(odeget(options,'NonNegative',[],'fast'));
if nonNegative
warning(message('MATLAB:ode15i:NonNegativeIgnored'));
end
% Handle the output
if nargout > 0
outputFcn = odeget(options,'OutputFcn',[],'fast');
else
outputFcn = odeget(options,'OutputFcn',@odeplot,'fast');
end
outputArgs = {};
if isempty(outputFcn)
haveOutputFcn = false;
else
haveOutputFcn = true;
outputs = odeget(options,'OutputSel',1:neq,'fast');
outputArgs = varargin;
end
refine = max(1,odeget(options,'Refine',1,'fast'));
if ntspan > 2
outputAt = 'RequestedPoints'; % output only at tspan points
elseif refine <= 1
outputAt = 'SolverSteps'; % computed points, no refinement
else
outputAt = 'RefinedSteps'; % computed points, with refinement
S = (1:refine-1) / refine;
end
printstats = strcmp(odeget(options,'Stats','off','fast'),'on');
% Handle the event function
[haveEventFcn,eventFcn,eventArgs,valt,teout,yeout,ieout] = ...
odeevents(FcnHandlesUsed,odeFcn,t0,y0,options,[{yp0(:)},varargin]);
if haveEventFcn
eventArgs = {eventArgs{2:end}}; % remove yp0
end
t = t0;
y = y0;
yp = yp0; % Assumes consistent initial slope is supplied.
% Initialize the partial derivatives
[Jac,dfdy,dfdyp,Jconstant,dfdy_options,dfdyp_options,nfcn] = ...
ode15ipdinit(odeFcn,t0,y0,yp0,f0,options,varargin);
npds = npds + 1;
nfevals = nfevals + nfcn;
PDscurrent = true;
maxk = odeget(options,'MaxOrder',5,'fast');
% Initialize method parameters for BDFs in Lagrangian form
% and constant step size. Column k corresponds to the formula
% of order k. lcf holds the leading coefficient, cf holds
% the rest.
lcf = [ 1 3/2 11/6 25/12 137/60 ];
cf = [ -1 -2 -3 -4 -5
0 1/2 3/2 3 5
0 0 -1/3 -4/3 -10/3
0 0 0 1/4 5/4
0 0 0 0 -1/5 ];
% derM(:,k) contains coefficients for calculating scaled
% derivative of order k using equally spaced mesh.
derM = [ 1 1 1 1 1 1
-1 -2 -3 -4 -5 -6
0 1 3 6 10 15
0 0 -1 -4 -10 -20
0 0 0 1 5 15
0 0 0 0 -1 -6
0 0 0 0 0 1 ];
maxit = 4;
% Adjust the warnings.
warnoffId = { 'MATLAB:singularMatrix', 'MATLAB:nearlySingularMatrix'};
for i = 1:length(warnoffId)
warnstat(i) = warning('query',warnoffId{i});
warnoff(i) = warnstat(i);
warnoff(i).state = 'off';
end
% hmin is a small number such that t + hmin is clearly different from t in
% the working precision, but with this definition, it is 0 if t = 0.
hmin = 16*eps*abs(t);
if isempty(htry)
% Compute an initial step size h using yp = y'(t0).
wt = max(abs(y),threshold);
rh = 1.25 * norm(yp ./ wt,inf) / sqrt(rtol);
absh = min(hmax, htspan);
if absh * rh > 1
absh = 1 / rh;
end
absh = max(absh, hmin);
else
absh = min(hmax, max(hmin, htry));
end
h = tdir * absh;
% Initialize. Set dummy value for klast to force the
% formation of iteration matrix.
k = 1;
klast = 0;
abshlast = absh;
raised_order = false;
% For j = 1:6
% t_{n+1-j} are stored in mesh(j)
% y_{n+1-j} are stored in meshsol(:,j)
mesh = zeros(1,maxk+2);
mesh(1) = t0;
meshsol = zeros(neq,maxk+2);
meshsol(:,1) = y0;
% Using the initial slope, create fictitious solution at t - h for
% starting the integration.
mesh(2) = t0 - h;
meshsol(:,2) = y0 - h*yp0;
nconh = 1;
% Allocate memory if we're generating output.
nout = 0;
tout = []; yout = [];
if nargout > 0
if output_sol
chunk = min(max(100,50*refine), refine+floor((2^11)/neq));
tout = zeros(1,chunk);
yout = zeros(neq,chunk);
kvec = zeros(1,chunk);
else
if ntspan > 2 % output only at tspan points
tout = zeros(1,ntspan);
yout = zeros(neq,ntspan);
else % alloc in chunks
chunk = min(max(100,50*refine), refine+floor((2^13)/neq));
tout = zeros(1,chunk);
yout = zeros(neq,chunk);
end
end
nout = 1;
tout(nout) = t;
yout(:,nout) = y;
end
% Initialize the output function.
if haveOutputFcn
feval(outputFcn,[t tfinal],y(outputs),'init',outputArgs{:});
end
% THE MAIN LOOP
done = false;
while ~done
hmin = 16*eps*abs(t);
absh = min(hmax, max(hmin, absh));
% Stretch the step if within 10% of tfinal-t.
if 1.1*absh >= abs(tfinal - t)
h = tfinal - t;
absh = abs(h);
done = true;
end
% LOOP FOR ADVANCING ONE STEP.
nfails = 0;
while true % Evaluate the formula.
gotynew = false; % is ynew evaluated yet?
invwt = 1 ./ max(abs(y),threshold);
while ~gotynew
h = tdir * absh;
tnew = t + h;
if done
tnew = tfinal; % Hit end point exactly.
end
h = tnew - t; % Purify h.
if (absh ~= abshlast) || (k ~= klast)
if absh ~= abshlast
nconh = 0;
end
Miter = dfdy + (lcf(k)/h)*dfdyp;
if issparse(Miter)
[L,U,P,Q,R] = lu(Miter);
else
[L,U,p] = lu(Miter,'vector');
end
ndecomps = ndecomps + 1;
havrate = false;
rate = 1; % Dummy value for test.
end
% Predict the solution and its derivative at tnew.
c = weights(mesh(1:k+1),tnew,1);
ynew = meshsol(:,1:k+1) * c(:,1);
ypnew = meshsol(:,1:k+1) * c(:,2);
ypred = ynew;
minnrm = 100*eps*norm(ypred .* invwt,inf);
% Compute local truncation error constant.
erconst = - 1/(k+1);
for j = 2:k
erconst = erconst - ...
cf(j,k)*prod(((t - (j-1)*h) - mesh(1:k+1)) ./ (h * (1:k+1)));
end
erconst = abs(erconst);
% Iterate with simplified Newton method.
tooslow = false;
for iter = 1:maxit
rhs = - feval(odeFcn,tnew,ynew,ypnew,odeArgs{:});
[lastmsg,lastid] = lastwarn('');
warning(warnoff);
if issparse(Miter)
del = Q * (U \ (L \ (P * (R \ rhs))));
else
del = U \ (L \ rhs(p));
end
warning(warnstat);
% If no new warnings or a muted warning, restore previous lastwarn.
[msg,msgid] = lastwarn;
if isempty(msg) || any(strcmp(msgid,warnoffId))
lastwarn(lastmsg,lastid);
end
newnrm = norm(del .* invwt,inf);
ynew = ynew + del;
ypnew = ypnew + (lcf(k)/h) * del;
if iter == 1
if newnrm <= minnrm
gotynew = true;
break;
end
savnrm = newnrm;
else
rate = (newnrm/savnrm)^(1/(iter-1));
havrate = true;
if rate > 0.9
tooslow = true;
break;
end
end
if havrate && ((newnrm * rate/(1 - rate)) <= 0.33*rtol)
gotynew = true;
break;
elseif iter == maxit
tooslow = true;
break;
end
end % end of Newton loop
nfevals = nfevals + iter;
nsolves = nsolves + iter;
if tooslow
nfailed = nfailed + 1;
abshlast = absh;
klast = k;
% Speed up the iteration by forming new linearization or reducing h.
if ~PDscurrent % always current if Jconstant
if isempty(dfdy_options) && isempty(dfdyp_options)
f = []; % No numerical approximations formed
else
f = feval(odeFcn,t,y,yp,odeArgs{:});
nfevals = nfevals + 1;
end
[dfdy,dfdyp,dfdy_options,dfdyp_options,NF] = ...
ode15ipdupdate(Jac,odeFcn,t,y,yp,f,dfdy,dfdyp,dfdy_options,dfdyp_options,odeArgs);
npds = npds + 1;
nfevals = nfevals + NF;
PDscurrent = true;
% Set a dummy value of klast to force formation of iteration matrix.
klast = 0;
elseif absh <= hmin
warning(message('MATLAB:ode15i:IntegrationTolNotMet', sprintf( '%e', t ), sprintf( '%e', hmin )));
solver_output = odefinalize(solver_name, sol,...
outputFcn, outputArgs,...
printstats, [nsteps, nfailed, nfevals,...
npds, ndecomps, nsolves],...
nout, tout, yout,...
haveEventFcn, teout, yeout, ieout,...
{kvec,yp});
if nargout > 0
varargout = solver_output;
end
return;
else
absh = 0.25 * absh;
done = false;
end
end
end % end of while loop for getting ynew
% Using the tentative solution, approximate scaled derivative
% used to estimate the error of the step.
sderkp1 = norm((ynew - ypred) .* invwt,inf) * ...
abs(prod((absh * (1:k+1)) ./ (tnew - mesh(1:k+1))));
erropt = sderkp1 / (k+1); % Error assuming constant step size.
err = sderkp1 * erconst; % Error accounting for irregular mesh.
% Approximate directly derivatives needed to consider lowering the
% order. Multiply by a power of h to get scaled derivative.
kopt = k;
if k > 1
if nconh >= k
sderk = norm(([ynew,meshsol(:,1:k)] * derM(1:k+1,k)) .* invwt,inf);
else
c = weights([tnew,mesh(1:k)],tnew,k);
sderk = norm(([ynew,meshsol(:,1:k)] * c(:,k+1)) .* invwt,inf) * absh^k;
end
if k == 2
if sderk <= 0.5*sderkp1;
kopt = k - 1;
erropt = sderk / k;
end
else
if nconh >= k-1
sderkm1 = norm(([ynew,meshsol(:,1:k-1)] * derM(1:k,k-1)) .* invwt,inf);
else
c = weights([tnew mesh(1:k-1)],tnew,k-1);
sderkm1 = norm(([ynew,meshsol(:,1:k-1)] * c(:,k)) .* invwt,inf) * absh^(k-1);
end
if max(sderkm1,sderk) <= sderkp1
kopt = k - 1;
erropt = sderk / k;
end
end
end
if err > rtol % Failed step
nfailed = nfailed + 1;
if absh <= hmin
warning(message('MATLAB:ode15i:IntegrationTolNotMet', sprintf( '%e', t ), sprintf( '%e', hmin )));
solver_output = odefinalize(solver_name, sol,...
outputFcn, outputArgs,...
printstats, [nsteps, nfailed, nfevals,...
npds, ndecomps, nsolves],...
nout, tout, yout,...
haveEventFcn, teout, yeout, ieout,...
{kvec,yp});
if nargout > 0
varargout = solver_output;
end
return;
end
abshlast = absh;
klast = k;
nfails = nfails + 1;
switch nfails
case 1
absh = absh * min(0.9,max(0.25, 0.9*(0.5*rtol/erropt)^(1/(kopt+1))));
case 2
absh = absh * 0.25;
otherwise
kopt = 1;
absh = absh * 0.25;
end
absh = max(absh,hmin);
if absh < abshlast
done = false;
end
k = kopt;
else % Successful step
break;
end
end % while true
nsteps = nsteps + 1;
if haveEventFcn
events_args = {eventFcn,tnew,ynew,ypnew,mesh(1:k),meshsol(:,1:k),eventArgs{:}};
[te,ye,ie,valt,stop] = odezero(@ntrp15i,@events_aux,events_args,valt,...
t,y,tnew,ynew,t0,mesh(1:k),meshsol(:,1:k));
if ~isempty(te)
if output_sol || (nargout > 2)
teout = [teout, te];
yeout = [yeout, ye];
ieout = [ieout, ie];
end
if stop % Stop on a terminal event.
% Adjust the interpolation data to [t te(end)].
tnew = te(end);
ynew = ye(:,end);
done = true;
end
end
end
if output_sol
nout = nout + 1;
if nout > length(tout)
tout = [tout, zeros(1,chunk)]; % requires chunk >= refine
yout = [yout, zeros(neq,chunk)];
kvec = [kvec, zeros(1,chunk)];
end
tout(nout) = tnew;
yout(:,nout) = ynew;
kvec(nout) = k;
end
if output_ty || haveOutputFcn
switch outputAt
case 'SolverSteps' % computed points, no refinement
nout_new = 1;
tout_new = tnew;
yout_new = ynew;
case 'RefinedSteps' % computed points, with refinement
tref = t + (tnew-t)*S;
nout_new = refine;
tout_new = [tref, tnew];
yout_new = [ntrp15i(tref,[],[],tnew,ynew,mesh(1:k),meshsol(:,1:k)), ynew];
case 'RequestedPoints' % output only at tspan points
nout_new = 0;
tout_new = [];
yout_new = [];
while next <= ntspan
if tdir * (tnew - tspan(next)) < 0
if haveEventFcn && stop % output tstop,ystop
nout_new = nout_new + 1;
tout_new = [tout_new, tnew];
yout_new = [yout_new, ynew];
end
break;
end
nout_new = nout_new + 1;
tout_new = [tout_new, tspan(next)];
if tspan(next) == tnew
yout_new = [yout_new, ynew];
else
yout_new = [yout_new, ntrp15i(tspan(next),[],[],tnew,ynew,...
mesh(1:k),meshsol(:,1:k))];
end
next = next + 1;
end
end
if nout_new > 0
if output_ty
oldnout = nout;
nout = nout + nout_new;
if nout > length(tout)
tout = [tout, zeros(1,chunk)]; % requires chunk >= refine
yout = [yout, zeros(neq,chunk)];
end
idx = oldnout+1:nout;
tout(idx) = tout_new;
yout(:,idx) = yout_new;
end
if haveOutputFcn
stop = feval(outputFcn,tout_new,yout_new(outputs,:),'',outputArgs{:});
if stop
done = true;
end
end
end
end
if done
break
end
% Advance the integration one step.
t = tnew;
y = ynew;
yp = ypnew;
mesh = [t mesh(1:end-1)];
meshsol = [y meshsol(:,1:end-1)];
PDscurrent = Jconstant;
klast = k;
abshlast = absh;
nconh = min(nconh+1,maxk+2);
% Estimate the scaled derivative of order k+2 if
% * at constant step size,
% * have not already decided to reduce the order,
% * not already at the maximum order, and
% * did not raise order on last step.
if (nconh >= k + 2) && ~(kopt < k) && ~(k == maxk) && ~raised_order
sderkp2 = norm((meshsol(:,1:k+3) * derM(1:k+3,k+2)) .* invwt,inf);
if (k > 1) && (sderk <= min(sderkp1,sderkp2))
kopt = k - 1;
erropt = sderk / k;
elseif ((k == 1) && (sderkp2 < 0.5*sderkp1)) || ...
((k > 1) && (sderkp2 < sderkp1))
kopt = k + 1;
erropt = sderkp2 / (k + 2);
end
end
temp = (erropt/(0.5*rtol))^(1/(kopt+1)); % hopt = absh/temp
if temp <= 1/2
absh = absh * 2;
elseif temp > 1
absh = absh * max(0.5,min(0.9,1/temp));
end
raised_order = kopt > k;
k = kopt;
end % while ~done
solver_output = odefinalize(solver_name, sol,...
outputFcn, outputArgs,...
printstats, [nsteps, nfailed, nfevals,...
npds, ndecomps, nsolves],...
nout, tout, yout,...
haveEventFcn, teout, yeout, ieout,...
{kvec,yp});
if nargout > 0
varargout = solver_output;
end
% -----------------------------
function c = weights(x,xi,maxder)
% Compute Lagrangian interpolation coeffients c for the value at xi
% of a polynomial interpolating at distinct nodes x(1),...,x(N) and
% derivatives of the polynomial of orders 0,...,maxder. c(j,d+1)
% is the coefficient of the function value corresponding to x(j) when
% computing the derivative of order d. Note that maxder <= N-1.
%
% This program is based on the Fortran code WEIGHTS1 of B. Fornberg,
% A Practical Guide to to Pseudospectral Methods, Cambridge University
% Press, 1995.
n = length(x) - 1;
c = zeros(n+1,maxder+1);
c(1,1) = 1;
tmp1 = 1;
tmp4 = x(1) - xi;
for i = 1:n
mn = min(i,maxder);
tmp2 = 1;
tmp5 = tmp4;
tmp4 = x(i+1) - xi;
for j = 0:i-1
tmp3 = x(i+1) - x(j+1);
tmp2 = tmp2*tmp3;
if j == i-1
for k = mn:-1:1
c(i+1,k+1) = tmp1*(k*c(i,k) - tmp5*c(i,k+1))/tmp2;
end
c(i+1,1) = - tmp1*tmp5*c(i,1)/tmp2;
end
for k = mn:-1:1
c(j+1,k+1) = (tmp4*c(j+1,k+1) - k*c(j+1,k))/tmp3;
end
c(j+1,1) = tmp4*c(j+1,1)/tmp3;
end
tmp1 = tmp2;
end
% --------------------------------------------------------------------------
function [vtry,isterminal,direction] = events_aux(ttry,ytry,eventfun,tnew,...
ynew,ypnew,mesh_k,meshsol_k,...
varargin)
% Auxiliary function used to detect events.
if ttry == tnew
yptry = ypnew;
else
[~,yptry] = ntrp15i(ttry,[],[],tnew,ynew,mesh_k,meshsol_k);
end
[vtry,isterminal,direction] = eventfun(ttry,ytry,yptry,varargin{:});