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
Y
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
Y
1 -12.30
2 0.95
