rka.m

matlab的数学物理方程数值算法源程序。这是"Numerical Methods for Physics"第二版的matlab源程序。

function [xSmall, t, tau] = rka(x,t,tau,err,derivsRK,param)
% Adaptive Runge-Kutta routine
% Inputs
%   x          Current value of the dependent variable
%   t          Independent variable (usually time)
%   tau        Step size (usually time step)
%   err        Desired fractional local truncation error
%   derivsRK   Right hand side of the ODE; derivsRK is the
%              name of the function which returns dx/dt
%              Calling format derivsRK(x,t,param).
%   param      Extra parameters passed to derivsRK
% Outputs
%   xSmall     New value of the dependent variable
%   t          New value of the independent variable
%   tau        Suggested step size for next call to rka

%* Set initial variables
tSave = t;  xSave = x;    % Save initial values
safe1 = .9;  safe2 = 4.;  % Safety factors

%* Loop over maximum number of attempts to satisfy error bound
maxTry = 100;  
for iTry=1:maxTry
	
  %* Take the two small time steps
  half_tau = 0.5 * tau;
  xTemp = rk4(xSave,tSave,half_tau,derivsRK,param);
  t = tSave + half_tau;
  xSmall = rk4(xTemp,t,half_tau,derivsRK,param);
  
  %* Take the single big time step
  t = tSave + tau;
  xBig = rk4(xSave,tSave,tau,derivsRK,param);
  
  %* Compute the estimated truncation error
  scale = err * (abs(xSmall) + abs(xBig))/2.;
  xDiff = xSmall - xBig;
  errorRatio = max( abs(xDiff)./(scale + eps) );
  
  %* Estimate new tau value (including safety factors)
  tau_old = tau;
  tau = safe1*tau_old*errorRatio^(-0.20);
  tau = max(tau,tau_old/safe2);
  tau = min(tau,safe2*tau_old);
  
  %* If error is acceptable, return computed values
  if (errorRatio < 1)  return;  end 
end

%* Issue error message if error bound never satisfied
error('ERROR: Adaptive Runge-Kutta routine failed');
return;