D201: First-order Differential Equations (Gragg-Bulirsch-Stoer)

Author(s): K.S. Kölbig	Library: MATHLIB
Submitter:	Submitted: 15.02.1989
Language: Fortran	Revised: 01.12.1994

Subroutine subprograms RDEQBS and DDEQBS advance the solution of the system of simultaneous first-order differential equations

from a specified value to a specified value of the independent variable x.

On computers other than CDC and Cray, only the double-precision version DDEQBS is available. On CDC and Cray computers, only the single-precision version RDEQBS is available.

Structure:

SUBROUTINE subprograms
User Entry Names : RDEQBS, DDEQBS
Obsolete User Entry Names: DEQBS RDEQBS
Files Referenced : Unit 6
External References: MTLMTR, ABEND, user-supplied SUBROUTINE subprogram

Usage:

For (type REAL), (type DOUBLE PRECISION),

    CALL tDEQBS(N,X1,X2,Y,H0,EPS,SUB,W)

N

(INTEGER) Number n of equations.

X1

(type according to t) Initial value of the independent variable x.

X2

(type according to t) Final value of the independent variable x.

Y

(type according to t) One-dimensional array of length N. On entry, , must contain . On exit, , contains approximate values .

H0

(type according to t) On entry, H0 must contain the proposed initial step-length . On exit, H0 contains the last computed step-length (See also Method and Notes).

EPS

(type according to t) The requested absolute accuracy . (EPS should not be smaller than approximately times the machine precision).

SUB

Name of a user-supplied SUBROUTINE subprogram, declared EXTERNAL in the calling program.

W

(type according to t) Array containing at least 36*N elements required as working-space.

The user-supplied subroutine SUB should be of the form

    SUBROUTINE SUB(X,Y,F)

where the variable X and the one-dimensional arrays Y(*) and F(*) are of type t. This subroutine must set

Method:

For the first integration step, starting at , the step-length h is chosen to be the smallest of the numbers for which not more than 9 stages of internal extrapolation yield an estimated error less than . This procedure is repeated until is reached. (For details, see Ref. 1).

Error handling:

Error D201.1: If the requestec accuracy cannot obtained, a message is written on Unit 6, unless subroutine MTLSET (N002) has been called.
For , or or , control is returned to the calling program without any change in Y.

Notes:

For well-conditioned systems of equations any reasonable value of the initial step length may be chosen. For ill-conditioned systems, the initial value of may be important, and tests with different values are advised. An inappropriate choice may lead to wrong results in such cases.

This subroutines is based on an Algol60 procedure given in Ref. 1. The adaption for integration over a given interval (not only over one step) is due to G. Janin.

References:

1. R. Bulirsch and J. Stoer, Numerical treatment of ordinary differential equations by extrapolation methods, Numer. Math. 8 (1966) 1-13.