stiff.m

matlab算法集 matlab算法集

function [t,X,e,k] = stiff (x0,t0,t1,m,tol,q,f,fx,ft) 
%---------------------------------------------------------------------- 
% Usage:       [t,X,e,k] = stiff (x0,t0,t1,m,tol,q,f,fx,ft);
%
% Description: Use the semi-implicit Bulirsch-Stoer variable order
%              extrapotion method to solve a stiff n-dimensional
%              system of ordinary differential equations:
%
%                 dx(t)/dt = f(t,x(t))   ,   x(t0) = x0
%
%              at m step equally spaced over the interval [t0,t1].
%                   
% Inputs:      x0  = n by 1 initial conditon vector
%              t0  = initial time
%              t1  = final time
%              m   = number of solution points (m >= 2)
%                
%              tol = upper bound on local truncation error (tol >= 0)
%              q   = maximum number of scalar function evaluations (q >= 1) 
%              f   = string containing name of function which defines the
%                    system of equations to be solved.
%              fx  = string containing name of optional function that
%                    defines the Jacobian matrix J = df(t,x)/dx
%              ft  = string containing name of optional function that
%                    defines df(t,x)/dt.  The forms of the user-supplied
%                    functions are as follows:
%
%                        function dx   = f(t,x)
%                        function dfdx = fx(t,x)
%                        function dfdt = ft(t,x)
%
%                    When f is called, it should evaluate the right-hand
%                    side of the system of differential equations and
%                    return the result in the n by 1 vector dx. When fx
%                    is called, it should evaluate the n by n matrix of
%                    partial derivatives of f(t,x) with respect to x. That
%                    is, dfdx(k,j) = df(k)/dx(j) for 1 <= k,j <= n. When
%                    ft is called, it should evaluate the n by 1 vector of
%                    partial derivatives of f(t,x) with respect to t. That
%                    is, dfdt(k) = df(k)/dt for 1 <= k <= n. The functions
%                    fx and ft are optional in the sense that stiff can be
%                    called with the arguments fx and/or ft set to the
%                    empty string, ''. When this is done, the partial 
%                    derivatives are approximated numerically by stiff 
%                    using using central differences. Thus the user does
%                    not have to supply the functions in this case. 
%
% Outputs:     t  = m by 1 vector containing solution times
%              X  = m by n matrix containing solution trajectory. 
%                   Tke kth row of X is the solution at time t(k).
%              e  = maximum local truncation error over [t0,t1]     
%              k  = number of scalar function evaluations.  If k < q,
%                   then the following error criterion was satisfied 
%                   by adjusting the extrapolation level.  Here ||E||
%                   denotes the infinity norm of the estimated local
%                   truncation error:
%
%                       ||E|| < max(||x||,1)*tol  
%---------------------------------------------------------------------- 

% Check calling arguments

   chkvec (x0,1,'stiff');
   m = args (m,2,m,4,'stiff');
   tol = args (tol,0,tol,5,'stiff');
   q = args (q,1,q,6,'stiff');
   chkfun (feval(f,t0,x0),7,'stiff')
   funft = getfun (ft);
   if funft
      chkfun (feval(ft,t0,x0),9,'stiff')
   end

% Initialize

   k = 0;  		
   n = length(x0);
   x = zeros (n,1);
   x = x0;
   X = zeros (m,n);
   X(1,:) = x0';
   t = linspace (t0,t1,m)'; 
   e = 0;
   
% Solve the system 

   hwbar = waitbar (0,'Solving Ordinary Differential Equations: stiff');
   for i = 2 : m
      waitbar ((i-2)/(m-2))
      [x,E,r] = stiff0 (x,t(i-1),t(i),tol,q,f,fx,ft);
      e = max(e,E);
      X(i,:) = x';
      k = k + r;
      if k >= q
         fprintf ('\nMaximum of %g function evaluations in stiff\n',k);
         break
      end
   end
   close(hwbar)