odecfl3.m

level set matlab code

function [ t, y, schemeData ] = ...
                            odeCFL3(schemeFunc, tspan, y0, options, schemeData)
% odeCFL3: integrate a CFL constrained ODE (eg a PDE by method of lines).
%
% [ t, y, schemeData ] = odeCFL3(schemeFunc, tspan, y0, options, schemeData)
%
% Integrates a system forward in time by CFL constrained timesteps
%   using a third order Total Variation Diminishing (TVD) Runge-Kutta
%   (RK) scheme.  Details can be found in O&F chapter 3.
%
% parameters:
%   schemeFunc	 Function handle to a CFL constrained ODE system
%                  (typically an approximation to an HJ term, see below).
%   tspan        Range of time over which to integrate (see below).
%   y0           Initial condition vector 
%                  (typically the data array in vector form).
%   options      An option structure generated by odeCFLset 
%                  (use [] as a placeholder if necessary).
%   schemeData   Structure passed through to schemeFunc.
%
%
%   t            Output time(s) (see below).
%   y            Output state (see below).
%   schemeData   Output version of schemeData (see below).
%
% A CFL constrained ODE system is described by a function with prototype
%
%        [ ydot, stepBound ] = schemeFunc(t, y, schemeData)
%
%   where t is the current time, y the current state vector and
%   schemeData is passed directly through.  The output stepBound
%   is the maximum allowed time step that will be taken by this function
%   (typically the option parameter factorCFL will choose a smaller step size).
%
% The time interval tspan may be given as 
%   1) A two entry vector [ t0 tf ], in which case the output will 
%      be scalar t = tf and a row vector y = y(tf).
%   2) A vector with three or more entries, in which case the output will
%      be column vector t = tspan and each row of y will be the solution
%      at one of the times in tspan.  Unlike Matlab's ode suite routines,
%      this version just repeatedly calls version (1), so it is not
%      particularly efficient.
%
% Note that using this routine for integrating HJ PDEs will usually
%   require that the data array be turned into a vector before the call
%   and reshaped into an array after the call.  Option (2) for tspan should 
%   not be used in this case because of the excessive memory requirements
%   for storing solutions at multiple timesteps.
%
% The output version of schemeData will normally be identical to the input
%   version, and therefore can be ignored.  However, if a PostTimestep 
%   routine is used (see odeCFLset) then schemeData may be modified during 
%   integration, and the version of schemeData at tf is returned in this
%   output argument.

% Copyright 2004 Ian M. Mitchell ([e-mail]).
% This software is used, copied and distributed under the licensing 
%   agreement contained in the file LICENSE in the top directory of 
%   the distribution.
%
% Ian Mitchell, 5/14/03.
% Calling parameters modified to more closely match Matlab's ODE suite
%   Ian Mitchell, 2/14/04.

  %---------------------------------------------------------------------------
  % How close (relative) do we need to be to the final time?
  small = 100 * eps;
  
  %---------------------------------------------------------------------------
  % Make sure we have the default options settings
  if((nargin < 4) | isempty(options))
    options = odeCFLset;
  end
  
  %---------------------------------------------------------------------------
  numT = length(tspan);
  
  %---------------------------------------------------------------------------
  % If we were asked to integrate forward to a final time.
  if(numT == 2)

    t = tspan(1);
    y = y0;
    steps = 0;
    startTime = cputime;
    
    while(tspan(2) - t > small * abs(tspan(2)))

      % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      % First substep: Forward Euler from t_n to t_{n+1}.
      
      % Approximate the derivative and CFL restriction.
      [ ydot, stepBound ] = feval(schemeFunc, t, y, schemeData);

      % CFL bound on timestep, but not beyond the final time.
      %   We'll use this fixed timestep for both substeps.
      deltaT = min([ options.factorCFL * stepBound; ...
                     tspan(2) - t; options.maxStep ]);
      
      % Take the first substep.
      y1 = y + deltaT * ydot;
      t1 = t + deltaT;

      % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      % Second substep: Forward Euler from t_{n+1} to t_{n+2}.

      % Approximate the derivative; we do not need the CFL information.
      ydot = feval(schemeFunc, t1, y1, schemeData);

      % We have already computed a deltaT.
      
      % Take the second substep.
      y2 = y1 + deltaT * ydot;
      t2 = t1 + deltaT;

      % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      % Combine t_n and t_{n+2} to get approximation at t_{n+1/2}
      yHalf = 0.25 * (3 * y + y2);
      tHalf = 0.25 * (3 * t + t2);

      % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      % Third substep: Forward Euler from t_{n+1/2} to t_{n+3/2}.

      % Approximate the derivative; we do not need the CFL information.
      ydot = feval(schemeFunc, tHalf, yHalf, schemeData);

      % We have already computed a deltaT.
      
      % Take the second substep.
      yThreeHalf = yHalf + deltaT * ydot;
      tThreeHalf = tHalf + deltaT;

      % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      % Combine t_n and t_{n+3/2} to get third order approximation of t_{n+1}.
      y = (1/3) * (y + 2 * yThreeHalf);
      t = (1/3) * (t + 2 * tThreeHalf);

      steps = steps + 1;

      % If there is a post-timestep routine, call it.
      if(~isempty(options.postTimestep))
        [ y schemeData ] = feval(options.postTimestep, t, y, schemeData);
      end

      % If we are in single step mode, then do not repeat.
      if(strcmp(options.singleStep, 'on'))
        break;
      end
    
    end
    
    endTime = cputime;

    if(strcmp(options.stats, 'on'))
      fprintf('\t%d steps in %g seconds from %g to %g\n', ...
              steps, endTime - startTime, tspan(1), t);
    end

  %---------------------------------------------------------------------------
  % If we were asked for the solution at multiple timesteps,
  %   just use recursion over each pair of timesteps.
  elseif(numT > 2)
    t = reshape(tspan, numT, 1);
    y = zeros(numT, length(y0));
    y(1,:) = y0';
    for i = 2 : numT
      [ garbage, yout ] = odeCFL3(schemeFunc, [ t(i-1), t(i) ], y(i-1,:)', ...
                                  schemeData, options)
      y(i,:) = yout';
    end
    
  %---------------------------------------------------------------------------
  else
    error('tspan must contain at least two entries');
  end