2. Axisymmetric virtual-casing principle

2.2. Virtual-casing principle for axisymmetric plasmas

Consider an axisymmetric plasma confined by external coils in equilibrium. Let $\boldsymbol {B}$ denote the total magnetic field. The poloidal field $\boldsymbol {B}^{{\rm pol}} = \boldsymbol {B} - B_{\phi } \boldsymbol {e}_{\phi }$ at any location is the sum of the poloidal field $\boldsymbol {B}_{{\rm ext}}^{{\rm pol}}$, due to the external coils, and of the poloidal field $\boldsymbol {B}_V^{{\rm pol}}$, due to the plasma current. The field $\boldsymbol {B}_V^{{\rm pol}}$ is given for all $\boldsymbol {x} \in \mathbb {R}^3$ by the Biot–Savart law:

(2.1)\begin{equation} \boldsymbol{B}_V^{{\rm pol}}(\boldsymbol{x})=\frac{\mu_{0}}{4{\rm \pi}} \iiint_{\varOmega}J^{{\rm tor}}_{V}(\boldsymbol{y}) \boldsymbol{e}_\phi(\phi(\boldsymbol{y})) \times \frac{\boldsymbol{x} - \boldsymbol{y}}{\| \boldsymbol{x} - \boldsymbol{y}\|^3} \mathrm{d} \boldsymbol{y}, \end{equation}

where $\mu _{0}$ is the permeability of free space, the integration volume $\varOmega$ is the plasma domain and $J^{{\rm tor}}_{V}$ is the toroidal current density in the plasma. The vector $\boldsymbol {e}_{\phi }(\phi (\boldsymbol {y}))$ is the unit vector in the toroidal direction at $\boldsymbol {y}$, where the argument $\phi (\boldsymbol {y})$ is the toroidal angle of the point $\boldsymbol {y}$.

Equation (2.1) is a volume integral, which is expensive to evaluate numerically. The virtual-casing principle gives a formula for $\boldsymbol {B}^{{\rm pol}}_V$ that depends only on the full field $\boldsymbol {B}^{{\rm pol}}$ at the plasma boundary and only requires the evaluation of a surface integral (i.e. line integral for axisymmetric domains) (Shafranov & Zakharov Reference Shafranov and Zakharov1972; Zakharov Reference Zakharov1973; Hanson Reference Hanson2015). Specifically, the virtual-casing principle states that if $\varGamma$ is the flux surface bounding the plasma, then $\boldsymbol {B}_V^{{\rm pol}}$ can be written in terms of a field generated by the toroidal surface current $\boldsymbol {J}^{{\rm tor}}_S$ such that $\mu _0 \boldsymbol {J}^{{\rm tor}}_S = -\boldsymbol {n} \times \boldsymbol {B}^{{\rm pol}}$, where $\boldsymbol {n}$ is the outward unit normal vector to $\varGamma$. This relation is given according to Hanson (Reference Hanson2015):

(2.2)\begin{equation} \boldsymbol{B}^{{\rm pol}}_V(\boldsymbol{x}) =\frac 1{4{\rm \pi}} \iint_{\varGamma} \left[\frac{(\boldsymbol{n}(\boldsymbol{y}) \times \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y})) \times (\boldsymbol{x} - \boldsymbol{y})}{\| \boldsymbol{x} - \boldsymbol{y}\|^3}\right] \mathrm{d} \varGamma(\boldsymbol{y}) +\begin{cases} \boldsymbol{B}^{{\rm pol}}(\boldsymbol{x}) & \boldsymbol{x} \in \varOmega, \\ \boldsymbol{B}^{{\rm pol}}(\boldsymbol{x})/2 & \boldsymbol{x} \in \varGamma, \\ \boldsymbol{0} & \boldsymbol{x} \notin \bar{\varOmega}. \end{cases} \end{equation}

In this letter, we are interested in the most common application of the virtual-casing principle, in which we want to compute the magnetic field due to the plasma currents on the plasma boundary, i.e. for $\boldsymbol {x}\in \varGamma$. That motivates why we focus on the poloidal magnetic field. Indeed, for axisymmetric plasma boundaries and magnetic fields, the toroidal magnetic field on the plasma boundary is entirely generated by external coils. To see this, let $C_{r,z}$ be any circle of radius $r$ in a plane perpendicular to the $z$-axis, at altitude $z$, that is tangent to the axisymmetric plasma boundary. Applying Stokes’ theorem along this circle, we can write

(2.3)\begin{equation} \oint_{C_{r,z}}B_{\phi}(r,z)r\,\mathrm{d}\phi = 2{\rm \pi} r B_{\phi}(r,z) = \mu_{0}\iint_{D_{r,z}} \boldsymbol{J} \boldsymbol{\cdot} \mathrm{d}\boldsymbol{S} = \mu_{0} I_{\mathrm{coils}}, \end{equation}

where $D_{r,z}$ is the flat disc perpendicular to the $z$-axis bounded by the circle $C_{r,z}$, and $I_{\mathrm {coils}}$ is the total coil-induced current. In the last equality, we have used the fact that the plasma currents are divergence free. We obtain the well-known expression for the vacuum toroidal field, $B_{\phi }(r)=\mu _{0}I_{\mathrm {coils}}/(2{\rm \pi} r)$, which holds all the way to the plasma boundary and confirms that for axisymmetric plasmas, plasma currents do not contribute to the toroidal magnetic on the plasma boundary.

2.3. Reduction of the virtual-casing principle to line integrals

Since we are considering an axisymmetric surface $\varGamma$, the surface integral in (2.2) can be expressed in terms of a line integral by integrating over the toroidal angle analytically.

The poloidal magnetic field ${\boldsymbol {B}}^{{\rm pol}}$ at any point $\boldsymbol {y}(\phi, t) \in \varGamma$ can be expressed in terms of its poloidal flux function $\psi (r,z)$ and the parametrization $(\phi,t) \mapsto \boldsymbol {y}(\phi, t)$ by (Freidberg Reference Freidberg2014)

(2.4)\begin{align} \boldsymbol{B}^{{\rm pol}}(\boldsymbol{y}(\phi,t)) & = \boldsymbol\nabla \psi(r(t),z(t)) \times \boldsymbol\nabla \phi \nonumber\\ & =-\frac 1{r(t)} \frac{\partial \psi}{\partial z}(r(t), z(t)) \boldsymbol{e}_r(\phi) + \frac 1{r(t)} \frac{\partial \psi}{\partial r}(r(t), z(t)) \boldsymbol{e}_z. \end{align}

Inserting this expression in (2.2) and integrating it analytically with respect to the toroidal angle, we may evaluate $\boldsymbol {B}_V^{{\rm pol}}$ directly through the remaining univariate integral. Presently, we consider evaluating on the plasma boundary $\varGamma$. Without loss of generality, we assume that the evaluation point lies on the cross-section of the plasma boundary with azimuthal angle $\phi =0$, which has cylindrical coordinate unit vectors that align with the rectangular unit vectors $\{\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z\}$. We denote such an evaluation point as $(R,0,Z)$ in cylindrical coordinates.

Upon integrating the toroidal angle, the univariate expression for $\boldsymbol {B}_V^{{\rm pol}}$ on the plasma boundary is

(2.5)\begin{equation} \boldsymbol{B}_V^{{\rm pol}}(R,Z)= \frac{1}{4{\rm \pi}} \int_{0}^{L} \boldsymbol{f}(t)\,\mathrm{d} t + \frac{1}{2} \boldsymbol{B}^{{\rm pol}}(R,Z),\quad (R,0,Z) \in \varGamma, \end{equation}

where (Toler et al. Reference Toler, Cerfon and Malhotra2023)

(2.6)\begin{align} \boldsymbol{f}(t) & = \frac {2} {r(t) \sqrt{\alpha + \beta}} \left[ \frac{\partial \psi}{\partial z} r'(t) - \frac{\partial \psi}{\partial r} z'(t) \right] \left\{\frac {Z - z(t)}{R}\left(-K(k^2)+ \frac{\alpha}{\alpha - \beta} E(k^2)\right) \right.\boldsymbol{e}_x \nonumber\\ & \quad \left.+ \left(K(k^2)+ \frac{r(t)^2 - R^2 - (Z-z(t))^2}{\alpha - \beta} E(k^2)\right)\boldsymbol{e}_z\right\}\!. \end{align}

Here, the functions $K$ and $E$ are the complete elliptic integrals of the first and second kind, primes denote derivatives with respect to the parametrization variable $t$, and we have introduced the quantities

(2.7)\begin{equation} \begin{cases} \alpha = \alpha(t; R,Z) = R^2 + r(t)^2 + (Z-z(t))^2, \\ \beta = \beta(t; R,Z) = 2 R r(t), \\ k^2 = k(t; R,Z)^2 = \dfrac{2\beta}{\alpha+\beta}. \end{cases} \end{equation}

Since we are evaluating on the plasma boundary, the evaluation point can be identified by its generating curve parameter $t$, say $t_0$, so $r(t_0)=R$ and $z(t_0)=Z$. There are two sources of singularity in $\boldsymbol {f}(t)$ near $t_0$. The mapping $t \mapsto K(k(t)^2)$ is $O(\log |t-t_0|)$ near $t_0$ (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014; Toler et al. Reference Toler, Cerfon and Malhotra2023), and the errors from Taylor series for $r(t)$ and $z(t)$ yield $1/(\alpha -\beta ) = O(1/(t-t_0)^2)$.

3. High-order principal value quadrature for the axisymmetric virtual-casing principle

Previously (Toler et al. Reference Toler, Cerfon and Malhotra2023), we showed that the integral in (2.5) is not convergent. That means that this integral must be interpreted in the Cauchy principal value sense for all components of $\boldsymbol {B}_V^{{\rm pol}}(R,Z)$ in (2.5) to be well defined. Since the observation point $(R,Z)$ lies on the plasma boundary, we may identify this point with a value of the boundary parameter $t$, say $t_0$.

To evaluate all the components of $\boldsymbol {B}_V^{{\rm pol}}(R,Z)$ via (2.5), we must evaluate the integral in this expression with the following substitution:

(3.1)\begin{equation} \int_{0}^{L} \boldsymbol{f}(t)\,\mathrm{d} t\quad \xrightarrow{\qquad} \quad \lim_{\epsilon \to 0^+} \left(\int_{0}^{t_0-\epsilon} \boldsymbol{f}(t)\,\mathrm{d} t + \int_{t_0+\epsilon}^{L} \boldsymbol{f}(t)\,\mathrm{d} t \right)\!, \end{equation}

where $\boldsymbol {f}$ is given in (2.6). (If $t_0=0$, the integration interval within the limit should be $[\epsilon, L-\epsilon ]$.) Let us discuss two quadrature rules to do so: the alternating trapezoidal rule and the periodic Kapur–Rokhlin rule (Kapur & Rokhlin Reference Kapur and Rokhlin1997; Toler et al. Reference Toler, Cerfon and Malhotra2023).

We begin by considering the alternating trapezoidal rule. We recall that the standard $N$-point periodic trapezoidal rule for the integral $\int _{0}^{L} g(t)\,\mathrm{d} t$ for an $L$-periodic integrand $g$ is defined by equispaced quadrature nodes with uniform weights $h=L/N$. The alternating trapezoidal rule for the case when $g$ has a singularity at $t_0$ is a specific case of the trapezoidal rule, where the quadrature nodes $\{s_i\}_{i=1}^{N}$ are placed to straddle $t_0$. That is, the nodes are $s_i = t_0 + (i - \frac {1}{2}) h$ for $i=1, \ldots, N$. The alternating trapezoidal rule for some types of Cauchy principal value integrals is supported by mathematical theory. In particular, it is known to converge spectrally accurately for principal value integrals of the form $\int _{-1}^1 (q(t)/t)\,\mathrm{d} t$ when $q$ is an analytic function (Kress & Martensen Reference Kress and Martensen1970; Trefethen & Weideman Reference Trefethen and Weideman2014). Additionally, if $q \in C^{P}[0,L]$ but is not analytic, then the quadrature converges with order $P$ (Sidi & Israeli Reference Sidi and Israeli1988). By staggering quadrature nodes with equal weights on each side of the singularity, the alternating trapezoidal rule captures the cancellation that is necessary for a well-defined Cauchy principal value. For our application, we note that the asymptotic behaviour of the integrand is not $q(t)/t$, but $q_1(t)/t + q_2(t) \log |t| + q_3(t)$ for bounded functions $q_1, q_2, q_3$. As a result of the logarithmically singular term, we do not expect high-order convergence of the quadrature error. Rather, we expect low-order convergence, as we will verify numerically in the next section.

In contrast, the Kapur–Rokhlin rule is precisely designed to integrate logarithmic singularities with high-order accuracy (Toler et al. Reference Toler, Cerfon and Malhotra2023). We recall that the $(N-1)$-point rule for the periodic setting of interest in our context can be written as

(3.2)\begin{equation} \int_{t_0-L/2}^{t_0+L/2} g(t)\,\mathrm{d} t = \sum_{\substack{j=-(N-1)/2 \\ j \ne 0}}^{(N-1)/2} w_j g_j + O(h^n), \end{equation}

where $g$ is logarithmically singular at $t_0$, with

(3.3)\begin{equation} w_j =\begin{cases} (1+\gamma_j + \gamma_{-j}) h, & 1 \le |j| \le n, \\ h, & \text{otherwise}, \end{cases} \end{equation}

and $g_j \equiv g(t_0 + j h)$. The convergence order $n$ may be chosen by the user and the weights $\{ \gamma _j \}$ may be read from a table of values (Kapur & Rokhlin Reference Kapur and Rokhlin1997). The spacing $h=L/N$ between quadrature points is uniform (except for the doubled space between nodes $-1$ and $1$), and we are using notation so that the meaning of $h$ is the same for both the alternating trapezoidal rule and the Kapur–Rokhlin rule. In (3.2), we assume that $N$ is odd to simplify notation. If $N$ is even, the points may be indexed by $j=\pm 1, \ldots, \pm (N/2-1), N/2$.

Like the alternating trapezoidal rule, the periodic Kapur–Rokhlin rule staggers nodes around the singularity with equal weights, so it also captures the principal value cancellation when integrating $\boldsymbol {f}$. Indeed, the weights in (3.3) satisfy $w_{j} = w_{-j}$. However, the quadrature error converges much faster than for the alternating trapezoidal rule, since the quadrature rule was also designed for the logarithmic singularity of $\boldsymbol {f}$. We next discuss practical performance of these numerical methods in the context of a virtual-casing principle calculation.

4. Numerical results

To compare and verify the numerical methods for calculating the principal value of the integral in (2.5), we use the same reference data as Toler et al. (Reference Toler, Cerfon and Malhotra2023) on the same geometry. Specifically, the axisymmetric plasma boundary is given by the level set ${\{ \psi =0 \}}$ of the poloidal flux function given by

(4.1)\begin{equation} \psi(r,z) = \frac{\kappa F_B}{2 R_0^3 q_0} \left[\frac 14 (r^2 - R_0^2)^2 + \frac 1{\kappa^2} r^2 z^2 - a^2 R_0^2 \right]\!, \end{equation}

which solves the Grad–Shafranov equation with the Solov'ev profiles $\mu _{0}p(\psi )=-[F_{B}(\kappa + 1/\kappa )/( R_0^{3}q_0)]\psi$ and $F(\psi )=F_B$, where $p(\psi )$ is the plasma pressure profile, and $F(\psi )=rB_{\phi }$, with $B_\phi$ the toroidal magnetic field (Lütjens, Bondeson & Sauter Reference Lütjens, Bondeson and Sauter1996; Lee & Cerfon Reference Lee and Cerfon2015). The parameters $R_0$ and $q_0$ may be interpreted as the major radius and safety factor at the magnetic axis, and $\kappa$ and $a$ as the elongation and minor radius of the plasma boundary, respectively. Throughout, we use the fusion relevant values $F_B = R_0 = q_0 = 1$, and $\kappa = 1.7$ and $a = 1/3$. The level set $\{ \psi = 0 \}$ may be parametrized by the functions (Lütjens et al. Reference Lütjens, Bondeson and Sauter1996; Lee & Cerfon Reference Lee and Cerfon2015)

(4.2a,b)\begin{equation} (r(t))^2 = R_0^2 + 2 a R_0 \cos t\quad \text{and}\quad z(t) = \kappa a \frac{R_0}{r(t)} \sin t \end{equation}

for $t \in [0, L]$ with $L=2{\rm \pi}$.

We take as ground truth the plasma-induced magnetic field from a high-resolution computation by the approach of Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019). The method used in that implementation is different from the one-dimensional integral approach presented here in several ways, making it appropriate for our verification. Specifically, the code presented by Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019) views the plasma equilibrium as a fully three-dimensional equilibrium and does not assume axisymmetry. Furthermore, Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019) obtain $\boldsymbol {B}^{{\rm pol}}_V$ on the plasma boundary directly by evaluating the virtual-casing principle as a surface integral, as opposed to first reducing to a univariate integral by axisymmetry. Finally, the Cauchy principal value of the surface integral is numerically evaluated via a partition of unity scheme to handle the singularity of the integrand. The reference solution was computed on a fine grid with 1,580 toroidal angle values and 1,200 poloidal points on each cross-section. In a self-convergence comparison with a solution from finer resolution, these data were found to be accurate to 11 digits of precision. Figure 1 illustrates the geometry and reference solution for this problem. This approach by Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019) is more robust by treating fully three-dimensional volumes with two-dimensional surfaces, but axisymmetric specializations like those we present here require discretization of one less dimension and result in faster runtimes.

Figure 1. Reference solution for the Solov'ev profiles of § 4. Colours represent $\| \boldsymbol {B}\|$ on the plasma surface.

With these reference values, we shall assess the numerical accuracy of the alternating trapezoidal rule and the Kapur–Rokhlin rule, which we discussed in the previous section. For each quadrature scheme, we compute both the radial component $\boldsymbol {B}_{V,R}^{{\rm pol}}$ of $\boldsymbol {B}_V^{{\rm pol}}$ and the vertical component $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ at 1,200 points on the plasma boundary which are equispaced in the generating curve parameter $t$. We measure the relative quadrature error for each component as the maximum of the pointwise absolute difference between the computed field component and the reference values. We normalize both errors by the maximum absolute value of the 2,400 reference values for $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$.

Figure 2 illustrates that the alternating trapezoidal rule indeed converges to the same values found with the scheme presented by Malhotra et al. (Reference Malhotra, Cerfon, O'Neil and Toler2019), but does so with a low order of convergence. The convergence rates are $O(h^3)$ for the radial component $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $O(h)$ for the vertical component $\boldsymbol {B}_{V,Z}^{{\rm pol}}$. The radial component converges quicker because its integrand has more regularity due to the following asymptotic behaviour. Recall from § 2.3 that $t \mapsto K(k(t)^2)$ is $O(\log |t-t_0|)$ near $t_0$. Hence, in the integrand $\boldsymbol {f}$ of (2.6), the product $(Z-z(t))K(k^2)$ has only a removable discontinuity at $t_0$ and not a logarithmic singularity. The only singularity in the radial component is the principal value singularity of $(Z-z(t))/(\alpha -\beta )$. No such product dampens the singularity of $K$ in the vertical component. Nevertheless, the convergence rates for both components of $\boldsymbol {B}_{V}^{{\rm pol}}$ are better than those guaranteed by the quadrature theory we discussed in § 3.

Figure 2. Maximum pointwise error using the alternating trapezoidal rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ via principal value.

We next consider the periodic Kapur–Rokhlin rule for the principal value quadrature. We view each component of the integrand $\boldsymbol {f}$ as the sum of a logarithmically singular term and a remainder term:

(4.3)\begin{equation} \begin{cases} \boldsymbol{e}_x \boldsymbol{\cdot} \boldsymbol{f}(t) = p_R(t) \log|t-t_0| + q_R(t), \\ \boldsymbol{e}_z \boldsymbol{\cdot} \boldsymbol{f}(t) = p_Z(t) \log|t-t_0| + q_Z(t). \end{cases} \end{equation}

In particular, we have

(4.4)\begin{equation} p_R(t) \log|t-t_0| = \frac {2} {r(t) \sqrt{\alpha + \beta}} \left[ \frac{\partial \psi}{\partial z} r'(t) - \frac{\partial \psi}{\partial r}z'(t) \right] \frac {Z - z(t)}{R} (-K(k^2)) \end{equation}

and

(4.5)\begin{equation} p_Z(t) \log|t-t_0| = \frac {2} {r(t) \sqrt{\alpha + \beta}} \left[\frac{\partial \psi}{\partial z} r'(t) - \frac{\partial \psi}{\partial r} z'(t) \right]K(k^2), \end{equation}

and we set $q_R$ and $q_Z$ to be the remaining corresponding terms in (2.6). In figures 3, 4 and 5, we show the results of applying the second-, sixth- and tenth-order Kapur–Rokhlin rules, respectively. These are the theoretical convergence rates guaranteed for smooth coefficient functions $p_R, p_Z, q_R, q_Z$. In this example, however, we see convergence rates that deviate slightly from these values. The sixth-order rule converges at the expected rate of $O(h^6)$. However, the second-order rule converges faster than predicted at rates of approximately $O(h^{2.71})$ and $O(h^{2.52})$ for the radial and vertical components, respectively. Additionally, the tenth-order rule converges at approximate rates of $O(h^{8.74})$ and $O(h^{8.73})$. We compute these approximate rates by performing an ordinary least squares line of best fit to the data measurements in the asymptotic regime of convergence.

Figure 3. Maximum pointwise error using the second-order Kapur–Rokhlin rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ via principal value.

Figure 4. Maximum pointwise error using the sixth-order Kapur–Rokhlin rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ via principal value. We compare the Kapur–Rokhlin (KR) error with the alternating trapezoidal (AT) error of figure 2.

Figure 5. Maximum pointwise error using the tenth-order Kapur–Rokhlin rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$ via principal value. We compare the Kapur–Rokhlin (KR) error with the alternating trapezoidal (AT) error of figure 2.

Notably, the tenth-order scheme converges slightly slower than predicted. We offer two potential explanations for this phenomenon. First, the smoothness assumptions of Kapur–Rokhlin are violated, as we explain presently. Similar to the alternating trapezoidal rule, the Kapur–Rokhlin rule requires that all of $\{p_R, p_Z, q_R, q_Z\}$ are sufficiently smooth to guarantee the theoretical convergence rate (Kapur & Rokhlin Reference Kapur and Rokhlin1997). However, both $q_R$ and $q_Z$ include the elliptic integral $E$. The mapping $t \mapsto E(k(t)^2)$ is continuous, but its first derivative is not. This is readily seen by the identity of Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2014) that

(4.6)\begin{equation} \frac{\partial E(k^2)}{\partial k^2} = \frac{E(k^2) - K(k^2)}{2 k^2}. \end{equation}

The first derivative of $t \mapsto E(k(t)^2)$, and hence the first derivatives of $q_R$ and $q_Z$, are logarithmically singular at $t_0$ due to the appearance of $K$. Second, very high orders of the Kapur–Rokhlin rule have weight corrections $\{\gamma _j\}$ which are large in magnitude and are sign-indefinite. This introduces numerical instability as a price to pay for the higher convergence order (Martinsson Reference Martinsson2019).

We hypothesize that the instability from the weight corrections also explains the observation that the quadrature error saturates at approximately nine digits of accuracy in figures 4 and 5. It is numerically common for algorithmic instabilities to raise the error floor from machine precision. The instability is known to worsen as the theoretical convergence order increases, so we see this phenomenon earlier and more markedly for the tenth-order rule compared with the sixth-order rule. In contemporary plasma physics applications relevant to this work, however, accuracy at the level of eight or nine digits is typically not a limiting factor in overall calculations, so we do not explore the error saturation further.

Though the second- and sixth-order rules have agreed with or exceeded the theoretical baselines in our experiment, we emphasize that we cannot guarantee such performance in general. The behaviour of all of the Kapur–Rokhlin schemes in this setting are less predictable, but nevertheless promising in practice.

When we directly compare the performance of the Kapur–Rokhlin scheme and the alternating trapezoidal scheme in figures 4 and 5, we see that both methods are competitive for computing $\boldsymbol {B}^{{\rm pol}}_{V,R}$, whose virtual-casing integrand contains a less severe singularity. The Kapur–Rokhlin rule is asymptotically faster to converge, but the alternating trapezoidal rule has better accuracy when the number of quadrature points is small. However, the Kapur–Rokhlin rule more clearly outperforms the alternating trapezoidal rule for computing $\boldsymbol {B}^{{\rm pol}}_{V,Z}$, whose virtual-casing integrand contains the more severe singularity. In this case, the alternating trapezoidal rule only has higher accuracy when the number of quadrature points is very limited: less than approximately 50.

Finally, we consider a third quadrature scheme which directly combines the alternating trapezoidal rule and the Kapur–Rokhlin rule. Returning to the decomposition (4.3), we integrate the logarithmically singular terms with $p_R$ and $p_Z$ of (4.4) and (4.5) via the periodic Kapur–Rokhlin rule, and we integrate the remaining terms $q_R$ and $q_Z$ via the alternating trapezoidal rule. One might imagine that this scheme would outperform both the pure Kapur–Rokhlin rule and the pure alternating trapezoidal rule, since this combined scheme uses each constituent quadrature on the term for which it is specialized. However, figure 6 reveals that the lack of smoothness in the integrand again stalls convergence. For few quadrature nodes, the Kapur–Rokhlin corrections considerably improve convergence, but for more than approximately 100 nodes, the slower convergence of the alternating trapezoidal rule dominates. We observe that the limiting asymptotic convergence rate appears to be third-order convergence in this experiment for both the radial and vertical components of $\boldsymbol {B}_{V}^{{\rm pol}}$. This improves the asymptotic rate from the pure alternating trapezoidal rule for $\boldsymbol {B}_{V,Z}^{{\rm pol}}$, indicating that the alternating trapezoidal rule originally achieved only $O(h)$ convergence because it ignored the logarithmic singularity.

Figure 6. Maximum pointwise error using a combination of the tenth-order Kapur–Rokhlin rule and the alternating trapezoidal rule to compute $\boldsymbol {B}_{V,R}^{{\rm pol}}$ and $\boldsymbol {B}_{V,Z}^{{\rm pol}}$. We compare the error with the alternating trapezoidal (AT) error of figure 2.

We additionally note that the asymptotic rates for the second- and tenth-order Kapur–Rokhlin rules in figures 2 and 4 experimentally appear to agree with rates of $O(h^3 \log h)$ and $O(h^9 \log h)$, respectively, though this is our hypothesis and we are not aware of quadrature theory that would imply the presence of a $\log h$ factor in the convergence rate in this context. However, convergence orders with $\log h$ factors are known to appear when ill-suited quadrature rules, like the alternating trapezoidal rule, are used to integrate a function with a purely logarithmic singularity (Martinsson Reference Martinsson2019). We display this feature for the tenth-order rule in figure 6.

The pure Kapur–Rokhlin method remains the most efficient quadrature scheme discussed for this principle value application. It empirically displays high-order convergence and accuracy close to the theorized rates. The tenth-order rule achieves a maximal accuracy of approximately nine digits with approximately 400 quadrature nodes.