stiff0.m

matlab算法集 matlab算法集

function [x,e,s] = stiff0 (x,t0,t1,tol,q,f,fx,ft,dm)
%------------------------------------------------------------------
% Usage:       [x,e,s] = stiff0 (x,t0,t1,tol,q,f,fx,ft,dm);
% 
% Description: Use the semi-implicit Bulirsch-Stoer variable order
%              extrapolation method to solve a stiff n-dimensional
%              system of ordinary differential equations:
%
%                 dx(t)/dt = f(t,x(t))   
%
%              over the time interval [t0,t1].
%                   
% Inputs:       x   = n by 1 initial conditon vector, x = x(t0)
%               t0  = initial time
%               t1  = final time
%               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 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 stiff0 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 stiff0
%                     using using central differences. Thus the user does
%                     not have to supply the functions in this case. 
%              dm   = optional display mode.  If present, intermediate
%                     results are displayed.   
%
% Outputs:     x  = n by 1 vector containing estimate of solution:
%                   x = x(t1)
%              e = maximum local truncation error over [t0,t1]     
%              s = number of scalar function evaluations.  If it is less
%                  than the user-specified maximum, q, then the following
%                  error criterion was satisfied by adjusting the
%                  extrapolation level.  Here ||E|| denotes theinfinity
%                  norm of the estimated local truncation error:
%
%                      ||E|| < max(||x||,1)*tol  
%------------------------------------------------------------------

% Initialize 

   chkvec (x,1,'stiff0');
   tol = args(tol,0,tol,4,'stiff0');
   q = args(q,0,q,5,'stiff0');
   chkfun (feval(f,t0,x),6,'stiff0');

   display = nargin > 8;
   s = 0;
   k = 0;			% starting level 
   levels = 8;		        % maximum level 
   r = 2^levels+1;  	        % maximum steps 
   n = length(x);
   I = eye(n,n);
   Ak  = zeros (n,levels);		% row k of A 
   Ak1 = zeros (n,levels);		% row k-1 of A 
   Q = zeros (r,n);			% midpoint estimates 	
   dx = zeros (n,1);			% derivative vector 
   ev = zeros (n,1);			% error vector 
   v = zeros (n,1);			% linear system solution  
   b = zeros (n,1);			% right-hand side 
   J = zeros (n,n);			% Jacobian matrix  
   Q(1,:) = x';
   e = max(norm(x,inf),1)*tol + 1;
   funfx = getfun (fx);
   funft = getfun (ft);
   if funft
      chkfun (feval(ft,t0,x),8,'stiff0');
   end

   while (k == 1) | ((e > max(norm(x,inf),1)*tol) & ...
         (k < levels) & (s < q)); 

% Apply semi-implicit midpoint method 
      
      k = k + 1;
      m = 2^k;
      h = (t1 - t0)/m;

% First step 

      dx = feval(f,t0,Q(1,:));
      s = s + n;
      if ~funfx
         J = fxstiff (t0,Q(1,:)',f);
         s = s + 2*n*n; 
      else 
         J = feval(fx,t0,Q(1,:)');
         s = s + n*n;
      end
      if ~funft
         b = ftstiff (t0,Q(1,:)',f);
         s = s + 2*n; 
      else 
         b = feval(ft,t0,Q(1,:)');
         s = s + n; 
      end 
      for i = 1 : n
         b(i) = h*(h*b(i) + dx(i));
      end
      J = I - h*J; 
      v = ludec (J,b); 
      for i = 1 : n
         Q(2,i) = Q(1,i) + v(i);
      end
      
% Midpoint steps 

      for j = 1 : m-1
         t = t0 + j*h;   
         b = feval(f,t,Q(j+1,:)');
         s = s + n;
         if ~funfx' 
            J = fxstiff (t,Q(j+1,:)',f);
            s = s + 2*n*n; 
         else  	
            J = feval(fx,t,Q(j+1,:)');
            s = s + n*n; 
         end
         for i = 1 : n
            for p = 1 : n
              b(i) = b(i) - J(i,p)*(Q(j+1,p) - Q(j,p));
            end
            b(i) = b(i)*2*h; 
         end
         J = I - h*J;
         v = ludec (J,b); 
         for i = 1 : n
            Q(j+2,i) = Q(j,i) + v(i); 
         end
      end
      
% Last step 

      dx = feval(f,t1,Q(m+1,:));
      if ~funfx
         J = fxstiff (t1,Q(m+1,:)',f);
         s = s + 2*n*n; 
      else 	
         J = feval(fx,t1,Q(m+1,:)');
         s = s + n*n; 
      end
      if ~funft 
         b = ftstiff (t0,Q(1,:)',f);
         s = s + 2*n; 
      else 
         b = feval(ft,t0,Q(1,:)');
         s = s + n; 
      end 
      for i = 1 : n
         b(i) = h*(h*b(i) + dx(i));
      end
      J = I - h*J;
      v = ludec (J,b);
      Ak1 = Ak;
      for i = 1 : n
         Ak(i,1) = (Q(m,i) + Q(m+1,i) + v(i))/2;
      end
         
% Update row k of matrix of extrapolated estimates 

      for j = 2 : k
         r0 = 4^(j-1);
         for i = 1 : n
            Ak(i,j) = (r0*Ak(i,j-1) - Ak1(i,j-1))/(r0 - 1);
            ev(i) = (Ak(i,j-1) - Ak1(i,j-1))/(r0 - 1); 
         end
         if display
            if j == 2
               fprintf ('\nInterval: [%.7g, %.7g]',t0,t1);
            end   
            fprintf ('\n %g %g',k,j);
            for i = 1 : n
               fprintf (' %10.3e',ev(i));
            end
            if j == k 
               fprintf ('\n');
               wait;
            end    
         end
      end
      e = norm (ev,inf);
      for i = 1 : n
         x(i) = Ak(i,k);
      end

   end
   s = min(s,q);