adaptrk45.m

用matlab写的一些数值算法

function  [x,y,neval] = adaptrk45(f, ab, y0, h0, tol)
% Runge-Kutta with adaptive step length control
% Input
% f     handle to function that defines the rhs in  y' = f(x,y)
% ab    range of integration is [a, b],  a=ab(1), b=ab(2)
% y0    vector of initial values
% h0    initial step length
% tol   tolerance for step length control
% Output
% x     vector of x-values, where the solution is computed
% y     matrix with approximations of the solution
% neval number of function evaluations

% Version 11.12.2003.  INCBOX

% Runge-Kutta-Fehlberg45 constants.  Alpha in A, beta in B
A = [0 1/4 3/8 12/13 1 1/2];
B = [0 0 0 0 0 0
     1/4 0 0 0 0 0
     3/32 9/32 0 0 0 0
     1932/2197 -7200/2197 7296/2197 0 0 0
     439/216 -8 3680/513 -845/4104 0 0
     -8/27 2 -3544/2565 1859/4104 -11/40 0]';
% Weights   
w4 = [25/216 0 1408/2565 2197/4104 -1/5 0]'; 
w5 = [16/135 0 6656/12825 28561/56430 -9/50 2/55]'; 
dw = w4 - w5;
% Initialize
a = ab(1);  b = ab(2);
m = length(y0);
x = a;  y = y0(:);  neval = 0;              % initialize output
n = 0;  xn = a;  yn = y0;  h = min(h0, b-a);
K = zeros(m,6);                         % for storing k1,...,k6
while  xn < b
  xn1 = min(xn + h, b);            % possible adjustment at end
  h = xn1 - xn;  
  for  j = 1 : 6
    K(:,j) = feval(f, xn+A(j)*h, yn+h*K*B(:,j));
  end
  neval = neval + 6;
  y5 = yn + h*K*w5;                        % RKF5 approximation
  d = h*norm(K*dw);                           % estimated error
  if  d <= tol                                % accept the step
    n = n+1;  xn = xn1;  yn = y5;
    x(n+1) = xn;   y(:,n+1) = yn;              % save in output
  end  
  h = h * min(0.8*(tol/d)^0.2, 5);         % Update step length
end
x = x';  y = y';                    % return in standard format