predcorr.m

matlab算法集 matlab算法集

function [X,E] = predcorr(X,t0,t1,f)
%------------------------------------------------------------------
% Usage:       X = predcorr(X,t0,t1,f)
% 
% Description: Perform one step of Adam's fourth-order predictor-corrector
%              method to solve an n-dimensional system of ordinary
%              differential equations:
%
%                  dx(t)/dt = f[t,x(t)]
%
% Inputs:      X  = 4 by n matrix X containing past solutions where
%                   h = t1 - t0:
%                    
%                       ---------------------------
%                       Row  On entry    On exit
%                       ---------------------------
%                        1     x(t0)       x(t0+h)
%                        2     x(t0-h)     x(t0)
%                        3     x(t0-2h)    x(t0-h)
%                        4     x(t0-3h)    x(t0-2h)
%                       ---------------------------
%
%              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:     X  = 4 by n matrix X updated to include current and
%                   past solutions (see above). 
%------------------------------------------------------------------
 
   chkfun(feval(f,t0,X(1,:)'),4,'predcorr');
   h = t1-t0;
   predictor = [55 -59 37 -9]';
   corrector = [9 19 -5 1]';

% Initialize 

   n = size(X,2);
   dx = zeros (n,1);
   xp = zeros (n,1);
   xc = zeros (n,1);
   
% Compute fourth-order Adams-Bashforth predictor: xp /

   xp = X(1,:)';
   for i = 1 : 4
      t = t0 - (i-1)*h; 	
      dx = feval(f,t,X(i,:)');
      for j = 1 : n
         xp(j) = xp(j) + h*predictor(i)*dx(j)/24;
      end
   end
   
% Compute fourth-order Adams-Moulton corrector: xc 

   xc = X(1,:)';
   for i = 1 : 4
      t = t0 - (i-1)*h;  
      if i == 1   	
         dx = feval(f,t+h,xp);
      else
         dx = feval(f,t+h,X(i-1,:)');
      end
      for j = 1 : n
         xc(j) = xc(j) + h*corrector(i)*dx(j)/24;
      end
   end
   
% Compute local truncation error 

   for j = 1 : n
      dx(j) = 19*(xp(j) - xc(j))/270;
   end
   E = norm(dx,inf);   

% Update X 

   for i = 4 : -1 : 2
      for j = 1 : n
         X(i,j) = X(i-1,j);
      end
   end 
   for j = 1 : n
      X(1,j) = xc(j) + dx(j);
   end