Beyond the Runge-Kutta-Wentzel-Kramers-Brillouin method

TL;DR

This paper investigates advanced numerical methods based on Magnus expansion and WKB for solving highly oscillatory differential equations, demonstrating improved speed and exploring new variants like Jordan-Magnus, with applications to cosmological models.

Contribution

It introduces higher-dimensional Magnus expansion-based methods, compares them with WKB, and proposes a new Jordan-Magnus approach for better performance in oscillatory systems.

Findings

Magnus expansion methods increase computation speed by ~50%

WKB and Magnus are distinct mathematical approaches

Jordan-Magnus shows WKB-like performance in 2D systems

Abstract

We explore higher-dimensional generalizations of the Runge-Kutta-Wentzel-Kramers-Brillouin method for integrating coupled systems of first-order ordinary differential equations with highly oscillatory solutions. Such methods could improve the performance and adaptability of the codes which are used to compute numerical solutions to the Einstein-Boltzmann equations. We test Magnus expansion-based methods on the Einstein-Boltzmann equations for a simple universe model dominated by photons with a small amount of cold dark matter. The Magnus expansion methods achieve an increase in run speed of about 50% compared to a standard Runge-Kutta integration method. A comparison of approximate solutions derived from the Magnus expansion and the Wentzel-Kramers-Brillouin (WKB) method implies the two are distinct mathematical approaches. Simple Magnus expansion solutions show inferior long range accuracy compared to WKB. However we also demonstrate how one can improve on the standard Magnus approach to obtain a new "Jordan-Magnus" method. This has a WKB-like performance on simple two-dimensional systems, although its higher-dimensional generalization remains elusive.

Tables

Table 1. Table 1: Ratio of integration times vs. RKF4(5) time for the different plots. Averaged over 20 runs.

Method	Airy eq.	Burst eq. ($n=40$)
First Magnus	$\mathrm{\hspace{0.33em}0.115}\pm 0.004$	$0.78\pm 0.02$
Magnus 4∘ GL	$\mathrm{\hspace{0.33em}0.244}\pm 0.004$	$0.81\pm 0.01$
Second Magnus	$\mathrm{\hspace{0.33em}0.253}\pm 0.003$	$0.503\pm 0.008$
RKWKB	$\mathrm{\hspace{0.33em}0.0029}\pm 0.0003$	$0.0065\pm 0.0004$
Jordan-Magnus	$\mathrm{\hspace{0.33em}0.050}\pm 0.007$	$0.246\pm 0.028$

Table 2. Table 2: Ratio of integration times for different methods vs. RKF4(5) for the photon & CDM system. Parameter settings and initial conditions are as in Figs.   8 and  9 . Averaged over 20 runs.

Method	Ratio
Magnus 6∘ GL	$\mathrm{\hspace{0.33em}1.58}\pm 0.06$
First Magnus	$\mathrm{\hspace{0.33em}0.71}\pm 0.03$
Magnus 4∘ GL	$\mathrm{\hspace{0.33em}0.62}\pm 0.03$

Equations

$$\dot{\mathbf{x}}(t)=\text{A}(t)\mathbf{x},$$

$$\ddot{x}=-{\omega }^2(t)x.$$

$$x=\exp\Big(\sum_{n=0}^{\infty}{S_{n}}(t)\Big{missing}).$$

$$S_0$$

$$S_1$$

$$2{\dot{S}}_0{\dot{S}}_n+{\ddot{S}}_{n-1}+$$

$$x(t)=\frac{A}{\sqrt{\omega }}{e}^{i{\int }_0^t\omega (t_1)\text{differential}{t}_1}+\frac{B}{\sqrt{\omega }}{e}^{-i{\int }_0^t\omega (t_1)\text{differential}{t}_1}.$$

$$\ddot{x}+{\omega }^{2}x=\pm i\dot{\omega }(t)x,$$

$$\bm{x}(t)=\exp\Big(\sum_{n=1}^{\infty}\textup{{\textOmega}}_{n}(t)\Big{missing})\bm{x}(0),$$

$${\mathbf{\Omega }}_1(t)$$

$${\mathbf{\Omega }}_2(t)$$

$$\displaystyle\begin{split}\bm{\Omega}_{3}(t)&=\frac{1}{6}\int_{0}^{t}\differential t_{1}\int_{0}^{t_{1}}\differential t_{2}\int_{0}^{t_{2}}\Big{(}\commutator{\textup{{A}}(t_{1})}{\commutator{\textup{{A}}(t_{2})}{\textup{{A}}(t_{3}}}\\ &+\commutator{\textup{{A}}(t_{3})}{\commutator{\textup{{A}}(t_{2})}{\textup{{A}}(t_{1}}}\Big{)}\differential t_{3}.\end{split}$$

$${\int }_{t_0}^T{\Vert \text{A}(s)\Vert }_2\text{differential}s<\pi ,$$

$${\Vert \text{A}(s)\Vert }_2={\text{max}}_{i=1,...,n}\sqrt{{\lambda }_i({\text{A}}^T\text{A})},$$

$$\mathbf{x}(t)=\text{M}(t_0,t)\mathbf{x}(t_0).$$

$${\mathbf{x}}_{n+1}=\text{M}(t_n,t_{n+1}){\mathbf{x}}_n.$$

$$\text{M}(t_0,t)=1+(t-t_0)\text{A}(t_0).$$

$$\text{A}(t)=\sum _{j=0}^{\infty }{a}_j(t-t_{1/2}{)}^j\text{where}{a}_j=\frac{1}{j!}{\text{derivative}[n]\text{A}t|}_{t=t_{1/2}}.$$

$${\text{textOmega}}^{[4]}={\alpha }_1-\frac{1}{12}\text{commutator}{\alpha }_1{\alpha }_2+\mathcal{O}(h^5),$$

$$\text{differential}{s}^2=a^2(\eta )[(1+2\Psi )\text{differential}{\eta }^2-(1-2\Phi ){\delta }_{ij}\text{differential}{x}^i\text{differential}{x}^j].$$

$$\ddot{\Phi }+\frac{4}{\eta }\dot{\Phi }+\frac{1}{3}{k}^2\Phi =0.$$

$$\Phi =A_k\frac{j_1(x)}{x}+B_k\frac{y_1(x)}{x},$$

$$\begin{array}{r}{j}_1(x)=\frac{\sin (x)}{x^2}-\frac{\cos (x)}{x},\\ y_1(x)=-\frac{\cos (x)}{x^2}-\frac{\sin (x)}{x}.\end{array}$$

$${\delta }_r=4{\Theta }_0=-\frac{2}{3}(k\eta {)}^2\Phi -2\eta \dot{\Phi }-2\Phi ,$$

$$\text{setcounter}MaxMatrixCols14\left[\begin{array}{cccccccccccccc}-2\mathcal{H}{\Omega }_{\gamma }& -k& \frac{12{\mathcal{H}}^3{\Omega }_{\gamma }}{k^2}& 0& 0& 0& -2\mathcal{H}{\Omega }_{\nu }& 0& \frac{12{\mathcal{H}}^3{\Omega }_{\nu }}{k^2}& -\frac{1}{2}\mathcal{H}{\Omega }_{\text{dm}}& 0& -\frac{1}{2}\mathcal{H}{\Omega }_b& 0& -\frac{k^2}{3\mathcal{H}}-\mathcal{H}\\ \frac{1}{3}k& \dot{\tau }& -\frac{4{\Omega }_{\gamma }{\mathcal{H}}^2}{k}-\frac{2}{3}k& 0& 0& 0& 0& 0& -\frac{4{\mathcal{H}}^2{\Omega }_{\nu }}{k}& 0& 0& 0& -\frac{1}{3}i\dot{\tau }& \frac{1}{3}k\\ 0& \frac{2}{5}k& \frac{21}{20}\dot{\tau }& \frac{1}{20}\dot{\tau }& 0& \frac{1}{20}\dot{\tau }& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -\frac{1}{2}\dot{\tau }& \frac{1}{2}\dot{\tau }& -k& -\frac{1}{2}\dot{\tau }& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{3}k& \dot{\tau }& -\frac{2}{3}k& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -\frac{1}{20}\dot{\tau }& -\frac{1}{20}\dot{\tau }& \frac{2}{5}k& \frac{19}{20}\dot{\tau }& 0& 0& 0& 0& 0& 0& 0& 0\\ -2\mathcal{H}{\Omega }_{\gamma }& 0& \frac{12{\mathcal{H}}^3{\Omega }_{\gamma }}{k^2}& 0& 0& 0& -2\mathcal{H}{\Omega }_{\nu }& -k& \frac{12{\mathcal{H}}^3{\Omega }_{\nu }}{k^2}& -\frac{1}{2}\mathcal{H}{\Omega }_{\text{dm}}& 0& -\frac{1}{2}\mathcal{H}{\Omega }_b& 0& -\frac{k^2}{3\mathcal{H}}-\mathcal{H}\\ 0& 0& -\frac{4{\mathcal{H}}^2{\Omega }_{\gamma }}{k}& 0& 0& 0& \frac{1}{3}k& 0& -\frac{4{\Omega }_{\nu }{\mathcal{H}}^2}{k}-\frac{2}{3}k& 0& 0& 0& 0& \frac{1}{3}k\\ 0& 0& 0& 0& 0& 0& 0& \frac{2}{5}k& 0& 0& 0& 0& 0& 0\\ -6\mathcal{H}{\Omega }_{\gamma }& 0& \frac{36{\mathcal{H}}^3{\Omega }_{\gamma }}{k^2}& 0& 0& 0& -6\mathcal{H}{\Omega }_{\nu }& 0& \frac{36{\mathcal{H}}^3{\Omega }_{\nu }}{k^2}& -\frac{3}{2}\mathcal{H}{\Omega }_{\text{dm}}& 0& -\frac{3}{2}\mathcal{H}{\Omega }_b& 0& -\frac{k^2}{\mathcal{H}}-3\mathcal{H}\\ 0& 0& \frac{12i{\mathcal{H}}^2{\Omega }_{\gamma }}{k}& 0& 0& 0& 0& 0& \frac{12i{\mathcal{H}}^2{\Omega }_{\nu }}{k}& 0& -H& 0& 0& -ik\\ -6\mathcal{H}{\Omega }_{\gamma }& 0& \frac{36{\mathcal{H}}^3{\Omega }_{\gamma }}{k^2}& 0& 0& 0& -6\mathcal{H}{\Omega }_{\nu }& 0& \frac{36{\mathcal{H}}^3{\Omega }_{\nu }}{k^2}& -\frac{3}{2}\mathcal{H}{\Omega }_{\text{dm}}& 0& -\frac{3}{2}\mathcal{H}{\Omega }_b& -ik& -\frac{k^2}{\mathcal{H}}-3\mathcal{H}\\ 0& \frac{4i\dot{\tau }{\Omega }_{\gamma }}{{\Omega }_b}& \frac{12i{\mathcal{H}}^2{\Omega }_{\gamma }}{k}& 0& 0& 0& 0& 0& \frac{12i{\mathcal{H}}^2{\Omega }_{\nu }}{k}& 0& 0& 0& \frac{4\dot{\tau }{\Omega }_{\gamma }}{3{\Omega }_b}-H& -ik\\ -2\mathcal{H}{\Omega }_{\gamma }& 0& \frac{12{\mathcal{H}}^3{\Omega }_{\gamma }}{k^2}& 0& 0& 0& -2\mathcal{H}{\Omega }_{\nu }& 0& \frac{12{\mathcal{H}}^3{\Omega }_{\nu }}{k^2}& -\frac{1}{2}\mathcal{H}{\Omega }_{\text{dm}}& 0& -\frac{1}{2}\mathcal{H}{\Omega }_b& 0& -\frac{k^2}{3\mathcal{H}}-\mathcal{H}\end{array}\right]$$

$$\text{vspace}0.2cm\left[\begin{array}{c}{\dot{\Theta }}_0\\ {\dot{\Theta }}_1\\ \dot{\delta }\\ \dot{v}\\ \dot{\Phi }\end{array}\right]=\left[\begin{array}{ccccc}-2\mathcal{H}{\Omega }_{\gamma }& -k& -\frac{1}{2}\mathcal{H}{\Omega }_{\text{dm}}& 0& -\frac{k^2}{3\mathcal{H}}-\mathcal{H}\\ \frac{k}{3}& 0& 0& 0& \frac{k}{3}\\ -6\mathcal{H}{\Omega }_{\gamma }& 0& -\frac{3}{2}\mathcal{H}{\Omega }_{\text{dm}}& -ik& -\frac{k^2}{\mathcal{H}}-3\mathcal{H}\\ 0& 0& 0& -H& -ik\\ -2\mathcal{H}{\Omega }_{\gamma }& 0& -\frac{1}{2}\mathcal{H}{\Omega }_{\text{dm}}& 0& -\frac{k^2}{3\mathcal{H}}-\mathcal{H}\end{array}\right]{\left[\begin{array}{c}{\Theta }_0\\ {\Theta }_1\\ \delta \\ v\\ \Phi \end{array}\right]}_{.}$$

$$\text{A}(k,\eta )={\left[\begin{array}{ccccc}-\frac{1}{\eta }& -k& -\frac{1}{4}& 0& -\frac{\eta k^2}{3}-\frac{1}{\eta }\\ \frac{k}{3}& 0& 0& 0& \frac{k}{3}\\ -\frac{3}{\eta }& 0& -\frac{3}{4}& -ik& -\eta k^2-\frac{3}{\eta }\\ 0& 0& 0& -\frac{1}{{\eta }^2}& -ik\\ -\frac{1}{\eta }& 0& -\frac{1}{4}& 0& -\frac{\eta k^2}{3}-\frac{1}{\eta }\end{array}\right]}_{.}\text{vspace}2.0cm$$

$$\mathbf{x}={\left[\begin{array}{c}x\\ \dot{x}\end{array}\right]}_{,}\text{A}={\left[\begin{array}{cc}0& 1\\ -{\omega }^2& 0\end{array}\right]}_{.}$$

$${\text{textOmega}}_1$$

$${\text{textOmega}}_2$$

$$a=\frac{1}{2}t{\int }_0^t{\omega }^2\text{differential}{t}_1-{\int }_0^t\text{differential}{t}_1{\int }_0^{t_1}{\omega }^2\text{differential}{t}_2=\frac{1}{2}{\int }_0^t{t}_1^2\text{derivative}t⟨{\omega }^2⟩\text{differential}{t}_1.$$

$$\text{M}={\left[\begin{array}{cc}\cos (\sqrt{⟨{\omega }^2⟩}t)& \frac{1}{\sqrt{⟨{\omega }^2⟩}}\sin (\sqrt{⟨{\omega }^2⟩}t)\\ \sqrt{⟨{\omega }^2⟩}\sin (\sqrt{⟨{\omega }^2⟩}t)& \cos (\sqrt{⟨{\omega }^2⟩}t)\end{array}\right]}_{.}$$

$$x(t)=A\cos (\sqrt{⟨{\omega }^2⟩}t)+B\frac{1}{\sqrt{⟨{\omega }^2⟩}}\sin (\sqrt{⟨{\omega }^2⟩}t).$$

Full text

††preprint: APS/123-QED

I Introduction

In the new era of high precision cosmology, one of the most important tools in the cosmologist’s toolbox is the suite of Boltzmann codes. These methods, including class [1], camb [2], and pycosmo [3] provide numerical solutions to the Einstein-Boltzmann equations, which describe the evolution of linear perturbations in the Universe. The dominant computational cost for these codes is integrating coupled systems of linear first order ordinary differential equations (ODEs) with highly oscillatory solutions [1, 2]. Currently the codes use several numerical and semianalytic approximations to get around these bottlenecks. However these approximations are designed for specific cosmological theories. A more generalized theory-independent method would make it easier to switch between different models and test new extensions to the current theory. New numerical integration methods could also reduce the run-time of the code. While the time needed to run a single simulation is currently only a few seconds, the many thousands of simulation runs needed for a typical Monte Carlo analysis add up to a substantial computational cost.

Handley, Lasenby and Hobson proposed in [4], and developed with Agocs in [5], a new numerical method based on a combination of the standard Runge-Kutta (RK) approach and the Wentzel-Kramers-Brillouin (WKB) approximation, dubbed the “Runge-Kutta-Wentzel-Kramers-Brillouin” (RKWKB) method. In highly oscillatory regions of the solution using the WKB approximation allows integration steps which span many oscillations, giving faster performance than standard RK methods. However RKWKB is limited to essentially one-dimensional systems [4, 5].

We can write a multidimensional coupled system of linear ODEs in the form of (1.1),

[TABLE]

where $\bm{x}$ is a vector of dynamical variables, A is a time-dependent matrix, and an overdot denotes differentiation with respect to time.

In the mathematics field of numerical analysis there has been substantial work examining the formal solution of (1.1), the Magnus expansion [6]. A number of authors have also used the Magnus expansion and related techniques to develop numerical methods to integrate highly oscillatory differential equations [6], for example Iserles [7, 8].

These Magnus-based numerical integrators have been extensively applied to problems in quantum and atomic physics [9, 10], as the time dependent Schrodinger equation can be written in the form of (1.1) if we use a basis of wave functions [11]. Some alternative methods for solving coupled differential equations have also been explored in the context of primordial cosmology [12]. However, Magnus-based methods remain unexplored with respect to the Einstein-Boltzmann equations and cosmological systems in general. Furthermore, to our knowledge, no previous authors have tested these methods with adaptive step size control.

In this paper we examine the application of Magnus-based techniques to the Einstein-Boltzmann equations. However, we also achieve results which apply to coupled systems of first order ODEs more generally, and hence have relevance to quantum systems, along with problems in many other areas of physics and applied mathematics [4].

II Background

II.1 The WKB approximation

The WKB method is an example of multiple scale analysis [13] applied to oscillatory systems. It is most commonly applied to equations of the form of (2.1) [14],

[TABLE]

The derivation begins by considering a general solution of the form

[TABLE]

Using the method of separation of scales we can obtain expressions for $S_{n}$ [14]

[TABLE]

Truncating with only the $S_{0},S_{1}$ terms gives the classic WKB approximation:

[TABLE]

If we include only the $S_{0}$ term we obtain

[TABLE]

implying we require $|\dot{\omega}(t)|\ll|\omega^{2}(t)|$ for the approximation to be reasonable.

II.2 The Magnus expansion

The Magnus expansion [6] provides a general solution to equations of the type (1.1). As with WKB, it starts by by considering a solution in the form of the exponential of an infinite series (2.8) [6],

[TABLE]

the first terms of which are

[TABLE]

Further terms can be generated by a recursive procedure [6] and involve an increasing number of nested integrals, commutators, and powers of A.

For the $\bm{\Omega}$ series to converge, A must be sufficiently small in some sense. It can be shown [6] that the Magnus series is absolutely convergent for a general complex matrix A, for $t_{0}\leq t\leq T$, if

[TABLE]

where $\|\cdot\|_{2}$ denotes the two-norm of a matrix, which is [15]

[TABLE]

where $\lambda_{i}(\textup{{Q}})$ denotes the $i$th eigenvalue of matrix Q. However, the actual convergence domain can be larger than that given by (2.12) [6].

II.3 Stepping numerical integration methods

A general stepping numerical integration method proceeds by using an approximate solution to the equation to advance from one integration point $\bm{x}_{n}$ at time $t_{n}$ to the next $\bm{x}_{n+1}$ at time $t_{n+1}$. For a system of first order linear ODEs such as (1.1) we can write the approximate solution for initial conditions $\bm{x}(t_{0})$ in terms of a solution matrix $\textup{{M}}(t_{0},t)$ [16],

[TABLE]

Then the general stepping procedure is

[TABLE]

Runge-Kutta methods work by using truncated Taylor series solutions for $\textup{{M}}(t_{0},t)$. For example, Euler’s method, the simplest Runge-Kutta method [4], has

[TABLE]

The RKWKB method (to second order in the WKB series) is obtained by replacing the Taylor series solution with the two independent WKB solutions in (2.6). As the WKB solutions typically remain good approximations to the true solution over many oscillations, unlike truncated Taylor series solutions, the RKWKB method is able to take much larger steps than Runge-Kutta methods, which in general reduces the required run-time [5]. Handley et al. [4] and Agocs et al. [5] also enable the program to switch to a standard Runge-Kutta-Fehlberg method in slowly oscillating regions; however, in this paper we use RKWKB to refer to the simpler nonswitching approach.

We can also obtain stepping numerical integration methods of the form of (2.15) from the Magnus expansion. If we can compute the analytic forms of $\textup{{\textOmega}}_{i}$ then we can use those directly. In that case the \textOmega series truncated at $\textup{{\textOmega}}=\bm{\Omega}_{1}+\cdots+\bm{\Omega}_{i}$ where $i=2s-1$ or $2s-2$ for integer $s$ gives a method order $h^{2s+1}$ in the step size $h$ [6]. However, if $\textup{{\textOmega}}_{i}$ cannot be computed analytically, we can use an approach from Blanes et al. [6] Sec. 5.4. The A is Taylor expanded about the midpoint $t_{1/2}=t_{n}+h/2$ of each step, where $h=t_{n+1}-t_{n}$ is the step size.

[TABLE]

One can then use the defining recursion relation for the Magnus expansion to express \textOmega to the desired order in $h$. For example to order four we have

[TABLE]

where $\alpha_{j}=h^{j}a_{j-1}$. The $\alpha_{j}$’s can then be approximated using various quadrature rules (including Gauss-Legendre, Simpson, Newton-Cotes) and linear combinations of A evaluated at different points in the interval $[t_{n},t_{n+1}]$.

II.4 The Einstein-Boltzmann equations

In this section we describe the background theory and the Einstein-Boltzmann equations that will be used in Sec. V.

II.4.1 Overview

We seek to express the Einstein-Boltzmann equations in the form of (1.1) [17]. Working in the conformal Newtonian gauge in flat spacetime, the metric is given by

[TABLE]

where the scalar perturbations are $\Psi,\Phi$, and $\eta$ is conformal time [18]. As we work at early times we neglect dark energy and treat the neutrinos as relativistic. We denote the photon and neutrino temperature perturbations as $\Theta$ and $\mathcal{N}$ respectively. The analogous perturbation $\Theta_{P}$ describes the strength of the photon polarization. We treat the cold dark matter (CDM) and baryonic matter as nonrelativistic and describe them with overdensities $\delta,\delta_{b}$ and peculiar velocities $v,v_{b}$.

To achieve a linear set of equations we expand the perturbations in spatial Fourier modes, introducing wave number $k$, and we expand the relativistic perturbations in multipole moments $\Theta_{l}(k,\eta)$ for $l=0,1,2,\ldots,\infty$. The resulting $14\times 14$ $\textup{{A}}(k,\eta)$ matrix for $l\leq 2$ is shown in Fig. 1.

The conformal Hubble rate is $\mathcal{H}=\dot{a}/a$, where an overdot is used to denote the derivative with respect to conformal time. We denote the nonconformal cosmic Hubble rate as $H=\mathcal{H}/a$. The $\dot{\tau}=-n_{e}{\sigma}_{T}a$ is the conformal time derivative of the optical depth, $\tau$, with $n_{e}$ the mean electron density and ${\sigma}_{T}$ the Thompson cross section. Mean densities are described using parameters of the form $\Omega=\rho/\rho_{\textup{crit}}$ where $\rho$ is the mean density of a component and $\rho_{\textup{crit}}$ is the critical density $=\frac{3H^{2}}{8\pi G}$. These parameters are denoted $\Omega_{\gamma},\Omega_{\nu},\Omega_{\textup{dm}}$ and $\Omega_{b}$ for the photons, neutrinos, CDM and baryonic matter respectively.

II.4.2 Photon and CDM universe

As a test case we consider a simplified universe with only photons and cold dark matter. By neglecting the baryonic matter we can avoid the complexities of the Coulomb interaction between the baryons and the photons. We also neglect the photon polarization and all harmonics $l>1$, approximating the photons as a fluid [17]. In this limit the gravitational potential $\Phi=\Psi$. By making these approximations we obtain a five-dimensional system with $\bm{x}=[\Theta_{0},\Theta_{1},\delta,v,\Phi]^{T}$. The reduced dimensions makes this much more amenable than the 14-dimensional system, while still preserving some of the key features of the early universe. By testing on this simplified system we provide a proof-of-principle study for application to the full Einstein-Boltzman equations.

For adiabatic perturbations and radiation domination we obtain a second order evolution equation for the potential [17, 19]

[TABLE]

This has a general solution (2.21) [17, 19]

[TABLE]

where $x=k\eta/\sqrt{3}$ and $j_{1},y_{1}$ are first order spherical Bessel functions of the first and second kinds respectively. These take the form

[TABLE]

As radiation dominates the energy density we also find

[TABLE]

where $\delta_{r}$ is the overdensity for radiation. Hence for $k\eta\gg 1$ we have $\Theta_{0}\propto-(k\eta)^{2}\Phi\sim A_{k}\cos{(k\eta/\sqrt{3})}+B_{k}\sin{(k\eta/\sqrt{3})}$. This gives the acoustic oscillations of the photon fluid at the sound speed $c_{s}=1/\sqrt{3}$ that one would expect. The solution in (2.21) suggests the potential also oscillates for $k\eta\gtrsim 1$, but with the amplitude decaying as $1/(k\eta)^{2}$.

To use the numerical methods we must first obtain the relevant A. Working from the $14\times 14$ matrix in Fig. 1 we can derive the reduced A matrix shown in (2.24),

[TABLE]

We set the conformal time at matter-radiation equality to be $\eta_{eq}=1$, set $a(\eta_{eq})=1$ and assume $\Omega_{\gamma}(\eta_{eq})=\Omega_{\textup{dm}}(\eta_{eq})=1/2$. This gives

[TABLE]

III Comparing WKB and the Magnus expansion

III.1 Comparing analytic solutions

To investigate the relationship between the WKB approximation and the Magnus expansion we compared the lowest order approximate analytic solutions to (2.1).

Equation (2.1) can be written in the form of (1.1) by setting

[TABLE]

For the Magnus expansion, the first two \textOmega terms are given by

[TABLE]

where

[TABLE]

Including only $\textup{{\textOmega}}_{1}$ we obtain $\bm{x}(t)=\textup{{M}}(t)\bm{x}(0)$, where $\textup{{M}}=\exp(\textup{{\textOmega}}_{1}(t))$ given by (3.5),

[TABLE]

We can then write $x(t)$ as

[TABLE]

The frequency is approximated as $\omega\sim\sqrt{\langle\omega^{2}\rangle}=\sqrt{\frac{1}{t}\int_{0}^{t}\omega^{2}\differential t_{1}}$. We also get some amplitude correction from the $B$ or $M_{12}$ term. We denote this solution “first Magnus.”

Including both $\textup{{\textOmega}}_{1}$ and $\textup{{\textOmega}}_{2}$ we obtain

[TABLE]

In (3.7) the frequency is approximated as $\hat{\omega}$,

[TABLE]

where $\dot{\langle\omega^{2}\rangle}$ denotes the time derivative of $\langle\omega^{2}\rangle$. Then we obtain (3.9) and (3.10),

[TABLE]

[TABLE]

with corrections to both the phase and the amplitude, where $\hat{\phi}$ is a time dependent phase term. We denote this solution “second Magnus.” For our matrix in (3.1) $\left\|\textup{{A}}(s)\right\|_{2}=\textup{max}(\omega^{2},1)$. So (2.12) implies a convergence domain of $T-t_{0}<\pi,\int_{t_{0}}^{T}\omega^{2}\differential t<\pi$. One can also see from (3.8) that at sufficiently large times $\hat{\omega}$ becomes imaginary, giving exponentially growing or shrinking solutions.

For the WKB approximation in Eq. (2.6) the $n=0$ term approximates the frequency as

[TABLE]

while the $n=1$ term provides a correction to the amplitude. We can write (2.6) in the form $\bm{x}(t)=\textup{{M}}(t_{0},t)\bm{x}(t_{0})$. Setting $t_{0}=0$ for simplicity we obtain

[TABLE]

where C=\big{(}\frac{\omega(0)}{\omega}\big{)}^{\frac{1}{2}}\cos(\langle\omega\rangle t) and S=\big{(}\frac{\omega(0)}{\omega}\big{)}^{\frac{1}{2}}\sin(\langle\omega\rangle t) are the two independent solutions. Applying the stepping procedure in (2.15) to this M gives you the RKWKB method (to second order) [4].

We can thus observe that the Magnus expansion and the WKB method take fundamentally different approaches for solving the same problem. The classic WKB approximation approximates the oscillation frequency as $\langle\omega\rangle$, while the first Magnus solution uses the root mean square, $\sqrt{\langle\omega^{2}\rangle}$. The terms of the WKB series involve successively higher order derivatives of $\omega$, while the successive Magnus expansion terms involve successively higher order integrals.

The WKB approximation is specialized to highly oscillatory solutions, and fails when $\omega$ goes to zero. The Magnus expansion can cope with nonoscillatory solutions and is less specialized, which comes at the cost of inferior long range accuracy compared to WKB, as we shall see in Sec. III.2.

III.2 Testing on trial equations

These solutions were plotted for two trial equations taken from Handley et al. [4], namely the Airy equation and the “burst” equation, which both have analytic exact solutions. The Airy equation takes the form

[TABLE]

and the burst equation is given by

[TABLE]

where $n$ is an integer.

The analytic WKB approximate solution (2.6) and the first and second Magnus solutions (3.6) and (3.9) for the Airy equation (3.13) are plotted in Fig. 2.

For the upper plot, with $t_{0}=1$, the WKB solution remains a very good approximation to the true solution; however, the first and second Magnus solutions perform poorly. The first Magnus reproduces the increase in oscillation frequency, but does not reproduce the changes in amplitude. The second Magnus solution breaks down completely around $t-t_{0}\sim\pi$ showing the exponential growth predicted in Sec. III.1.

If we shift the starting point to $t=11$, as in the bottom plot in Fig. 2, the first Magnus solution follows the true solution more closely. However unlike the WKB solution it fails to follow the frequency of the true oscillation at large times. The second Magnus solution still breaks down as predicted.

Next we consider the burst equation (3.14). Figure 3 shows the analytic WKB, first Magnus and second Magnus approximate solutions from $t_{0}=-10$ up to $t=+10$ with $n=40$. While WKB provides an excellent match to the true solution throughout the whole range, both Magnus solutions give very poor approximations. The first Magnus solution fails to follow the changes in amplitude or frequency, and the second Magnus solution once again diverges after a relatively short period.

Figure 4 shows the result if we try to apply the WKB and Magnus solutions to only the small, central, highly oscillatory region. This region has a smaller $\dot{\omega}/\omega$ ratio, so it should be easier to model. Even then the WKB solution dramatically outperforms both Magnus solutions.

III.2.1 A closer look at the Magnus error

To show why the first Magnus solution compares so poorly in terms of accuracy we can differentiate the Magnus solutions. Let $C_{M}$ and $S_{M}$ be the two independent first Magnus solutions,

[TABLE]

Let $1+\epsilon=\omega^{2}/\langle\omega^{2}\rangle$ and $\xi=\dot{\omega}/\omega$, and then

[TABLE]

[TABLE]

One can see that for $\epsilon=0,\xi=0$ we get back (2.1). As the errors can increase with $\epsilon$ as well as $\xi$, they are potentially much larger than for WKB.

IV The Jordan-Magnus method

The WKB approximation provides good long range accuracy but is limited to one-dimensional systems. In contrast the truncated Magnus expansion provides a matrix-based approach suitable for multidimensional systems, but with inferior long range accuracy as we saw in Sec. III. Ideally we want a way to combine the positive aspects of both WKB and Magnus, that is achieve WKB-like long range accuracy with a solution that works on multidimensional systems. In this section we introduce a new approach based on the Magnus expansion, dubbed the “Jordan-Magnus” method, which goes some way toward this goal.

If the A matrix were time-independent, we could solve the system by diagonalizing A to obtain the independent oscillating modes. Consider introducing a linear transformation to (1.1). Let $\bm{x}_{P}=\textup{{P}}^{-1}\bm{x}$, and let $\textup{{J}}=\textup{{P}}^{-1}\textup{{A}}\textup{{P}}$ for some transformation matrix $\textup{{P}}(t)$. Then

[TABLE]

where $\dot{(\textup{{P}}^{-1})}$ denotes the time derivative of $\textup{{P}}^{-1}$. We want to choose a P such that $\textup{{A}}_{P}$ is as diagonal as possible. Let $\textup{{A}}_{P}=\bm{\Lambda}+\textup{{K}}$ where $\bm{\Lambda}$ is diagonal and K has zeros on the diagonal.

Once again we seek a solution in the form $\bm{x}(t)=\textup{{M}}(t)\bm{x}(0)$. We can obtain an approximate solution by neglecting K,

[TABLE]

This avoids the potentially costly matrix exponential. Alternatively we can include K and obtain

[TABLE]

which is the first Magnus solution to (4.1).

An obvious choice of linear transformation is one that diagonalizes A. Now in general A may not have a diagonal form; however, all $N\times N$ matrices can be placed in Jordan normal form [20] (which reduces to diagonal form for diagonalizable matrices).

For \textup{{A}}=\big{[}\begin{smallmatrix}0&1\;\\ -\omega^{2}&0\;\end{smallmatrix}\big{]}, A is indeed diagonalizable, and we obtain

[TABLE]

We then find the $\bm{\Lambda}$ only solution (4.2) returns

[TABLE]

where $C$ and $S$ are the independent WKB solutions from Sec. II.1. Thus we have obtained a WKB-like framework from a Magnus-like approach.

The result is not identical to the WKB solution matrix (3.12). We have no $\dot{\omega}$ terms, and indeed because the Magnus expansion proceeds by integration as opposed to differentiation, extending the series does not introduce any further derivatives. This means that, while similar, the Jordan-Magnus approach remains distinct from WKB. In certain senses this lack of $\dot{\omega}$ terms may be an advantage as derivatives can be difficult to evaluate numerically [5].

If we try plotting the results of the Jordan-Magnus method for the Airy and burst equations (Figs. 5 and 6) we can see it dramatically outperforms the first and second Magnus solutions in accuracy, achieving results similar to the WKB solution.

However, this method also has some disadvantages over standard Magnus methods. It is ultimately reliant on being able to obtain $\textup{{P}},\bm{\Lambda}$, and K as functions of $t$. In addition, numerical Jordan decomposition is inherently unstable [20]. This implies we need to compute the analytic, symbolic, Jordan normal form of A to apply the method. While this is possible for the two-dimensional (1.1) system, it quickly becomes prohibitive for more complicated systems with more dimensions.

While there is no stable numerical algorithm for obtaining Jordan-normal form, algorithms do exist for diagonalization [21]. Furthermore, although a complex square matrix is not guaranteed to be diagonalizable, a very large proportion of such matrices are diagonalizable [22]. Hence diagonalization algorithms could form the basis for a more fully numerical Jordan-Magnus method, although at present efficient methods remain elusive.

V Results

A python code was developed [23] to implement linear numerical integration methods of the form of (2.15) using the WKB, truncated Magnus, and Jordan-Magnus approximate solutions. The numpy package [24] was used to perform the numerical calculations, and the sympy package [25] was used to perform the symbolic manipulation. For comparison a RKF4(5) method and a second order (in the WKB series) RKWKB method were also implemented. The method of adaptive step size control used is described in Appendix A.

Figure 7 shows the result of the numerical code for the different methods on the Airy and burst equations described in Sec. III.2. For both equations we can see that the RKF4(5) method requires the smallest step size, and the Magnus methods use roughly similar intermediate step sizes. As expected the RKWKB method gives the largest step sizes, quickly reaching the maximum allowed step size $h_{\textup{max}}$, although it also gives the largest relative error. The Jordan-Magnus method gives smaller step sizes than RKWKB, but larger step sizes compared to the other methods.

Each plot includes a note of $T$, the time taken for the integration for each method. This of course varies substantially depending on the machine and the code implementation. However the ratio between the run-times for the different methods should give a reasonable indication of the relative computational cost. Table 1 shows the averaged results for Fig. 7. We can see that the methods with the largest step size give the smallest time ratio and thus the smallest computational cost. The RKWKB method performs best and Jordan-Magnus second best. This result also supports the approach taken in the previous sections of seeking an approximate solution with good long range accuracy, as this allows large integration step sizes.

We note that the relative errors shown in Fig. 7 are not representative of what could be achieved by an optimized code (for example see Agocs et al. [5] for a practical implementation of RKWKB to fourth order in the WKB series). Nonetheless, the results concerning the relative step sizes and relative run-times achieved should still apply.

V.1 Application to the Einstein-Boltzmann equations: Photon and CDM universe

Finally we applied the methods and numerical integration code to the Einstein-Boltzmann equations for a simple radiation dominated universe with only photons and CDM, as described in Sec. II.4.2.

Unfortunately it was not possible to use the new Jordan-Magnus method on this system due to the difficulty of obtaining a symbolic, analytic Jordan normal form. However we can test the analytic first order Magnus solution, fourth and sixth order numerical Magnus solutions, and the RKF4(5) method. The results for $\Theta_{0}$ and $\Phi$ are shown in Figs. 8 and 9.

Fitting a solution of the form $A_{0}\cos{(k\eta}/\sqrt{3})+B_{0}\sin{(k\eta}/\sqrt{3})$, where $A_{0}$ and $B_{0}$ are constants, to the results gives good agreement for Fig. 8 for all but the earliest times as expected.

For Fig. 9, fitting the adiabatic solution for $\Phi$ (2.21) to the initial conditions gives good agreement only for very early times up to about $\eta\approx 0.015$. However, a separate adiabatic solution (2.21) fitted to the RKF4(5) data for $\eta>0.015$ gives good agreement for $\eta\gtrsim 0.03$ onwards. Hence while the results do not fit a single adiabatic solution they can be well described by two adiabatic regimes with a transition region between $0.015\lesssim\eta\lesssim 0.03$.

The plots of $\log_{10}(h)$ show RKF4(5) gives the smallest step size, the sixth order GL Magnus method has intermediate step sizes, and the first Magnus and fourth order GL Magnus methods give the largest step sizes.

As in Sec. V we can compare the average ratio of run-times (Table 2). This shows that the sixth order GL Magnus method actually underperforms RKF4(5), with the smaller number of steps not compensating for the larger time per step. However, the analytic first Magnus and fourth order GL Magnus methods both give improved performance compared to RKF4(5).

This shows that despite their inferior performance compared to WKB in Sec. III, Magnus expansion based methods can still outperform RKF4(5) on a real cosmological higher-dimensional system.

VI Conclusion

In this work we sought to explore methods for efficiently integrating coupled systems of first order ODEs with highly oscillatory solutions, which would be a higher-dimensional extension of the RKWKB method, with particular applications to cosmological systems of interest.

The work of previous authors [7, 8, 4, 5] suggests methods based on the Magnus expansion could be candidates for a multidimensional equivalent of RKWKB. However comparison of the analytic Magnus expansion and WKB approximate solutions implies the two methods are fundamentally distinct mathematical approaches. For equations of the form of (2.1) WKB approximates the frequency as $\langle\omega\rangle$, whereas the first Magnus solution approximates it as $\sqrt{\langle\omega^{2}\rangle}$. The WKB solution series involves successively higher derivatives of $\omega$ while the Magnus expansion involves successively higher order integrals. Testing on trial equations of the form of (2.1), the WKB solution give much higher accuracy. This indicates that the Magnus expansion and WKB based methods are orthogonal lines of research; however, the relationship between the two may benefit from more formal mathematical inquiry.

The simplest Magnus expansion based approximate solutions were found to gave inferior long range accuracy compared to WKB. By introducing a linear transformation we can convert the A matrix to a more diagonal form, obtaining a new “Jordan-Magnus” method. This gave excellent long range accuracy on simple two-dimensional problems which closely resembled the WKB solutions. However, currently the method relies on finding the symbolic Jordan normal form of A which is very difficult for high-dimensional cases. It could form the basis for a more fully numerical method which could be applied to arbitrarily high dimensions.

A python code was created to implement the corresponding numerical integration methods, incorporating an adaptive step size procedure and the RKWKB and RKF4(5) methods for comparison. The code leaves plenty of room for optimization, and one could certainly improve the performance by switching to a compiled language like c. On simple two-dimensional systems, simple Magnus expansion based numerical methods ran significantly faster than RKF4(5), with the Jordan-Magnus method running faster still, although not achieving the speed of RKWKB.

Finally we applied this code and these methods to a simplified set of Einstein-Boltzmann equations corresponding to photon and CDM oscillations in the radiation era. The results again demonstrate that Magnus expansion based methods performed better than RKF4(5) by giving a speed improvement of about $50\%$. However, the improvement was not as impressive as for the simpler systems, and we could not implement the Jordan-Magnus method. For large times the numerical results fit the expected solutions for adiabatic perturbations. However, we could not fit a single adiabatic solution to the numerical results for the potential $\Phi$. Instead the results suggest a transition between two different adiabatic regimes. This could indicate a breakdown of the adiabatic approximation for that system.

Future work would involve testing these Magnus expansion based methods on other multidimensional systems, including the $14\times 14$ Einstein-Boltzmann matrix for $l\leq 2$. If successful, one could extend up to high $l$ values and see how one might include these methods in future Boltzmann solver codes.

These results demonstrate that multidimensional numerical integration methods based on the Magnus expansion show promise for improving the performance of not only cosmological Boltzmann solvers but also many similar problems in other areas of physics.

Acknowledgements

J.B. thanks the University of Cambridge Physics Part III project scheme, which formed the basis for this work. W.H. was supported by a Gonville & Caius research fellowship. We are also grateful for the advice of Fruzsina Agocs and Lukas Hergt of the Kavli Institute for Cosmology, Cambridge.

Appendix A Procedure for adaptive step size control

We implemented the adaptive step size procedure laid out in Press and Flannery (1988) [26]. To estimate the error, and thus the correct step size for each step, we compute the estimate of $\bm{x}_{n+1}$ from one step of size $h$, and a separate estimate using two steps size $h/2$, denoting this $\bm{x}^{*}_{n+1}$. The error is then $\Delta_{i}=|x_{i\;{n+1}}^{*}-x_{i\;n+1}|$ for each component $i$. We define a maximum tolerable error as $\Delta^{\textup{max}}_{i}=\epsilon(a_{\textup{tol}}+r_{\textup{tol}}|x_{i\;n+1}^{*}|)$, where $a_{\textup{tol}}$ and $r_{\textup{tol}}$ parametrize the absolute and relative error tolerance, and $\epsilon$ is a small scale factor.

The procedure is then as follows:

Calculate $\bm{x}^{*}_{n+1}$ and $\bm{x}_{n+1}$ for step size $h$. 2. 2.

Calculate $\Delta_{i}$ and $\Delta^{\textup{max}}_{i}$ and the error ratio $R=\sqrt{\frac{1}{N}\sum_{i}^{N}(\Delta_{i}/\Delta^{\textup{max}}_{i})^{2}}$. 3. 3.

If $R\leq 1$, compute a new $h$ as $h_{\textup{new}}=hSR^{-1/(\nu+1)}$ where $S$ = safety factor $\sim 0.99$, $\nu$ = order of method.

Go to 5. 4. 4.

Else if $R>1$, compute a new $h$ as $h_{\textup{new}}=hSR^{-1/\nu}$.

If $h_{\textup{new}}<h_{\textup{min}}$, set $h_{\textup{new}}=h_{\textup{min}}$ and go to 5.

Otherwise go to 1. 5. 5.

If $h_{\textup{new}}<0.2h$, set $h_{\textup{new}}=0.2h$.

If $h_{\textup{new}}>h_{\textup{max}}$, set $h_{\textup{new}}=h_{\textup{max}}$.

If $h_{\textup{new}}>10h$, set $h_{\textup{new}}=10h$. 6. 6.

Use current $\bm{x}^{*}_{n+1}$ as the estimate. Advance time. Set $h=h_{\textup{new}}$. Go to next step.

The RKWKB method used the adaptive step size procedure outlined above, setting $\nu=1$. The RKF4(5) was implemented following [26]. The approach to step size control is the same, except the error is calculated using the difference between the fourth and fifth order RK methods.

The $a_{\textup{tol}}$ and $r_{\textup{tol}}$ values for the RKF4(5) and Magnus methods were chosen using trial and error to achieve the largest step sizes for the target error $\epsilon$ when applied to the burst equation as in Fig. 7. For the RKWKB method the values for the Magnus methods were used.