naginterfaces.library.ode.bvp_​shoot_​genpar_​algeq

naginterfaces.library.ode.bvp_shoot_genpar_algeq(p, n, n1, pe, pf, dp, swp, icount, brkpts, bc, fcn, ymax, monlev, comm, e=None, eqn=None, constr=None, monit=None, prsol=None, data=None, io_manager=None)

bvp_shoot_genpar_algeq solves a two-point boundary value problem for a system of first-order ordinary differential equations with boundary conditions, combined with additional algebraic equations. It uses initial value techniques and a modified Newton iteration in a shooting and matching method.

For full information please refer to the NAG Library document for d02sa

https://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/d02/d02saf.html

Parameters

pfloat, array-like, shape $\left(m\right)$

$p\left[\text{i}\right]$ must be set to an estimate of the $\text{i}$th argument, $p_{\text{i}}$, for $\text{i}=1,2,\dots ,m$.

nint

$n$, the total number of differential equations.

n1int

$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 $n-n1$ equations are used as the driving equations.

pefloat, array-like, shape $\left(m\right)$

$pe[\text{i}-1]$, for $\text{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.

pffloat, array-like, shape $\left(m\right)$

$pf[\text{i}-1]$, for $\text{i}=1,2,\dots ,m$, should be set to a ‘floor’ value in the convergence test on the $i$th argument $p_i$. If $pf[i-1]\le 0.0$ on entry then it is set to the small positive value $\sqrt{ϵ}$ (where $ϵ$ 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 ($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 pe[i-1]\times max(\left|p_i\right|,pf[i-1\left]\right)\text{,}i=1,2,\dots ,m\text{,}$$

where $p_i$ is the current value of the $i$th argument.

The values $pf[i-1]$ are also used in determining the Newton iterates as discussed in Notes, see equation (6).

dpfloat, array-like, shape $\left(m\right)$

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 $dp\left[i\right]=0.0$ on entry, an estimate is made internally by setting

$$dp\left[i\right]=\sqrt{ϵ}\times max\left(pf\right[i],\left|p_i\right|)\text{,}$$

where $p_i$ is the initial value of the argument supplied by you and $ϵ$ 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 $\text{i}$, for $\text{i}=1,2,\dots ,m$, $p_{\text{i}}$ is perturbed to $p_{\text{i}}+dp\left[\text{i}\right]$ and the $\text{i}$th column of $J$ is estimated as

$$(r(p_{\text{i}}+dp[\text{i}\left]\right)-r\left(p_{\text{i}}\right))/dp\left[\text{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_{\text{i}}$ to $p_{\text{i}}-dp\left[\text{i}\right]$ and if this also fails then central differences are used with $p_{\text{i}}$ perturbed to $p_{\text{i}}+10.0\times dp\left[\text{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_{\text{i}}$, for $\text{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 $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$).

swpfloat, array-like, shape $(\text{npoint},6)$

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

$swp[i-1,2]$ must contain a lower bound for the modulus of the step size on the $\text{i}$th sub-interval $\left[x\right[\text{i}-1],x[\text{i}+1-1\left]\right]$, for $\text{i}=1,2,\dots ,\text{npoint}-1$.

If $swp[i-1,2]=0.0$ on entry, a very small default value is used.

By setting $swp[i-1,2]>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.

$swp[\text{i}-1,1]$ must contain an upper bound on the modulus of the step size to be used in the integration on $\left[x\right[\text{i}-1],x[\text{i}+1-1\left]\right]$, for $\text{i}=1,2,\dots ,\text{npoint}-1$.

If $swp[i-1,1]=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.

icountint

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

brkptscallable x = brkpts(npoint, p, data=None)

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

Parameters

npointint

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

pfloat, ndarray, shape $\left(m\right)$

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

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

xfloat, array-like, shape $\left(npoint\right)$

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

$$a=x\left[0\right]<x\left[1\right]<\cdots <x[npoint-1]=b$$

or

$$a=x\left[0\right]>x\left[1\right]>\cdots >x[npoint-1]=b\text{.}$$

bccallable (g1, g2) = bc(p, n, data=None)

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

Parameters

pfloat, ndarray, shape $\left(m\right)$

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

nint

$n$, the number of differential equations.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

g1float, array-like, shape $\left(n\right)$

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

g2float, array-like, shape $\left(n\right)$

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

fcncallable f = fcn(x, y, p, i, data=None)

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

Parameters

xfloat

$x$, the value of the argument.

yfloat, ndarray, shape $\left(n\right)$

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

pfloat, ndarray, shape $\left(m\right)$

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

iint

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

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

ffloat, array-like, shape $\left(n\right)$

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

ymaxfloat

A non-negative value which is used as a bound on all values ${\parallel y\left(x\right)\parallel }_{\infty }$ where $y\left(x\right)$ is the solution at any point $x$ between $x\left[0\right]$ and $x[\text{npoint}-1]$ 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 $ymax=0$ on entry, no bound on the solution $y$ is used; that is the integrations proceed without any checking on the size of ${\parallel y\parallel }_{\infty }$.

monlevint

Setting $monlev=0$ disables monitoring of the pseudo-Newton iteration. Setting $monlev=1$ enables this monitoring.

commdict, communication object, modified in place

Communication structure.

On initial entry: need not be set.

eNone or float, array-like, shape $\left(n\right)$, optional

Note: if this argument is None then a default value will be used, determined as follows: ${10}^{-5}$.

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

$$\left|{\text{err}}_i\right|\le e[i-1]\times max\{\sqrt{ϵ},\left|y_i\right|\}\text{,}$$

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

eqnNone or callable e = eqn(q, p, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

$eqn$ is used to describe the additional algebraic equations to be solved in the determination of the parameters, $p_{\text{i}}$, for $\text{i}=1,2,\dots ,m$.

If there are no additional algebraic equations (i.e., $m=n_1$) then $eqn$ is never called and may be None.

Parameters

qint

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

pfloat, ndarray, shape $\left(m\right)$

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

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

efloat, array-like, shape $\left(q\right)$

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.

constrNone or callable retval = constr(p, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

$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$. If None is supplied the function will act as if $constr$ returns the value $True$ at all times, and in the first instance this is recommended as the actual argument.

Parameters

pfloat, ndarray, shape $\left(m\right)$

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

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

retvalbool

Must be $True$ if the constraints are satisfied by the parameters $p$.

monitNone or callable monit(istate, iflag, ifail1, p, f, pnorm, pnorm1, eps, d, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

$monit$ enables you to monitor the values of various quantities during the calculation. If $monlev=0$ on entry no monitoring is performed and this argument is never referenced. Otherwise, if $monlev=1$ on entry, you may supply None, to use NAG-defined default monitoring, or you may supply your own custom monitoring function. The function is called by bvp_shoot_genpar_algeq after every calculation of the norm ${\Vert d^{-1}{\tilde{J}}^{+}r\Vert }_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.

Parameters

istateint

The state of the Newton iteration.

$istate=0$

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

$istate=1$ to $5$

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

$istate=6$

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

$istate=-6$ to $-1$

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

$istate=-7$

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

iflagint

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

ifail1int

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

pfloat, ndarray, shape $\left(m\right)$

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

ffloat, ndarray, shape $\left(m\right)$

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

pnormfloat

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

pnorm1float

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

epsfloat

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 ${ϵ}^{3/8}$ or greater (where $ϵ$ 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.

dfloat, ndarray, shape $\left(m\right)$

$J$, the singular values of the current modified Jacobian matrix. If $d[m-1]$ is small relative to $d\left[0\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$, $brkpts$ or $eqn$). Alternatively, the system of differential equations may be very ill-conditioned when viewed as an initial value problem, in which case bvp_shoot_genpar_algeq is unsuitable. This may also be indicated by some singular values being very large. These values of $d[\text{i}-1]$, for $\text{i}=1,2,\dots ,m$, should not be changed.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

prsolNone or callable z = prsol(z, y, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

$prsol$ can be used to obtain values of the solution $y$ at a selected point $z$ by integration across the final range $\left[x\right[0],x[npoint-1\left]\right]$.

If no output is required None can be used as the actual argument.

Parameters

zfloat

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

yfloat, ndarray, shape $\left(n\right)$

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

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

zfloat

The next point at which output is required. The new point must be nearer $x[npoint-1]$ than the old.

If $z$ is set to a point outside $\left[x\right[0],x[npoint-1\left]\right]$ the process stops and control returns from bvp_shoot_genpar_algeq to the (sub)program from which bvp_shoot_genpar_algeq is called.

Otherwise the next call to $prsol$ is made by bvp_shoot_genpar_algeq at the point $z$, with solution values $y_1,y_2,\dots ,y_n$ at $z$ contained in $y$.

If $z$ is set to $x[npoint-1]$ exactly, the final call to $prsol$ is made with $y_1,y_2,\dots ,y_n$ as values of the solution at $x[npoint-1]$ produced by the integration.

In general the solution values obtained at $x[npoint-1]$ 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 $x\left[0\right],x\left[1\right],\dots ,x[npoint-1]$ are available in the locations $swp[0,3],swp[1,3],\dots ,swp[npoint-1,3]$ at all times.

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

pfloat, ndarray, shape $\left(m\right)$

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

pffloat, ndarray, shape $\left(m\right)$

The values actually used.

dpfloat, ndarray, shape $\left(m\right)$

The values actually used.

swpfloat, ndarray, shape $(\text{npoint},6)$

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

$swp[\text{i}-1,1]$, for $\text{i}=1,2,\dots ,\text{npoint}-1$, is usually unchanged.

If the maximum step size $swp[i-1,1]$ is so small or the length of the range $\left[x\right[i-1],x[i+1-1\left]\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 $e[i-1]$ but by these other factors, $swp[\text{npoint}-1,1]$ is set to the floating-point value of $i$ if the problem last occurred in $\left[x\right[i-1],x[i+1-1\left]\right]$.

Any results obtained when this value is returned as nonzero should be viewed with caution.

$swp[\text{i}-1,2]$, for $\text{i}=1,2,\dots ,\text{npoint}-1$, are unchanged.

If an error exit with $errno$ = 4, 5 or 6 (see Exceptions) occurs on the integration made from $x[i-1]$ to $x\left[i\right]$ the floating-point value of $i$ is returned in $swp[\text{npoint}-1,0]$.

The actual point $x\in \left[x\right[i-1],x[i+1-1\left]\right]$ where the error occurred is returned in $swp[0,4]$ (see also the specification of $w$).

The floating-point value of $\text{npoint}$ is returned in $swp[\text{npoint}-1,0]$ if the error exit is caused by a call to $bc$.

If an error exit occurs when estimating the Jacobian matrix ($errno$ = 7, 8, 9, 10, 11 or 12, see Exceptions) and if argument $p_i$ was the cause of the failure then on exit $swp[\text{npoint}-1,0]$ contains the floating-point value of $i$.

$swp[\text{i}-1,3]$ contains the point $x[\text{i}-1]$, for $\text{i}=1,2,\dots ,\text{npoint}$, used at the solution $p$ or at the final values of $p$ if an error occurred.

$swp$ is also partly used as workspace.

wfloat, ndarray, shape $(max(n,m),3\times m+12+max(11,m))$

In the case of an error exit of the type where the point of failure is returned in $swp[0,4]$, the solution at this point of failure is returned in $w[\text{i}-1,0]$, for $\text{i}=1,2,\dots ,n$.

Otherwise $w$ is used for workspace.

Raises

NagValueError

(errno $1$)

On entry, $monlev=⟨value⟩$.

Constraint: $monlev=0$ or $1$.

(errno $1$)

On entry, an element of the parameter convergence control array is zero or negative.

(errno $1$)

On entry an element of the local error control array is zero or negative.

(errno $1$)

On entry, $ymax=⟨value⟩$.

Constraint: $ymax\ge 0.0$.

(errno $1$)

On entry, $\text{npoint}=⟨value⟩$.

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

(errno $1$)

On entry, $icount=⟨value⟩$.

Constraint: $icount\ge 0$.

(errno $1$)

On entry, $m=⟨value⟩$ and $n1=⟨value⟩$.

Constraint: $m\ge n1$.

(errno $1$)

On entry, $n=⟨value⟩$ and $n1=⟨value⟩$.

Constraint: $n\ge n1$.

(errno $1$)

On entry, $n1=⟨value⟩$.

Constraint: $n1\ge 1$.

(errno $2$)

The constraints have been violated by the initial parameters.

(errno $3$)

The sequence of break-points is not strictly monotonic in their initial specification.

(errno $4$)

In the integration from first to last break-point with initial or final parameters, the step size was reduced too far for the integration to proceed. Consider reversing the order of break-points or relaxing the local error control. If this error exit still results, either this function is not a suitable method for solving the problem, or the initial choice of parameters is very poor.

(errno $5$)

In the integration from first to last break-point with initial or final parameters, a suitable initial step could not be found on one of the intervals. Consider reversing the order of break-points or relaxing the local error control. If this error exit still results, either this function is not suitable for solving the problem, or the initial choice of parameters is very poor.

(errno $6$)

During integration the solution exceeded $ymax$ in magnitude, where, on entry, $ymax=⟨value⟩$.

(errno $7$)

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

(errno $8$)

On perturbing the parameters when calculating the initial approximation to the Jacobian, the monotonicity condition on break-points is violated.

(errno $9$)

An initial step-length could be found for integration to proceed with the current parameters.

(errno $10$)

The step-length required to calculate the Jacobian to sufficient accuracy is too small.

(errno $11$)

On calculating the initial approximation to the Jacobian, the solution to the system of differential equations exceeded $ymax$ in magnitude, where, on entry, $ymax=⟨value⟩$.

(errno $12$)

The Jacobian has an insignificant column. Make sure that the solution vector depends on all the parameters.

(errno $13$)

An internal singular value decomposition has failed. This error can be avoided by changing the initial parameter estimates.

(errno $14$)

The Newton iteration has failed to converge. This can indicate a poor initial choice of parameters or a very difficult problem. Consider varying elements of the parameter convergence control if the residuals are small; otherwise vary initial parameter estimates.

(errno $15$)

The number of iterations permitted by $icount$ has been exceeded where this was set to a positive value on entry. $icount=⟨value⟩$.

(errno $16$)

Internal error in calculating Jacobian. Please contact NAG.

(errno $17$)

Internal error in calculating residual. Please contact NAG.

(errno $18$)

Internal error in integration. Please contact NAG.

(errno $19$)

Internal error in Newton method. Please contact NAG.

Notes

No equivalent traditional C interface for this routine exists in the NAG Library.

bvp_shoot_genpar_algeq solves a two-point boundary value problem for a system of $n$ first-order ordinary differential equations with separated boundary conditions by determining certain unknown arguments $p_1,p_2,\dots ,p_m$. (There may also be additional algebraic equations to be solved in the determination of the arguments and, if so, these equations are defined by $eqn$.) The arguments may be, but need not be, boundary values; they may include eigenvalues, arguments in the coefficients of the differential equations, coefficients in series expansions or asymptotic expansions for boundary values, the length of the range of definition of the system of differential equations, etc.

It is assumed that we have a system of $n$ differential equations of the form

$$y^{\prime }=f(x,y,p)\text{,}$$

where $p={(p_1,p_2,\dots ,p_m)}^T$ is the vector of arguments, and that the derivative $f$ is evaluated by $fcn$. Also, $n_1$ of the equations are assumed to depend on $p$. For $n_1<n$ the $n-n_1$ equations of the system are not involved in the matching process. These are the driving equations; they should be independent of $p$ and of the solution of the other $n_1$ equations. In numbering the equations in $fcn$ and $bc$ the driving equations must be put first (as they naturally occur in most applications). The range of definition [$a,b$] of the differential equations is defined by $brkpts$ and may depend on the arguments $p_1,p_2,\dots ,p_m$ (that is, on $p$). $brkpts$ must define the points $x_1,x_2,\dots ,x_{\text{npoint}}$, $\text{npoint}\ge 2$, which must satisfy

$$a=x_1<x_2<\cdots <x_{\text{npoint}}=b$$

(or a similar relationship with all the inequalities reversed).

If $\text{npoint}>2$ the points $x_1,x_2,\dots ,x_{\text{npoint}}$ can be used to break up the range of definition. Integration is restarted at each of these points. This means that the differential equations (1) can be defined differently in each sub-interval $[x_{\text{i}},x_{\text{i}+1}]$, for $\text{i}=1,2,\dots ,\text{npoint}-1$. Also, since initial and maximum integration step sizes can be supplied on each sub-interval (via the array $swp$), you can indicate parts of the range $[a,b]$ where the solution $y\left(x\right)$ may be difficult to obtain accurately and can take appropriate action.

The boundary conditions may also depend on the arguments and are applied at $a=x_1$ and $b=x_{\text{npoint}}$. They are defined (in $bc$) in the form

$$y\left(a\right)=g_1\left(p\right),>y\left(b\right)=g_2\left(p\right)\text{.}$$

The boundary value problem is solved by determining the unknown arguments $p$ by a shooting and matching technique. The differential equations are always integrated from $a$ to $b$ with initial values $y\left(a\right)=g_1\left(p\right)$. The solution vector thus obtained at $x=b$ is subtracted from the vector $g_2\left(p\right)$ to give the $n_1$ residuals $r_1\left(p\right)$, ignoring the first $n-n_1$, driving equations. Because the direction of integration is always from $a$ to $b$, it is unnecessary, in $bc$, to supply values for the first $n-n_1$ boundary values at $b$, that is the first $n-n_1$ components of $g_2$ in (3). For $n_1<m$ then $r_1\left(p\right)$. Together with the $m-n_1$ equations defined by $eqn$,

$$r_2\left(p\right)=0\text{,}$$

these give a vector of residuals $r$, which at the solution, $p$, must satisfy

$$\begin{array}{c}\left(\begin{array}{c}{r}_1\left(p\right)\\ r_2\left(p\right)\end{array}\right)=0\text{.}\end{array}$$

These equations are solved by a pseudo-Newton iteration which uses a modified singular value decomposition of $J=\frac{\partial r}{\partial p}$ when solving the linear equations which arise. The Jacobian $J$ used in Newton’s method is obtained by numerical differentiation. The arguments at each Newton iteration are accepted only if the norm ${\Vert D^{-1}{\tilde{J}}^{+}r\Vert }_2$ is much reduced from its previous value. Here ${\tilde{J}}^{+}$ is the pseudo-inverse, calculated from the singular value decomposition, of a modified version of the Jacobian $J$ ($J^{+}$ is actually the inverse of the Jacobian in well-conditioned cases). $D$ is a diagonal matrix with

$$d_{ii}=max(\left|p_i\right|,pf\left[i\right])$$

where $pf$ is an array of floor values.

See Deuflhard (1974) for further details of the variants of Newton’s method used, Gay (1976) for the modification of the singular value decomposition and Gladwell (1979) for an overview of the method used.

Two facilities are provided to prevent the pseudo-Newton iteration running into difficulty. First, you are permitted to specify constraints on the values of the arguments $p$ via a $constr$. These constraints are only used to prevent the Newton iteration using values for $p$ which would violate them; that is, they are not used to determine the values of $p$. Secondly, you are permitted to specify a maximum value $y_{max}$ for ${\parallel y\left(x\right)\parallel }_{\infty }$ at all points in the range $[a,b]$. It is intended that this facility be used to prevent machine ‘overflow’ in the integrations of equation (1) due to poor choices of the arguments $p$ which might arise during the Newton iteration. When using this facility, it is presumed that you have an estimate of the likely size of ${\parallel y\left(x\right)\parallel }_{\infty }$ at all points $x\in [a,b]$. $y_{max}$ should then be chosen rather larger (say by a factor of $10$) than this estimate.

You are strongly advised to set $monlev$ to monitor the progress of the pseudo-Newton iteration. You can output the solution of the problem $y\left(x\right)$ by supplying a suitable $prsol$.

bvp_shoot_genpar_algeq is designed to try all possible options before admitting failure and returning to you. Provided the function can start the Newton iteration from the initial point $p$ it will exhaust all the options available to it (though you can override this by specifying a maximum number of iterations to be taken). The fact that all its options have been exhausted is the only error exit from the iteration. Other error exits are possible, however, whilst setting up the Newton iteration and when computing the final solution.

If you require more background information about the solution of boundary value problems by shooting methods you are recommended to read the appropriate modules of Hall and Watt (1976), and for a detailed description of bvp_shoot_genpar_algeq Gladwell (1979) is recommended.

References

Deuflhard, P, 1974, A modified Newton method for the solution of ill-conditioned systems of nonlinear equations with application to multiple shooting, Numer. Math. (22), 289–315

Gay, D, 1976, On modifying singular values to solve possibly singular systems of nonlinear equations, Working Paper 125, Computer Research Centre, National Bureau for Economics and Management Science, Cambridge, MA

Gladwell, I, 1979, The development of the boundary value codes in the ordinary differential equations module of the NAG Library, Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science, (eds B Childs, M Scott, J W Daniel, E Denman and P Nelson) (76), Springer–Verlag

Hall, G and Watt, J M (ed.), 1976, Modern Numerical Methods for Ordinary Differential Equations, Clarendon Press, Oxford