D202 First-order Differential Equations (Runge--Kutta--Merson)

Subroutine subprograms RDEQMR and DDEQMR 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. The integration step-length is automatically adjusted to keep the estimated error per step less than a specified value.

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 : RDEQMR, DDEQMR
Obsolete User Entry Names: DEQMR $\equiv$ RDEQMR
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 tDEQMR(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 6*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:

The method is a modification by Merson of the Runge--Kutta method. The initial value of the step-length h is taken to be the first of the numbers $h$₀, h₀/2, h₀/4, ... for which the estimated relative error at the end of the step is less than $\epsilon$ . At each susequent step, an estimate of the integration error for that step (proportional to $h5$ ) is computed. If the corresponding relative error exceeds $\epsilon$ , the current step-length is halfed; if it is less than $\epsilon /32$ the step-length is doubled. This process is continued until $x$₂ is reached. (For details, see Ref. 1).

Error handling:

Error D202.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.

References:

1. G.N. Lance, Numerical methods for high-speed computers, (Iliffe & Sons, London 1960) 56

D203

K.S. Kölbig