Numerical Methods and Special Functions
Source:vignettes/numerical-foundations/numerical-foundations.Rmd
numerical-foundations.RmdIntroduction
Canonical scattering models replace a field defined throughout space with a series of geometry-matched eigenfunctions. This is the principal numerical advantage of modal scattering theory (Waterman 2009). It does not eliminate numerical approximation. The infinite series must be truncated, special functions must be evaluated in finite precision, overlap integrals may require quadrature, and the resulting coefficient systems may be ill-conditioned.
This page develops those numerical ingredients and the other integration, projection, and summation methods used by the package. Model-specific derivations remain on their theory pages.
See Notation and Symbols for the shared definitions of k, ka, J_m, j_n, P_n^m, S_{mn}, and R_{mn}^{(q)}.
Geometry-matched modal bases
Let \{\Psi_\nu\} be a basis adapted to the target boundary. A separated field is represented formally as:
p(\boldsymbol{x}) = \sum_{\nu=0}^{\infty} a_\nu\Psi_\nu(\boldsymbol{x}).
For a sphere, cylinder, or spheroid, the boundary is a coordinate surface. Orthogonality then decouples the boundary conditions by mode or reduces them to small blocks. A penetrable spheroid is an important exception. Its interior and exterior angular bases generally have different spheroidal parameters, so projection between those bases produces a coupled system.
The basis is exact for the stated canonical geometry. A numerical calculation is still finite because it retains only a subset of the coefficients a_\nu. The errors of concern are representation error from an inappropriate geometry, truncation error from too few modes, evaluation error in the special functions, quadrature error, and error amplified by an ill-conditioned solve.
Cylindrical and spherical functions
Bessel and Hankel functions
Separation of the Helmholtz equation in cylindrical coordinates gives the cylindrical Bessel equation, while spherical separation gives its spherical counterpart (Folver and Maximon 2026). For order m and degree n, respectively, the equations are:
\begin{aligned} z^2 y'' + z y' + \left(z^2-m^2\right)y &= 0, \\ z^2 y'' + 2z y' + \left[z^2-n(n+1)\right]y &= 0. \end{aligned}
The spherical functions can be evaluated through half-integer cylindrical functions. For example:
j_n(z) = \sqrt{\frac{\pi}{2z}}J_{n+1/2}(z), \qquad y_n(z) = \sqrt{\frac{\pi}{2z}}Y_{n+1/2}(z).
Regular incident or interior fields use J_m or j_n. Outgoing scattered fields use the corresponding Hankel combinations:
H_m^{(1)}(z) = J_m(z) + iY_m(z), \qquad h_n^{(1)}(z) = j_n(z) + iy_n(z).
Boundary conditions also require derivatives of these functions. At high order or near zeros, ratios of a function and its derivative can lose accuracy through underflow, overflow, or cancellation. Stable recurrences and scaled representations are therefore as important as the defining equations (Abramowitz and Stegun 1964).
Legendre functions and recurrences
The angular factor in spherical separation is an associated Legendre function. It satisfies:
\frac{d}{d\mu}\left[(1-\mu^2)\frac{dP_n^m}{d\mu}\right] + \left[n(n+1)-\frac{m^2}{1-\mu^2}\right]P_n^m = 0.
For the axisymmetric case, the Legendre polynomials obey the orthogonality relation:
\int_{-1}^{1}P_n(\mu)P_\ell(\mu)\,d\mu = \frac{2}{2n+1}\delta_{n\ell}.
This relation isolates spherical modal coefficients after the boundary conditions are projected onto P_n. Associated functions provide the corresponding structure when azimuthal order is nonzero (Dunster 2026).
Integer-degree Legendre values are generated efficiently through the three-term recurrence:
(n+1)P_{n+1}(\mu) = (2n+1)\mu P_n(\mu)-nP_{n-1}(\mu).
The package’s general Legendre utilities also support cases outside the integer modal path. Those branches use hypergeometric series, numerical integration, or finite-difference derivatives as needed. Values near branch points or singular endpoints deserve additional scrutiny because those general evaluators do not have the same numerical behavior as the integer recurrences.
Prolate spheroidal wave functions
Prolate spheroidal coordinates use an angular coordinate -1\leq\eta\leq 1, a radial coordinate \xi\geq 1, and a spheroidal parameter c. Separation produces angular and radial functions with the same order m, degree n, and eigenvalue \lambda_{mn}(c) (Flammer 1957).
The angular function satisfies:
\frac{d}{d\eta}\left[(1-\eta^2)\frac{dS_{mn}}{d\eta}\right] + \left[\lambda_{mn}(c)-c^2\eta^2- \frac{m^2}{1-\eta^2}\right]S_{mn} = 0.
The radial function satisfies:
\frac{d}{d\xi}\left[(\xi^2-1)\frac{dR_{mn}}{d\xi}\right] - \left[\lambda_{mn}(c)-c^2\xi^2+ \frac{m^2}{\xi^2-1}\right]R_{mn} = 0.
When c\rightarrow 0, the angular eigenvalue approaches n(n+1) and the angular functions approach associated Legendre functions. At nonzero c, both the eigenvalue and the basis depend on the acoustic size. This dependence is one reason spheroidal evaluation is more demanding than evaluation of a fixed Legendre basis (Volkmer 2026).
For fixed m and c, the angular functions are orthogonal under the normalization used by the surrounding derivation. Writing the squared norm as N_{mn}(c), the relation is:
\int_{-1}^{1}S_{mn}(c,\eta)S_{m\ell}(c,\eta)\,d\eta = N_{mn}(c)\delta_{n\ell}.
The first radial kind is regular and the second is the independent singular solution. The third and fourth kinds are formed as:
R_{mn}^{(3)} = R_{mn}^{(1)} + iR_{mn}^{(2)}, \qquad R_{mn}^{(4)} = R_{mn}^{(1)} - iR_{mn}^{(2)}.
The package uses the third kind for an outgoing spheroidal field. The PSMS theory page gives the coordinate definitions, field expansions, normalization factors, and boundary projections.
Numerical implementation
The Smn() and Rmn() functions validate
inputs and dispatch to compiled code. They are not themselves the
spheroidal algorithm. The numerical methods come from the Van Buren and
Boisvert profcn program, including the Bouwkamp eigenvalue
method and alternative expansions for radial functions over different
parameter ranges (Van Buren and Boisvert
2002). The treatment of the second radial kind is described
separately by Van Buren and Boisvert (2004).
The actual special-function implementation used by
acousticTS is the vendored and package-adapted profcn
source in src/prolate_swf.f90. It calculates
eigenvalues, angular functions, radial functions of the first and second
kinds, derivatives, normalization quantities, and accuracy indicators.
The same file is compiled in both double and quadruple precision, as
shown in src/Makevars.
The rest of the source tree connects that kernel to R:
-
R/spheroidal.Rprovides user-facing validation and dispatch. -
src/psms.cppdefines the R-to-C++ entry points. -
src/psms_smn.hextracts angular values and derivatives. -
src/psms_rmn.hextracts radial kinds and forms the complex incoming and outgoing functions. -
src/psms_support.hhandles batched calls and base-10 exponent rescaling.
The vendored kernel was adapted for the package, including batching
and cached support calculations. Its numerical ancestry and the original
program remain available in the upstream
prolate_swf repository. For an audit of the
special-function mathematics, begin with
src/prolate_swf.f90, then use the C++ headers to trace how
its outputs enter the package models.
Many backend results are stored as a mantissa and a base-10 exponent. The C++ support layer reconstructs a returned component as:
x = \widehat{x}\,10^{e}.
Keeping \widehat{x} and e separate during the difficult part of the calculation extends the usable dynamic range. Batching consecutive degrees for a fixed order also avoids repeating the expensive eigenvalue and expansion setup for each (m,n) pair.
Modal truncation and convergence
An infinite modal result has the form:
F = \sum_{n=0}^{\infty} a_n.
The computed result retains a finite upper order N:
F_N = \sum_{n=0}^{N} a_n, \qquad E_N = F-F_N = \sum_{n=N+1}^{\infty}a_n.
The tail E_N is not normally known. Acoustic size supplies a starting scale because progressively higher orders become important as ka increases. The current package defaults are model-specific:
| Modal calculation | Starting or hard upper order |
|---|---|
| Spherical fluid and related boundaries | N=\operatorname{round}(ka+20) |
| Elastic spherical shell | N=\operatorname{round}(ka)+10 |
| Viscous-layer elastic spherical shell | N=\max\{2,\operatorname{round}(ka)+10\} |
| Finite cylinder | N=\lceil ka\rceil+10 |
| Solid elastic sphere | N_0=\operatorname{round}(ka)+10 before adaptive extension |
| Prolate spheroid | m_{\max}=\lceil 2ka\rceil, n_{\max}=m_{\max}+\lceil\chi_1/2\rceil |
Here a is the relevant exterior radius and \chi_1=k_1q is the exterior spheroidal parameter based on semifocal length q. These formulas are initial resolution choices, not mathematical error bounds. Material contrast, incidence angle, resonances, and cancellation can all require additional resolution.
When two truncations can be compared, a useful convergence diagnostic is:
\varepsilon_N = \frac{|F_{N+\Delta N}-F_N|} {\max\!\left(|F_{N+\Delta N}|,F_{\mathrm{scale}}\right)}.
The positive scale F_{\mathrm{scale}} prevents division by a value near a physical null. A small \varepsilon_N supports convergence with respect to that refinement. It does not establish that the physical model or geometry is appropriate.
Adaptive summation applies the same idea to individual modal terms. The solid elastic-sphere calculation begins from its usual ka-based order and extends the series until its tail term is sufficiently small. The spheroidal solver can also suppress negligible modal tails and select quadrature order from the retained angular content. Fixed controls remain useful for reproducibility and for explicit convergence studies.
Quadrature and continuous integrals
Gauss-Legendre overlap quadrature
Interior and exterior spheroidal angular functions have different parameters when the wave speeds differ. Their cross-products are not diagonal under the ordinary orthogonality relation, so the penetrable problem contains overlap integrals such as:
K_{n\ell}^{(m)} = \int_{-1}^{1} S_{mn}(c_1,\eta)S_{m\ell}(c_2,\eta)\,d\eta.
The package evaluates these finite-interval integrals with Gauss-Legendre quadrature. For Q nodes, the approximation is:
\int_{-1}^{1}f(\eta)\,d\eta \approx \sum_{q=1}^{Q}w_q f(\eta_q).
The nodes are the roots of P_Q, and their weights are:
P_Q(\eta_q)=0, \qquad w_q=\frac{2} {(1-\eta_q^2)\left[P_Q'(\eta_q)\right]^2}.
The compiled quadrature routine obtains the nodes by Newton iteration starting from Chebyshev-like estimates. Gauss-Legendre quadrature is exact for polynomials through degree 2Q-1, but an overlap of oscillatory spheroidal functions is not a polynomial. Quadrature order must therefore increase with the angular content that the modal system retains (Temme 2026).
Adaptive integration of oscillatory amplitudes
Weak-scattering and curved-ray formulations contain continuous line integrals whose phase can oscillate rapidly. Their common numerical form is:
F(k)=C(k)\int_a^b A(s,k)e^{i\Phi(s,k)}\,ds.
The compiled DWBA path linearly interpolates the stored centerline
and radius within each body segment, then evaluates the real and
imaginary parts with the adaptive QAGS routine from QUADPACK. The
curved-cylinder Fresnel calculation uses the same real-and-imaginary
decomposition through R’s integrate() interface. QAGS
combines Gauss-Kronrod error estimates with interval subdivision and
extrapolation (Piessens et al. 1983).
Relative tolerance, absolute tolerance, and the subdivision limit control the integration error. Geometry resolution is a separate issue. Tightening the quadrature tolerance cannot recover curvature or radius variation that was lost when the body was represented by too few segments.
Projection and linear systems
Transition-matrix projection
A transition matrix maps incident coefficients to scattered coefficients (Waterman 1969). In block form, the map is:
\boldsymbol{b}^{(m)} =\boldsymbol{T}^{(m)}\boldsymbol{a}^{(m)}.
For spherical-coordinate branches, the package evaluates the boundary operator on Gauss-Legendre surface nodes and projects the residual back onto the retained angular basis. A projected boundary equation has the form:
\left\langle \Psi_\mu, \mathcal{B}\!\left[p_{\mathrm{inc}}+p_{\mathrm{scat}}\right] \right\rangle_Q=0, \qquad \mu=1,\ldots,N.
The subscript Q indicates that the surface inner product is evaluated by quadrature. Increasing basis order without enough surface nodes can therefore produce a larger algebraic system that represents the boundary no more accurately. The current finite-cylinder branch uses a cylindrical modal basis rather than forcing its sidewall-endcap corner into this spherical surface projection.
Dense boundary systems
After truncation and projection, a coupled boundary problem has the matrix form:
\boldsymbol{K}\boldsymbol{a}=\boldsymbol{b}.
This structure occurs in penetrable spheroidal kernels, projected T-matrix blocks, and the modewise viscous-layer shell calculation. The latter solves a 6\times6 monopole system and a 10\times10 system for each higher spherical order. Direct solves are used first where appropriate, with more robust solves or a pseudoinverse used as fallbacks.
If columns of \boldsymbol{K} become nearly dependent, small evaluation or quadrature errors can produce large changes in \boldsymbol{a}. With the singular value decomposition, the system matrix is written as:
\boldsymbol{K} =\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{*}.
For the double-precision dense spheroidal kernel solve, the package constructs a pseudoinverse by retaining singular values above the cutoff:
\tau = N\,\sigma_{\max}\,\epsilon_{\mathrm{mach}}, \qquad \boldsymbol{K}^{+} =\boldsymbol{V}\boldsymbol{\Sigma}_{\tau}^{+}\boldsymbol{U}^{*}.
Here N is the block dimension and singular values no larger than \tau are assigned a zero reciprocal. This prevents the solver from amplifying numerically unresolved directions without bound. The package relies on LAPACK and Armadillo for its standard dense decompositions (Anderson et al. 1999).
A common summary of sensitivity is the two-norm condition number:
\kappa_2(\boldsymbol{K}) =\frac{\sigma_{\max}}{\sigma_{\min}}.
When \sigma_{\min} is small, increasing precision may reduce roundoff, but it cannot repair an under-resolved truncation or quadrature rule. The quadruple-precision spheroidal path therefore combines higher precision with a complete-pivoting linear solve and iterative refinement. It is a numerical resolution option for difficult parameter ranges, not a different scattering model (Higham 2002).
Summation, segmentation, and stochastic convergence
Cancellation and compensated sums
Modal far fields and segmented-body approximations add complex terms with different phases. If large terms nearly cancel, ordinary left-to-right summation can discard low-order digits. The spheroidal far-field code uses a Kahan-style compensated update. For a running sum s_j and compensation c_j, one step is:
\begin{aligned} y_j &= a_j-c_j, \\ t_j &= s_j+y_j, \\ c_{j+1} &= (t_j-s_j)-y_j, \\ s_{j+1} &= t_j. \end{aligned}
Compensation reduces rounding error in the retained sum. It does not make an insufficient modal cutoff converge (Higham 2002).
Segmented geometries
Several elongated-body models replace a continuous body with axial segments. Their coherent amplitude has the discrete form:
F_N(k)=\sum_{s=1}^{N}F_s(k).
The segment count controls how well the stored centerline, radius, curvature, and phase are resolved. This affects DWBA-based models, phase-compensated variants, and the composite fish-body calculation. A segment-refinement study should keep the physical body fixed, resample it more finely, and compare the complex amplitude before converting it to target strength. Agreement only in dB can conceal compensating phase errors near a null.
Stochastic realizations
The stochastic DWBA adds random phase perturbations to segment contributions and averages linear backscatter over M realizations (Demer and Conti 2003). For realization values X_r, the Monte Carlo mean and its estimated standard error are:
\overline{X}_M=\frac{1}{M}\sum_{r=1}^{M}X_r, \qquad \operatorname{SE}(\overline{X}_M) =\frac{s_X}{\sqrt{M}}.
The familiar M^{-1/2} decrease is statistical rather than deterministic. Reproducibility requires recording the random seed and the number of realizations. Numerical convergence also requires checking segment resolution because increasing M cannot correct a poorly resolved body.
Interpreting numerical convergence
A numerical result should be stable under moderate, targeted refinement. The relevant check depends on the calculation:
- increase the modal upper order and compare the complex amplitude or target strength,
- increase quadrature order when modal overlaps or surface projections are present,
- tighten adaptive integration tolerances for continuous line integrals,
- refine segmented geometry without changing the underlying target,
- increase stochastic realizations and report the seed for stochastic models, and
- compare double and quadruple precision when the spheroidal system is sensitive.
Agreement under one refinement tests only that numerical choice. Agreement between two truncations does not validate material properties, boundary conditions, or the use of a canonical geometry. Conversely, an oscillatory spectrum is not automatically a numerical artifact. Resonances and interference are physical when their locations and amplitudes persist under the relevant refinements.