rkfstep.m

matlab算法集 matlab算法集

function [x1,r] = rkfstep (x0,t0,t1,f)
%----------------------------------------------------------------
% Usage:       [x1,r] = rkfstep (x0,t0,t1,f);
%
% Description: Perform one step of the fifth order Runge-Kutta Fehlberg
%              method to solve an n-dimensional system of ordinary
%              differential equations:
%
%              dx(t)/dt = f(t,x(t)),        x(t0) = x0 
%
% Inputs       x0 = n by 1 initial condition vector, x0 = x(t0)
%              t0 = initial time
%              t1 = final time
%              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:    x1 = n by 1 vector containing estimate of solution: 
%                  x1 = x(t1)
%             r  = the infinity norm of the local truncation error E
%                  at time t1. 
%----------------------------------------------------------------
                    
   chkvec (x0,1,'rkfstep');
   chkfun (feval(f,t0,x0),4,'rkfstep');
   h = t1 - t0;
   n = length (x0);
   K = zeros (n,6);
   x = zeros (n,1);
   e = zeros (n,1);
   x1 = zeros (n,1);
       
% Evaluate derivatives at points inside [t,t+h] 

   K(:,1) = feval(f,t0,x0);
   for i = 1 : n
      d = x0(i)+ h*K(i,1)/4;
      x(i) = d;
   end

   K(:,2) = feval(f,t0+h/4,x);
   for i = 1 : n
      d = x0(i)+ h*(3*K(i,1) + 9*K(i,2))/32.;
      x(i) = d;
   end

   K(:,3) = feval(f,t0+3*h/8,x);     
   for i = 1 : n
      d = x0(i)+ h*(1932*K(i,1) - 7200*K(i,2) + 7296*K(i,3))/2197;
      x(i) = d;
   end

   K(:,4) = feval(f,t0+12*h/13,x);
   for i = 1 : n
      d = x0(i)+ h*(439*K(i,1)/216 - 8*K(i,2) + 3680*K(i,3)/513 ...
          - 845*K(i,4)/4104);
      x(i) = d;
   end

   K(:,5) = feval(f,t1,x);
   for i = 1 : n
      d = x0(i)+ h*(-8*K(i,1)/27 + 2*K(i,2) - 3544*K(i,3)/2565 ...
          + 1859*K(i,4)/4104 - 11*K(i,5)/40);
      x(i) = d;
   end

% Compute x1 = x(t+h) 

   K(:,6) = feval(f,t0+h/2,x);
   for i = 1 : n
      d = x0(i)+ h*(16*K(i,1)/135 + 6656*K(i,3)/12825 ...
          + 28561*K(i,4)/56430 - 9*K(i,5)/50 + 2*K(i,6)/55);
      x1(i) = d;
   end

% Compute error E 

   for i = 1 : n 
      d = h*(K(i,1)/360 - 128*K(i,3)/4275 - 2197*K(i,4)/75240 ...
          + K(i,5)/50 + 2*K(i,6)/55);
      e(i) = d;
   end
   r = norm (e,inf);