IVOAM

Solves an initial-value problem for a system of ordinary differential equations of order one or two using a variable order Adams method.

Required Arguments

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

IDO | State |
---|---|

1 | Initial entry input value. |

2 | Normal re-entry input value. On output, if IDO = 2 then the integration is finished. If the integrator is called with a new value for TEND, the integration continues. If the integrator is called with TEND unchanged, an error message is issued. |

3 | Input value to use on final call to release workspace. |

>3 | Output value that indicates that a fatal error has occurred. |

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 only 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 (IDO, T, Y, HIDRVS [, …]), where

The usage is CALL FCN (IDO, T, Y, HIDRVS [, …]), where

Required Arguments

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

This flag corresponds to the IDO argument described above. If FCN has complicated subexpressions, which depend only weakly or not at all on Y then these subexpressions need only be computed when IDO = 1 and their values then reused when IDO = 2.

This flag corresponds to the IDO argument described above. If FCN has complicated subexpressions, which depend only weakly or not at all on Y then these subexpressions need only be computed when IDO = 1 and their values then reused when IDO = 2.

T — Independent variable, t. (Input)

Y — Array of length k containing the dependent variable values, y, and first derivatives, if any. k will be the sum of the orders of the equations in the system of equations to solve. (Input)

HIDRVS — Array of length n = NEQ, where n is the number of equations in the system to solve, containing the values of the highest order derivatives evaluated at (t, y). (Output)

IVOAM uses size(HIDRVS) to set the default value of NEQ unless the optional argument NEQ is used.

Optional Arguments

FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied subroutine. (Input/Output)

For a detailed description of this argument see FCN_DATA below.

FCN must be declared EXTERNAL in the calling program.

T — Independent variable, t. (Input/Output)

On input, T contains the initial independent variable value. On output, T is replaced by TEND unless error conditions arise. See IDO for details. (Input/Output)

On input, T contains the initial independent variable value. On output, T is replaced by TEND unless error conditions arise. See IDO for details. (Input/Output)

TEND — Value of t = tend where the solution is required. (Input)

Y — Array of length k containing the dependent variables, y(t), and first derivatives, if any. (Input/Output)

k will be the sum of the orders of the equations in the system of equations to solve. On input, Y contains the initial values, y(t0) and y’(t0) (if needed). On output, Y contains the approximate solution, y(t). For example, for a system of first order equations, Y(i) is the i-th dependent variable. For a system of second order equations, Y(2i-1) is the i-th dependent variable and Y(2i) is the derivative of the i-th dependent variable. For systems of equations in which one or more equations is of order 2, optional argument KORDER must be used to denote the order of each equation so that the derivatives in Y can be identified. By default it is assumed that all equations are of order 1 and Y contains only dependent variables.

k will be the sum of the orders of the equations in the system of equations to solve. On input, Y contains the initial values, y(t0) and y’(t0) (if needed). On output, Y contains the approximate solution, y(t). For example, for a system of first order equations, Y(i) is the i-th dependent variable. For a system of second order equations, Y(2i-1) is the i-th dependent variable and Y(2i) is the derivative of the i-th dependent variable. For systems of equations in which one or more equations is of order 2, optional argument KORDER must be used to denote the order of each equation so that the derivatives in Y can be identified. By default it is assumed that all equations are of order 1 and Y contains only dependent variables.

HIDRVS — Array of length n = NEQ, where n is the number of equations in the system to solve, containing the highest order derivatives at the point Y. (Output)

IVOAM uses size(HIDRVS) to set the default value of NEQ unless the optional argument NEQ is used.

IVOAM uses size(HIDRVS) to set the default value of NEQ unless the optional argument NEQ is used.

Optional Arguments

NEQ — Number of differential equations in the system of equations to solve. (Input)

Default: NEQ = size (HIDRVS).

Default: NEQ = size (HIDRVS).

KORDER — An array of length NEQ specifying the orders of the equations in the system of equations to solve. The elements of KORDER can be 1 or 2. KORDER must be used with argument Y to define systems of mixed or higher order. (Input)

Default: KORDER = (1,1,1 …,1).

Default: KORDER = (1,1,1 …,1).

EQNERR — An array of length NEQ specifying the error tolerance for each equation. (Input)

Let e(i) be the error tolerance for equation i. Then

Let e(i) be the error tolerance for equation i. Then

Value | Explanation |
---|---|

e(i) > 0 | Implies an absolute error tolerance of e(i) is to be used for equation i. |

e(i) = 0 | Implies that the default absolute error tolerance (defined below) is to be used for equation i. |

e(i) < 0 | Implies a relative error test is to be performed for equation i. In this case, the base error tolerance used will be |e(i)| and the relative error factor used will be (15/16 * ∣e(i)∣). Thus the actual absolute error tolerance used will be ∣e(i)∣ * (15/16 * ∣e(i)∣). |

Default: An absolute error tolerance of 1.E-5 is used for single precision and 1.D-10 for double precision for all equations.

HINC — Factor used for increasing the stepsize. (Input)

One should set HINC such that 9/8 <= HINC <= 4.

Default: HINC = 2.0.

One should set HINC such that 9/8 <= HINC <= 4.

Default: HINC = 2.0.

HDEC — Factor used for decreasing the stepsize. (Input)

One should set HDEC such that 1/4 <= HDEC <= 7/8.

Default: HDEC = 0.5.

One should set HDEC such that 1/4 <= HDEC <= 7/8.

Default: HDEC = 0.5.

HMIN — Absolute value of the minimum stepsize permitted. (Input)

Default: HMIN = 10.0/amach(2) for single precision and 10.0/dmach(2) for double precision.

Default: HMIN = 10.0/amach(2) for single precision and 10.0/dmach(2) for double precision.

HMAX — Absolute value of the maximum stepsize permitted. (Input)

Default: HMAX = amach(2) for single precision and dmach(2) for double precision.

Default: HMAX = amach(2) for single precision and dmach(2) for double precision.

FCN_DATA – A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied subroutine. (Input/Output)

The derived type, s_fcn_data, is defined as:

The derived type, s_fcn_data, is defined as:

type s_fcn_data

real(kind(1e0)), pointer, dimension(:) :: rdata

integer, pointer, dimension(:) :: idata

end type

in module mp_types. The double precision counterpart to s_fcn_data is named d_fcn_data. The user must include a use mp_types statement in the calling program to define this derived type.

Note that if this optional argument is present then FCN_DATA must also be defined as an optional argument in the user-supplied subroutine.

Fortran 90 Interface

Generic: CALL IVOAM (IDO, FCN, T, TEND, Y, HIDRVS [, …])

Specific: The specific interface names are S_IVOAM and D_IVOAM.

Description

Routine IVOAM is based on the JPL Library routine SIVA. IVOAM uses a variable order Adams method to solve the initial value problem

or more generally

where y is the vector

zi(k) is the kth derivative of zi with respect to t, di is the order of the ith differential equation, and η is a vector with the same dimension as y.

Note that the systems of equations solved by IVOAM can be of order one, order two, or mixed order one and two.

Comments

Informational errors

Type | Code | Description |
---|---|---|

3 | 1 | The requested error tolerance is too small. |

3 | 2 | The stepsize has been reduced too rapidly. The integrator is going to do a restart. |

Examples

Example 1

In this example a system of two equations of order two is solved.

The initial conditions are

Since the system is of order two, optional argument KORDER must be used to specify the orders of the equations. Also, because the system is of order two, Y(1) contains the first dependent variable, Y(2) contains the derivative of the first dependent variable, Y(3) contains the second dependent variable, and Y(4) contains the derivative of the second dependent variable.

USE IVOAM_INT

USE UMACH_INT

USE CONST_INT

IMPLICIT NONE

INTEGER IDO, IEND, NOUT, KORDER(2)

REAL T, TEND, Y(4), HIDRVS(2), DELTA

EXTERNAL FCN

! Initialize

IDO = 1

T = 0.0

Y(1) = 1.0

Y(2) = 0.0

Y(3) = 0.0

Y(4) = 1.0

KORDER = 2

! Write title

CALL UMACH (2, NOUT)

WRITE (NOUT,99997)

! Integrate ODE

IEND = 0

DELTA = CONST('PI')

DELTA = 2.0*DELTA

DO

IEND = IEND + 1

TEND = T + DELTA

IF(TEND .GT. 20.0) TEND = 20.0

CALL IVOAM (IDO, FCN, T, TEND, Y, HIDRVS, KORDER=KORDER)

IF (IEND .LE. 4) THEN

WRITE (NOUT,99998) T, Y(1), Y(2), HIDRVS(1)

WRITE (NOUT,99999) Y(3), Y(4), HIDRVS(2)

! Finish up

IF (IEND .EQ. 4) IDO = 3

CYCLE

END IF

EXIT

END DO

99997 FORMAT (11X, 'T', 12X, 'Y1/Y2', 9X, 'Y1P/Y2P', 7X, 'Y1PP/Y2PP')

99998 FORMAT (4F15.4)

99999 FORMAT (15X, 3F15.4)

END

SUBROUTINE FCN (IDO, T, Y, HIDRVS)

INTEGER IDO

REAL T, Y(*), HIDRVS(*)

REAL TP

TP = Y(1)*Y(1) + Y(3)*Y(3)

TP = 1.0E0/(TP*SQRT(TP))

HIDRVS(1) = -Y(1)*TP

HIDRVS(2) = -Y(3)*TP

RETURN

END

Output

T Y1/Y2 Y1P/Y2P Y1PP/Y2PP

6.2832 1.0000 -0.0000 -1.0000

0.0000 1.0000 0.0000

12.5664 1.0000 -0.0000 -1.0000

0.0000 1.0000 -0.0000

18.8496 1.0000 -0.0000 -1.0000

0.0000 1.0000 -0.0000

20.0000 0.4081 -0.9129 -0.4081

0.9129 0.4081 -0.9129

Example 2

This contrived example illustrates how to use IVOAM to solve a system of equations of mixed order.

The height, y(t), of an object of mass m above the surface of the Earth can be modelled using Newton's second law as:

or

where -mg is the downward force of gravity and -ky′ is the force due to air resistance, in a direction opposing the velocity. If the object is a meteor, the mass, m, and air resistance, k, will decrease as the meteor burns up in the atmosphere. The mass is proportional to r3 (r = radius) and the air resistance, presumably dependent on the surface area, may be assumed to be proportional to r2, so that k/m = k0/r. The rate at which the meteor's radius decreases as it burns up may depend on r, on the velocity y′, and, since the density of the atmosphere depends on y, on y itself. However, we will construct a very simple model where the rate is just proportional to the square of the velocity,

We solve (1) and (2), with k0 = 0.005, c0 = 10-8, g = 9.8 and initial conditions y(0) = 100,000 meters, y'(0) = ‑1000 meters/second, r(0) = 1 meter.

USE IVOAM_INT

USE UMACH_INT

IMPLICIT NONE

INTEGER IDO, IEND, NOUT, KORDER(2)

REAL T, TEND, Y(3), HIDRVS(2), DELTA, EQNERR(2)

EXTERNAL FCN

! Initialize

IDO = 1

T = 0.0

Y(1) = 100000.0

Y(2) = -1000.0

Y(3) = 1.0

KORDER(1) = 2

KORDER(2) = 1

EQNERR = .003

! Write title

CALL UMACH (2, NOUT)

WRITE (NOUT,99997)

! Integrate ODE

IEND = 0

DELTA = 10.0

DO

IEND = IEND + 1

TEND = T + DELTA

IF(TEND .GT. 50.0) TEND = 50.0

CALL IVOAM (IDO, FCN, T, TEND, Y, HIDRVS, &

KORDER=KORDER, EQNERR=EQNERR)

IF (IEND .LE. 5) THEN

WRITE (NOUT,99998) T, Y(1), Y(2), HIDRVS(1)

WRITE (NOUT,99999) Y(3), HIDRVS(2)

! Finish up

IF (IEND .EQ. 5) IDO = 3

CYCLE

END IF

EXIT

END DO

99997 FORMAT (11X, 'T', 10X, 'Y1/Y2', 11X, 'Y1P', 11X, 'Y1PP/Y2PP')

99998 FORMAT (4F15.4)

99999 FORMAT (2(15X, F15.4))

END

SUBROUTINE FCN (IDO, T, Y, HIDRVS)

INTEGER IDO

REAL T, Y(*), HIDRVS(*)

HIDRVS(1) = -9.8 - .005/y(3)*y(2)

HIDRVS(2) = -1.0E-8 * y(2)*y(2)

RETURN

END

Output

T Y1/Y2 Y1P Y1PP/Y2PP

10.0000 89773.0391 -1044.0096 -3.9701

0.8954 -0.0109

20.0000 79150.9844 -1078.6334 -2.9083

0.7826 -0.0116

30.0000 68240.9453 -1101.0380 -1.5031

0.6635 -0.0121

40.0000 57184.9062 -1106.9635 0.4253

0.5413 -0.0121

50.0000 46178.1367 -1089.8292 3.1700

0.4201 -0.0119