NAG Library Function Document

1Purpose

nag_ode_ivp_adams_roots (d02qfc) is a function 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.

2Specification

 #include #include
 void (*fcn)(Integer neqf, double x, const double y[], double f[], Nag_User *comm),
double *t, double y[], double tout,
 double (*g)(Integer neqf, double x, const double y[], const double yp[], Integer k, Nag_User *comm),
Nag_User *comm, Nag_ODE_Adams *opt, NagError *fail)

3Description

Given the initial values $x,{y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$ the function integrates a non-stiff system of first order ordinary differential equations of the type, ${y}_{\mathit{i}}^{\prime }={f}_{\mathit{i}}\left(x,{y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{neqf}}$, from $x={\mathbf{t}}$ to $x={\mathbf{tout}}$ using a variable-order variable-step Adams' method. The system is defined by fcn, which evaluates ${f}_{i}$ in terms of $x$ and ${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$, and ${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$ are supplied at $x={\mathbf{t}}$. The function is capable of finding roots (values of $x$) of prescribed event functions of the form
 $g j x,y, y ′ = 0 , j = 1 , 2 , … , neqg .$
(See nag_ode_ivp_adams_setup (d02qwc) for the specification of 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 function 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. The event functions are defined by g, which evaluates ${g}_{j}$ in terms of $x,{y}_{1},\dots ,{y}_{{\mathbf{neqf}}}$ and ${y}_{1}^{\prime },\dots ,{y}_{{\mathbf{neqf}}}^{\prime }$. In one-step mode the function 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 need to select the mode of operation, the error control, the root-finding technique and various integration inputs with a prior call of the setup function nag_ode_ivp_adams_setup (d02qwc).
For a description of the practical implementation of an Adams' formula see Shampine and Gordon (1975) and Shampine and Watts (1979).
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

5Arguments

1:    $\mathbf{neqf}$IntegerInput
On entry: the number of differential equations.
Constraint: ${\mathbf{neqf}}\ge 1$.
2:    $\mathbf{fcn}$function, supplied by the userExternal Function
fcn must evaluate the functions ${f}_{i}$ (that is the first derivatives ${y}_{i}^{\prime }$) for given values of its arguments $x,{y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$.
The specification of fcn is:
 void fcn (Integer neqf, double x, const double y[], double f[], Nag_User *comm)
1:    $\mathbf{neqf}$IntegerInput
On entry: the number of differential equations.
2:    $\mathbf{x}$doubleInput
On entry: the current value of the argument $x$.
3:    $\mathbf{y}\left[{\mathbf{neqf}}\right]$const doubleInput
On entry: ${\mathbf{y}}\left[\mathit{i}-1\right]$ contains the current value of the argument ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neqf}}$.
4:    $\mathbf{f}\left[{\mathbf{neqf}}\right]$doubleOutput
On exit: ${\mathbf{f}}\left[\mathit{i}-1\right]$ must contain the value of ${f}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neqf}}$.
5:    $\mathbf{comm}$Nag_User *
Pointer to a structure of type Nag_User with the following member:
pPointer
On entry/exit: the pointer $\mathbf{comm}\mathbf{\to }\mathbf{p}$ should be cast to the required type, e.g., struct user *s = (struct user *)comm → p, to obtain the original object's address with appropriate type.
Note: fcn should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by nag_ode_ivp_adams_roots (d02qfc). If your code inadvertently does return any NaNs or infinities, nag_ode_ivp_adams_roots (d02qfc) is likely to produce unexpected results.
3:    $\mathbf{t}$double *Input/Output
On entry: after a call to nag_ode_ivp_adams_setup (d02qwc) with ${\mathbf{state}}=\mathrm{Nag_NewStart}$ (i.e., an initial entry), t must be set to the initial value of the independent variable $x$.
On 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.
4:    $\mathbf{y}\left[{\mathbf{neqf}}\right]$doubleInput/Output
On entry: the initial values of the solution ${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$.
On 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.
5:    $\mathbf{tout}$doubleInput
On 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 ${\mathbf{tout}}={\mathbf{t}}$ on exit, tout must be reset beyond t in the direction of integration, before any continuation call.
6:    $\mathbf{g}$function, supplied by the userExternal Function
g must evaluate a given component of $g\left(x,y,{y}^{\prime }\right)$ at a specified point.
If root-finding is not required the actual argument for g must be the NAG defined null double function pointer NULLDFN.
The specification of g is:
 double g (Integer neqf, double x, const double y[], const double yp[], Integer k, Nag_User *comm)
1:    $\mathbf{neqf}$IntegerInput
On entry: the number of differential equations.
2:    $\mathbf{x}$doubleInput
On entry: the current value of the independent variable.
3:    $\mathbf{y}\left[{\mathbf{neqf}}\right]$const doubleInput
On entry: the current values of the dependent variables.
4:    $\mathbf{yp}\left[{\mathbf{neqf}}\right]$const doubleInput
On entry: the current values of the derivatives of the dependent variables.
5:    $\mathbf{k}$IntegerInput
On entry: the component of $g$ which must be evaluated.
6:    $\mathbf{comm}$Nag_User *
Pointer to a structure of type Nag_User with the following member:
pPointer
On entry/exit: the pointer $\mathbf{comm}\mathbf{\to }\mathbf{p}$ should be cast to the required type, e.g., struct user *s = (struct user *)comm → p, to obtain the original object's address with appropriate type.
Note: g should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by nag_ode_ivp_adams_roots (d02qfc). If your code inadvertently does return any NaNs or infinities, nag_ode_ivp_adams_roots (d02qfc) is likely to produce unexpected results.
7:    $\mathbf{comm}$Nag_User *
Pointer to a structure of type Nag_User with the following member:
pPointer
On entry/exit: the pointer $\mathbf{comm}\mathbf{\to }\mathbf{p}$, of type Pointer, allows you to communicate information to and from fcn and g. An object of the required type should be declared, e.g., a structure, and its address assigned to the pointer $\mathbf{comm}\mathbf{\to }\mathbf{p}$ by means of a cast to Pointer in the calling program. E.g. comm.p = (Pointer)&s.
8:    $\mathbf{opt}$Nag_ODE_Adams *
Pointer to a structure of type Nag_ODE_Adams as initialized by the setup function nag_ode_ivp_adams_setup (d02qwc) with the following members:
rootNag_BooleanOutput
On exit: if root-finding was required (${\mathbf{neqg}}>0$ in a call to the setup function nag_ode_ivp_adams_setup (d02qwc) ), then $\mathbf{root}$ specifies whether or not the output value of the argument t is a root of one of the event functions. If $\mathbf{root}=\mathrm{Nag_FALSE}$, then no root was detected, whereas $\mathbf{root}=\mathrm{Nag_TRUE}$ indicates a root.
If root-finding was not required (${\mathbf{neqg}}=0$) then on exit $\mathbf{root}=\mathrm{Nag_FALSE}$.
If $\mathbf{root}=\mathrm{Nag_FALSE}$, then $\mathbf{opt}\mathbf{\to }\mathbf{index}$, $\mathbf{opt}\mathbf{\to }\mathbf{type}$, $\mathbf{opt}\mathbf{\to }\mathbf{events}$ and $\mathbf{opt}\mathbf{\to }\mathbf{resids}$ are indeterminate.
indexIntegerOutput
On exit: the index $k$ of the event equation ${g}_{k}\left(x,y,{y}^{\prime }\right)=0$ for which the root has been detected.
typeIntegerOutput
On exit: information about the root detected for the event equation defined by $\mathbf{opt}\mathbf{\to }\mathbf{index}$. The possible values of $\mathbf{type}$ with their interpretations are as follows:
• If $\mathbf{type}=1$, a simple root, or lack of distinguishing information available.
• If $\mathbf{type}=2$, a root of even multiplicity is believed to have been detected, that is no change in sign of the event function was found.
• If $\mathbf{type}=3$, a high order root of odd multiplicity.
• If $\mathbf{type}=4$, a possible root, but due to high multiplicity or a clustering of roots accurate evaluation of the event function was prohibited by round-off error and/or cancellation.
In general, the accuracy of the root is less reliable for values of $\mathbf{type}>1$.
eventsInteger *Output
On exit: array pointer containing information about the $k$th event function on a very small interval containing the root, t. All roots lying in this interval are considered indistinguishable numerically and therefore should be regarded as defining a root at t. The possible values of $\mathbf{events}\left[j-1\right]$, $j=1,2,\dots ,$neqg, with their interpretations are as follows:
• $\mathbf{events}\left[j-1\right]=0$, the $j$th event function did not have a root;
• $\mathbf{events}\left[j-1\right]=-1$, the $j$th event function changed sign from positive to negative about a root, in the direction of integration;
• $\mathbf{events}\left[j-1\right]=1$, the $j$th event function changed sign from negative to positive about a root, in the direction of integration;
• $\mathbf{events}\left[j-1\right]=2$, a root was identified, but no change in sign was observed.
residsdoubleOutput
On exit: array pointer, $\mathbf{opt}\mathbf{\to }\mathbf{resids}\left[j-1\right]$, $j=1,2,\dots ,$neqg, contains value of the $j$th event function computed at the root, t.
ypdoubleOutput
On exit: array pointer to the approximate derivative of the solution component ${y}_{i}$ at the output value of t. These values are obtained by the evaluation of ${y}^{\prime }=f\left(x,y\right)$ except when the output value of the argument t is tout and $\mathbf{opt}\mathbf{\to }\mathbf{tcurr}\ne {\mathbf{tout}}$, in which case they are obtained by interpolation.
tcurrdoubleOutput
On exit: the value of the independent variable which the integrator has actually reached. $\mathbf{tcurr}$ will always be at least as far as the output value of the argument t in the direction of integration, but may be further.
hlastdoubleOutput
On exit: the last successful step size used in the integration.
hnextdoubleOutput
On exit: the next step size which the integration would attempt.
ord_lastIntegerOutput
On exit: the order of the method last used (successfully) in the integration.
ord_nextIntegerOutput
On exit: the order of the method which the integration would attempt on the next step.
nsuccessIntegerOutput
On exit: the number of integration steps attempted that have been successful since the start of the current problem.
nfailIntegerOutput
On exit: the number of integration steps attempted that have failed since the start of the current problem.
tolfacdoubleOutput
On exit: a tolerance scale factor, $\mathbf{tolfac}\ge 1.0$, returned when nag_ode_ivp_adams_roots (d02qfc) exits with ${\mathbf{fail}}\mathbf{.}\mathbf{code}={\mathbf{NE_ODE_TOL}}$. If rtol and atol are uniformly scaled up by a factor of $\mathbf{tolfac}$ and nag_ode_ivp_adams_setup (d02qwc) is called, the next call to nag_ode_ivp_adams_roots (d02qfc) is deemed likely to succeed.
9:    $\mathbf{fail}$NagError *Input/Output
The NAG error argument (see Section 3.7 in How to Use the NAG Library and its Documentation).

6Error Indicators and Warnings

NE_DIRECTION_CHANGE
The value of tout, $〈\mathit{\text{value}}〉$, indicates a change in the integration direction. This is not permitted on a continuation call.
NE_MAX_STEP
The maximum number of steps have been attempted.
If integration is to be continued then the function may be called again and a further max_step steps will be attempted (see nag_ode_ivp_adams_setup (d02qwc) for details of max_step ).
NE_NEQF
The value of neqf supplied is not the same as that given to the setup function nag_ode_ivp_adams_setup (d02qwc). ${\mathbf{neqf}}=〈\mathit{\text{value}}〉$ but the value given to
nag_ode_ivp_adams_setup (d02qwc) was $〈\mathit{\text{value}}〉$.
NE_NO_G_FUN
Root finding has been requested by setting ${\mathbf{neqg}}>0$, ${\mathbf{neqg}}=〈\mathit{\text{value}}〉$, but argument g is a null function.
NE_NO_SETUP
The setup function nag_ode_ivp_adams_setup (d02qwc) has not been called.
NE_ODE_TOL
The error tolerances are too stringent. rtol and atol should be scaled up by the factor $\mathbf{opt}\mathbf{\to }\mathbf{tolfac}$ and the integration function re-entered. $\mathbf{opt}\mathbf{\to }\mathbf{tolfac}=〈\mathit{\text{value}}〉$ (see Section 9).
NE_SETUP_ERROR
The call to setup function nag_ode_ivp_adams_setup (d02qwc) produced an error.
NE_SINGULAR_POINT
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 function again.
NE_STIFF_PROBLEM
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 fail and calling the function again.
NE_T_CHANGED
The value of t has been changed from $〈\mathit{\text{value}}〉$ to $〈\mathit{\text{value}}〉$. This is not permitted on a continuation call.
NE_T_SAME_TOUT
On entry, ${\mathbf{tout}}={\mathbf{t}}$, t is $〈\mathit{\text{value}}〉$.
NE_TOUT_TCRIT
${\mathbf{tout}}=〈\mathit{\text{value}}〉$ but crit was set Nag_TRUE in setup call and integration cannot be attempted beyond ${\mathbf{tcrit}}=〈\mathit{\text{value}}〉$.
NE_WEIGHT_ZERO
An error weight has become zero during the integration, see d02qwc document; ${\mathbf{atol}}\left[〈\mathit{\text{value}}〉\right]$ was set to 0.0 but ${\mathbf{y}}\left[〈\mathit{\text{value}}〉\right]$ is now 0.0. Integration successful as far as ${\mathbf{t}}=〈\mathit{\text{value}}〉$.
The value of the array index is returned in ${\mathbf{fail}}\mathbf{.}\mathbf{errnum}$.

7Accuracy

The accuracy of integration is determined by the arguments vectol, rtol and atol in a prior call to nag_ode_ivp_adams_setup (d02qwc). 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 nag_ode_ivp_adams_roots (d02qfc) 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, atol and rtol in the function document for nag_ode_ivp_adams_setup (d02qwc).
The accuracy of any roots located will depend on the accuracy of integration and may also be restricted by the numerical properties of $g\left(x,y,{y}^{\prime }\right)$. When evaluating $g$ you should try to write the code so that unnecessary cancellation errors will be avoided.

8Parallelism and Performance

If the function fails with ${\mathbf{fail}}\mathbf{.}\mathbf{code}={\mathbf{NE_ODE_TOL}}$, then the combination of atol and rtol may be so small that a solution cannot be obtained, in which case the function should be called again using larger values for rtol and/or atol when calling the setup function nag_ode_ivp_adams_setup (d02qwc). 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. The function 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 function 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 ${\mathbf{fail}}\mathbf{.}\mathbf{code}={\mathbf{NE_ODE_TOL}}$, but usually an error exit with ${\mathbf{fail}}\mathbf{.}\mathbf{code}={\mathbf{NE_MAX_STEP}}$ or NE_STIFF_PROBLEM. The Adams' methods are not efficient in such cases. A high proportion of failed steps (see argument $\mathbf{opt}\mathbf{\to }\mathbf{nfail}$) may indicate stiffness but there may be other reasons for this phenomenon.
nag_ode_ivp_adams_roots (d02qfc) can be used for producing results at short intervals (for example, for graph plotting); you should set ${\mathbf{crit}}=\mathrm{Nag_TRUE}$ and tcrit to the last output point required in a prior call to nag_ode_ivp_adams_setup (d02qwc) and then set tout appropriately for each output point in turn in the call to nag_ode_ivp_adams_roots (d02qfc).
The structure opt will contain pointers which have been allocated memory by calls to nag_ode_ivp_adams_setup (d02qwc). This allocated memory is then accessed by nag_ode_ivp_adams_roots (d02qfc) and, if required, nag_ode_ivp_adams_interp (d02qzc). When all calls to these functions have been completed the function nag_ode_ivp_adams_free (d02qyc) may be called to free memory allocated to the structure.

10Example

We solve the equation
 $y ′′ = -y , y 0 = 0 , y ′ 0 = 1$
reposed as
 $y 1 ′ = y 2 y 2 ′ = - y 1$
over the range $\left[0,10.0\right]$ with initial conditions ${y}_{1}=0.0$ and ${y}_{2}=1.0$ using vector error control (${\mathbf{vectol}}=\mathrm{Nag_TRUE}$) and computation of the solution at ${\mathbf{tout}}=10.0$ with ${\mathbf{tcrit}}=10.0$ (${\mathbf{crit}}=\mathrm{Nag_TRUE}$). Also, we use nag_ode_ivp_adams_roots (d02qfc) to locate the positions where ${y}_{1}=0.0$ or where the first component has a turning point, that is ${y}_{1}^{\prime }=0.0$.

10.1Program Text

Program Text (d02qfce.c)

None.

10.3Program Results

Program Results (d02qfce.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017