Shooting for Boundary-Value Problems

Shooting replaces a boundary-value problem by an initial-value integration and a root search over missing initial data. It is attractive because it reuses mature ODE solvers, but it inherits the growth and decay of the underlying fundamental solutions. A boundary residual can be tiny even when the required initial amplitude has lost all reliable digits. The method must therefore be judged by sensitivity and interior residuals, not by endpoint agreement alone.

Required background. Ordinary Differential Equations supplies initial-value flow maps, boundary-value problems, and root finding.

Helpful background. Formulations of General Relativity provides constrained relativistic systems that reduce to radial boundary problems, while Partial Differential Equations and Green Functions supplies boundary conditions and separated operators.

From a boundary condition to a root

For a first-order system $y'=f(x,y,p)$ with left boundary $B_L(y(a),p)=0$ and right boundary $B_R(y(b),p)=0$, parameterize the admissible left data by a vector $s$. Integrate to $b$ and define

$$F(s)=B_R\bigl(y(b;s),p\bigr).$$

Single shooting solves $F(s)=0$. Newton’s method needs the sensitivity matrix $\partial y(b;s)/\partial s$, obtained by finite differences, automatic differentiation, or the variational equation. The last option evolves

$$S'=f_yS, \qquad S(a)=\frac{\partial y(a;s)}{\partial s}.$$

Here $p$ is fixed. If $p=p(s)$ is also unknown, the sensitivity equation gains the forcing $f_p\,\partial p/\partial s$, and differentiating $F$ also includes the explicit parameter dependence of $B_R$. Count one additional boundary or normalization condition for each new parameter. Single shooting and its sensitivity equations are developed in Ascher, Mattheij, and Russell 1995, Chapter 4 § 4.1, pp. 132–135, while the conditioning distinction used here is treated in Chapter 3 § 3.4, pp. 110–126. The sensitivity also makes missing modes visible: if the parameterization of left data cannot span the boundary residual, $F'(s)$ is rank deficient.

A well-conditioned boundary problem with unstable shooting

Consider

$$y''=\lambda^2y, \qquad y(0)=1, \qquad y(1)=e^{-\lambda},$$

with exact solution

$$y(x)=e^{-\lambda x}, \qquad \lambda\gt0.$$

Writing $y'(0)=s$ gives

$$y(x;s)=\cosh(\lambda x)+\frac{s}{\lambda}\sinh(\lambda x), \qquad s_*=-\lambda.$$

At the correct slope the growing exponentials cancel. A slope perturbation $\delta s$ instead produces

$$\delta y(x)=\frac{\delta s}{\lambda}\sinh(\lambda x), \qquad F'(s)=\frac{\sinh\lambda}{\lambda}.$$

At $\lambda=50$, $F'(s)\simeq5.18\times10^{19}$. Even a tiny absolute error in the initial slope or an integration step can seed the growing solution. Evaluate the exact reference as $e^{-\lambda x}$, not as the difference of two large hyperbolic terms.

This is instability of the shooting representation, not exponential ill-conditioning of the boundary-value problem in the maximum norm. With endpoint perturbations $\delta a,\delta b$, the exact solution perturbation is

$$\delta y(x)=\delta a\frac{\sinh(\lambda(1-x))}{\sinh\lambda} +\delta b\frac{\sinh(\lambda x)}{\sinh\lambda}.$$

Both weights are nonnegative and their sum is at most one, so $\lVert\delta y\rVert_\infty\le\max(\lvert\delta a\rvert,\lvert\delta b\rvert)$. A large absolute derivative of the shooting map is also not, by itself, the condition number of its scalar inverse. The lost cancellation and directional error growth are what this example tests.

Run single shooting with an adaptive high-order integrator at $\mathrm{rtol}=10^{-8},10^{-12},10^{-16}$ and componentwise absolute tolerances $(\mathrm{atol}_y,\mathrm{atol}_{y'})=\mathrm{rtol}(10^{-6},10^{-6}\lambda)$, using 53-bit, 100-bit, and 200-bit arithmetic. The $10^{-16}$ request is a saturation test in binary64, not an expectation that it deliver sixteen reliable digits. These are absolute-error experiments on a profile of maximum value one; they do not request relative accuracy in the exponentially small right endpoint. Report the absolute endpoint residual, the maximum solution error against the exact formula on an independent 2001-point grid, the differential residual $y''-\lambda^2y$, the sensitivity $F'(s)$, and stable digits of $s$. Increasing only the integrator tolerance may not recover digits once arithmetic precision or unstable-mode contamination dominates.

Multiple shooting and collocation

Multiple shooting divides $[0,1]$ into subintervals, treats each initial state as an unknown, and imposes continuity at interfaces. Each growth factor is then roughly $e^{\lambda\Delta x}$ rather than $e^\lambda$. The nonlinear system is larger but better scaled. Single and multiple shooting are developed in Ascher, Mattheij, and Russell 1995, Chapter 4 §§ 4.1–4.3, pp. 132–154; the relation between boundary-value conditioning and multiple-shooting stability is derived independently in Mattheij 1982, §§ 3–5, pp. 966–974. Collocation avoids directional integration altogether by enforcing the ODE throughout the domain; its construction and comparison with shooting are treated in Ascher, Mattheij, and Russell 1995, Chapter 5 § 5.4, pp. 218–225, and §§ 5.6.2–5.6.4, pp. 247–265.

Compare single shooting, ten-interval multiple shooting, and collocation on the exact benchmark. Hold the solution norm and boundary tolerance definitions fixed. Refine the number of shooting intervals or collocation nodes and vary arithmetic precision independently. A credible result exhibits decreasing interior error and residual, stable left slope, and agreement among methods. Endpoint residual alone is weak: Newton can force $y(1)=1$ while the interior profile is inaccurate.

Parameters, resonances, and branches

In relativistic applications the root variable may include a frequency or central value. Complex resonance shooting is particularly sensitive because outgoing and ingoing asymptotic solutions can differ exponentially. State the $e^{-i\omega t}$ convention, complex sheet, endpoint approximation, and normalization. Follow a branch continuously from a known parameter value rather than selecting the nearest root independently at each step.

The validity regime excludes endpoints where the asymptotic basis is inaccurate and parameter values where the chosen left-data basis becomes singular. Vary both endpoints and compare with a collocation or spectral calculation. A root that moves under modest endpoint changes is not established.

Reproducible execution and failure criteria

Pin the ODE library and version, method, tolerances, arithmetic precision, Jacobian construction, root solver, interval partition, scaling, and exact input equations. A clean rerun emits the initial slope, endpoint and interior residuals, condition estimate, work counters, refinement and precision sweeps, and comparison with the analytic solution. No absent program is implied to have performed this sequence.

Reject shooting when Newton converges to a wrong branch, the sensitivity matrix is singular, stable digits decrease under a tighter tolerance, or agreement appears only at endpoints. Other failure modes are overflow in a growing basis, loss of a decaying mode, finite-difference sensitivities dominated by cancellation, and a normalization that becomes zero. Multiple shooting reduces growth per interval; it does not repair an ill-posed continuum boundary condition.

Check your understanding

For $\lambda=50$, a single-shooting run reports $\lvert y(1)-e^{-50}\rvert\lt10^{-13}$ but has a maximum interior error of $10^{-3}$. Which diagnostic decides whether the solution is acceptable?

Solution

The endpoint residual does not decide it; the run fails the interior error and differential-residual checks. Increase precision, inspect the sensitivity, and compare multiple shooting or collocation on a refinement sequence. Accept only when the interior profile converges and the recovered initial slope has stable digits as well as satisfying both boundaries.

References

• Ascher, U. M., Mattheij, R. M. M., and Russell, R. D. (1995), Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, SIAM. DOI.
• Mattheij, R. M. M. (1982), “The Conditioning of Linear Boundary Value Problems,” SIAM Journal on Numerical Analysis 19, 963–978. DOI.

Further reading

• Hairer, E., Nørsett, S. P., and Wanner, G. (1993), Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd revised edition, Springer. DOI.