IVMRK

Solves an initial-value problem yʹ = f(t, y) for ordinary differential equations using Runge-Kutta pairs of various orders.

Required Arguments

IDO — Flag indicating the state of the computation. (Input/Output)

IDO	State
1	Initial entry
2	Normal re-entry
3	Final call to release workspace
4	Return after a step
5	Return for function evaluation (reverse communication)

Normally, the initial call is made with IDO = 1. The routine then sets IDO = 2, and this value is used for all but the last call that is made with IDO = 3. This final call is used to release workspace, which was automatically allocated by the initial call with IDO = 1.

FCN — User-supplied subroutine to evaluate functions. The usage is
CALL FCN (N, T, Y, YPRIME), where

N — Number of equations. (Input)

T — Independent variable. (Input)

Y — Array of size N containing the dependent variable values, y. (Input)

YPRIME — Array of size N containing the values of the vector yʹ evaluated at (t, y). (Output)

FCN must be declared EXTERNAL in the calling program.

T — Independent variable. (Input/Output)
On input, T contains the initial value. On output, T is replaced by TEND unless error conditions have occurred.

TEND — Value of t where the solution is required. (Input)
The value of TEND may be less than the initial value of t.

Y — Array of size N of dependent variables. (Input/Output)
On input, Y contains the initial values. On output, Y contains the approximate solution.

YPRIME — Array of size N containing the values of the vector y' evaluated at (t, y). (Output)

Optional Arguments

N — Number of differential equations. (Input)
Default: N= size (Y,1).

FORTRAN 90 Interface

Generic: CALL IVMRK (IDO, FCN, T, TEND, Y, YPRIME [, …])

Specific: The specific interface names are S_IVMRK and D_IVMRK.

FORTRAN 77 Interface

Single: CALL IVMRK (IDO, N, FCN, T, TEND, Y, YPRIME)

Double: The double precision name is DIVMRK.

Description

Routine IVMRK finds an approximation to the solution of a system of first-order differential equations of the form yʹ = f(t, y) with given initial data. Relative local error is controlled according to a user-supplied tolerance. For added efficiency, three Runge-Kutta formula pairs, of orders 3, 5, and 8, are available.

Optionally, the values of the vector yʹ can be passed to IVMRK by reverse communication, avoiding the user-supplied subroutine FCN. Reverse communication is especially useful in applications that have complicated algorithmic requirement for the evaluations of f(t, y). Another option allows assessment of the global error in the integration.

The routine IVMRK is based on the codes contained in RKSUITE, developed by R. W. Brankin, I. Gladwell, and L. F. Shampine (1991).

Comments

1. Workspace may be explicitly provided, if desired, by use of I2MRK/DI2MRK. The reference is:

CALL I2MRK (IDO, N, FCN, T, TEND, Y, YPRIME, TOL, THRES, PARAM, YMAX, RMSERR, WORK, IWORK)

The additional arguments are as follows:

TOL — Tolerance for error control. (Input)

THRES — Array of size N. (Input)
THRES(I) is a threshold for solution component Y(I). It is chosen so that the value of Y(L) is not important when Y(L) is smaller in magnitude than THRES(L). THRES(L) must be greater than or equal to sqrt(amach(4)).

PARAM — A floating-point array of size 50 containing optional parameters. (Input/Output)
If a parameter is zero, then a default value is used. These default values are given below. The following parameters must be set by the user:

	PARAM	Definition
1	HINIT	Initial value of the step size. Must be chosen such that 0.01 ≥ HINIT ≥ 10.0 amach(4). Default: automatic selection of stepsize
2	METHOD	1 - use the (2, 3) pair 2 - use the (4, 5) pair 3 - use the (7, 8) pair. Default: METHOD = 1 if 1.e-2 ≥ tol > 1.e-4 METHOD = 2 if 1.e-4 ≥ tol > 1.e-6 METHOD = 3 if 1.e-6 ≥ tol
3	ERREST	ERREST = 1 attempts to assess the true error, the difference between the numerical solution and the true solution. The cost of this is roughly twice the cost of the integration itself with METHOD = 2 or METHOD = 3, and three times with METHOD = 1. Default: ERREST = 0.
4	INTRP	If nonzero, then return the IDO = 4 before each step. See Comment 3. Default: 0.
5	RCSTAT	If nonzero, then reverse communication is used to get derivative information. See Comment 4. Default: 0.
6 - 30		Not used.
	The following entries are set by the program:
31	HTRIAL	Current trial step size.
32	NSTEP	Number of steps taken.
33	NFCN	Number of function evaluations.
34	ERRMAX	The maximum approximate weighted true error taken over all solution components and all steps from T
35	TERRMX	First value of the independent variable where an approximate true error attains the maximum value ERRMAX.

YMAX — Array of size N, where YMAX(L) is the largest value of ABS(Y(L)) computed at any step in the integration so far.

RMSERR — Array of size N where RMSERR(L) approximates the RMS average of the true error of the numerical solution for the L-th solution component, L = 1, …, N. The average is taken over all steps from T through the current integration point. RMSERR is accessed and set only if PARAM(3) = 1.

WORK — Floating point work array of size 39N using the working precision. The contents of WORK must not be changed from the first call with IDO = 1 until after the final call with IDO = 3.

IWORK — Length of array work. (Input)

2. Informational errors

Type	Code	Description
4	1	It does not appear possible to achieve the accuracy specified by TOL and THRES(*) using the current precision and METHOD. A larger value for METHOD, if possible, will permit greater accuracy with this precision. The integration must be restarted.
4	2	The global error assessment may not be reliable beyond the current integration point T. This may occur because either too little or too much accuracy has been requested or because f(t, y) is not smooth enough for values of t just past TEND and current values of the solution y. This return does not mean that you cannot integrate past TEND, rather that you cannot do it with PARAM(3) = 1.

3. If PARAM(4) is nonzero, the subroutine returns with IDO = 4 and will resume calculation at the point of interruption if re-entered with IDO = 4. Some parameters the user might want to examine are IDO, HTRIAL, NSTEP, NFCN, T, and Y. The array Y contains the newly computed trial value for y(t), accepted or not.

4. If PARAM(5) is nonzero, the subroutine will return with IDO = 5. At this time, evaluate the derivatives at T, place the result in YPRIME, and call IVMRK again. The dummy function I40RK/DI40RK may be used in place of FCN.

Examples

Example 1

This example integrates the small system (A.2.B2) from the test set of Enright and Pryce (1987):

USE IVMRK_INT

USE WRRRN_INT

IMPLICIT NONE

INTEGER N

PARAMETER (N=3)

! Specifications for local variables

INTEGER IDO

REAL T, TEND, Y(N), YPRIME(N)

EXTERNAL FCN

! Set initial conditions

T = 0.0

TEND = 20.0

Y(1) = 2.0

Y(2) = 0.0

Y(3) = 1.0

IDO = 1

CALL IVMRK (IDO, FCN, T, TEND, Y, YPRIME)

! Final call to release workspace

IDO = 3

CALL IVMRK (IDO, FCN, T, TEND, Y, YPRIME)

CALL WRRRN ('Y', Y)

END

SUBROUTINE FCN (N, T, Y, YPRIME)

! Specifications for arguments

INTEGER N

REAL T, Y(*), YPRIME(*)

YPRIME(1) = -Y(1) + Y(2)

YPRIME(2) = Y(1) - 2.0*Y(2) + Y(3)

YPRIME(3) = Y(2) - Y(3)

RETURN

END

Output

1 1.000

2 1.000

3 1.000

Example 2

This problem is the same mildly stiff problem (A.1.F2) from the test set of Enright and Pryce as Example 2 for IVPRK.

Although not a stiff solver, one notes the greater efficiency of IVMRK over IVPRK, in terms of derivative evaluations. Reverse communication is also used in this example. Users will find this feature particularly helpful if their derivative evaluation scheme is difficult to isolate in a separate subroutine.

USE I2MRK_INT

USE UMACH_INT

USE AMACH_INT

IMPLICIT NONE

INTEGER N

PARAMETER (N=2)

! Specifications for local variables

INTEGER IDO, ISTEP, LWORK, NOUT

REAL PARAM(50), PREC, RMSERR(N), T, TEND, THRES(N), TOL, &

WORK(1000), Y(N), YMAX(N), YPRIME(N)

REAL AK1, AK2, AK3

SAVE AK1, AK2, AK3

! Specifications for intrinsics

INTRINSIC SQRT

REAL SQRT

! Specifications for subroutines

EXTERNAL I40RK

! Specifications for functions

DATA AK1, AK2, AK3/294.0, 3.0, 0.01020408/

CALL UMACH (2, NOUT)

! Set initial conditions

T = 0.0

Y(1) = 1.0

Y(2) = 0.0

! Set tolerance for error control,

! threshold vector and parameter

! vector

TOL = .001

PREC = AMACH(4)

THRES = SQRT (PREC)

PARAM = 0.0E0

LWORK = 1000

! Turn on derivative evaluation by

! reverse communication

PARAM(5) = 1

IDO = 1

ISTEP = 24

! Print header

WRITE (NOUT,99998)

10 CONTINUE

TEND = ISTEP

CALL I2MRK (IDO, N, I40RK, T, TEND, Y, YPRIME, TOL, THRES, PARAM,&

YMAX, RMSERR, WORK, LWORK)

IF (IDO .EQ. 5) THEN

! Evaluate derivatives

YPRIME(1) = -Y(1) - Y(1)*Y(2) + AK1*Y(2)

YPRIME(2) = -AK2*Y(2) + AK3*(1.0-Y(2))*Y(1)

GO TO 10

ELSE IF (ISTEP .LE. 240) THEN

! Integrate to 10 equally spaced points

WRITE (NOUT,'(I6,3F12.3)') ISTEP/24, T, Y

IF (ISTEP .EQ. 240) IDO = 3

ISTEP = ISTEP + 24

GO TO 10

END IF

! Show number of derivative evaluations

WRITE (NOUT,99999) PARAM(33)

99998 FORMAT (3X, 'ISTEP', 5X, 'TIME', 9X, 'Y1', 10X, 'Y2')

99999 FORMAT (/, 4X, 'NUMBER OF DERIVATIVE EVALUATIONS WITH IVMRK =', &

F6.0)

END

! DUMMY FUNCTION TO TAKE THE PLACE OF DERIVATIVE EVALUATOR

SUBROUTINE I40RK (N, T, Y, YPRIME)

INTEGER N

REAL T, y(*), YPRIME(*)

RETURN

END

Output

ISTEP TIME Y1 Y2

1 24.000 0.688 0.002

2 48.000 0.634 0.002

3 72.000 0.589 0.002

4 96.000 0.549 0.002

5 120.000 0.514 0.002

6 144.000 0.484 0.002

7 168.000 0.457 0.002

8 192.000 0.433 0.001

9 216.000 0.411 0.001

10 240.000 0.391 0.001

NUMBER OF DERIVATIVE EVALUATIONS WITH IVMRK = 1375.

Example 3

This example demonstrates how exceptions may be handled. The problem is from Enright and Pryce (A.2.F1), and has discontinuities. We choose this problem to force a failure in the global error estimation scheme, which requires some smoothness in y. We also request an initial relative error tolerance which happens to be unsuitably small in this precision.

If the integration fails because of problems in global error assessment, the assessment option is turned off, and the integration is restarted. If the integration fails because the requested accuracy is not achievable, the tolerance is increased, and global error assessment is requested. The reason error assessment is turned on is that prior assessment failures may have been due more in part to an overly stringent tolerance than lack of smoothness in the derivatives.

When the integration is successful, the example prints the final relative error tolerance, and indicates whether or not global error estimation was possible.

USE IMSL_LIBRARIES

IMPLICIT NONE

INTEGER N

PARAMETER (N=2)

! Specifications for local variables

INTEGER IDO, LWORK, NOUT

REAL PARAM(50), PREC, RMSERR(N), T, TEND, THRES(N), TOL,&

WORK(100), Y(N), YMAX(N), YPRIME(N)

! Specifications for intrinsics

INTRINSIC SQRT

REAL SQRT

! Specifications for subroutines

! Specifications for functions

EXTERNAL FCN

CALL UMACH (2, NOUT)

! Turn off stopping for FATAL errors

CALL ERSET (4, -1, 0)

! Initialize input, turn on global

! error assessment

LWORK = 100

PREC = AMACH(4)

TOL = SQRT(PREC)

PARAM = 0.0E01

THRES = TOL

TEND = 20.0E0

PARAM(3) = 1

10 CONTINUE

! Set initial values

T = 0.0E0

Y(1) = 0.0E0

Y(2) = 0.0E0

IDO = 1

CALL I2MRK (IDO, N, FCN, T, TEND, Y, YPRIME, TOL, THRES, PARAM,&

YMAX, RMSERR, WORK, LWORK)

IF (IERCD() .EQ. 32) THEN

! Unable to achieve requested

! accuracy, so increase tolerance.

! Activate global error assessment

TOL = 10.0*TOL

PARAM(3) = 1

WRITE (NOUT,99995) TOL

GO TO 10

ELSE IF (IERCD() .EQ. 34) THEN

! Global error assessment has failed,

! cannot continue from this point,

! so restart integration

WRITE (NOUT,99996)

PARAM(3) = 0

GO TO 10

END IF

! Final call to release workspace

IDO = 3

CALL I2MRK (IDO, N, FCN, T, TEND, Y, YPRIME, TOL, THRES, PARAM,&

YMAX, RMSERR, WORK, LWORK)

! Summarize status

WRITE (NOUT,99997) TOL

IF (PARAM(3) .EQ. 1) THEN

WRITE (NOUT,99998)

ELSE

WRITE (NOUT,99999)

END IF

CALL WRRRN ('Y', Y)

99995 FORMAT (/, 'CHANGING TOLERANCE TO ', E9.3, ' AND RESTARTING ...'&

, /, 'ALSO (RE)ENABLING GLOBAL ERROR ASSESSMENT', /)

99996 FORMAT (/, 'DISABLING GLOBAL ERROR ASSESSMENT AND RESTARTING ...'&

, /)

99997 FORMAT (/, 72('-'), //, 'SOLUTION OBTAINED WITH TOLERANCE = ',&

E9.3)

99998 FORMAT ('GLOBAL ERROR ASSESSMENT IS AVAILABLE')

99999 FORMAT ('GLOBAL ERROR ASSESSMENT IS NOT AVAILABLE')

END

SUBROUTINE FCN (N, T, Y, YPRIME)

USE CONST_INT

! Specifications for arguments

INTEGER N

REAL T, Y(*), YPRIME(*)

! Specifications for local variables

REAL A

REAL PI

LOGICAL FIRST

SAVE FIRST, PI

! Specifications for intrinsics

INTRINSIC INT, MOD

INTEGER INT, MOD

! Specifications for functions

DATA FIRST/.TRUE./

IF (FIRST) THEN

PI = CONST('PI')

FIRST = .FALSE.

END IF

A = 0.1E0

YPRIME(1) = Y(2)

IF (MOD(INT(T),2) .EQ. 0) THEN

YPRIME(2) = 2.0E0*A*Y(2) - (PI*PI+A*A)*Y(1) + 1.0E0

ELSE

YPRIME(2) = 2.0E0*A*Y(2) - (PI*PI+A*A)*Y(1) - 1.0E0

END IF

RETURN

END

Output

*** FATAL ERROR 34 from i2mrk. The global error assessment may not

*** be reliable for T past 9.994749E-01. The integration is

*** being terminated.

DISABLING GLOBAL ERROR ASSESSMENT AND RESTARTING ...

*** FATAL ERROR 32 from i2mrk. In order to satisfy the error

*** requirement I6MRK would have to use a step size of

*** 3.647129E- 06 at TNOW = 9.999932E-01. This is too small

*** for the current precision.

CHANGING TOLERANCE TO 0.345E-02 AND RESTARTING ...

ALSO (RE)ENABLING GLOBAL ERROR ASSESSMENT

*** FATAL ERROR 34 from i2mrk. The global error assessment may

*** not be reliable for T past 9.986024E-01. The integration

*** is being terminated.

DISABLING GLOBAL ERROR ASSESSMENT AND RESTARTING ...

------------------------------------------------------------------------

SOLUTION OBTAINED WITH TOLERANCE = 0.345E-02

GLOBAL ERROR ASSESSMENT IS NOT AVAILABLE

1 -12.30

2 0.95