NAG Toolbox: nag_ode_bvp_shoot_genpar_algeq (d02sa)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{p}\left(\mathbf{m}\right)$ – double array

$\mathbf{p}\left(\mathit{i}\right)$ must be set to an estimate of the $\mathit{i}$th argument, $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.

2:     $\mathrm{n1}$ – int64int32nag_int scalar

$n_1$, the number of differential equations active in the matching process. The active equations must be placed last in the numbering in fcn and bc. The first $\mathbf{n}-\mathbf{n1}$ equations are used as the driving equations.

Constraint: $\mathbf{n1}\le \mathbf{n}$, $\mathbf{n1}\le \mathbf{m}$ and $\mathbf{n1}>0$.

3:     $\mathrm{pe}\left(\mathbf{m}\right)$ – double array

$\mathbf{pe}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,m$, must be set to a positive value for use in the convergence test in the $i$th argument $p_i$. See the description of pf for further details.

Constraint: $\mathbf{pe}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,m$.

4:     $\mathrm{pf}\left(\mathbf{m}\right)$ – double array

$\mathbf{pf}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,m$, should be set to a ‘floor’ value in the convergence test on the $i$th argument $p_i$. If $\mathbf{pf}\left(i\right)\le 0.0$ on entry then it is set to the small positive value $\sqrt{\epsilon }$ (where $\epsilon$ may in most cases be considered to be machine precision); otherwise it is used unchanged.

The Newton iteration is presumed to have converged if a full Newton step is taken ($\mathbf{istate}=1$ in the specification of monit), the singular values of the Jacobian are not being significantly perturbed (also see monit) and if the Newton correction $C_i$ satisfies

$$\left|C_i\right|\le \mathbf{pe}\left(i\right)\times \max \left(\left|p_i\right|,\mathbf{pf}\left(i\right)\right)\text{, \hspace{1em}}i=1,2,\dots ,m\text{,}$$

where $p_i$ is the current value of the $i$th argument. The values $\mathbf{pf}\left(i\right)$ are also used in determining the Newton iterates as discussed in Description, see equation (6).

5:     $\mathrm{e}\left(\mathbf{n}\right)$ – double array

Default: $\mathbf{e}\left(i\right)={10}^{-5}$

Values for use in controlling the local error in the integration of the differential equations. If ${\mathit{err}}_{\mathit{i}}$ is an estimate of the local error in $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$, then

$$\left|{\mathit{err}}_i\right|\le \mathbf{e}\left(i\right)\times \max \left\{\sqrt{\epsilon },\left|y_i\right|\right\}\text{,}$$

where $\epsilon$ may in most cases be considered to be machine precision.

Constraint: $\mathbf{e}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,\mathbf{n}$.

6:     $\mathrm{dp}\left(\mathbf{m}\right)$ – double array

A value to be used in perturbing the argument $p_i$ in the numerical differentiation to estimate the Jacobian used in Newton's method. If $\mathbf{dp}\left(i\right)=0.0$ on entry, an estimate is made internally by setting

$$\mathbf{dp}\left(i\right)=\sqrt{\epsilon }\times \max \left(\mathbf{pf}\left(i\right),\left|p_i\right|\right)\text{,}$$

(7)

where $p_i$ is the initial value of the argument supplied by you and $\epsilon$ may in most cases be considered to be machine precision. The estimate of the Jacobian, $J$, is made using forward differences, that is for each $\mathit{i}$, for $\mathit{i}=1,2,\dots ,m$, $p_{\mathit{i}}$ is perturbed to $p_{\mathit{i}}+\mathbf{dp}\left(\mathit{i}\right)$ and the $\mathit{i}$th column of $J$ is estimated as

$$\left(r\left(p_{\mathit{i}}+\mathbf{dp}\left(\mathit{i}\right)\right)-r\left(p_{\mathit{i}}\right)\right)/\mathbf{dp}\left(\mathit{i}\right)$$

where the other components of $p$ are unchanged (see (3) for the notation used). If this fails to produce a Jacobian with significant columns, backward differences are tried by perturbing $p_{\mathit{i}}$ to $p_{\mathit{i}}-\mathbf{dp}\left(\mathit{i}\right)$ and if this also fails then central differences are used with $p_{\mathit{i}}$ perturbed to $p_{\mathit{i}}+10.0\times \mathbf{dp}\left(\mathit{i}\right)$. If this also fails then the calculation of the Jacobian is abandoned. If the Jacobian has not previously been calculated then an error exit is taken. If an earlier estimate of the Jacobian is available then the current argument set, $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$, is abandoned in favour of the last argument set from which useful progress was made and the singular values of the Jacobian used at the point are modified before proceeding with the Newton iteration. You are recommended to use the default value $\mathbf{dp}\left(i\right)=0.0$ unless you have prior knowledge of a better choice. If any of the perturbations described are likely to lead to an unfortunate set of argument values then you should use constr to prevent such perturbations (all changes of arguments are checked by a call to constr).

7:     $\mathrm{swp}\left(\mathit{ldswp},6\right)$ – double array

ldswp, the first dimension of the array, must satisfy the constraint $\mathit{ldswp}\ge \mathbf{npoint}$.

$\mathbf{swp}\left(i,1\right)$ must contain an estimate for an initial step size for integration across the $i$th sub-interval $\left[\mathbf{x}\left(\mathit{i}\right),\mathbf{x}\left(\mathit{i}+1\right)\right]$, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}-1$, (see range). $\mathbf{swp}\left(i,1\right)$ should have the same sign as $\mathbf{x}\left(i+1\right)-\mathbf{x}\left(i\right)$ if it is nonzero. If $\mathbf{swp}\left(i,1\right)=0.0$, on entry, a default value for the initial step size is calculated internally. This is the recommended mode of entry.

$\mathbf{swp}\left(i,3\right)$ must contain a lower bound for the modulus of the step size on the $\mathit{i}$th sub-interval $\left[\mathbf{x}\left(\mathit{i}\right),\mathbf{x}\left(\mathit{i}+1\right)\right]$, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}-1$. If $\mathbf{swp}\left(i,3\right)=0.0$ on entry, a very small default value is used. By setting $\mathbf{swp}\left(i,3\right)>0.0$ but smaller than the expected step sizes (assuming you have some insight into the likely step sizes) expensive integrations with arguments $p$ far from the solution can be avoided.

$\mathbf{swp}\left(\mathit{i},2\right)$ must contain an upper bound on the modulus of the step size to be used in the integration on $\left[\mathbf{x}\left(\mathit{i}\right),\mathbf{x}\left(\mathit{i}+1\right)\right]$, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}-1$. If $\mathbf{swp}\left(i,2\right)=0.0$ on entry no bound is assumed. This is the recommended mode of entry unless the solution is expected to have important features which might be ‘missed’ in the integration if the step size were permitted to be chosen freely.

8:     $\mathrm{icount}$ – int64int32nag_int scalar

An upper bound on the number of Newton iterations. If $\mathbf{icount}=0$ on entry, no check on the number of iterations is made (this is the recommended mode of entry).

Constraint: $\mathbf{icount}\ge 0$.

9:     $\mathrm{range}$ – function handle or string containing name of m-file

range must specify the break-points $x_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}$, which may depend on the arguments $p_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,m$.

[x] = range(npoint, p, m)

Input Parameters

1:     $\mathrm{npoint}$ – int64int32nag_int scalar

$2$ plus the number of break-points in $\left(a,b\right)$.

2:     $\mathrm{p}\left(\mathbf{m}\right)$ – double array

The current estimate of the $\mathit{i}$th argument, for $\mathit{i}=1,2,\dots ,m$.

3:     $\mathrm{m}$ – int64int32nag_int scalar

$m$, the number of arguments.

Output Parameters

1:     $\mathrm{x}\left(\mathbf{npoint}\right)$ – double array

The $\mathit{i}$th break-point, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}$. The sequence $\left(\mathbf{x}\left(i\right)\right)$ must be strictly monotonic, that is either

$$a=\mathbf{x}\left(1\right)<\mathbf{x}\left(2\right)<\cdots <\mathbf{x}\left(\mathbf{npoint}\right)=b$$

or

$$a=\mathbf{x}\left(1\right)>\mathbf{x}\left(2\right)>\cdots >\mathbf{x}\left(\mathbf{npoint}\right)=b\text{.}$$

10:   $\mathrm{bc}$ – function handle or string containing name of m-file

bc must place in g1 and g2 the boundary conditions at $a$ and $b$ respectively.

[g1, g2] = bc(p, m, n)

Input Parameters

1:     $\mathrm{p}\left(\mathbf{m}\right)$ – double array

An estimate of the $\mathit{i}$th argument, $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.

2:     $\mathrm{m}$ – int64int32nag_int scalar

$m$, the number of arguments.

3:     $\mathrm{n}$ – int64int32nag_int scalar

$n$, the number of differential equations.

Output Parameters

1:     $\mathrm{g1}\left(\mathbf{n}\right)$ – double array

The value of $y_{\mathit{i}}\left(a\right)$, for $\mathit{i}=1,2,\dots ,n$, (where this may be a known value or a function of the parameters $p_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,m$).

2:     $\mathrm{g2}\left(\mathbf{n}\right)$ – double array

The value of $y_{\mathit{i}}\left(b\right)$, for $\mathit{i}=1,2,\dots ,n$, (where these may be known values or functions of the parameters $p_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,m$). If $n>n_1$, so that there are some driving equations, then the first $n-n_1$ values of g2 need not be set since they are never used.

11:   $\mathrm{fcn}$ – function handle or string containing name of m-file

fcn must evaluate the functions $f_{\mathit{i}}$ (i.e., the derivatives $y_{\mathit{i}}^{\prime }$), for $\mathit{i}=1,2,\dots ,n$.

[f] = fcn(x, y, n, p, m, ii)

Input Parameters

1:     $\mathrm{x}$ – double scalar

$x$, the value of the argument.

2:     $\mathrm{y}\left(\mathbf{n}\right)$ – double array

$y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, the value of the argument.

3:     $\mathrm{n}$ – int64int32nag_int scalar

$\mathit{n}$, the number of equations.

4:     $\mathrm{p}\left(\mathbf{m}\right)$ – double array

The current estimate of the $\mathit{i}$th argument $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.

5:     $\mathrm{m}$ – int64int32nag_int scalar

$m$, the number of arguments.

6:     $\mathrm{ii}$ – int64int32nag_int scalar

Specifies the sub-interval $\left[x_i,x_{i+1}\right]$ on which the derivatives are to be evaluated.

Output Parameters

1:     $\mathrm{f}\left(\mathbf{n}\right)$ – double array

The derivative of $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$, evaluated at $x$. $\mathbf{f}\left(i\right)$ may depend upon the parameters $p_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,m$. If there are any driving equations (see Description) then these must be numbered first in the ordering of the components of f.

12:   $\mathrm{eqn}$ – function handle or string containing name of m-file

eqn is used to describe the additional algebraic equations to be solved in the determination of the parameters, $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$. If there are no additional algebraic equations (i.e., $m=n_1$) then eqn is never called and the string nag_ode_bvp_shoot_genpar_algeq_dummy_eqn (d02hbz) should be used as the actual argument.

[e] = eqn(q, p, m)

Input Parameters

1:     $\mathrm{q}$ – int64int32nag_int scalar

The number of algebraic equations, $m-n_1$.

2:     $\mathrm{p}\left(\mathbf{m}\right)$ – double array

The current estimate of the $\mathit{i}$th argument $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.

3:     $\mathrm{m}$ – int64int32nag_int scalar

$m$, the number of arguments.

Output Parameters

1:     $\mathrm{e}\left(\mathbf{q}\right)$ – double array

The vector of residuals, $r_2\left(p\right)$, that is the amount by which the current estimates of the arguments fail to satisfy the algebraic equations.

13:   $\mathrm{constr}$ – function handle or string containing name of m-file

constr is used to prevent the pseudo-Newton iteration running into difficulty. constr should return the value true if the constraints are satisfied by the parameters $p_1,p_2,\dots ,p_m$. Otherwise constr should return the value false. Usually the dummy function nag_ode_bvp_shoot_genpar_algeq_dummy_constr (d02hby), which returns the value true at all times, will suffice and in the first instance this is recommended as the actual argument.

[result] = constr(p, m)

Input Parameters

1:     $\mathrm{p}\left(\mathbf{m}\right)$ – double array

An estimate of the $\mathit{i}$th argument, $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.

2:     $\mathrm{m}$ – int64int32nag_int scalar

$m$, the number of arguments.

Output Parameters

1:     $\mathrm{result}$ – logical scalar

$\mathbf{result}=\mathit{true}$ if the constraints are satisfied by the parameters p.

14:   $\mathrm{ymax}$ – double scalar

A non-negative value which is used as a bound on all values ${‖y\left(x\right)‖}_{\infty }$ where $y\left(x\right)$ is the solution at any point $x$ between $\mathbf{x}\left(1\right)$ and $\mathbf{x}\left(\mathbf{npoint}\right)$ for the current arguments $p_1,p_2,\dots ,p_m$. If this bound is exceeded the integration is terminated and the current arguments are rejected. Such a rejection will result in an error exit if it prevents the initial residual or Jacobian, or the final solution, being calculated. If $\mathbf{ymax}=0$ on entry, no bound on the solution $y$ is used; that is the integrations proceed without any checking on the size of ${‖y‖}_{\infty }$.

15:   $\mathrm{monit}$ – function handle or string containing name of m-file

monit enables you to monitor the values of various quantities during the calculation. It is called by nag_ode_bvp_shoot_genpar_algeq (d02sa) after every calculation of the norm ${‖{\mathbf{d}}^{-1}{\tilde{J}}^{+}r‖}_2$ which determines the strategy of the Newton method, every time there is an internal error exit leading to a change of strategy, and before an error exit when calculating the initial Jacobian. string nag_ode_bvp_shoot_genpar_algeq_sample_monit (d02hbx)nag_file_set_unit_advisory (x04ab)

monit(istate, iflag, ifail1, p, m, f, pnorm, pnorm1, eps, d)

Input Parameters

1:     $\mathrm{istate}$ – int64int32nag_int scalar

The state of the Newton iteration.

$\mathbf{istate}=0$

The calculation of the residual, Jacobian and ${‖{\mathbf{d}}^{-1}{\tilde{J}}^{+}r‖}_2$ are taking place.

$\mathbf{istate}=1$ to $5$

During the Newton iteration a factor of $2^{\left(-\mathbf{istate}+1\right)}$ of the Newton step is being used to try to reduce the norm.

$\mathbf{istate}=6$

The current Newton step has been rejected and the Jacobian is being re-calculated.

$\mathbf{istate}=-6$ to $-1$

An internal error exit has caused the rejection of the current set of argument values, $p$. $-\mathbf{istate}$ is the value which istate would have taken if the error had not occurred.

$\mathbf{istate}=-7$

An internal error exit has occurred when calculating the initial Jacobian.

2:     $\mathrm{iflag}$ – int64int32nag_int scalar

Whether or not the Jacobian being used has been calculated at the beginning of the current iteration. If the Jacobian has been updated then $\mathbf{iflag}=1$; otherwise $\mathbf{iflag}=2$. The Jacobian is only calculated when convergence to the current argument values has been slow.

3:     $\mathrm{ifail1}$ – int64int32nag_int scalar

If $-6\le \mathbf{istate}\le -1$, ifail1 specifies the ifail error number that would be produced were control returned to you. ifail1 is unspecified for values of istate outside this range.

4:     $\mathrm{p}\left(\mathbf{m}\right)$ – double array

The current estimate of the $\mathit{i}$th argument $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.

5:     $\mathrm{m}$ – int64int32nag_int scalar

$m$, the number of arguments.

6:     $\mathrm{f}\left(\mathbf{m}\right)$ – double array

$r$, the residual corresponding to the current argument values, provided $1\le \mathbf{istate}\le 5$ or $\mathbf{istate}=-7$. f is unspecified for other values of istate.

7:     $\mathrm{pnorm}$ – double scalar

A quantity against which all reductions in norm are currently measured.

8:     $\mathrm{pnorm1}$ – double scalar

$p$, the norm of the current arguments. It is set for $1\le \mathbf{istate}\le 5$ and is undefined for other values of istate.

9:     $\mathrm{eps}$ – double scalar

Gives some indication of the convergence rate. It is the current singular value modification factor (see Gay (1976)). It is zero initially and whenever convergence is proceeding steadily. eps is ${\epsilon }^{3/8}$ or greater (where $\epsilon$ may in most cases be considered machine precision) when the singular values of $J$ are approximately zero or when convergence is not being achieved. The larger the value of eps the worse the convergence rate. When eps becomes too large the Newton iteration is terminated.

10:   $\mathrm{d}\left(\mathbf{m}\right)$ – double array

$J$, the singular values of the current modified Jacobian matrix. If $\mathbf{d}\left(m\right)$ is small relative to $\mathbf{d}\left(1\right)$ for a number of Jacobians corresponding to different argument values then the computed results should be viewed with suspicion. It could be that the matching equations do not depend significantly on some argument (which could be due to a programming error in fcn, bc, range or eqn). Alternatively, the system of differential equations may be very ill-conditioned when viewed as an initial value problem, in which case nag_ode_bvp_shoot_genpar_algeq (d02sa) is unsuitable. This may also be indicated by some singular values being very large. These values of $\mathbf{d}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,m$, should not be changed.

16:   $\mathrm{prsol}$ – function handle or string containing name of m-file

prsol can be used to obtain values of the solution $y$ at a selected point $z$ by integration across the final range $\left[\mathbf{x}\left(1\right),\mathbf{x}\left(\mathbf{npoint}\right)\right]$. If no output is required nag_ode_bvp_shoot_genpar_algeq_dummy_prsol (d02hbw) can be used as the actual argument.

[z] = prsol(z, y, n)

Input Parameters

1:     $\mathrm{z}$ – double scalar

Contains $x_1$ on the first call. On subsequent calls z contains its previous output value.

2:     $\mathrm{y}\left(\mathbf{n}\right)$ – double array

The solution value $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$, at $z$.

3:     $\mathrm{n}$ – int64int32nag_int scalar

$\mathit{n}$, the total number of differential equations.

Output Parameters

1:     $\mathrm{z}$ – double scalar

The next point at which output is required. The new point must be nearer $\mathbf{x}\left(\mathbf{npoint}\right)$ than the old.

If z is set to a point outside $\left[\mathbf{x}\left(1\right),\mathbf{x}\left(\mathbf{npoint}\right)\right]$ the process stops and control returns from nag_ode_bvp_shoot_genpar_algeq (d02sa) to the (sub)program from which nag_ode_bvp_shoot_genpar_algeq (d02sa) is called. Otherwise the next call to prsol is made by nag_ode_bvp_shoot_genpar_algeq (d02sa) at the point z, with solution values $y_1,y_2,\dots ,y_n$ at z contained in y. If z is set to $\mathbf{x}\left(\mathbf{npoint}\right)$ exactly, the final call to prsol is made with $y_1,y_2,\dots ,y_n$ as values of the solution at $\mathbf{x}\left(\mathbf{npoint}\right)$ produced by the integration. In general the solution values obtained at $\mathbf{x}\left(\mathbf{npoint}\right)$ from prsol will differ from the values obtained at this point by a call to bc. The difference between the two solutions is the residual $r$. You are reminded that the points $\mathbf{x}\left(1\right),\mathbf{x}\left(2\right),\dots ,\mathbf{x}\left(\mathbf{npoint}\right)$ are available in the locations $\mathbf{swp}\left(1,4\right),\mathbf{swp}\left(2,4\right),\dots ,\mathbf{swp}\left(\mathbf{npoint},4\right)$ at all times.

Optional Input Parameters

1:     $\mathrm{m}$ – int64int32nag_int scalar

Default: the dimension of the arrays p, pe, pf, dp. (An error is raised if these dimensions are not equal.)

$m$, the number of arguments.

Constraint: $\mathbf{m}>0$.

2:     $\mathrm{n}$ – int64int32nag_int scalar

Default: the dimension of the array e.

$n$, the total number of differential equations.

Constraint: $\mathbf{n}>0$.

3:     $\mathrm{npoint}$ – int64int32nag_int scalar

Default: the first dimension of the array swp.

2 plus the number of break-points in the range of definition of the system of differential equations (1).

Constraint: $\mathbf{npoint}\ge 2$.

Output Parameters

1:     $\mathrm{p}\left(\mathbf{m}\right)$ – double array

The corrected value for the $i$th argument, unless an error has occurred, when it contains the last calculated value of the argument.

2:     $\mathrm{pf}\left(\mathbf{m}\right)$ – double array

The values actually used.

3:     $\mathrm{dp}\left(\mathbf{m}\right)$ – double array

The values actually used.

4:     $\mathrm{swp}\left(\mathit{ldswp},6\right)$ – double array

$\mathbf{swp}\left(\mathit{i},1\right)$ contains the initial step size used on the last integration on $\left[\mathbf{x}\left(\mathit{i}\right),\mathbf{x}\left(\mathit{i}+1\right)\right]$, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}-1$, (excluding integrations during the calculation of the Jacobian).

$\mathbf{swp}\left(\mathit{i},2\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}-1$, is usually unchanged. If the maximum step size $\mathbf{swp}\left(i,2\right)$ is so small or the length of the range $\left[\mathbf{x}\left(i\right),\mathbf{x}\left(i+1\right)\right]$ is so short that on the last integration the step size was not controlled in the main by the size of the error tolerances $\mathbf{e}\left(i\right)$ but by these other factors, then $\mathbf{swp}\left(\mathbf{npoint},2\right)$ is set to the floating-point value of $i$ if the problem last occurred in $\left[\mathbf{x}\left(i\right),\mathbf{x}\left(i+1\right)\right]$. Any results obtained when this value is returned as nonzero should be viewed with caution.

$\mathbf{swp}\left(\mathit{i},3\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}-1$, are unchanged.

If an error exit with $\mathbf{ifail}=\mathbf{4}$, $\mathbf{5}$ or $\mathbf{6}$ (see Error Indicators and Warnings) occurs on the integration made from $\mathbf{x}\left(i\right)$ to $\mathbf{x}\left(i+1\right)$ the floating-point value of $i$ is returned in $\mathbf{swp}\left(\mathbf{npoint},1\right)$. The actual point $x\in \left[\mathbf{x}\left(i\right),\mathbf{x}\left(i+1\right)\right]$ where the error occurred is returned in $\mathbf{swp}\left(1,5\right)$ (see also the specification of w). The floating-point value of npoint is returned in $\mathbf{swp}\left(\mathbf{npoint},1\right)$ if the error exit is caused by a call to bc.

If an error exit occurs when estimating the Jacobian matrix ($\mathbf{ifail}=\mathbf{7}$, $\mathbf{8}$, $\mathbf{9}$, $\mathbf{10}$, $\mathbf{11}$ or $\mathbf{12}$, see Error Indicators and Warnings) and if argument $p_i$ was the cause of the failure then on exit $\mathbf{swp}\left(\mathbf{npoint},1\right)$ contains the floating-point value of $i$.

$\mathbf{swp}\left(\mathit{i},4\right)$ contains the point $\mathbf{x}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{npoint}$, used at the solution $p$ or at the final values of $p$ if an error occurred.

swp is also partly used as workspace.

5:     $\mathrm{w}\left(\mathit{ldw},\mathit{sdw}\right)$ – double array

$\mathit{sdw}=3\times \mathbf{m}+12+\max \left(11,\mathbf{m}\right)$.

In the case of an error exit of the type where the point of failure is returned in $\mathbf{swp}\left(1,5\right)$, the solution at this point of failure is returned in $\mathbf{w}\left(\mathit{i},1\right)$, for $\mathit{i}=1,2,\dots ,n$.

Otherwise w is used for workspace.

6:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

   $\mathbf{ifail}=1$

One or more of the arguments n, n1, m, ldswp, npoint, icount, ldw, sdw, e, pe or ymax is incorrectly set.

   $\mathbf{ifail}=2$

The constraints have been violated by the initial arguments.

   $\mathbf{ifail}=3$

The condition $\mathbf{x}\left(1\right)<\mathbf{x}\left(2\right)<\cdots <\mathbf{x}\left(\mathbf{npoint}\right)$ (or $\mathbf{x}\left(1\right)>\mathbf{x}\left(2\right)>\cdots >\mathbf{x}\left(\mathbf{npoint}\right)$) has been violated on a call to range with the initial arguments.

   $\mathbf{ifail}=4$

In the integration from $\mathbf{x}\left(1\right)$ to $\mathbf{x}\left(\mathbf{npoint}\right)$ with the initial or the final arguments, the step size was reduced too far for the integration to proceed. Consider reversing the order of the points $\mathbf{x}\left(1\right),\mathbf{x}\left(2\right),\dots ,\mathbf{x}\left(\mathbf{npoint}\right)$. If this error exit still results, it is likely that nag_ode_bvp_shoot_genpar_algeq (d02sa) is not a suitable method for solving the problem, or the initial choice of arguments is very poor, or the accuracy requirement specified by $\mathbf{e}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,n$, is too stringent.

   $\mathbf{ifail}=5$

In the integration from $\mathbf{x}\left(1\right)$ to $\mathbf{x}\left(\mathbf{npoint}\right)$ with the initial or final arguments, an initial step could not be found to start the integration on one of the intervals $\mathbf{x}\left(i\right)$ to $\mathbf{x}\left(i+1\right)$. Consider reversing the order of the points. If this error exit still results it is likely that nag_ode_bvp_shoot_genpar_algeq (d02sa) is not a suitable function for solving the problem, or the initial choice of arguments is very poor, or the accuracy requirement specified by $\mathbf{e}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,n$, is much too stringent.

   $\mathbf{ifail}=6$

In the integration from $\mathbf{x}\left(1\right)$ to $\mathbf{x}\left(\mathbf{npoint}\right)$ with the initial or final arguments, the solution exceeded ymax in magnitude (when $\mathbf{ymax}>0.0$). It is likely that the initial choice of arguments was very poor or ymax was incorrectly set.

   $\mathbf{ifail}=7$

On calculating the initial approximation to the Jacobian, the constraints were violated.

   $\mathbf{ifail}=8$

On perturbing the arguments when calculating the initial approximation to the Jacobian, the condition $\mathbf{x}\left(1\right)<\mathbf{x}\left(2\right)<\cdots <\mathbf{x}\left(\mathbf{npoint}\right)$ (or $\mathbf{x}\left(1\right)>\mathbf{x}\left(2\right)>\cdots >\mathbf{x}\left(\mathbf{npoint}\right)$) is violated.

   $\mathbf{ifail}=9$

On calculating the initial approximation to the Jacobian, the integration step size was reduced too far to make further progress (see $\mathbf{ifail}=\mathbf{4}$).

   $\mathbf{ifail}=10$

On calculating the initial approximation to the Jacobian, the initial integration step size on some interval was too small (see $\mathbf{ifail}=\mathbf{5}$).

   $\mathbf{ifail}=11$

On calculating the initial approximation to the Jacobian, the solution of the system of differential equations exceeded ymax in magnitude (when $\mathbf{ymax}>0.0$).

   $\mathbf{ifail}=12$

On calculating the initial approximation to the Jacobian, a column of the Jacobian is found to be insignificant. This could be due to an element $\mathbf{dp}\left(i\right)$ being too small (but nonzero) or the solution having no dependence on one of the arguments (a programming error).

   $\mathbf{ifail}=13$

After calculating the initial approximation to the Jacobian, the calculation of its singular value decomposition failed. It is likely that the error will never occur as it is usually associated with the Jacobian having multiple singular values. To remedy the error it should only be necessary to change the initial arguments. If the error persists it is likely that the problem has not been correctly formulated.

   $\mathbf{ifail}=14$

The Newton iteration has failed to converge after exercising all its options. You are strongly recommended to monitor the progress of the iteration via monit. There are many possible reasons for the iteration not converging. Amongst the most likely are:

(a)	there is no solution;
(b)	the initial arguments are too far away from the correct arguments;
(c)	the problem is too ill-conditioned as an initial value problem for Newton's method to choose suitable corrections;
(d)	the accuracy requirements for convergence are too restrictive, that is some of the components of pe (and maybe pf) are too small – in this case the final value of this norm output via monit will usually be very small; or
(e)	the initial arguments are so close to the solution arguments $p$ that the Newton iteration cannot find improved arguments. The norm output by monit should be very small.

   $\mathbf{ifail}=15$

The number of iterations permitted by icount has been exceeded (in the case when $\mathbf{icount}>0$ on entry).

   $\mathbf{ifail}=16$

   $\mathbf{ifail}=17$

   $\mathbf{ifail}=18$

   $\mathbf{ifail}=19$

These indicate that there has been a serious error in an internal call. Check all function calls and array dimensions. Seek expert help.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: d02sa_example

function d02sa_example


fprintf('d02sa example results\n\n');

% For communication with prsol.
global ykeep ncall xkeep;

% Initialize variables and arrays.
p   = [1.2; 0.032; 2.5;0.2];
n1  = int64(3);
pe  = [0.001; 0.001; 0.001; 0.001];
pf  = [1e-06; 1e-06; 1e-06; 1e-06];
e   = [1e-05; 1e-05; 1e-05];
dp  = [0; 0; 0; 0];
swp = zeros(3, 6);
icount = int64(0);
ymax   = 0;
ncall  = 0;
ykeep  = zeros(1,3);
xkeep  = zeros(1,1);

fprintf('   x        y(1)        y(2)        y(3)\n')

[p, pf, dp, swp, w, ifail] = ...
  d02sa(...
        p, n1, pe, pf, e, dp, swp, icount, @range, @bc, @fcn, ...
        @eqn, @constr, ymax, @monit, @prsol);

% Plot results.
fig1 = figure;
display_plot(xkeep,ykeep)


function z = prsol(z, y, n)
  % Obtains values for solution and stores them for plotting.

  % For communication with main routine.
  global ykeep ncall xkeep;
  ncall = ncall+1;
  fprintf(' %6.4f',z);
  xkeep(ncall,1) = z;
  fprintf('%10.4f  ', y(1:n));
  fprintf('\n');
  z = z+0.5;
  if (abs(z-5) < 0.25)
    z = 5;
  end
  ykeep(ncall,:) = y;

function x = range(npoint, p, m)
  % Specify the break points.
  x = zeros(npoint,1);
  x(1) = 0;
  x(2) = p(3);
  x(3) = 5;

function [g1, g2] = bc(p, m, n)
  % Specify boundary conditions.
  g1 = [0, 0.5, p(1)];
  g2 = [0, 0.45, -1.2];

function f = fcn(x, y, n, p, m, interval)
  % Evaluate derivative vector.
  f = zeros(n,1);
  f(1) = tan(y(3));
  if (interval == 1)
    f(2) = -0.032d0*tan(y(3))/y(2) - 0.02d0*y(2)/cos(y(3));
    f(3) = -0.032d0/y(2)^2;
  else
    f(2) = -p(2)*tan(y(3))/y(2) - p(4)*y(2)/cos(y(3));
    f(3) = -p(2)/y(2)^2;
  end

function eq = eqn(q,p,m)
  % Specify additional algebraic equations.
  eq = zeros(q,1);
  eq(1) = 0.02 - p(4) - 1.0d-5*p(3);

function result = constr(p,n)
  % Check constraints.
  result = true;
  for i=1:n
    if p(i) < 0
      result = false;
    end
  end
  if p(3) > 5
    result = false;
  end

function monit(istate, iflag, ifail1, p, m, f, pnorm, pnorm1, eps, d)
  % Allows monitoring of parameter values during solution (not used here).

function display_plot(xkeep,ykeep)
  % Plot the results.
  plot(xkeep, ykeep(:,1), '-+', ...
       xkeep, ykeep(:,2), '--x', ...
       xkeep, ykeep(:,3), ':*');
  % Add title.
  title('Projectile in Two Media Problem');
  % Label the axes.
  xlabel('x');
  ylabel('Solution');
  % Add a legend.
  legend('height','velocity','angle','Location','SouthWest');

  % Now get the range of the x and y values, and use them to make nice axes.
  xmin = min(xkeep);
  xmax = max(xkeep);
  ymin = round(min(min(ykeep))-1);
  ymax = round(max(max(ykeep))+1);
  axis([xmin xmax ymin ymax]);

  stem(xkeep, ykeep, ymin, ymax, 'determined interface');

function stem(x, y, ymin, ymax, label)
  % Find the x bin that the crossover point between y2 and y3 lies in.
  for i = 1:length(x)
    if y(i,2) > y(i,3)
      break
    end
  end
  % Use interpolation to find the crossover point more accurately.
  m(1:2) = (y(i,2:3) - y(i-1,2:3))/(x(i)-x(i-1));
  c(1:2) = y(i,2:3) - m(1:2)*x(i);
  xx = (c(2)-c(1))/(m(1)-m(2));
  hold on
  % Draw line, and add label.
  plot([xx,xx],[ymin,ymax],'k-s');
  text('Position', [xx-0.1,0.5], 'String', label, 'Rotation', 90.0)

d02sa example results

   x        y(1)        y(2)        y(3)
 0.0000    0.0000      0.5000      1.1753  
 0.5000    1.0881      0.4127      1.0977  
 1.0000    1.9501      0.3310      0.9802  
 1.5000    2.5768      0.2582      0.7918  
 2.0000    2.9606      0.2019      0.4796  
 2.5000    3.0958      0.1773      0.0245  
 3.0000    2.9861      0.1935     -0.4353  
 3.5000    2.6289      0.2409     -0.7679  
 4.0000    2.0181      0.3047     -0.9767  
 4.5000    1.1454      0.3759     -1.1099  
 5.0000   -0.0000      0.4500     -1.2000