Numerical integrators based on the Magnus expansion for nonlinear dynamical systems

Explicit numerical integration algorithms up to order four based on the Magnus expansion for nonlinear non-autonomous ordinary differential equations are presented and tested on problems possessing qualitative (very often, geometric) features that is convenient to preserve under numerical discretization. The range of applications covers augmented dynamical systems, highly oscillatory problems and nonlinear non-autonomous partial differential equations of evolution previously discretized in space.


Introduction
Differential equations of the form where both Y and A are n×n matrices, are relevant in applications, especially when they are formulated in a particular matrix Lie group G and A belongs to the corresponding Lie algebra g. In that case, it can be shown that the solution Y (t) belongs to the given Lie group for any time t ≥ 0 (if the equation satisfies the usual regularity requirements). Relevant examples include the special orthogonal group G = SO(n) in rigid body mechanics and the symplectic group Sp(n) in Hamiltonian dynamics (in which case n is an even number) [1]. The fact that the evolution of Y is formulated in a Lie algebra is the mathematical counterpart of certain symmetries the physical system presents. Examples in point are the symplectic character of the Hamiltonian flow and the unitarity of the evolution operator in quantum mechanics.
The linear case, i.e., when Eq. (1) reads is also very relevant in applications, especially in problems related with the time-dependent Schrödinger equations, such as the average Hamiltonian theory and quantum computing. In these fields, the Magnus expansion plays a significant role, both as a theoretical tool and as a numerical integrator. In essence, the solution of Eq. (2) is expressed as the exponential of a certain infinite series each term containing iterated integrals of nested commutators of A evaluated at different times. Thus, for the first terms one has Ω 1 (t) = t 0 A(t 1 )dt 1 , where [A, B] ≡ AB − BA denotes the usual commutator.
In the particular case where A(t) is periodic in Eq. (2), i.e., when there exists some T > 0 such that A(t + T ) = A(t), the so-called Floquet-Magnus expansion has been proposed in [2] and subsequently applied in [3] in the context of solid state nuclear magnetic resonance problems, also comparing this approach with previous treatments based on average Hamiltonian theory and static perturbation theory.
There are numerous problems where Eqs. (1) and (2) have to be solved numerically, but in such a way that the qualitative properties of the exact solution are also shared by the approximations thus constructed. In particular, if Eq.
(2) describes a quantum system, it is required that any approximation to Y (t) be unitary. Otherwise, the probability conservation property might be compromised. Numerical integrators providing approximate solutions to Eqs.
(1) and (2) belonging to the same Lie group G are called Lie-group methods, and they have received a great attention in the literature during the last two decades (see, e.g., [1,4,5] and references therein).
Another area where using numerical integration methods preserving the Lie group structure is of paramount importance is in the context of the augmented dynamical systems technique. This was first proposed by Liu in [6] to integrate numerically any system of k differential equations x = f (t, x), x ∈ R k , and later applied to a variety of settings, ranging from stiff differential equations to boundary value problems and systems with constraints [7][8][9][10][11].
Essentially, the idea consists in integrating not only the state variables x, but also the norm x . To proceed, the (k + 1)-dimensional vector y = (x, x ) T is introduced, so that it obeys the differential equation It is assumed that f satisfies a Lipschitz condition and x > 0, so that the system is well defined. Otherwise, a convenient reparameterization renders the problem into this form. Although the dimension has increased by one, the augmented dynamical system (5) has an inherent symmetry, namely A ∈ so(k, 1), the Lie algebra of the Lorentz group SO(k, 1) formed by all matrices X with unit determinant such that XJX T = J, with whereas A verifies AJ + JA T = 0. Notice that y T Jy = x · x − x 2 = 0, and thus this cone condition has to be exactly preserved by any numerical integrator applied to Eq. (5) to guarantee consistency. In this way, it is argued in [6] that the so-called group preserving schemes (GPS) conserve the group structure (and thus the cone condition). In consequence, they have the same asymptotic behavior as the original continuous system and do not induce spurious solutions or ghost fixed points. The simplest method in this setting is a first-order Lie-Euler discretization, either based on the matrix exponential y n+1 = exp(hA(t n , y n ))y n or the Cayley transformation and/or Padé approximants [1,6].
As mentioned before, in many situations one has to use numerical techniques to construct approximate solutions of Eqs. (1) and (2) on a given time interval in such a way that the approximations thus obtained still possess as many qualitative properties of the exact solution as possible. Lie-group methods lead to approximations Y k ≈ Y (t k ) on a grid t 0 = 0 < t 1 < · · · < t n still belonging to the Lie group G. By contrast, general-purpose integrators such as explicit Runge-Kutta or extrapolation methods, when applied to Eqs. (1) and (2), typically furnish accurate solutions that do not stay in G, and thus they do not reproduce the correct qualitative behavior of the system, especially for long time integrations.
There are several possible strategies to integrate numerically equation (1) in such a way that the numerical approximation still belongs to G. For example, in the Runge-Kutta-Munthe-Kaas (RKMK) class of methods the problem is lifted from G to the underlying Lie algebra g with the inverse exponential map, then the associated differential equation in g is solved with a Runge-Kutta method and finally this solution is mapped back to G [1,12,13].
In the linear case, numerical methods based on the Magnus expansion have shown to be more efficient than other standard integrators whereas preserving the Lie group structure of the problem [14,15]. Given the excellent behavior shown by these schemes, it is hardly surprising that several attempts to generalize them to the nonlinear problem (1) have been conducted in the past. Thus, in [16], schemes of order p can be derived by applying the idea of collocation to the differential equation defined in g. It turns out, however, that the methods obtained are implicit, although by relaxing the collocation conditions it is still possible to get explicit methods of order up to three.
In contrast, the approach followed in [17] leads in a natural manner to explicit integrators based on the Magnus expansion which are expressed in terms of arbitrary quadrature rules. This feature allows for a great flexibility of the procedure, in the sense that many different classes of problems modelled by Eq.
(1) can be addressed essentially within the same framework just by introducing different discretizations to the integrals appearing in the formalism. Thus, one could choose a Newton-Cotes quadrature rule when A is only known at equispaced points [18], a Gaussian rule to rise the order with the minimum number of evaluations [18], a Filon-type quadrature for highly oscillatory problems [19], etc.
Whereas the Magnus expansion for the nonlinear system (1) was explicitly derived in [17] and some particular numerical integrators were presented, it is our belief that a more exhaustive analysis of the technique and in particular a sound assessment of the schemes is still necessary. Thus, in this paper we show how to adapt the technique to different types of systems, ranging from (non-)autonomous ordinary differential equations evolving in Lie groups to highly oscillatory problems to non-autonomous partial differential equations of evolution type discretized in space. In this respect, we also show how this adaptivity can be further extended by introducing variable step size techniques, thus increasing significantly the efficiency of the proposed algorithms. More specifically, (i) we show that the nonlinear Magnus expansion applied in combination with the augmented dynamical systems technique produces higher order GPS, (ii) we introduce two variants of the procedure to deal with highly oscillatory problems and (iii) we show that the technique is also able to produce accurate results when applied to reaction-diffusion equations previously discretized in space.
2 A review of numerical integrators based on the Magnus expansion

Schemes of orders 2, 3 and 4
The numerical methods we consider in this paper were first proposed in [17] and are based on a generalization of the Magnus expansion for the solution of Eq. (1). In particular, it is shown that Y (t) = exp(Ω(t))Y 0 , where successive approximations of the matrix function Ω(t) are obtained by iteration, taking Ω [0] (t) ≡ 0. In particular, it is shown that, if terms up to k = m − 2 are only kept in Eq. (8), then Equivalently, the Taylor  Thus, an approximation of order four is achieved by computing up to m = 4, etc. As a matter of fact, only appropriate estimates of the integrals are indeed necessary to get consistent approximations up to the order considered, as we will show in the sequel.
If instead of the matrix equation (1), one considers the differential equation the same procedure is obviously valid, except that now the solution is obtained as To construct a numerical integrator of order p on a time interval [t 0 , t f ] for equation (10) one takes some grid 0 = t 0 < t 1 < · · · < t N = t f with associated time increments h n = t n+1 − t n for 0 ≤ n ≤ N − 1, and then determines the recurrence where the dependence of Ω [p] (t) also on the initial time has been made explicit. In doing so, it is sufficient to approximate in a consistent way all the integrals appearing in Ω [1] , Ω [2] , . . . up to the order considered and finally to compute the action of exp(Ω [k] ) on vectors (if one is considering equation (10)) or matrices (for equation (1)). Special techniques (Chebyshev polynomials, Krylov methods, etc.) exist for the second case requiring typically less computational effort than the evaluation of the matrix exponential itself ( [20,21] and references therein).
Thus, in particular, we consider the following schemes order 2, 3 and 4.
Order two. According with the previous comments, only the first two iterations in Eq. (8) are required: If (as it is often the case), Ω [1] can be evaluated exactly, and the trapezoidal rule is applied to Ω [2] , one has If this is not possible, we can always take the approximation u = h n A(t n , y n ) in Eq. (13). This scheme will be referred to in the sequel as M2.
Order three. Several schemes are possible depending on the particular quadrature rule used to approximate the integrals in Ω [1] , Ω [2] and Ω [3] . We have taken the 2-points Gauss-Legendre quadrature and the Simpson rule [18]. Since methods based on the later provide better results for the examples considered, we collect it here (M3 scheme): y n+1 = e u 5 y n Order four. The same quadratures can be used to achieve order four by considering now up to Ω [4] in Eq. (8). If the Simpson rule is applied one ends up with exactly the same sequence as in Eq. (14) with the additional computations This feature opens a new possibility, namely the straightforward implementation of a variable step size strategy. Notice that, in the 4th-order scheme, one determines two numerical solutions at t n+1 of orders three and four, respectively:ŷ n+1 = e u 5 y n and y n+1 = e v y n .
Then the quantity is used to estimate the local error. If E r computed at t n+1 is below some prescribed tolerance tol, then the step from t n to t n+1 is accepted and then we proceed to compute the approximation to the solution at t n+2 . If E r > tol instead, then the approximation at t n+1 is rejected and a smaller step is chosen to compute a new approximation at t n+1 . In either case, it turns out that the new step size can be taken as [22] h new = s h c tol E r where h c denotes the current value of the step size and s is a "safety factor" (typically, a number close to 1, i.e., 0.8 or 0.9) chosen to decrease the probability of a rejection at the next step. With this error estimate and the resulting new step size, one typically advances with the higher-order result y n+1 [23]. Moreover, in our case the estimate E r is obtained for free, sinceŷ n+1 has to be computed anyway as an intermediate stage in the calculation of y n+1 . The resulting scheme is called M4.
In practice, we scale the i-th component of y n+1 −ŷ n+1 by a factor d i = |(y n+1 ) i | to work with relative errors and additional parameters are introduced in order not to increase nor to decrease the step size too fast [23]. As for the initial step size, a possible choice is just h 0 = tol/2.
It is worth remarking that, although Magnus integrators have been originally devised for explicitly time dependent problems, it is worth noticing that the previous algorithms also provide numerical approximations in the autonomous setting, i.e., for nonlinear equations of the form This is in contrast with the linear case, since in that case the Magnus expansion just reproduces the exact solution.

Numerical examples
To illustrate the main features of the previous schemes, we consider next two simple examples.
Example 1. The motion of a free rigid body is described by the standard Euler equations [24]: where I 1 , I 2 , I 3 > 0 are the three principal moments of inertia of the body, and Π = (Π 1 , Π 2 , Π 3 ) is the angular momentum in the moving frame. As is well known, this system possesses two quadratic invariants, namely the length of the constant angular momentum in the orthogonal body frame, and the energy of the system. When the previous schemes based on the Magnus expansion are applied to Eq. (20), C is exactly preserved by construction (up to round-off), whereas keeping track of the error in H along the integration may help to give an assessment on the performance of the different integrators.
For our experiment we take I 1 = 3, I 2 = 2, I 1 = 3/2, initial conditions Π(0) = (1, 1, 1) and integrate until t f = 100 with different Magnus methods and different step sizes and/or tolerances. Then we compute the relative error in energy at t f and plot this error in terms of the required computational cost. This can be estimated as the total number of evaluations of the matrix A plus the number of exponentials multiplied by vectors and the commutators required by each method, according with the following table: Method Order A eval Exp-vec Comm One could think that our estimate of the computational cost of the different methods is rather crude. It nevertheless captures quite accurately the relative performance of the schemes: if one considers instead the CPU time required by each integration, a similar relative configuration of the respective curves is obtained.
Example 2. The next example is designed to illustrate the use of the Magnus integrators in conjunction with the augmented dynamical systems technique.
To this end, we take the equation [25]: with exact solution Here f (t, x) = (x 2 , −x 1 − x 2 2 + ln(t)) T , and so we apply Magnus methods to equation (5), where now A is a 3 × 3 matrix. Of course, for solving this particular problem other methods can be more efficient, but we stress that our purpose here is to show that the Magnus integrators can be successfully applied in augmented dynamical systems.
We compute the numerical solution in the interval t ∈ [1, 101] with the first order scheme (7) and methods M2, M3 and M4, then determine the error as the Euclidean norm of the difference between x e and the numerical solution at t f = 101 for several step sizes and finally we plot this error as a function of the computational cost of each algorithm estimated as in Example 1. Thus we get Figure 2. It clearly illustrates the advantage of using a higher order integration scheme in problems where Liu's approach is applied: they not only preserve the inherent symmetry of the augmented dynamical system, but they do so with a much reduced computational cost for a given accuracy. We see then that, by applying Magnus expansion on augmented systems results in higher order GPS. These schemes could be applied, e.g., to get higher order approximations in combination with radial basis functions (RBFs) for PDEs [26,27].
It is worth noticing that M2, M3 and M4 have essentially the same computational cost as the first-order scheme (7), but greater accuracy is always achieved. So we recommend replacing Eq. (7) by M2, M3 or M4 whenever the augmented dynamical system technique is employed.

Generalization to highly oscillation problems
During the last decade so much effort has been devoted to the design and analysis of numerical integrators especially adapted to solve different classes of ordinary differential equations exhibiting high oscillation. For example, we can cite trigonometric integrators and their modifications for solving the highly oscillatory second order autonomous system q = −Aq + g(q) with a symmetric semi-definite matrix A [28,29], adiabatic integrators for linear systems with time-dependent skew-Hermitian matrices [4], heterogeneous multiscale methods for mechanical systems subjected to fast vibrations [30], etc., for which special analytical techniques have been developed, such as the modulated Fourier expansions [4] and the averaging based on word series [31]. In all cases, special integrators are constructed (or adapted) to particular classes of problems exhibiting high oscillation.
There are, however, highly oscillatory time-dependent systems which are not easily formulated in such a way that these special integrators can be advantageously applied, but still can be described by Eq. (10). This is the case, in particular, of classical mechanical systems with time-dependent frequencies or nonlinear equations of the form y = ε −1 A(t, y)y, where ε is a small parameter. For these systems we show next two different procedures to get numerical approximations based on the Magnus expansion depending on the the particular matrix A(t, y).
First approach. The first technique constitutes a generalization of the so-called modified Magnus methods proposed by Iserles [32,33] and further extended in [17] to the scalar equation y + a(t, y, y )y = 0. Here we consider a generic problem (10) with a high oscillatory n × n matrix A.
To begin with, assume that y n ≈ y(t n ) has been computed and the goal is to determine the approximation y n+1 at t n+1 = t n + h n . We first factorize the solution at intermediate times as By inserting Eq. (25) into y = A(t, y)y, we obtain the differential equation satisfied by z(τ ), namely where B(τ, z) = e −τ A(tn,yn) A t n + τ, e τ A(tn,yn) z − A(t n , y n ) e τ A(tn,yn) .
The new variable z(τ ) may be seen thus as a correction to the approximation provided by the first term of the Magnus expansion discretized with Euler's method, and satisfies a differential equation of the same type as y(t), but with some important differences: (i) B A by construction, and (ii) their entries are highly oscillatory since they are obtained from the entries of exp(τ A). Therefore, if the Magnus expansion is applied to Eq. (26), approximations with a significantly smaller error are expected, even if the same quadratures are used.
To be more specific, consider the first Magnus iteration applied to Eq. (26) with a fixed step size h, and denote for simplicity Ω [1] (h) ≡ Ω [1] (t n+1 , t n ), where G 1 (τ ) = e −τ A(tn,yn) A t n + τ, e τ A(tn,yn) y n − A(t n , y n ) e τ A(tn,yn) , and discretize this integral with the Simpson rule, so that since G 1 (0) = 0. The resulting scheme y n+1 = e hA(tn,yn) e u(h) y n , although simple, already provides a reasonable description of the system, as is illustrated by the following example.
Example 3. It is defined by with g(t) = t cos(t 2 ) (−4t + cos 2 (t 2 )) − 2 sin(t 2 )). The function g is chosen so that the exact solution when x 1 (0) = 1, x 2 (0) = 0 reads clearly exhibiting an increasingly oscillatory character as time goes on. System (32) can be recast in the form Eq. (10) by introducing the vector y = (x 1 , x 2 , x 3 = 1) T and the augmented matrix We integrate this problem in the interval t ∈ [0, 20] with the modified scheme (31) (which we call M1mf ), M2, M3 and M4, determining the relative error in the solution as a function of the CPU time needed for the computation.
The results are depicted in Figure 3.
We see that the modified scheme (31) with a fixed step size provides a similar accuracy as the much more elaborated 4th-order method with variable step size implementation. Notice that (modified) trigonometric integrators cannot applied to this system in a straightforward way.
Second approach. In the special, but important, situation when A(t, y) has a purely imaginary spectrum and can be diagonalized for fixed values t and y, yet another specially adapted procedure can be applied. In that case we can write so that The first term of the Magnus expansion applied to Eq. (26) reads thus Notice that the entries of the matrix in the previous integral are given by As before, the next stage in the procedure consists in applying a conveniently chosen quadrature rule to Eq. (38). Given the highly oscillatory character of Eq. (39), it seems appropriate to consider a Filon-type method [19]. Since the entries in Eq. (38) are of the form h 0 g(s) e iωs ds, we propose the quadrature rule where the weights b i are determined in such a way that Eq. (40) is exact for g(τ ) = 1, τ, τ 2 : Notice that when ω → 0 and thus the Simpson rule is obtained in this limit. In our case, since M (0) = 0, the approximation reads where now ω k ≡ λ − λ k . Finally, y n+1 = e hA(tn,yn) e u(h) y n .
Example 4. With we fix x 1 (0) = 1, x 2 (0) = 0 and compute the norm of error in t = 20. As the exact solution we take, for convenience, the output generated by the function NDSolve of Mathematica with a very stringent tolerance, although any other procedure with the same requirements could aslo be used (e.g., the DOPR853 routine designed by Hairer and Wanner [23]). We can then apply Eq. (43) but also the first approach (31). The corresponding results are collected in Figure  4. We observe that the procedure (31) (line M1mf ) gives the best results, whereas Eq. (43) (line M1ms) is as efficient as the much elaborated scheme M4. Notice that both M1mf and M1ms do not require the computation of commutators, and thus both are particularly favorable when the dimension of the problem is not so small.
4 Application to non-autonomous partial differential equations of evolution type The Magnus expansion can also be applied for the integration in time of certain linear time-dependent partial differential equations (PDEs) of evolution. A common approach in this setting is to first discretize in space, for which procedures such as finite differences, pseudo-spectral methods based on collocation with trigonometric polynomials, etc. exist. The resulting (large) system matrix of linear ordinary differential equation inherits fundamental properties of the corresponding PDE. For instance, if the original PDE corresponds to the Schrödinger equation of quantum mechanics, then the coefficient matrix of the resulting ordinary differential equation is skew-Hermitian. Although interpolatory Magnus integrators present some drawbacks in this setting due to the presence of commutators, it is possible to introduce some modifications and design efficient new classes of commutator-free schemes [34].
Motivated by the results achieved for the linear case, we try to generalize the treatment to reaction-diffusion equations. These appear in models describing the time evolution of chemical or biological species in a medium such as water or air. Another well-known example is the evolution of the the population of one or several species distributed in space under the action of two concurrent phenomena: reaction between species and diffusion which makes the species spread out in space. From a mathematical point of view, they belong to the class of semi-linear parabolic partial differential equations of the form where each component of u(x, t) ∈ R d represents the concentration (or the population) of one species, C is a diagonal matrix containing the diffusion rates and F models local interactions between species. The case F (u, t) = r(t)u(1 − u) corresponds to Fisher's (also called Kolmogorov-Petrovsky-Piskunov) equation, and is used to describe the spreading of biological populations.
For simplicity, we only consider one species and one space dimension, i.e., with periodic boundary conditions in the space domain [0, 1]. We analyze, in particular, a periodic time dependence in the reaction term, and initial condition u 0 (x) = sin(2πx). After discretization in space with finite differences we arrive at the differential equation where U = (u 1 , . . . , u N ) T ∈ R N , the second derivative has been approximated by the matrix D of size N × N given by and F (U, t) is now defined by We choose N = 100, β = ω = 1 and compute the relative error (in the 2norm) at the final time t = 1 by applying the previous algorithms based on the Magnus expansion, in particular schemes M2 and M4. We also compare with the well known 2nd-order Strang splitting scheme. This requires taking t as a new coordinate x t and consider instead the enlarged system The solution of each subsystem at time t = h is given by ϕ [1] h : , and ϕ [2] h : respectively, with c ≡ β cos(ω x t (0)). The Strang splitting applied to this system results in the following algorithm for the step t n → t n+1 , providing an approx- imation U n+1 ≈ U (t n + h): Here again we compute the reference solution by applying the function NDSolve of Mathematica with a very stringent tolerance. The corresponding results are depicted in Figure 5. Notice that the scheme M4 is more efficient that the Strang splitting SS. It is worth noticing that, although one could construct higher order methods by composing the 2nd-order Strang splitting, the resulting schemes involve necessarily some negative coefficients, and thus present severe instabilities. Only splitting or composition methods with complex coefficients are appropriate [5]. This is not the case of methods based on the Magnus expansion.
Nonlinear differential equations of the form Eq. (1) possessing specific qualitative properties appear frequently in applications. Very often these properties are associated with a special structure of the matrix A(t, y), namely A belongs to some Lie algebra. In particular, if A is skew-symmetric, then the norm of the vector solution y(t) is preserved along the evolution, a feature that seems convenient to retain also when the differential equation is approximately solved by numerical schemes. In contrast with other numerical integrators like explicit Runge-Kutta and multi-step methods, the explicit schemes based on the Magnus expansion considered here do preserve by construction these features, thus providing the corresponding approximations a qualitatively superior performance. The methods are formulated as the exponential of a conveniently chosen approximation of an integral involving the matrix A evaluated at different times and commutators for order higher than two.
Although the methods are originally designed to deal with the explicitly time-dependent nonlinear equation (1), they prove to be also competitive in the treatment of autonomous problems with exactly the same formulation, as shown by the examples considered here. This is an attractive feature of the schemes in comparison with other structure-preserving integrators such as splitting methods: although in principle they can be applied to nonautonomous problems, removing previously the time dependency by considering time as a new coordinate is typically required, and this may deteriorate their overall efficiency [4].
Using the Magnus expansion in combination with the augmented dynamical systems technique produces higher order GPS, so that we can apply the variable-step methods presented here instead of the basic first-order GPS for the same family of problems. In addition, highly oscillatory problems formulated as Eq. (1) constitute a natural area of application of these schemes. With just a previous linear transformation, M1 already provides good results for general matrices A. If A, on the other hand, has a purely imaginary spectrum and can be diagonalized, the formalism can be adapted and special quadrature rules can be implemented.
Magnus-based schemes can also be applied to the time integration of reactiondiffusion equations previously discretized in space. We have illustrated this feature by applying several Magnus integrators to the well known Fisher equation and we plan to extend the treatment to other classes of explicitly timedependent nonlinear PDEs, such as the Gross-Pitaevskii equation. In that case, however, and as in linear problems, the presence of commutators in the algorithms might be problematic, and so one possibility might be to generalize the treatment carried out in [34] to the nonlinear setting.
Analytic approximations based on the Fer expansion have also been considered in the linear case (Eq. (2) and applied to different physical problems [15,35]. In contrast to the Magnus expansion, the Fer approach consists in writing the solution of Eq. (2 as an infinite product of exponentials of matrices, each one containing higher order corrections. It is even possible to construct numerical integrators by truncating appropriately the infinite series appearing in each exponential, and replacing the integrals by suitable quadratures, just as in the case of Magnus integrators [36]. It is then natural to consider their generalization to nonlinear problems, just as we have carried out here for the Magnus expansion, and this may constitute an interesting line of research to pursue.