bsext0.m

matlab算法集 matlab算法集

function [x,e,p] = bsext0 (x,t0,t1,tol,q,f,dm)
%------------------------------------------------------------------
% Usage:       [x,e,p] = bsext0 (x,t0,t1,tol,q,f,dm);
%
% Description: Use the Bulirsch-Stoer variable order extrapotion method
%              to solve an 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. The form of f is:
%
%                    function dx = f(t,x)
%
%                    When f is called, it must evaluate the right-hand
%                    side of the system of differential equations and
%                    return the result in the n by 1 vector dx.
%             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]     
%             p = number of scalar function evaluations.  If p < 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  
%------------------------------------------------------------------

% Initialize 

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

   k = 0;
   levels = 8;		              % maximum level 
   r =2^levels + 1;	              % maximum steps 
   p = 0;		
   n = length(x);
   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 
   Q(1,:) = x';
   e = max(norm(x,inf),1)*tol + 1;

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

% Apply modified midpoint method 

      k = k + 1;
      m = 2^k;
      h = (t1 - t0)/m;

% First step 
      
      dx = feval(f,t0,Q(1,:));
      p = p + n;
      for i = 1 : n
         Q(2,i) = Q(1,i) + h*dx(i);
      end

% Midpoint steps 

      for j = 1 : m-1
         t = t0 + j*h;      
         dx = feval(f,t,Q(j+1,:));
         p = p + n;
         for i = 1 : n
            Q(j+2,i) = Q(j,i) + 2*h*dx(i); 
         end
      end

% Last step 

      Ak1 = Ak;      
      dx = feval(f,t1,Q(m+1,:));
      p = p + n;
      for i = 1 : n
        Ak(i,1) = (Q(m+1,i) + Q(m,i) + h*dx(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
            fprintf ('\n %g %g',k,j);
            for i = 1 : n
               fprintf (' %10.3e',ev(i));
            end  
         end
      end 
         if display
            fprintf ('\n');
            wait;
         end
      e = norm (ev,inf);
%      show ('e',e)   
      for i = 1 : n
         x(i) = Ak(i,k);
      end
      
%      show ('Ak',Ak')

   end
   p = min(p,q);