rkf0.m

matlab算法集 matlab算法集

function [x,h,e,k] = rkf0 (x,t0,t1,h,tol,q,f)
%----------------------------------------------------------------
% Usage:       [x,h,e,k] = rkf0 (x,t0,t1,h,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))   
%
%              over the time interval [t0,t1].
%                   
% Inputs:      x   = n by 1 initial conditon vector, x = x(t0)
%              t0  = initial time
%              t1  = final time
%              h   = initial step size
%              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:     x = n by 1 vector containing estimate of solution:
%                  x = x(t1)
%              h = final step size
%              e = maximum local truncation error over [t0,t1]     
%              k = 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 step
%                  size. Here ||E|| denotes the infinity norm of the
%                  estimated local truncation error:
%
%                       ||E|| < max(||x||,1)*tol  
%----------------------------------------------------------------

% Initialize 

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

   k = 0;			% function evaluations 	
   gamma = 0.8;	% safety factor 
   t = t0;			% independent variable 
   em = sqrt(eps);
   n = length(x);   
   d = em*max(abs(t0),abs(t1));	       % smallest step 
   hmin = d;
   x1 = zeros (n,1);
   y = zeros (n,1);
   e = 0;

% Adjust step size 

   delta = t1 - t0;
   dt = delta;
   x1 = x;
   E0 = 0;
   E1 = E0 + 1; 
   Dt = hmin + 1;

% Solve the system 

   while  ((delta > 0) & (t < t1)) | ((delta <= 0) & (t > t1)) 

      E1 = E0 + 1;
      while (E1 > E0) & (abs(Dt) >= hmin) & (k < q)

         [y,E1] = rkfstep (x1,t,t+dt,f);         % actual error 
         k = k + 6*n; 
         E0 = max(norm(y,inf),1)*tol;		% desired error 
         if E1 <= E0 
            t = t + dt;
            x1 = y;
         end
         r = max(em*em,E1);
         d = gamma*dt*((E0/r)^0.2);	        % new step 
         Dt = d;
         h = Dt; 
         if Dt > 0      		        % final step ? 
            dt = min(Dt,t1-t);
         else
            dt = max(Dt,t1-t);
         end 

      end

      e = max(e,E1);
      if E1 > E0
         t = t + dt;
         x1 = y;
      end
      if abs(Dt) < hmin
         fprintf ('\nStep size was smaller than %g in rkf0\n',hmin);
      end
      
   end
   
% Finalize       

   k = min(q,k);
   x = y;