Errors

We have framed our discretization in terms of \(\epsilon\) but, in practice, a better choice would be \(\varepsilon = \Delta s/s_{\text{max}}\) since this more directly correlates with the processing and memory requirements. Therefore we convert a specified \(\varepsilon\) into a value of \(\epsilon\) by inverting \[\varepsilon = h\, \epsilon^{s_0}/\log \epsilon^{-s_0} \; ,\] which can either be done numerically or else using the approximation \[\epsilon^{s_0} \simeq -\frac{L_0}{6 L_1^3}\, \left[ 6 L_1^4 + 6 L_1^2 L_2(L_1 + 1) - 3 L_1 L_2 (L_2 - 2) + L_2 (2 L_2^2 - 9 L_2 + 6) \right] \; ,\] where \[L_0 \equiv \varepsilon/h \quad , \quad L_1 \equiv \log(L_0) \quad , \quad L_2 \equiv \log(-L_1) \; ,\] which is valid for \(L_0 \lesssim 0.35\). The discretizations in Fig. \ref{fig:Ggrid} are calculated with a fixed \(\varepsilon = 0.01\), and so have corresponding values of \(\epsilon\) that increase with \(\ell\).

In order to provide more direct control of the transform accuracy, we implement a driver routine that starts from some initial \(\varepsilon\) (nominally \(0.01\)) and then reduces its value by factors until¹ \[\left| f_{\ell}(r_j; \varepsilon) - f_{\ell}(r_j; 2\varepsilon)\right| < \max\left(e_{\text{abs},\ell}\cdot r_j^p, e_{\text{rel},\ell}\cdot\left|f_{\ell}(r_j; \varepsilon)\right|\right)\] for all \(j\), where the \(r_j\) are equally spaced over the range \(r_{\text{min}} \le r \le r_{\text{max}}\) and \(p\) allows for an \(r\)-weighing of the ‘constant’ error criterion (but we normally use \(p = 0\)). We abandon the halving of \(\varepsilon\) once it reaches some prescribed minimum value (nominally \(10^{-6}\)), to protect against pathological cases.

For the purposes of transforming a set of multipoles that are subsequently recombined using equation \ref{eqn:multipole2}, \ref{eqn:multipole3}, or \ref{eqn:multipole3mu}, we require that² \[e_{\text{rel},\ell} = \frac{e_{\text{rel}}}{M}\, \frac{\left|\sum_{\ell'} f_{\ell'}(r) L_{\ell'}(\mu_r)\right|} {\left|f_{\ell}(r) L_{\ell}(\mu_r)\right|} \quad , \quad e_{\text{abs},\ell} = \frac{e_{\text{abs}}}{M} \; ,\] for all \((r,\mu_r)\) on a grid covering the region of interest, where \(M\) is the number of multipoles included in the expansion.