d02qgf is a reverse communication routine for integrating a non-stiff system of first-order ordinary differential equations using a variable-order variable-step Adams' method. A root-finding facility is provided.

2 Specification

Fortran Interface

Subroutine d02qgf (

neqf, t, y, tout, neqg, root, irevcm, trvcm, yrvcm, yprvcm, grvcm, kgrvcm, rwork, lrwork, iwork, liwork, ifail)

Integer, Intent (In)	::	neqf, neqg, lrwork, liwork
Integer, Intent (Inout)	::	irevcm, kgrvcm, iwork(liwork), ifail
Integer, Intent (Out)	::	yrvcm, yprvcm
Real (Kind=nag_wp), Intent (In)	::	tout, grvcm
Real (Kind=nag_wp), Intent (Inout)	::	t, y(neqf), rwork(lrwork)
Real (Kind=nag_wp), Intent (Out)	::	trvcm
Logical, Intent (Out)	::	root

C Header Interface

#include <nag.h>

void

d02qgf_ (const Integer *neqf, double *t, double y[], const double *tout, const Integer *neqg, logical *root, Integer *irevcm, double *trvcm, Integer *yrvcm, Integer *yprvcm, const double *grvcm, Integer *kgrvcm, double rwork[], const Integer *lrwork, Integer iwork[], const Integer *liwork, Integer *ifail)

The routine may be called by the names d02qgf or nagf_ode_ivp_adams_roots_revcom.

3 Description

Given the initial values

x, y_{1}, y_{2}, \dots, y_{neqf}

d02qgf integrates a non-stiff system of first-order differential equations of the type

y_{i}^{'} = f_{i} (x, y_{1}, y_{2}, \dots, y_{neqf}), i = 1, 2, \dots, neqf,

from

x = t

x = tout

using a variable-order variable-step Adams' method. You define the system by reverse communication, evaluating

f_{i}

in terms of

x

and

y_{1}, y_{2}, \dots, y_{neqf}

, and

y_{1}, y_{2}, \dots, y_{neqf}

are supplied at

x = t

by d02qgf. The routine is capable of finding roots (values of

x

) of prescribed event functions of the form

g_{j} (x, y, y^{'}) = 0, j = 1, 2, \dots, neqg .

Each

g_{j}

is considered to be independent of the others so that roots are sought of each

g_{j}

individually. The root reported by the routine will be the first root encountered by any

g_{j}

. Two techniques for determining the presence of a root in an integration step are available: the sophisticated method described in Watts (1985) and a simplified method whereby sign changes in each

g_{j}

are looked for at the ends of each integration step. You also define each

g_{j}

by reverse communication. In one-step mode the routine returns an approximation to the solution at each integration point. In interval mode this value is returned at the end of the integration range. If a root is detected this approximation is given at the root. You select the mode of operation, the error control, the root-finding technique and various optional inputs by a prior call to the setup routine d02qwf.

For a description of the practical implementation of an Adams' formula see Shampine and Gordon (1975).

4 References

Shampine L F and Gordon M K (1975) Computer Solution of Ordinary Differential Equations – The Initial Value Problem W H Freeman & Co., San Francisco

Shampine L F and Watts H A (1979) DEPAC – design of a user oriented package of ODE solvers Report SAND79-2374 Sandia National Laboratory

Watts H A (1985) RDEAM – An Adams ODE code with root solving capability Report SAND85-1595 Sandia National Laboratory

5 Arguments

Note: this routine uses reverse communication. Its use involves an initial entry, intermediate exits and re-entries, and a final exit, as indicated by the argument irevcm. Between intermediate exits and re-entries, all arguments other than grvcm and rwork must remain unchanged.

1: $neqf$ – Integer Input

On initial entry: the number of first-order ordinary differential equations to be solved by d02qgf. It must contain the same value as the argument neqf used in the prior call to d02qwf.

Constraint:

neqf \geq 1

2: $t$ – Real (Kind=nag_wp) Input/Output

On initial entry: that is after a call to d02qwf with

statef ='S'

, t must be set to the initial value of the independent variable

x

On final exit: the value of

x

at which

y

has been computed. This may be an intermediate output point, a root, tout or a point at which an error has occurred. If the integration is to be continued, possibly with a new value for tout, t must not be changed.

3: $y (neqf)$ – Real (Kind=nag_wp) array Input/Output

On initial entry: the initial values of the solution

y_{1}, y_{2}, \dots, y_{neqf}

On final exit: the computed values of the solution at the exit value of t. If the integration is to be continued, possibly with a new value for tout, these values must not be changed.

4: $tout$ – Real (Kind=nag_wp) Input

On initial entry: the next value of

x

at which a computed solution is required. For the initial t, the input value of tout is used to determine the direction of integration. Integration is permitted in either direction. If

tout = t

on exit, tout must be reset beyond t in the direction of integration, before any continuation call.

5: $neqg$ – Integer Input

On initial entry: the number of event functions which you are defining for root-finding. If root-finding is not required the value for neqg must be

\leq 0

. Otherwise it must be the same value as the argument neqg used in the prior call to d02qwf.

6: $root$ – Logical Output

On final exit: if root-finding was required (

neqg > 0

on entry), root specifies whether or not the output value of the argument t is a root of one of the event functions. If

root = .FALSE.

, no root was detected, whereas

root = .TRUE.

indicates a root and you should make a call to d02qyf for further information.

If root-finding was not required (

neqg = 0

on entry),

root = .FALSE.

7: $irevcm$ – Integer Input/Output

On initial entry: must have the value

0

On intermediate exit: specifies what action you must take before re-entering d02qgf with irevcm unchanged.

$irevcm = 1$ , $2$ , $3$ , $4$ , $5$ , $6$ or $7$: Indicates that you must supply $y^{'} = f (x, y)$ , where $x$ is given by trvcm and $y_{i}$ is returned in $y (i)$ , for $i = 1, 2, \dots, neqf$ when $yrvcm = 0$ and $rwork (yrvcm + i - 1)$ , for $i = 1, 2, \dots, neqf$ when $yrvcm \neq 0$ . $y_{i}^{'}$ should be placed in location $rwork (yprvcm + i - 1)$ , for $i = 1, 2, \dots, neqf$ .
$irevcm = 8$: Indicates that the current step was not successful due to error test failure. The only information supplied to you on this return is the current value of the independent variable t, as given by trvcm. No values must be changed before re-entering d02qgf. This facility enables you to determine the number of unsuccessful steps.
$irevcm = 9$ , $10$ , $11$ or $12$: Indicates that you must supply $g_{k} (x, y, y^{'})$ , where $k$ is given by kgrvcm, $x$ is given by trvcm, $y_{i}$ is given by $y (i)$ and $y_{i}^{'}$ is given by $rwork (yprvcm - 1 + i)$ . The result $g_{k}$ should be placed in the variable grvcm.

On final exit: has the value

0

, which indicates that an output point or root has been reached or an error has occurred (see ifail).

Note: any values you return to d02qgf as part of the reverse communication procedure should not include floating-point NaN (Not a Number) or infinity values, since these are not handled by d02qgf. If your code does inadvertently return any NaNs or infinities, d02qgf is likely to produce unexpected results.

8: $trvcm$ – Real (Kind=nag_wp) Output

On intermediate exit: the current value of the independent variable.

9: $yrvcm$ – Integer Output

On intermediate exit: with

irevcm = 1

2

3

4

5

6

7

9

10

11

12

, yrvcm specifies the locations of the dependent variables

y

for use in evaluating the differential system or the event functions.

$yrvcm = 0$: $y_{i}$ is given by $y (i)$ , for $i = 1, 2, \dots, neqf$ .
$yrvcm \neq 0$: $y_{i}$ is given by $rwork (yrvcm + i - 1)$ , for $i = 1, 2, \dots, neqf$ .

10: $yprvcm$ – Integer Output

On intermediate exit: with

irevcm = 1

2

3

4

5

6

7

, yprvcm specifies the positions in rwork at which you should place the derivatives

y^{'}

y_{i}^{'}

should be placed in location

rwork (yprvcm + i - 1)

, for

i = 1, 2, \dots, neqf

With

irevcm = 9

10

11

12

, yprvcm specifies the locations of the derivatives

y^{'}

for use in evaluating the event functions.

y_{i}^{'}

is given by

rwork (yprvcm + i - 1)

, for

i = 1, 2, \dots, neqf

11: $grvcm$ – Real (Kind=nag_wp) Input

On initial entry: need not be set.

On intermediate re-entry: with

irevcm = 9

10

11

12

, grvcm must contain the value of

g_{k} (x, y, y^{'})

, where

k

is given by kgrvcm.

12: $kgrvcm$ – Integer Input/Output

On intermediate re-entry: with

irevcm = 9

10

11

12

, kgrvcm must remain unchanged from a previous call to d02qgf.

On intermediate exit: with

irevcm = 9

10

11

12

, kgrvcm specifies which event function

g_{k} (x, y, y^{'})

you must evaluate.

13: $rwork (lrwork)$ – Real (Kind=nag_wp) array Communication Array

This must be the same argument rwork as supplied to d02qwf. It is used to pass information from d02qwf to d02qgf, and from d02qgf to d02qxf, d02qyf and d02qzf. Therefore, the contents of this array must not be changed before the call to d02qgf or calling any of the routines d02qxf, d02qyf and d02qzf.

14: $lrwork$ – Integer Input

On initial entry: the dimension of the array rwork as declared in the (sub)program from which d02qgf is called.

This must be the same argument lrwork as supplied to d02qwf.

15: $iwork (liwork)$ – Integer array Communication Array

This must be the same argument iwork as supplied to d02qwf. It is used to pass information from d02qwf to d02qgf, and from d02qgf to d02qxf, d02qyf and d02qzf. Therefore, the contents of this array must not be changed before the call to d02qgf or calling any of the routines d02qxf, d02qyf and d02qzf.

16: $liwork$ – Integer Input

On initial entry: the dimension of the array iwork as declared in the (sub)program from which d02qgf is called.

This must be the same argument liwork as supplied to d02qwf.

17: $ifail$ – Integer Input/Output

On initial entry: ifail must be set to

0

- 1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

- 1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

- 1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

- 1

is recommended since useful values can be provided in some output arguments even when

ifail \neq 0

on exit. When the value $- 1$ or $1$ is used it is essential to test the value of ifail on exit.

On final exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$: On entry, irevcm had an illegal value, $⟨ value ⟩$ .

On entry, $liwork = ⟨ value ⟩$ and $liwork = ⟨ value ⟩$ in d02qwf.
Constraint: $liwork = liwork$ in d02qwf.

On entry, $lrwork = ⟨ value ⟩$ and $lrwork = ⟨ value ⟩$ in d02qwf.
Constraint: $lrwork = lrwork$ in d02qwf.

On entry, $neqf = ⟨ value ⟩$ and $neqf = ⟨ value ⟩$ in d02qwf.
Constraint: $neqf = neqf$ in d02qwf.

On entry, $neqg = ⟨ value ⟩$ and $neqg = 0$ in d02qwf.
Constraint: if $neqg > 0$ then $neqg > 0$ in d02qwf.

On entry, $neqg = ⟨ value ⟩$ and $neqg = 0$ in d02qwf.
Constraint: if $neqg \leq 0$ then $neqg = neqg$ in d02qwf.

On entry, $neqg = ⟨ value ⟩$ and $neqg = ⟨ value ⟩$ in d02qwf.
Constraint: if $neqg \leq 0$ then $neqg \leq 0$ in d02qwf.

On entry, $tout = ⟨ value ⟩$ but $crit = .TRUE.$ in d02qwf.
Integration cannot be attempted beyond $tcrit = ⟨ value ⟩$ .

On entry, $tout = t = ⟨ value ⟩$ .

The call to setup routine d02qwf produced an error.

The setup routine d02qwf has not been called.

The value of t has been changed from $⟨ value ⟩$ to $⟨ value ⟩$ .
This is not permitted on a continuation call.

The value of tout, $⟨ value ⟩$ , indicates a change in the integration direction. This is not permitted on a continuation call.

$ifail = 2$: The maximum number of steps has been attempted.
If integration is to be continued then the routine may be called again and a further maxstp in d02qwf steps will be attempted.

$ifail = 3$: The error tolerances are too stringent.

$ifail = 4$: Error weight $i$ has become zero during the integration.
$atol (i) = 0.0$ in d02qwf but $y (i)$ is now $0.0$ . Integration successful as far as t.

$ifail = 5$: The problem appears to be stiff.
(See the D02 Chapter Introduction for a discussion of the term ‘stiff’). Although it is inefficient to use this integrator to solve stiff problems, integration may be continued by resetting ifail and calling the routine again.

$ifail = 6$: A change in sign of an event function has been detected but the root-finding process appears to have converged to a singular point of t rather than a root. Integration may be continued by calling the routine again.

$ifail = 7$: Two successive errors detected at the current value of t, $⟨ value ⟩$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The accuracy of integration is determined by the arguments vectol, rtol and atol in a prior call to d02qwf. Note that only the local error at each step is controlled by these arguments. The error estimates obtained are not strict bounds but are usually reliable over one step. Over a number of steps the overall error may accumulate in various ways, depending on the properties of the differential equation system. The code is designed so that a reduction in the tolerances should lead to an approximately proportional reduction in the error. You are strongly recommended to call d02qgf with more than one set of tolerances and to compare the results obtained to estimate their accuracy.

The accuracy obtained depends on the type of error test used. If the solution oscillates around zero a relative error test should be avoided, whereas if the solution is exponentially increasing an absolute error test should not be used. If different accuracies are required for different components of the solution then a component-wise error test should be used. For a description of the error test see the specifications of the arguments vectol, rtol and atol in the routine document for d02qwf.

The accuracy of any roots located will depend on the accuracy of integration and may also be restricted by the numerical properties of

g (x, y, y^{'})

. When evaluating

g

you should try to write the code so that unnecessary cancellation errors will be avoided.

8 Parallelism and Performance

d02qgf is not thread safe and should not be called from a multithreaded user program. Please see Section 1 in FL Interface Multithreading for more information on thread safety.

d02qgf is not threaded in any implementation.

9 Further Comments

If d02qgf fails with

ifail = 3

then the combination of atol and rtol may be so small that a solution cannot be obtained, in which case the routine should be called again with larger values for rtol and/or atol (see d02qwf). If the accuracy requested is really needed then you should consider whether there is a more fundamental difficulty. For example:

(a)in the region of a singularity the solution components will usually be of a large magnitude. d02qgf could be used in one-step mode to monitor the size of the solution with the aim of trapping the solution before the singularity. In any case numerical integration cannot be continued through a singularity, and analytical treatment may be necessary;
(b)for ‘stiff’ equations, where the solution contains rapidly decaying components, the routine will require a very small step size to preserve stability. This will usually be exhibited by excessive computing time and sometimes an error exit with $ifail = 3$ , but usually an error exit with $ifail = 2$ or $5$ . The Adams' methods are not efficient in such cases and you should consider using a routine from the Sub-chapter D02M–N. A high proportion of failed steps (see argument nfail in d02qxf) may indicate stiffness but there may be other reasons for this phenomenon.

d02qgf can be used for producing results at short intervals (for example, for graph plotting); you should set

crit = .TRUE.

and tcrit to the last output point required in a prior call to d02qwf and then set tout appropriately for each output point in turn in the call to d02qgf.

10 Example

This example solves the following system (for a projectile)

\begin{array}{l} y^{'} = \tan ϕ \\ v^{'} = \frac{- 0.032 \tan ϕ}{v} - \frac{0.02 v}{\cos ϕ} \\ ϕ^{'} = \frac{- 0.032}{v^{2}} \end{array}

over an interval

[0.0, 10.0]

starting with values

y = 0.5

v = 0.5

and

ϕ = π / 5

using scalar error control (

vectol = .FALSE.

) until the first point where

y = 0.0

is encountered.

Also, d02qgf is used to produce output at intervals of

2.0

d02qg: FL CL CPP AD

NAG FL Interfaced02qgf (ivp_​adams_​roots_​revcom)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
d02qgf (ivp_adams_roots_revcom)