{
Ordinary differential equations
author: Nikolai Shokhirev http://www.shokhirev.com/nikolai.html)
created: 2002.02.02
last modified: 2004.12.24
ŠNikolai V. Shokhirev, 2001-2004 }
unit uODE;
interface
uses
uMatTypes, uDynArrays;
type
TDiffEqProc = procedure(const t: TFloat; const y, yp: IFArr1D);
{
Pascal version by Nikolai V. Shokhirev
PURPOSE:
Integrate a system of neqn first order differential equations by the Fehlberg
fourth-fifth order Runge-Kutta method.
Ref.
George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler.
"Computer Methods for Mathematical Computations". Prentice-Hall, 1977.
ORIGINAL AUTHORS:
RKF45 - H.A. Watts and L.F. Shampine (Sandia)
NOTES:
RKF45 is primarily designed to solve non-stiff and mildly stiff differential
equations when derivative evaluations are inexpensive.
Rkf45 should generally not be used when the user is demanding high accuracy.
The arrays yp, f1, f2, f3, f4, and f5 (of dimension at least neqn) and the
variables h, savre, savae, nfe, kop, init,jflag,and kflag are used internally
by the code and appear in the call list to eliminate local retention of
variables between calls. Accordingly, they should not be altered. Items of
possible interest are
yp - derivative of solution vector at t
h - an appropriate stepsize to be used for the next step
nfe- counter on the number of derivative function evaluations
ABSTRACT:
Subroutine Rkf45 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[i] 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 Double 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 Rkf45 again (and perhaps
define a new value for tout). Rkfs45 calls subroutine Fehl which computes an
approximate solution over one step.
Rkf45 uses the Runge-Kutta-Fehlberg(4,5) method described in the reference
E. Fehlberg: Low-order Classical Runge-Kutta formulas with Stepsize
Control, NASA TR R-315.
The performance of Rkf45 is illustrated in the reference
L.F.Shampine,H.A.Watts,S.Davenport: "Solving Non-stiff Ordinary
Differential Equations - The State of the Art.
Sandia Laboratories Report SAND75-0182, to appear in SIAM Review.
First call to Rkf45:
The user must provide storage in his calling program for the arrays in
the call list y, yp, f1, f2, f3, f4, f5 - 1D float arrays with the same
limits and dimension = neqn, supply procedure F(t,y,yp) and initialize the
following parameters:
y - vector of initial conditions
t - starting point of integration , must be a variable
tout - output point at which solution is desired.
t = tout is allowed on the first call only, in which case Rkf45 returns with
iflag=2 if continuation is possible.
relerr, abserr - relative and absolute local error tolerances which must be
non-negative. relerr must be a variable while abserr may be a constant.
The code should normally not be used with relative error control smaller
than about i.e-8 to avoid limiting precision difficulties the code
requires relerr to be larger than an internally computed relative error
parameter which is machine dependent. In particular, pure absolute error
is not permitted. If a smaller than allowable value of relerr is
attempted, rkf45 increases relerr appropriately and returns control to the
user before continuing the integration.
iflag - +1,-1 indicator to initialize the code for each new problem.
Normal input is +1. the user should set iflag = -1 only when one-step
integrator control is essential. In this case, Rkf45 attempts to advance
the solution a Double step in the direction of tout each time it is
called. Since this mode of operation results in extra computing overhead,
it should be avoided unless needed.
Output from Rkf45
y - solution at t
t - last point reached in integration.
iflag = 2 - integration reached tout. Indicates successful return and is the
normal mode for continuing integration.
= -2 - a Double successful step in the direction of tout has been taken.
Normal mode for continuing integration one step at a time.
= 3 - integration was not completed because relative error tolerance
was too small. relerr has been increased appropriately for
continuing.
= 4 - integration was not completed because more than 3000 derivative
evaluations were needed. This is approximately 500 steps.
= 5 - integration was not completed because solution vanished making
a pure relative error test impossible. Must use 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 requested accuracy could
not be achieved using smallest allowable stepsize. User must
increase the error tolerance before continued integration can
be attempted.
= 7 - it is likely that Rkf45 is inefficient for solving this problem.
Too much output is restricting the natural stepsize choice.
Use the one-step integrator mode.
= 8 - invalid input parameters. This indicator occurs if any of the
following is satisfied:
t = tout and iflag <> +1 or -1
relerr or abserr < 0.
iflag = 0 or < -2 or > 8
h, nfe, yp, f1, f2, f3, f4, f5 - information which is usually of no interest
to the user but necessary for subsequent calls.
- yp contain the first derivatives of the solution vector y at t.
- h contains the stepsize to be attempted on the next step.
- nfe contains the derivative evaluation counter.
Subsequent calls to Rkf45
Subroutine Rkf45 returns with all information needed to continue the
integration. If the integration reached tout, the user need only define
a new tout and call Rkf45 again. In the one-step integrator mode (iflag=-2)
the user must keep in mind that each step taken is in the direction of the
current tout. Upon reaching tout (indicated by changing iflag to 2),the
user must then define a new tout and reset iflag to -2 to continue in the
one-step integrator mode.
If the integration was not completed but the user still wants to continue
(iflag = 3,4 cases), he just calls Rkf45 again. With iflag = 3, the relerr
parameter has been adjusted appropriately for continuing the integration.
In the case of iflag=4 the function counter will be reset to 0 and
another 3000 function evaluations are allowed.
However,in the case iflag = 5, the user must first alter the error criterion
to use a positive value of abserr before integration can proceed. If he
does not, execution is terminated.
Also, in the case iflag = 6, it is necessary for the user to reset iflag
to 2 (or -2 when the one-step integration mode is being used) as well as
increasing either abserr,relerr or both before the integration can be
continued. If this is not done, execution will be terminated. The
occurrence of iflag=6 indicates a trouble spot (solution is changing
rapidly,singularity may be present) and it often is inadvisable to continue.
If (flag = 7 is encountered, the user should use the one-step integration mode
with the stepsize determined by the code or consider switching to the
Adams codes DE/STEP, INTRP. If the user insists upon continuing the
integration with Rkf45, he must reset iflag to 2 before calling Rkf45 again.
Otherwise, execution will be terminated.
If iflag = 8 is obtained, integration can not be continued unless the invalid
input parameters are corrected.
It should be noted that the arrays work,iwork contain information required
for subsequent integration. Accordingly, work and iwork should not be
altered.
This interfacing routine merely relieves the user of a long calling list
via the splitting apart of two working storage arrays.
INPUT
neqn - number of equations
y - solution vector at t
t - independent variable
tout - output point at which solution is desired
relerr - relative error tolerance
abserr - absolute error tolerance
OUTPUT
iflag - indicator for status of work
h - step size
savre;
savae;
nfe - number of function evaluations
kop;
IN, OUT
yp - derivatives
f1, f2, f3, f4, f5
init;
jflag;
kflag; }
procedure Rkfs45(F: TDiffEqProc; var t: TFloat; const tout, abserr: TFloat;
var relerr: TFloat; var iflag: TInt;
const y, yp, f1, f2, f3, f4, f5: IFArr1D;
var h, savre, savae: TFloat;
var nfe, kop, init, jflag, kflag: TInt);
{
Pascal version by Nikolai V. Shokhirev
Ref.
George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler.
"Computer Methods for Mathematical Computations". Prentice-Hall, 1977.
PURPOSE:
Integrate a system of neqn first-order ordinary differential equations
of the form dy/dt=F(t,y) [y is a vector], where the initial values of y
and the initial values of yp, the derivatives are specified at the starting
point t. Fehl advances the solution over the fixed step h and returns the
fifth order (sixth-order locally) solution approximation at t+h in array s.
neqn: TInt; // number of equations
INPUT
t - starting point
h - step size
y - array of length neqn; function at t
yp - array of length neqn; derivatives at t
OUTPUT
f1 - array of length neqn for internal use
f2 - array of length neqn for internal use
f3 - array of length neqn for internal use
f4 - array of length neqn for internal use
f5 - array of length neqn for internal use
s - array of length neqn; the results }
procedure Fehl(F: TDiffEqProc; const t, h: TFloat;
const y, yp, f1, f2, f3, f4, f5, s : IFarr1D);
implementation
uses
math;
initialization
end.
Back
Generated by Lore's Source to HTML Converter(http://www.newty.de/lsc/index.html)