function [ y, yp, t, flag ] = rkf45 ( f, neqn, y, yp, t, tout, relerr, ...
abserr, flag )
%*****************************************************************************80
%
%% rkf45 carries out the Runge-Kutta-Fehlberg method.
%
% Discussion:
%
% This routine is primarily designed to solve non-stiff and mildly stiff
% differential equations when derivative evaluations are inexpensive.
% It should generally not be used when the user is demanding
% high accuracy.
%
% This routine integrates a system of NEQN first-order ordinary differential
% equations of the form:
%
% dY(i)/dT = F(T,Y(1),Y(2),...,Y(NEQN))
%
% where the Y(1:NEQN) are given at T.
%
% Typically the subroutine is used to integrate from T to TOUT but it
% can be used as a one-step integrator to advance the solution a
% single step in the direction of TOUT. On return, the parameters in
% the call list are set for continuing the integration. The user has
% only to call again (and perhaps define a new value for TOUT).
%
% Before the first call, the user must
%
% * supply the subroutine F(T,Y,YP) to evaluate the right hand side;
% and declare F in an EXTERNAL statement;
%
% * initialize the parameters:
% NEQN, Y(1:NEQN), T, TOUT, RELERR, ABSERR, FLAG.
% In particular, T should initially be the starting point for integration,
% Y should be the value of the initial conditions, and FLAG should
% normally be +1.
%
% Normally, the user only sets the value of FLAG before the first call, and
% thereafter, the program manages the value. On the first call, FLAG should
% normally be +1 (or -1 for single step mode.) On normal return, FLAG will
% have been reset by the program to the value of 2 (or -2 in single
% step mode), and the user can continue to call the routine with that
% value of FLAG.
%
% (When the input magnitude of FLAG is 1, this indicates to the program
% that it is necessary to do some initialization work. An input magnitude
% of 2 lets the program know that that initialization can be skipped,
% and that useful information was computed earlier.)
%
% The routine returns with all the information needed to continue
% the integration. If the integration reached TOUT, the user need only
% define a new TOUT and call again. In the one-step integrator
% mode, returning with FLAG = -2, the user must keep in mind that
% each step taken is in the direction of the current TOUT. Upon
% reaching TOUT, indicated by the output value of FLAG switching to 2,
% the user must define a new TOUT and reset FLAG to -2 to continue
% in the one-step integrator mode.
%
% In some cases, an error or difficulty occurs during a call. In that case,
% the output value of FLAG is used to indicate that there is a problem
% that the user must address. These values include:
%
% * 3, integration was not completed because the input value of RELERR, the
% relative error tolerance, was too small. RELERR has been increased
% appropriately for continuing. If the user accepts the output value of
% RELERR, then simply reset FLAG to 2 and continue.
%
% * 4, integration was not completed because more than MAXNFE derivative
% evaluations were needed. This is approximately (MAXNFE/6) steps.
% The user may continue by simply calling again. The function counter
% will be reset to 0, and another MAXNFE function evaluations are allowed.
%
% * 5, integration was not completed because the solution vanished,
% making a pure relative error test impossible. The user must use
% a non-zero ABSERR to continue. Using the one-step integration mode
% for one step is a good way to proceed.
%
% * 6, integration was not completed because the requested accuracy
% could not be achieved, even using the smallest allowable stepsize.
% The user must increase the error tolerances ABSERR or RELERR before
% continuing. It is also necessary to reset FLAG to 2 (or -2 when
% the one-step integration mode is being used). The occurrence of
% FLAG = 6 indicates a trouble spot. The solution is changing
% rapidly, or a singularity may be present. It often is inadvisable
% to continue.
%
% * 7, it is likely that this routine is inefficient for solving
% this problem. Too much output is restricting the natural stepsize
% choice. The user should use the one-step integration mode with
% the stepsize determined by the code. If the user insists upon
% continuing the integration, reset FLAG to 2 before calling
% again. Otherwise, execution will be terminated.
%
% * 8, invalid input parameters, indicates one of the following:
% NEQN <= 0;
% T = TOUT and |FLAG| /= 1;
% RELERR < 0 or ABSERR < 0;
% FLAG == 0 or FLAG < -2 or 8 < FLAG.
%
% Licensing:
%
% This code is distributed under the MIT license.
%
% Modified:
%
% 05 April 2011
%
% Author:
%
% Original FORTRAN77 version by Herman Watts, Lawrence Shampine.
% MATLAB version by John Burkardt.
%
% Reference:
%
% Erwin Fehlberg,
% Low-order Classical Runge-Kutta Formulas with Stepsize Control,
% NASA Technical Report R-315, 1969.
%
% Lawrence Shampine, Herman Watts, S Davenport,
% Solving Non-stiff Ordinary Differential Equations - The State of the Art,
% SIAM Review,
% Volume 18, pages 376-411, 1976.
%
% Input:
%
% function F: a user-supplied function to evaluate the
% derivatives dYdT(T), of the form:
% function yp = f ( t, y )
%
% integer NEQN: the number of equations to be integrated.
%
% real Y(NEQN): the current solution vector at T.
%
% real YP(NEQN): the current value of the
% derivative of the dependent variable. The user should not set or alter
% this information!
%
% real T: the current value of the independent variable.
%
% real TOUT, the output point at which solution is
% desired. TOUT = T is allowed on the first call only, in which case
% the routine returns with FLAG = 2 if continuation is possible.
%
% real RELERR, ABSERR: the relative and absolute
% error tolerances for the local error test. At each step the code
% requires:
% abs ( local error ) <= RELERR * abs ( Y ) + ABSERR
% for each component of the local error and the solution vector Y.
% RELERR cannot be "too small". If the routine believes RELERR has been
% set too small, it will reset RELERR to an acceptable value and return
% immediately for user action.
%
% integer FLAG: chooses normal run or single step option.
% set FLAG to +1 for normal use, or to -1 for single step mode.
%
% Output:
%
% real Y(NEQN): the updated solution vector at T.
%
% real YP(NEQN): the updated derivative of the dependent variable.
%
% real T: the updated value of the independent variable.
%
% integer FLAG: indicator for status of integration.
% A value of 2 or -2 indicates normal progress,
% while any other value indicates a problem that should be addressed.
%
persistent abserr_save;
persistent flag_save;
persistent h;
persistent init;
persistent kflag;
persistent kop;
persistent nfe;
persistent relerr_save;
if ( isempty ( abserr_save ) )
abserr_save = -1.0;
h = -1.0;
flag_save = -1000;
init = -1000;
kflag = -1000;
kop = -1;
nfe = -1;
relerr_save = -1.0;
end
maxnfe = 3000;
remin = 1.0E-12;
%
% Check the input parameters.
%
if ( neqn < 1 )
flag = 8;
return
end
if ( relerr < 0.0 )
flag = 8;
return
end
if ( abserr < 0.0 )
flag = 8;
return
end
if ( flag == 0 || 8 < flag || flag < -2 )
flag = 8;
return
end
mflag = abs ( flag );
%
% Is this a continuation call?
%
if ( mflag ~= 1 )
if ( t == tout && kflag ~= 3 )
flag = 8;
return
end
if ( mflag == 2 )
if ( kflag == 3 )
flag = flag_save;
mflag = abs ( flag );
elseif ( init == 0 )
flag = flag_save;
elseif ( kflag == 4 )
nfe = 0;
elseif ( kflag == 5 && abserr == 0.0 )
fprintf ( 1, '\n' );
fprintf ( 1, 'RKF45 - Fatal error!\n' );
fprintf ( 1, ' KFLAG = 5 and ABSERR = 0.\n' );
error ( 'RKF45 - Fatal error!' );
elseif ( kflag == 6 && relerr <= relerr_save && abserr <= abserr_save )
fprintf ( 1, '\n' );
fprintf ( 1, 'RKF45 - Fatal error!\n' );
fprintf ( 1, ' KFLAG = 6 and\n' );
fprintf ( 1, ' RELERR <= RELERR_SAVE and\n' );
fprintf ( 1, ' ABSERR <= ABSERR_SAVE\n' );
error ( 'RKF45 - Fatal error!' );
end
%
% FLAG = 3, 4, 5, 6, 7 or 8.
%
else
if ( flag == 3 )
flag = flag_save;
if ( kflag == 3 )
mflag = abs ( flag );
end
elseif ( flag == 4 )
nfe = 0;
flag = flag_save;
if ( kflag == 3 )
mflag = abs ( flag );
end
elseif ( flag == 5 && 0.0 < abserr )
flag = flag_save;
if ( kflag == 3 )
mflag = abs ( flag );
end
%
% Integration cannot be continued because the user did not respond to
% the instructions pertaining to FLAG = 5, 6, 7 or 8.
%
else
fprintf ( 1, '\n' );
fprintf ( 1, 'RKF45 - Fatal error!\n' );
fprintf ( 1, ' Integration cannot continue because the user\n' );
fprintf ( 1, ' did not respond to the instructions pertaining\n' );
fprintf ( 1, ' to the output flag warning FLAG = 5, 6, 7 or 8.\n' );
error ( 'RKF45 - Fatal error!' );
end
end
end
%
% Save the input value of FLAG.
% Set the continuation flag KFLAG for subsequent input checking.
%
flag_save = flag;
kflag = 0;
%
% Save RELERR and ABSERR for checking input on subsequent calls.
%
relerr_save = relerr;
abserr_save = abserr;
%
% Restrict the relative error tolerance to be at least
%
% 2 * EPS + REMIN
%
% to avoid limiting precision difficulties arising from impossible
% accuracy requests.
%
relerr_min = 2.0 * eps + remin;
%
% Is the relative error tolerance too small?
%
if ( relerr < relerr_min )
relerr = relerr_min;
flag = 3;
kflag = 3;
return
end
dt = tout - t;
%
% Initialization:
%
% Set the initialization completion indicator, INIT;
% set the indicator for too many output points, KOP;
% evaluate the initial derivatives;
% set the counter for function evaluations, NFE;
% estimate the starting stepsize.
%
if ( mflag == 1 )
init = 0;
kop = 0;
yp = f ( t, y );
nfe = 1;
if ( t == tout )
flag = 2;
return
end
end
if ( init == 0 )
init = 1;
h = abs ( dt );
toln = 0.0;
for k = 1 : neqn
tol = relerr * abs ( y(k) ) + abserr;
if ( 0.0 < tol )
toln = tol;
ypk = abs ( yp(k) );
if ( tol < ypk * h^5 )
h = ( tol / ypk )^0.2;
end
end
end
if ( toln <= 0.0 )
h = 0.0;
end
h = max ( h, 26.0 * eps * max ( abs ( t ), abs ( dt ) ) );
if ( flag < 0 )
flag_save = - 2;
else
flag_save = + 2;
end
end
%
% Set the stepsize for integration in the direction from T to TOUT.
%
if ( dt < 0.0 )
h = - abs ( h );
else
h = + abs ( h );
end
%
% Test to see if too may output points are being requested.
%
if ( 2.0 * abs ( dt ) <= abs ( h ) )
kop = kop + 1;
end
%
% Unnecessary frequency of output.
%
if ( kop == 100 )
kop = 0;
flag = 7;
return
end
%
% If we are too close to the output point, then simply extrapolate and return.
%
if ( abs ( dt ) <= 26.0 * eps * abs ( t ) )
t = tout;
y(1:neqn) = y(1:neqn) + dt * yp(1:neqn);
yp = f ( t, y );
nfe = nfe + 1;
flag = 2;
return
end
%
% Initialize the output point indicator.
%
output = 0;
%
% To avoid premature underflow in the error tolerance function,
% scale the error tolerances.
%
scale = 2.0 / relerr;
ae = scale * abserr;
%
% Step by step integration.
%
while ( true )
hfaild = 0;
%
% Set the smallest allowable stepsize.
%
hmin = 26.0 * eps * abs ( t );
%
% Adjust the stepsize if necessary to hit the output point.
%
% Look ahead two steps to avoid drastic changes in the stepsize and
% thus lessen the impact of output points on the code.
%
dt = tout - t;
if ( 2.0 * abs ( h ) <= abs ( dt ) )
else
%
% Will the next successful step complete the integration to the output point?
%
if ( abs ( dt ) <= abs ( h ) )
output = 1;
h = dt;
else
h = 0.5 * dt;
end
end
%
% Here begins the core integrator for taking a single step.
%
% The tolerances have been scaled to avoid premature underflow in
% computing the error tolerance function ET.
% To avoid problems with zero crossings, relative error is measured
% using the average of the magnitudes of the solution at the
% beginning and end of a step.
% The error estimate formula has been grouped to control loss of
% significance.
%
% To distinguish the various arguments, H is not permitted
% to become smaller than 26 units of roundoff in T.
% Practical limits on the change in the stepsize are enforced to
% smooth the stepsize selection process and to avoid excessive
% chattering on problems having discontinuities.
% To prevent unnecessary failures, the code uses 9/10 the stepsize
% it estimates will succeed.
%
% After a step failure, the stepsize is not allowed to increase for
% the next attempted step. This makes the code more efficient on
% problems having discontinuities and more effective in general
% since local extrapolation is being used and extra caution seems
% warranted.
%
% Test the number of derivative function evaluations.
% If okay, try to advance the integration from T to T+H.
%
while ( true )
%
% Have we done too much work?
%
if ( maxnfe < nfe )
flag = 4;
kflag = 4;
return
end
%
% Advance an approximate solution over one step of length H.
%
[ f1, f2, f3, f4, f5, f6 ] = fehl ( f, neqn, y, t, h, yp );
f1 = f6;
nfe = nfe + 5;
%
% Compute and test allowable tolerances versus local error estimates
% and remove scaling of tolerances. The relative error is
% measured with respect to the average of the magnitudes of the
% solution at the beginning and end of the step.
%
eeoet = 0.0;
for k = 1 : neqn
et = abs ( y(k) ) + abs ( f1(k) ) + ae;
if ( et <= 0.0 )
flag = 5;
return
end
ee = abs ...
( ( -2090.0 * yp(k) ...
+ ( 21970.0 * f3(k) - 15048.0 * f4(k) ) ...
) ...
+ ( 22528.0 * f2(k) - 27360.0 * f5(k) ) ...
);
eeoet = max ( eeoet, ee / et );
end
esttol = abs ( h ) * eeoet * scale / 752400.0;
if ( esttol <= 1.0 )
break
end
%
% Unsuccessful step. Reduce the stepsize, try again.
% The decrease is limited to a factor of 1/10.
%
hfaild = 1;
output = 0;
if ( esttol < 59049.0 )
s = 0.9 / esttol^0.2;
else
s = 0.1;
end
h = s * h;
if ( abs ( h ) < hmin )
flag = 6;
kflag = 6;
return
end
end
%
% We exited the loop because we took a successful step.
% Store the solution for T+H, and evaluate the derivative there.
%
t = t + h;
y(1:neqn) = f1(1:neqn);
yp = f ( t, y );
nfe = nfe + 1;
%
% Choose the next stepsize. The increase is limited to a factor of 5.
% If the step failed, the next stepsize is not allowed to increase.
%
if ( 0.0001889568 < esttol )
s = 0.9 / esttol^0.2;
else
s = 5.0;
end
if ( hfaild )
s = min ( s, 1.0 );
end
if ( h < 0.0 )
h = - max ( s * abs ( h ), hmin );
else
h = + max ( s * abs ( h ), hmin );
end
%
% End of core integrator
%
% Should we take another step?
%
if ( output )
t = tout;
flag = 2;
return
end
if ( flag <= 0 )
break
end
end
%
% One step integration mode.
%
flag = - 2;
return
end