rkf.m

matlab算法集 matlab算法集

function [t,X,e,k] = rkf (x0,t0,t1,m,tol,q,f)
%---------------------------------------------------------------------- 
% Usage:      [t,X,e,k] = rkf (x0,t0,t1,m,tol,q,f);
%
% Description: Use the fifth order Runge-Kutta Fehlberg method with
%              automatic step size control to solve an n-dimensional
%              system of ordinary differential equations:
%
%                  dx(t)/dt = f(t,x(t))   ,   x(t0) = x0
%
%              at m steps equally spaced over the time interval [t0,t1].
%
% Inputs:      x   = n by 1 initial conditon vector, x = x(t0)
%              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. 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.
%
% 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 step size. Here ||E|| denotes the
%                  infinity norm of the estimated local truncation 
%                  error:
%
%                      ||E|| < max(||x||,1)*tol  
%---------------------------------------------------------------------- 


% Initialize 

   chkvec (x0,1,'rkf');
   m   = args (m,2,m,4,'rkf');
   tol = args (tol,0,tol,5,'rkf');
   q   = args (q,1,q,6,'rkf');
   chkfun (feval(f,t0,x0),7,'rkf')

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


   hwbar = waitbar(0,'Solving Ordinary Differential Equations: rkf');
   for i = 2 : m
      waitbar ((i-1)/(m-1));
      [x,h,E,r] = rkf0 (x,t(i-1),t(i),h,tol,q,f);
      k = k + r;
      e = max(e,E);
      if k >= q
         fprintf ('\nA maximum of %g steps taken in rkf.\n',k);
         break
      end
      X(i,:) = x';
   end
   close (hwbar)