Formal Solutions for Polarized Radiative Transfer. II. High-order Methods

The best known high-order scheme of this class is probably the classical Runge–Kutta 4 method (RK4) and its application to the formal solution for polarized light was already proposed by Landi Degl'Innocenti (1976). The method is described by

$$\begin{eqnarray}&&{{\boldsymbol{I}}}_{k+1}={{\boldsymbol{I}}}_{k}+\displaystyle \frac{1}{6}[{{\boldsymbol{k}}}_{1}+2{{\boldsymbol{k}}}_{2}+2{{\boldsymbol{k}}}_{3}+{{\boldsymbol{k}}}_{4}],\end{eqnarray} \tag{ 3 }$$

where the four stage values are given by

$$\begin{eqnarray*}{{\boldsymbol{k}}}_{1} & = & {\rm{\Delta }}{s}_{k}{\boldsymbol{F}}({s}_{k},{{\boldsymbol{I}}}_{k}),\\ {{\boldsymbol{k}}}_{2} & = & {\rm{\Delta }}{s}_{k}{\boldsymbol{F}}({s}_{k}+{\rm{\Delta }}{s}_{k}/2,{{\boldsymbol{I}}}_{k}+{{\boldsymbol{k}}}_{1}/2),\\ {{\boldsymbol{k}}}_{3} & = & {\rm{\Delta }}{s}_{k}{\boldsymbol{F}}({s}_{k}+{\rm{\Delta }}{s}_{k}/2,{{\boldsymbol{I}}}_{k}+{{\boldsymbol{k}}}_{2}/2),\\ {{\boldsymbol{k}}}_{4} & = & {\rm{\Delta }}{s}_{k}{\boldsymbol{F}}({s}_{k+1},{{\boldsymbol{I}}}_{k}+{{\boldsymbol{k}}}_{3}).\end{eqnarray*}$$

The right side of Equation (3) does not contain the term ${{\boldsymbol{I}}}_{k+1}$ and RK4 is therefore classified as an explicit method.

RK4 is well known for its fourth-order accuracy, as indicated in Table 1, and its convergence analysis is easily found in the literature (e.g., Deuflhard & Bornemann 2002; Frank & Leimkuhler 2012). However, the quantities ${\bf{K}}$ and ${\boldsymbol{\epsilon }}$ at the intermediate point ${s}_{k}+{\rm{\Delta }}{s}_{k}/2$ must be properly provided through interpolation. In order to maintain fourth-order accuracy, one needs an interpolation of degree $q\geqslant 3$, i.e., at least cubic. A lower-order interpolation would decrease the order of accuracy: a parabolic interpolation results in a third-order accurate method and a linear interpolation provides a second-order accurate method.

An alternative strategy for providing accurate absorption and emission quantities at the intermediate point is based on high-order interpolations of the atmospheric model parameters, such as the temperature, the microscopic and macroscopic velocities, and the strength and orientation of the magnetic field. The quantities ${\bf{K}}$ and ${\boldsymbol{\epsilon }}$ at ${s}_{k}+{\rm{\Delta }}{s}_{k}/2$ are then evaluated through the interpolated thermodynamic parameters.

As explained in Paper I, the stability of a numerical method is often deduced through the simple autonomous scalar initial value problem (IVP) given by

$$\begin{eqnarray}y^{\prime} (t) & = & \lambda y(t),\\ y(0) & = & {y}_{0},\end{eqnarray} \tag{ 4 }$$

with $\lambda \in {\mathbb{C}}$. The solution $y(t)={y}_{0}{e}^{\lambda t}$ converges to zero as $t\to \infty$ for $\mathrm{Re}(\lambda )\lt 0$. Defining $z=\lambda {\rm{\Delta }}t$, where ${\rm{\Delta }}t$ denotes the cell width, the RK4 method applied to the IVP (4) is recast into the form

$$\begin{eqnarray*}&&{y}_{k+1}={\phi }_{\mathrm{RK}4}(z){y}_{k},\end{eqnarray*}$$

where the stability function ${\phi }_{\mathrm{RK}4}$ reads

$$\begin{eqnarray*}&&{\phi }_{\mathrm{RK}4}(z)=1+z+\displaystyle \frac{{z}^{2}}{2}+\displaystyle \frac{{z}^{3}}{6}+\displaystyle \frac{{z}^{4}}{24}.\end{eqnarray*}$$

Stability is guaranteed by the condition $\parallel {\phi }_{\mathrm{RK}4}(z)\parallel \lt 1$. The stability region of RK4 presented in Figure 1(a) is clearly bounded. This indicates that the method could suffer from magnification of numerical errors for large z values. This problem is relevant for optically thick cells. In fact, the eigenvalues of the propagation operator $-{\bf{K}}$, which have real parts that are always negative, increase with the total absorption coefficient ${\eta }_{I}$ (Landi Degl'Innocenti & Landolfi 2004).

Zoom In Zoom Out Reset image size

In order to face this problem, a hybrid technique can be used. This strategy applies an A-stable method, e.g., the trapezoidal method, for optical thick cells, making use of RK4 elsewhere. One correctly argues that this hybrid technique possesses the lowest order of accuracy among the used methods, i.e., second-order accuracy if the trapezoidal method is chosen. However, the strong attenuation induced by optically thick cells could reduce the error propagation. In fact, this hybrid method maintains fourth-order convergence, as clearly shown in Figure 2 where it is labeled Runge–Kutta 4.

Zoom In Zoom Out Reset image size

Note that the assumption of a constant eigenvalue λ in Equation (4) is a limitation of this simplified stability analysis. In fact, variations of ${\bf{K}}$ along the integration path usually affect the stability region of the numerical method (see Paper I). Using an optical depth scale usually supports the assumption of a constant eigenvalue λ in Equation (4).

As mentioned above, the RK4 method is explicit: it avoids the additional solution of the 4 × 4 implicit linear system, which is required by implicit methods. This fact could significantly reduce the total amount of computational time required. When applied to the polarized formal solution, RK4 has often been classified as computationally costly, because of the very small step size required (Rees et al. 1989). In light of the previous stability analysis, one is led to believe that the requirement of small numerical cells is mainly due to the bounded stability region of the RK4 method and that the use of a hybrid strategy should overcome this problem.

The p-step Adams–Moulton's method approximates the integrand ${\boldsymbol{F}}$ in Equation (2) by the p-order Lagrange polynomial

$$\begin{eqnarray*}&&{\boldsymbol{F}}(s,{\boldsymbol{I}}(s))\approx \displaystyle \sum _{i=k-p+1}^{k+1}{{\boldsymbol{F}}}_{i}{{\ell }}_{i}(s)\ \mathrm{for}\ s\in [{s}_{k},{s}_{k+1}],\end{eqnarray*}$$

which matches the numerical values ${{\boldsymbol{F}}}_{i}=-{{\bf{K}}}_{i}{{\boldsymbol{I}}}_{i}+{{\boldsymbol{\epsilon }}}_{i}$ at positions s_i. The Lagrange basis polynomials ℓ_i given by

$$\begin{eqnarray*}&&{{\ell }}_{i}(s)=\displaystyle \prod _{\mathop{k-p+1\leqslant m\leqslant k+1}\limits_{m\ne i}}\displaystyle \frac{s-{s}_{m}}{{s}_{i}-{s}_{m}},\end{eqnarray*}$$

satisfy the relation ${{\ell }}_{i}({s}_{j})={\delta }_{{ij}}$, where the Kronecker delta ${\delta }_{{ij}}$ is defined by

$$\begin{eqnarray*}{\delta }_{{ij}}=\left\{\begin{array}{ll}1 & \ \mathrm{if}\ i=j,\\ 0 & \ \mathrm{if}\ i\ne j.\end{array}\right.\end{eqnarray*}$$

The integral in Equation (2) can then be solved by parts, yielding, after some algebra, a linear system of the form of Equation (5). The presence of the term ${{\boldsymbol{F}}}_{k+1}$ in the Lagrange interpolation provides ${\beta }_{k+1}\ne 0$, indicating the method to be implicit. By contrast, Adams–Bashford methods define the Lagrange interpolation through ${{\boldsymbol{F}}}_{k-p},\ldots ,\,{{\boldsymbol{F}}}_{k}$, providing a linear system with ${\beta }_{k+1}=0$, which is therefore explicit.

At this point, one has to note the critical difference between the Adams–Moulton family and the DELO family discussed in Paper I. In the former, the Lagrange polynomial approximates the integrand ${\boldsymbol{F}}$, while in the latter the interpolation is applied to the effective source function. Nonetheless, the two different strategies share similar convergence properties, because the local truncation error depends on the interpolation degree in both cases.

Although outside the assumption $p\geqslant 1$, it is habit to include the p = 0 case in the Adams–Moulton family. In this instance, the integrand in Equation (2) is approximated as ${\boldsymbol{F}}\approx {{\boldsymbol{F}}}_{k+1}$. Doing so, one obtains the common first-order accurate backward (or implicit) Euler method (e.g., Deuflhard & Bornemann 2002).

Now, if the first-order Lagrange interpolant is used to approximate the integrand ${\boldsymbol{F}}$ in Equation (2), one obtains the famous second-order accurate trapezoidal method as assessed in Paper I. Note that both the backward Euler method and the trapezoidal method are one-step methods, which also belong to the Runge–Kutta class (see Section 2).

If a parabolic Lagrange interpolation is performed through ${{\boldsymbol{F}}}_{k-1}$, ${{\boldsymbol{F}}}_{k}$, and ${{\boldsymbol{F}}}_{k+1}$, one obtains the implicit linear system given by

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{k+1}{{\boldsymbol{I}}}_{k+1}={{\boldsymbol{\Phi }}}_{k}{{\boldsymbol{I}}}_{k}+{{\boldsymbol{\Phi }}}_{k-1}{{\boldsymbol{I}}}_{k-1}+{{\boldsymbol{\Psi }}}_{k+1}+{{\boldsymbol{\Psi }}}_{k}+{{\boldsymbol{\Psi }}}_{k-1},\end{eqnarray} \tag{ 6 }$$

and the coefficients ${{\boldsymbol{\Phi }}}_{k-1}$, ${{\boldsymbol{\Phi }}}_{k}$, ${{\boldsymbol{\Phi }}}_{k+1}$, ${{\boldsymbol{\Psi }}}_{k-1}$, ${{\boldsymbol{\Psi }}}_{k}$, and ${{\boldsymbol{\Psi }}}_{k+1}$ are provided in Appendix A. The two-step numerical scheme described by Equation (6) is called the Adams–Moulton 3 method.

A cubic Lagrange interpolation through ${{\boldsymbol{F}}}_{k-2}$, ${{\boldsymbol{F}}}_{k-1}$, ${{\boldsymbol{F}}}_{k}$, and ${{\boldsymbol{F}}}_{k+1}$, provides the following implicit linear system:

$$\begin{eqnarray}{{\boldsymbol{\Phi }}}_{k+1}{{\boldsymbol{I}}}_{k+1} & = & {{\boldsymbol{\Phi }}}_{k}{{\boldsymbol{I}}}_{k}+{{\boldsymbol{\Phi }}}_{k-1}{{\boldsymbol{I}}}_{k-1}+{{\boldsymbol{\Phi }}}_{k-2}{{\boldsymbol{I}}}_{k-2}\\ & & +{{\boldsymbol{\Psi }}}_{k+1}+{{\boldsymbol{\Psi }}}_{k}+{{\boldsymbol{\Psi }}}_{k-1}+{{\boldsymbol{\Psi }}}_{k-2},\end{eqnarray} \tag{ 7 }$$

and the coefficients ${{\boldsymbol{\Phi }}}_{k-2}$, ${{\boldsymbol{\Phi }}}_{k-1}$, ${{\boldsymbol{\Phi }}}_{k}$, ${{\boldsymbol{\Phi }}}_{k+1}$, ${{\boldsymbol{\Psi }}}_{k-2}$, ${{\boldsymbol{\Psi }}}_{k-1}$, ${{\boldsymbol{\Psi }}}_{k}$, and ${{\boldsymbol{\Psi }}}_{k+1}$ are provided in Appendix A. The three-step numerical scheme described by Equation (7) is called the Adams–Moulton 4 method.

This family of formal solvers can be further expanded by just increasing the interpolation degree of the integrand ${\boldsymbol{F}}$. However, the complexity of the numerical methods would increase and the expressions for the ${\boldsymbol{\Phi }}$ and ${\boldsymbol{\Psi }}$ coefficients would become gradually more cumbersome.

The local truncation error of the Adams–Moulton methods is due to the fact that the integrand ${\boldsymbol{F}}$ in Equation (2) is approximated by a polynomial. A Lagrange polynomial of degree p is known to be $(p+1)\mathrm{th}$-order accurate. The resulting local truncation error satisfies

$$\begin{eqnarray*}&&{L}^{{\rm{A}}-{\rm{M}}}[p]\approx O({\rm{\Delta }}{s}^{p+2}),\end{eqnarray*}$$

indicating a p-step Adams–Moulton method as $(p+1)$-order accurate (see Paper I). Accordingly, the Adams–Moulton 3 method described by Equation (6) is third-order accurate, whereas the Adams–Moulton 4 method described by Equation (7) is fourth-order accurate, as summarized in Table 1, and a numerical confirmation is given by Figure 2.

The simple stability analysis performed above for one-step methods cannot be applied to this class, because of the multiterm contribution. Therefore, linear multistep methods require a more complex derivation of the stability region (Frank & Leimkuhler 2012).

The numerical scheme given by Equation (5) applied to the IVP (4) gives

$$\begin{eqnarray*}&&\displaystyle \sum _{i=k-p+1}^{k+1}{\alpha }_{i}{y}_{i}=\lambda {\rm{\Delta }}t\displaystyle \sum _{i=k-p+1}^{k+1}{\beta }_{i}{y}_{i},\end{eqnarray*}$$

and, defining $z=\lambda {\rm{\Delta }}t$, one gets

$$\begin{eqnarray*}&&\displaystyle \sum _{i=k-p+1}^{k+1}({\alpha }_{i}-z\,{\beta }_{i}){y}_{i}=0.\end{eqnarray*}$$

For any z, this is a linear difference equation with the characteristic polynomial

$$\begin{eqnarray}&&\displaystyle \sum _{i=k-p+1}^{k+1}({\alpha }_{i}-z\,{\beta }_{i}){\zeta }^{i}=0=\rho (\zeta )-z\,\sigma (\zeta ).\end{eqnarray} \tag{ 8 }$$

The stability region of a linear multistep method is the set of complex values z for which all roots ζ of the polynomial Equation (8) lie on the unit disk, i.e., $| \zeta | \leqslant 1$, and those with a modulus of one are simple. On the boundary of the stability region, precisely one root has a modulus of one. Therefore, an explicit representation for the boundary of the stability region is given by

$$\begin{eqnarray*}&&\delta S=\left\{z=\displaystyle \frac{\rho ({e}^{i\theta })}{\sigma ({e}^{i\theta })}\right\},\ \mathrm{for}\ \theta \in [-\pi ,\pi ].\end{eqnarray*}$$

The stability regions of the Adams–Moulton 3 and Adams–Moulton 4 methods are clearly bounded, as shown by Figures 1(c) and (d). As in the case of the RK4 method, stiffness could appear in optically thick cells, imposing a reduction of the cell width or a switch to A-stable methods to maintain convergence. Therefore, stability constraints are clearly a disadvantage when using high-order Adams–Moulton methods. In this sense, the second Dahlquist barrier clarifies the situation, stating that an A-stable linear multistep method has an order of accuracy of $p\leqslant 2$.

Although the use of optical depth usually supports the assumption of a constant eigenvalue λ in Equation (4), this assumption limits the validity of the stability analysis. An additional limitation is given by the fact that the stability analysis of linear multistep methods assumes a homogeneous discrete grid. Strongly variable meshes, such as logarithmically spaced grids, could therefore alter the stability conditions for linear multistep schemes.

As mentioned above, the Adams–Moulton methods are implicit and require the solution of a 4 × 4 implicit linear system. The similarity to the DELO methods suggests a similar computational cost.

Here, the cubic Hermite interpolation is chosen to approximate ${\boldsymbol{F}}$ in Equation (2). For notational simplicity one defines the normalized variable $t\in [0,1]$ as

$$\begin{eqnarray*}&&t=\displaystyle \frac{s-{s}_{k}}{{\rm{\Delta }}{s}_{k}},\ \mathrm{for}\ s\in [{s}_{k},{s}_{k+1}].\end{eqnarray*}$$

The cubic Hermite interpolation, approximating ${\boldsymbol{F}}$ inside the interval $[0,1]$, reads

$$\begin{eqnarray}{\boldsymbol{F}}(t) & \approx & {{\boldsymbol{F}}}_{k}\cdot (1-3{t}^{2}+2{t}^{3})+{{\boldsymbol{F}}}_{k}^{{\prime} }\cdot {\rm{\Delta }}{s}_{k}(t-2{t}^{2}+{t}^{3})\\ & & +{{\boldsymbol{F}}}_{k+1}\cdot (3{t}^{2}-2{t}^{3})+{{\boldsymbol{F}}}_{k+1}^{{\prime} }\cdot {\rm{\Delta }}{s}_{k}(-{t}^{2}+{t}^{3}).\end{eqnarray} \tag{ 9 }$$

In addition to the grid values ${{\boldsymbol{F}}}_{k}$ and ${{\boldsymbol{F}}}_{k+1}$, the first derivatives ${{\boldsymbol{F}}}_{k}^{{\prime} }$ and ${{\boldsymbol{F}}}_{k+1}^{{\prime} }$ are also specified and Equation (9) provides the unique third-degree polynomial that matches both node values and node first derivatives at s_k and ${s}_{k+1}$. Moreover, the first derivative of ${\boldsymbol{F}}$ satisfies

$$\begin{eqnarray*}{\boldsymbol{F}}^{\prime} (s) & = & -{\bf{K}}^{\prime} (s){\boldsymbol{I}}(s)-{\bf{K}}(s){\boldsymbol{I}}^{\prime} (s)+{\boldsymbol{\epsilon }}^{\prime} (s)\\ & = & [{\bf{K}}(s){\bf{K}}(s)-{\bf{K}}^{\prime} (s)]{\boldsymbol{I}}(s)-{\bf{K}}(s){\boldsymbol{\epsilon }}(s)+{\boldsymbol{\epsilon }}^{\prime} (s),\end{eqnarray*}$$

where Equation (1) is used to replace the Stokes vector first derivative ${\boldsymbol{I}}^{\prime}$. Inserting numerical approximations, the first derivatives ${{\boldsymbol{F}}}_{k}^{{\prime} }$ and ${{\boldsymbol{F}}}_{k+1}^{{\prime} }$ can be written as

$$\begin{eqnarray*}{{\boldsymbol{F}}}_{k}^{{\prime} } & = & [{{\bf{K}}}_{k}{{\bf{K}}}_{k}-{{\bf{K}}}_{k}^{{\prime} }]{{\boldsymbol{I}}}_{k}-{{\bf{K}}}_{k}{{\boldsymbol{\epsilon }}}_{k}+{{\boldsymbol{\epsilon }}}_{k}^{{\prime} },\\ {{\boldsymbol{F}}}_{k+1}^{{\prime} } & = & [{{\bf{K}}}_{k+1}{{\bf{K}}}_{k+1}-{{\bf{K}}}_{k+1}^{{\prime} }]{{\boldsymbol{I}}}_{k+1}-{{\bf{K}}}_{k+1}{{\boldsymbol{\epsilon }}}_{k+1}+{{\boldsymbol{\epsilon }}}_{k+1}^{{\prime} }.\end{eqnarray*}$$

Replacing the integrand ${\boldsymbol{F}}$ in Equation (2) with the cubic Hermite interpolant given by Equation (9), one evaluates the integral by parts. Making use of the previous identities for ${{\boldsymbol{F}}}_{k}^{{\prime} }$ and ${{\boldsymbol{F}}}_{k+1}^{{\prime} }$ and performing some algebra, one recovers the following implicit linear system

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{k+1}{{\boldsymbol{I}}}_{k+1}={{\boldsymbol{\Phi }}}_{k}{{\boldsymbol{I}}}_{k}+{{\boldsymbol{\Psi }}}_{k+1}+{{\boldsymbol{\Psi }}}_{k},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray*}{{\boldsymbol{\Phi }}}_{k} & = & 1-\displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\bf{K}}}_{k}+\displaystyle \frac{{\rm{\Delta }}{s}_{k}^{2}}{12}[{{\bf{K}}}_{k}{{\bf{K}}}_{k}-{{\bf{K}}}_{k}^{{\prime} }],\\ {{\boldsymbol{\Phi }}}_{k+1} & = & 1+\displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\bf{K}}}_{k+1}+\displaystyle \frac{{\rm{\Delta }}{s}_{k}^{2}}{12}[{{\bf{K}}}_{k+1}{{\bf{K}}}_{k+1}-{{\bf{K}}}_{k+1}^{{\prime} }],\\ {{\boldsymbol{\Psi }}}_{k} & = & \displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\boldsymbol{\epsilon }}}_{k}+\displaystyle \frac{{\rm{\Delta }}{s}_{k}^{2}}{12}[{{\boldsymbol{\epsilon }}}_{k}^{{\prime} }-{{\bf{K}}}_{k}{{\boldsymbol{\epsilon }}}_{k}],\\ {{\boldsymbol{\Psi }}}_{k+1} & = & \displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\boldsymbol{\epsilon }}}_{k+1}-\displaystyle \frac{{\rm{\Delta }}{s}_{k}^{2}}{12}[{{\boldsymbol{\epsilon }}}_{k+1}^{{\prime} }-{{\bf{K}}}_{k+1}{{\boldsymbol{\epsilon }}}_{k+1}].\end{eqnarray*}$$

The one-step numerical method described by Equation (10) corresponds exactly to the one proposed by Bellot Rubio et al. (1998), but here it is derived through a different strategy. Moreover, the first derivatives ${\bf{K}}^{\prime}$ and ${\boldsymbol{\epsilon }}^{\prime}$ are usually not provided by the problem and must be numerically approximated. The accuracy of the numerical derivatives could affect the order of accuracy of the entire method and Bellot Rubio et al. (1998) first opted for an expensive procedure based on cubic spline interpolation. When considering a physical quantity u, they also mentioned the possibility of calculating the numerical first derivative ${u}_{k}^{{\prime} }$ at the node s_k, assuming a parabolic dependence along ${s}_{k-1}$, s_k, and ${s}_{k+1}$. The explicit formula adapted to a non-uniform spatial grid reads

$$\begin{eqnarray}&&{u}_{k}^{{\prime} }={w}_{k-1}{u}_{k-1}+{w}_{k}{u}_{k}+{w}_{k+1}\,{u}_{k+1},\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray*}{w}_{k-1} & = & \displaystyle \frac{1}{{\rm{\Delta }}{s}_{k+1}+{\rm{\Delta }}{s}_{k}}-\displaystyle \frac{1}{{\rm{\Delta }}{s}_{k}},\\ {w}_{k} & = & \displaystyle \frac{1}{{\rm{\Delta }}{s}_{k}}-\displaystyle \frac{1}{{\rm{\Delta }}{s}_{k+1}},\\ {w}_{k+1} & = & \displaystyle \frac{1}{{\rm{\Delta }}{s}_{k+1}}-\displaystyle \frac{1}{{\rm{\Delta }}{s}_{k+1}+{\rm{\Delta }}{s}_{k}},\end{eqnarray*}$$

which is a second-order accurate approximation for the first derivatives. Fritsch & Butland (1984) proposed an alternative formula to recover second-order accurate first derivatives for producing monotone piecewise cubic Hermite interpolants.

The local truncation error is due to the fact that the integrand ${\boldsymbol{F}}$ is approximated by a cubic Hermite polynomial. One can show that the cubic Hermite interpolant is fourth-order accurate if the derivatives are at least third-order, third-order if the derivatives are second-order, and so on (Dougherty et al. 1989). Therefore, assuming derivatives are at least third-order accurate, the global error scales as $O({\rm{\Delta }}{s}^{4})$, indicating the cubic Hermitian method as fourth-order accurate (see Table 1). In confirmation of this, Bellot Rubio et al. (1998) perform an alternative convergence analysis, providing the same order of accuracy. The local truncation error analysis based on Taylor expansion proposed in Appendix B reveals that second-order accurate numerical derivatives, such as, for instance, the one given by Equation (11), are already sufficient to maintain the fourth-order accuracy of the cubic Hermitian method, as confirmed by Figure 3. However, De la Cruz Rodríguez & Piskunov (2013) indicate the cubic Hermitian method as third-order accurate. It can be surmised that this is due to the first-order accurate derivatives used in the method there, which are not sufficient to maintain fourth-order accuracy.

Zoom In Zoom Out Reset image size

The stability function of the cubic Hermitian method is easily deduced through the IVP (4) and reads

$$\begin{eqnarray*}&&{\phi }_{{\rm{H}}}(z)=\displaystyle \frac{1+z/2+{z}^{2}/12}{1-z/2+{z}^{2}/12},\end{eqnarray*}$$

with $z=\lambda {\rm{\Delta }}t$. Stability is then given by the condition $\parallel {\phi }_{{\rm{H}}}(z)\parallel \lt 1$. As displayed in Figure 1(b), the stability region contains the whole left side of the complex plane, indicating the cubic Hermitian method as A-stable. Paper I argues that this is an important feature to avoid numerical instability in the formal solution.

Once more, the assumption of a constant eigenvalue λ in Equation (4) is a limitation for the stability analysis because variations of ${\bf{K}}$ along the integration path affect the stability region of the cubic Hermitian method. The use of the optical depth scale usually supports the assumption of a constant eigenvalue λ in Equation (4).

The cubic Hermitian method is implicit and it requires the solution of the 4 × 4 linear system given by Equation (10). The additional matrix-by-matrix multiplications and the calculation of numerical derivatives increase the computational effort. However, Bellot Rubio et al. (1998) suggest that, when a certain accuracy is required, the high accuracy of the cubic Hermitian method allows one to use coarser spatial grids, reducing the total computational cost of the problem.

If the quadratic Bézier curve is used to approximate ${\boldsymbol{F}}$ inside the interval [0, 1], one gets

$$\begin{eqnarray}&&{\boldsymbol{F}}(t)\approx {{\boldsymbol{F}}}_{k}\cdot {(1-t)}^{2}+{\boldsymbol{C}}\cdot 2t(1-t)+{{\boldsymbol{F}}}_{k+1}\cdot {t}^{2}.\end{eqnarray} \tag{ 12 }$$

While the grid values ${{\boldsymbol{F}}}_{k}$ and ${{\boldsymbol{F}}}_{k+1}$ determine the start and the end points of the curve, the presence of the control point ${\boldsymbol{C}}$ allows one to shape the Bézier curve inside the interval. Auer (2003) proposed two different expressions for the control point ${\boldsymbol{C}}$, such that the quadratic Bézier curve corresponds to a quadratic Hermite polynomial, namely,

$$\begin{eqnarray*}{{\boldsymbol{C}}}^{(A)} & = & {{\boldsymbol{F}}}_{k}+\displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\boldsymbol{F}}}_{k}^{{\prime} },\\ {{\boldsymbol{C}}}^{(B)} & = & {{\boldsymbol{F}}}_{k+1}-\displaystyle \frac{{\rm{\Delta }}{s}_{k}}{2}{{\boldsymbol{F}}}_{k+1}^{{\prime} },\end{eqnarray*}$$

intending to maximize the accuracy of the interpolation. In fact, Auer (2003) points out that from the standpoints of continuity and accuracy, Hermite interpolation is preferred. Moreover, De la Cruz Rodríguez & Piskunov (2013) suggest that if both ${{\boldsymbol{C}}}^{(A)}$ and ${{\boldsymbol{C}}}^{(B)}$ can be computed, it is desirable to take the mean, i.e.,

$$\begin{eqnarray}&&{\boldsymbol{C}}=\displaystyle \frac{{{\boldsymbol{C}}}^{(A)}+{{\boldsymbol{C}}}^{(B)}}{2},\end{eqnarray} \tag{ 13 }$$

recovering a more symmetric interpolation. Replacing the integrand ${\boldsymbol{F}}$ in Equation (2) with the quadratic Bézier interpolant given by Equation (12) and inserting the symmetric control point given by Equation (13), one evaluates the integral by parts. Making use of the identities for ${{\boldsymbol{F}}}_{k}^{{\prime} }$ and ${{\boldsymbol{F}}}_{k+1}^{{\prime} }$ given in the previous section and performing some algebra, one recovers the implicit linear system described by Equation (10). Therefore, the obtained method corresponds exactly to the cubic Hermitian method, sharing its order of accuracy, stability, and computational cost properties.

This result intuitively explains the fourth-order accuracy obtained by the quadratic DELO-Bézier method described by De la Cruz Rodríguez & Piskunov (2013), for which the complexity of the method prevents an analytical prediction of the order of accuracy.

If the cubic Bézier curve is used to approximate ${\boldsymbol{F}}$ inside the interval $[0,1]$, one gets

$$\begin{eqnarray}{\boldsymbol{F}}(t) & \approx & {{\boldsymbol{F}}}_{k}\cdot {(1-t)}^{3}+{\tilde{{\boldsymbol{C}}}}_{k}\cdot 3t{(1-t)}^{2}\\ & & +{\tilde{{\boldsymbol{C}}}}_{k+1}\cdot 3{t}^{2}(1-t)+{{\boldsymbol{F}}}_{k+1}\cdot {t}^{3}.\end{eqnarray} \tag{ 14 }$$

The cubic Bézier curve is forced to be identical to the cubic Hermite polynomial given by Equation (9) by adopting the following control points:

$$\begin{eqnarray*}{\tilde{{\boldsymbol{C}}}}_{k} & = & {{\boldsymbol{F}}}_{k}+\displaystyle \frac{{\rm{\Delta }}{s}_{k}}{3}{{\boldsymbol{F}}}_{k}^{{\prime} },\\ {\tilde{{\boldsymbol{C}}}}_{k+1} & = & {{\boldsymbol{F}}}_{k+1}-\displaystyle \frac{{\rm{\Delta }}{s}_{k}}{3}{{\boldsymbol{F}}}_{k+1}^{{\prime} }.\end{eqnarray*}$$

Therefore, if the integrand ${\boldsymbol{F}}$ in Equation (2) is replaced by the cubic Bézier curve given by Equation (14) with the control points specified above, the resulting implicit linear system corresponds, once again, to the one given by Equation (10).