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

Subroutine subprograms RDEQBS and DDEQBS advance the solution of the system of $n\; \ge 1$ simultaneous first-order differential equations

$\{dy$_idx} = f_i(x,y₁,...,y_n), (i = 1,2,...,n),

from a specified value $x$₁ to a specified value $x$₂ 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 $\equiv$ RDEQBS
Files Referenced : Unit 6
External References: MTLMTR (N002), ABEND (Z035), user-supplied SUBROUTINE subprogram

Usage:

For $t=R$ (type REAL), $t=D$ (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 $x$₁ of the independent variable x.

X2

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

Y

(type according to t) One-dimensional array of length $\ge$ N. On entry, $Y(i),(i=1,...,N)$ , must contain $y$_i(x₁) . On exit, $Y(i),(i=1,...,N)$ , contains approximate values $y$_i(x₂) .

H0

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

EPS

(type according to t) The requested absolute accuracy $\epsilon$ . ( EPS should not be smaller than approximately $103$ 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.

    SUBROUTINE SUB(X,Y,F)

$F(I)=\; f$_I(X,Y(1),...,Y(N))(I = 1,2,...,N).

Method:

For the first integration step, starting at $x=x$₁ , the step-length h is chosen to be the smallest of the numbers $h$₀, h₀/2, h₀/4,... for which not more than 9 stages of internal extrapolation yield an estimated error less than $\epsilon$ . This procedure is repeated until $x=x$₂ 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 $X1=X2$ or $H0=0$ , 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 $h$₀ may be chosen. For ill-conditioned systems, the initial value of $h$₀ 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.

D202

K.S. Kölbig