Symplectic integrators for second-order linear non-autonomous equations

TL;DR

This paper introduces two new symplectic integrators derived from the Magnus expansion, tailored for second-order linear non-autonomous equations, with one family optimized for low to moderate dimensions and the other for large-scale systems.

Contribution

The paper presents two novel symplectic methods based on the Magnus expansion, designed specifically for different problem sizes of second-order linear non-autonomous equations.

Findings

Effective for low to moderate dimensions

Suitable for large systems like discretized wave equations

Numerical experiments demonstrate scheme advantages

Abstract

Two families of symplectic methods specially designed for second-order time-dependent linear systems are presented. Both are obtained from the Magnus expansion of the corresponding first-order equation, but otherwise they differ in significant aspects. The first family is addressed to problems with low to moderate dimension, whereas the second is more appropriate when the dimension is large, in particular when the system corresponds to a linear wave equation previously discretised in space. Several numerical experiments illustrate the main features of the new schemes.

Figures

Tables

Table 1. Table 1 : One- and two-exponential 4th- and 6th-order symplectic method, respectively, using the sixth-order Gauss–Legendre quadrature rule.

$c_1=\frac{1}{2}-\frac{\sqrt{15}}{10},c_2=\frac{1}{2},c_3=\frac{1}{2}+\frac{\sqrt{15}}{10}.$

Table 2. Table 2 : Cost of the methods in terms of matrix–matrix products, 𝒞 𝒞 \mathcal{C} , and matrix–vector products, 𝒱 𝒱 \mathcal{V} .

Method	$𝒞$	$𝒱$
${{\rm Y}}_1^{\left[4\right]}$	$17⁤\frac{1}{3}$	–
${{\rm Y}}_1^{[4,6]}$	4	–
${{\rm Y}}_1^{[4,8]}$	5	–
${\mathrm{RK}}_4^{\left[4\right]}$	8	4
${\mathrm{RKGL}}_2^{\left[4\right]}$	$4\times {\varrho }^{\left[4\right]}$	$2\times {\varrho }^{\left[4\right]}$
${\mathrm{RKNb}}_6^{\left[4\right]}$	12	6

Table 3. Table 3 : The orders of decomposition of the two best performing methods. The better one comes first.

$\omega$	1/125	1/25	1/5	1	5	25	125
	$\epsilon =1$
$p$ of ${{\rm Y}}_1^{[4,p]}$	6, 10	6, 10	6, 10	6, 10	12, 10	8, 12	8, 12
$q$ of ${{\rm Y}}_2^{[6,q]}$	10, 8	10, 8	10, 8	10, 8	12, 8	12, 8	12, 8
	$\epsilon =1/10$
$p$ of ${{\rm Y}}_1^{[4,p]}$	10, 6	10, 6	10, 6	10, 12	12, 8	8, 12	8, 12
$q$ of ${{\rm Y}}_2^{[6,q]}$	8, 10	8, 10	8, 10	8, 10	12, 8	12, 8	12, 8

Equations

$$x^{\prime \prime }(t)+M(t)x(t)=0,x(0)=x_0,x^{\prime }(0)=x_0^{\prime },$$

$${\partial }_t^{2}u(x,t)=f(x,t){\partial }_x^{2}u(x,t)+g(x,t)u(x,t),x\in \mathbb{R},t\ge 0,$$

$$z^{\prime }(t)=A(t)z(t),\text{mbox}withA(t)=\left(\begin{array}{cc}0& I\\ -M(t)& 0\end{array}\right),$$

$$\left(\begin{array}{c}q(t_n+h)\\ p(t_n+h)\end{array}\right)=e^{\Omega (t_n,h)}\left(\begin{array}{c}q(t_n)\\ p(t_n)\end{array}\right),$$

$${{\rm Y}}_1^{[4,q]}=\left(\begin{array}{cc}I& 0\\ D_3^{[4,q]}& I\end{array}\right)\left(\begin{array}{cc}0& D_2^{[4,q]}\\ I& 0\end{array}\right)\left(\begin{array}{cc}I& 0\\ D_1^{[4,q]}& I\end{array}\right),$$

$${\Psi }_m^{[p]}=\left(\begin{array}{cc}I& h{a}_{m+1}I\\ 0& I\end{array}\right)\left(\begin{array}{cc}I& 0\\ h{C}_m& I\end{array}\right)\left(\begin{array}{cc}I& h{a}_{m}I\\ 0& I\end{array}\right)\cdots \left(\begin{array}{cc}I& 0\\ h{C}_1& I\end{array}\right)\left(\begin{array}{cc}I& a_{1}I\\ 0& I\end{array}\right),$$

$$\Phi (t_n,h)=\exp \Omega (t_n,h),\Omega (t_n,h)=\sum _{k=1}^{\infty }{\Omega }_k(t_n,h),$$

$${\Omega }_1(t_n,h)=\underset{t_n}{\overset{t_n+h}{\int }}A(t_1)d{\tau }_1,{\Omega }_2(t_n,h)=\frac{1}{2}\underset{t_n}{\overset{t_n+h}{\int }}\underset{t_n}{\overset{{\tau }_1}{\int }}\left[A({\tau }_1),A({\tau }_2)\right]d{\tau }_{2}d{\tau }_1,$$

$$\frac{d\tilde{z}(t)}{dt}=\tilde{A}(t)\tilde{z}(t),\tilde{z}(t_n)=z(t_n),t\in [t_n,t_n+h],$$

$$\tilde{z}(t_n+h)-z(t_n+h)=\mathcal{O}(h^{2s+1}).$$

$$\widetilde{A}(t)=\sum_{i=1}^{s}{\cal L}_{i}\Big{(}\frac{t-t_{n}}{h}\Big{)}A_{i}=\frac{1}{h}\sum_{i=1}^{s}\left(\frac{\sigma}{h}\right)^{i-1}\alpha_{i},\qquad\quad\sigma\in\left[-\frac{h}{2},\frac{h}{2}\right],$$

$${\alpha }_{i+1}={h^{i+1}\frac{1}{i!}\frac{d^i\tilde{A}(t_{1/2}+\sigma )}{d{\sigma }^i}|}_{\sigma =0},i=0,\dots ,s-1,$$

$${\alpha }_1=\left(\begin{array}{cc}0& I\\ -{\mu }_1& 0\end{array}\right),{\alpha }_j=\left(\begin{array}{cc}0& 0\\ -{\mu }_j& 0\end{array}\right),j>1,$$

$${\mu }_{i+1}={h^{i+1}\frac{1}{i!}\frac{d^i\tilde{M}(t_{1/2}+\sigma )}{d{\sigma }^i}|}_{\sigma =0},i=0,\dots ,s-1.$$

$${\alpha }_1=h{A}_2,{\alpha }_2=\frac{\sqrt{15}h}{3}(A_3-A_1),{\alpha }_3=\frac{10h}{3}(A_3-2A_2+A_1),$$

$${\Omega }^{[6]}={\alpha }_1+\frac{1}{12}{\alpha }_3-\frac{1}{12}[12]+\frac{1}{240}[23]+\frac{1}{360}[113]-\frac{1}{240}[212]+\frac{1}{720}[1112],$$

$$[212]=\left(\begin{array}{cc}0& 0\\ 2{\mu }_2^2& 0\end{array}\right).$$

$$E_1=\exp \left(\begin{array}{cc}D& B\\ C& -D^T\end{array}\right),E_2=\exp \left(\begin{array}{cc}0& I\\ C& 0\end{array}\right),E_3=\exp \left(\begin{array}{cc}0& 0\\ C& 0\end{array}\right).$$

$${{\rm Y}}_1^{[4]}=\left(\begin{array}{cc}I& 0\\ h{C}_2^{[4]}& I\end{array}\right)\exp \left(h\left(\begin{array}{cc}0& I\\ D_1^{[4]}& 0\end{array}\right)\right)\left(\begin{array}{cc}I& 0\\ h{C}_1^{[4]}& I\end{array}\right)$$

$${{\rm Y}}_2^{[6]}=\left(\begin{array}{cc}I& 0\\ h{C}_2^{[6]}& I\end{array}\right)\exp \left(\frac{h}{2}\left(\begin{array}{cc}0& I\\ D_2^{[6]}& 0\end{array}\right)\right)\exp \left(\frac{h}{2}\left(\begin{array}{cc}0& I\\ D_1^{[6]}& 0\end{array}\right)\right)\left(\begin{array}{cc}I& 0\\ h{C}_1^{[6]}& I\end{array}\right),$$

$$\Phi (h)=\exp \left(h\left(\begin{array}{cc}0& I\\ C& 0\end{array}\right)\right)=\left(\begin{array}{cc}\cosh h\sqrt{C}& {\sqrt{C}}^{-1}\sinh h\sqrt{C}\\ \sqrt{C}\sinh h\sqrt{C}& \cosh h\sqrt{C}\end{array}\right).$$

$$\Phi (h)=\left(\begin{array}{cc}I& 0\\ R& I\end{array}\right)\left(\begin{array}{cc}I& Q\\ 0& I\end{array}\right)\left(\begin{array}{cc}I& 0\\ R& I\end{array}\right),$$

$$\begin{array}{rl}Q(C)& =\frac{\sinh h\sqrt{C}}{\sqrt{C}}=hI+\frac{C{h}^3}{6}+\frac{C^2{h}^5}{120}+\frac{C^3{h}^7}{5040}+\frac{C^4{h}^9}{362880}+\frac{C^5{h}^{11}}{39916800}+\frac{C^6{h}^{13}}{6227020800}+\mathcal{O}\left(h^{15}\right)\\ R(C)& =\sqrt{C}\tanh \left(\frac{h\sqrt{C}}{2}\right)=\frac{Ch}{2}-\frac{C^2{h}^3}{24}+\frac{C^3{h}^5}{240}-\frac{17C^4{h}^7}{40320}+\frac{31C^5{h}^9}{725760}-\frac{691C^6{h}^{11}}{159667200}+\mathcal{O}\left(h^{13}\right)\end{array}$$

$$M_1=M(t_n+c_{1}h),M_2=M(t_n+c_{2}h),M_3=M(t_n+c_{3}h),$$

$$K=M_1-M_3,L=-M_1+2M_2-M_3,F=h^2{K}^2.$$

$$\begin{array}[]{ll}C_{1}^{[4]}=\frac{\sqrt{15}}{36}K+\frac{5}{36}L&D_{1}^{[4]}=-M_{2}\\ C_{2}^{[4]}=-\frac{\sqrt{15}}{36}K+\frac{5}{36}L&\end{array}$$

$${{\rm Y}}_1^{[4,q]}=\left(\begin{array}{cc}I& 0\\ h{C}_2^{[4]}+R_1^{[4,q]}& I\end{array}\right)\left(\begin{array}{cc}I& Q_1^{[4,q+2]}\\ 0& I\end{array}\right)\left(\begin{array}{cc}I& 0\\ h{C}_1^{[4]}+R_1^{[4,q]}& I\end{array}\right).$$

$$\begin{array}[]{ll}C_{1}^{[6]}=-\frac{\sqrt{15}}{180}K+\frac{1}{18}L+\frac{1}{12960}F&D_{1}^{[6]}=-M_{2}-\frac{4}{3\sqrt{15}}K+\frac{1}{6}L\\ C_{2}^{[6]}=+\frac{\sqrt{15}}{180}K+\frac{1}{18}L+\frac{1}{12960}F&D_{2}^{[6]}=-M_{2}+\frac{4}{3\sqrt{15}}K+\frac{1}{6}L\end{array}$$

$${{\rm Y}}_2^{[6,q]}=\left(\begin{array}{cc}I& 0\\ h{C}_2^{[6]}+R_2^{[6,q]}& I\end{array}\right)\left(\begin{array}{cc}I& Q_2^{[6,q+2]}\\ 0& I\end{array}\right)\left(\begin{array}{cc}I& 0\\ R_2^{[6,q]}+R_1^{[6,q]}& I\end{array}\right)\left(\begin{array}{cc}I& Q_1^{[6,q+2]}\\ 0& I\end{array}\right)\left(\begin{array}{cc}I& 0\\ h{C}_1^{[6]}+R_1^{[6,q]}& I\end{array}\right).$$

$$\displaystyle A(t)=B(t)+D\qquad\mbox{ with }\qquad B(t)=\left(\begin{array}[]{cc}0&0\\ -M(t)&0\end{array}\right),\qquad D=\left(\begin{array}[]{cc}0&I\\ 0&0\end{array}\right),$$

Full text

††journal: J. Comput. Appl. Maths.

1 Introduction

The problem we address in this paper is the numerical integration of the second order time-dependent linear equation

[TABLE]

where $t\in\mathbb{R}$, $x(t)\in\mathbb{C}^{r}$ and primes denote time-derivatives. It has many practical applications, for example, in periodically variable systems such as quadrupole mass filter and quadrupole devices, [14, 24], microelectromechanical systems [26], in Bose–Einstein condensates, spatially linear electric fields, dynamic buckling of structures, electrons in crystal lattices, waves in periodic media, etc. (see [5, 19, 20, 22] and references therein). In these cases the matrix $M(t)$ is usually time periodic with period $T$ ((1) is then an example of a matrix Hill equation) and of moderately large size $r$. Parametric resonances can occur so that it is very important to know the stability regions in terms of the parameters of the system.

Another important example of this type is the linear time-dependent wave equation

[TABLE]

equipped with initial conditions $u(x,0)=u_{0}(x)$, and $u_{t}(x,0)=u^{\prime}_{0}(x)$. Once discretized in space in a bounded region, eq. 2 also leads to eq. 1.

When the dimension $r$ of eq. 1 is relatively small, the numerical computation of its fundamental matrix solution for one period (usually repeatedly many times for different values of the parameters of the system) is feasible and allows one to analyse the stability of the configuration. In this case numerical schemes that involve matrix–matrix products can be used. However, when eq. 1 results from a semidiscretized PDE like eq. 2 then $r\gg 1$, and so only numerical schemes involving matrix–vector products are suitable.

Taking these considerations into account, in this work we present two classes of numerical schemes: one that involves matrix–matrix products showing a high performance for the computation of the fundamental matrix solution for one period, and another class involving matrix–vector products, addressed for the numerical treatment of PDEs.

Equation (1) can be written as a first order system by introducing new variables $q=x,p=x^{\prime}$ as

[TABLE]

and $z=(q,p)^{T}$, $\,z_{0}=(q_{0},p_{0})^{T}\in\mathbb{C}^{2r}$. The solution of eq. 3 evolves through a linear transformation (evolution operator or fundamental matrix solution) given by $z(t)=\Phi(t,0)z(0)$. In the usual case in which $M(t)$ is a real and symmetric $r\times r$ matrix valued function ($M^{T}=M$) then $\Phi(t,0)$ is a symplectic transformation. The eigenvalues of $\Phi(t,0)$ occur in reciprocal pairs, say $\{\lambda,\lambda^{*},1/\lambda,1/\lambda^{*}\}$, where $\lambda^{*}$ denotes the complex conjugate of $\lambda$. As a result, for stable systems, all of the eigenvalues must lie on the unit circle (see for example [13]). This is a very important property that is not preserved in general by standard numerical integrators. Thus, one can be forced to use very small time steps to avoid the undesirable numerical instabilities or asymptotic stabilities arising from this fact, with the resulting degradation in the efficiency of the algorithms.

The solution of the non-autonomous eq. 3 cannot be written in general in closed form. Nevertheless, we can always formally write it as a single exponential by using the Magnus expansion [11, 18],

[TABLE]

where $\Omega(t_{n},\tau)$ satisfies a (highly nonlinear) differential equation. Although different approximations can be found in the literature (see e.g. [11] and references therein), they involve the computation of nested commutators of $A(t)$, evaluated at different times, and so for many problems may be computationally expensive111We must remark however that in some cases, like for the linear Schrödinger equation with time-dependent potential, the evaluation of the commutators can be efficiently carried using appropriate approach [6].

For this reason, here we consider two different strategies starting from the formal solution provided by the Magnus expansion. The first leads us to optimised methods specifically addressed to problems defined by (1) when matrix–matrix products are feasible; to these we call Magnus-decomposition methods. The second strategy allows us to build optimised methods when only matrix–vector products are suitable. These are referred to as Magnus–splitting methods.

In Magnus-decomposition methods, in particular, we approximate $\exp\Omega(t_{n},h)$ by a composition of simpler symplectic maps that for the problem at hand are considerably faster to compute. As an illustration, a 4th-order method within this family contains only one costly matrix exponential and this exponential is cheaply approximated by a symplectic approximation of order $q$ leading to an scheme given by the composition:

[TABLE]

where $D_{i}^{[4,q]},\ i=1,2,3$, are certain matrices to be defined later on.

On the other hand, in Magnus–splitting methods, we consider different decompositions that are suitable for low-oscillatory and high-dimensional problems, where only matrix–vector products are reasonable in the numerical scheme. Methods of order $p$ within this family are of the form

[TABLE]

where the coefficients $a_{i},b_{i}$ are obtained numerically by solving a set of non-linear polynomial equations and the matrices $C_{i}$ are linear combinations of $M(t)$ evaluated at certain quadrature points. The product $\Psi_{m}^{[p]}z_{0}$ is done with only $m$ products of matrices $C_{i}$ on a vector of dimension $r$.

Since all the schemes presented in this work are based on the Magnus expansion, they become explicit symplectic integrators when $M$ is real and symmetric and can be considered as geometric integrators [8, 15, 17, 25], showing a favourable behaviour on long-time simulations.

2 Magnus-based methods

2.1 Magnus expansion

As stated before, the Magnus expansion [18] expresses the solution of  eq. 3 on the time interval [$t_{n}$, $t_{n}+h$] in the form of a single exponential of an infinite series

[TABLE]

whose first terms are given by

[TABLE]

where $[\alpha,\beta]=\alpha\beta-\beta\alpha$ is the matrix commutator. When matrix $A(t)$ is given by (3), the exponent $\Omega$ and as well as any truncation of the series at order $p$, $\Omega^{[p]}$, belong to the symplectic Lie algebra, and thus symplecticity is automatically preserved,

Approximations of $\Omega$ in terms of $A(t)$ evaluated at the nodes of some quadrature rule can be obtained as follows. First we consider the polynomial $\widetilde{A}(t)$ of degree $s-1$ in $t$ that interpolates $A(t)$ on $[t_{n},t_{n}+h]$ at the points $t_{n}+c_{j}h,\ j=1,\ldots,s$, where $c_{i}$ are the nodes of the Gauss–Legendre quadrature rule of order $2s$. The perturbed problem reads

[TABLE]

where $z(t_{n})$ is the exact solution of (3) at $t_{n}$. From a direct application of the Alekseev–Gröbner lemma [27] (see also [16, 21, 23]), we have that

[TABLE]

Letting $t=t_{n}+\frac{h}{2}+\sigma$, we write the interpolation polynomial as

[TABLE]

with $A_{i}=A(t_{n}+c_{i}h)=\widetilde{A}(t_{n}+c_{i}h)$ and the usual Lagrange polynomials ${\cal L}_{i}(t)$. Notice that for our problem (3) we have

[TABLE]

where

[TABLE]

and

[TABLE]

In particular, for sixth-order methods ($s=3$) we have

[TABLE]

where $A_{i}=A(t_{n}+c_{i}h),\ c_{1}=\frac{5-\sqrt{15}}{10},\ c_{2}=\frac{1}{2},\ c_{3}=\frac{5+\sqrt{15}}{10}$. One can readily check that $\alpha_{1}={\cal O}(h),\alpha_{2}={\cal O}(h^{2}),\alpha_{3}={\cal O}(h^{3})$.

The integrals in the Magnus expansion for $\widetilde{A}(t)$ can be computed from (8) analytically. This results in an approximation of $z(t)$ expressed in terms of $\alpha_{1},\>\alpha_{2},\>\alpha_{3}$ up to order $2s$. Specifically, a sixth-order approximation $\Omega^{[6]}=\Omega+\mathcal{O}(h^{7})$ results from

[TABLE]

where $[ij\ldots kl]$ represents the nested commutator $[\alpha_{i},[\alpha_{j},[\ldots,[\alpha_{k},\alpha_{l}]\ldots]]]$.

At this point, we can proceed in different ways and we present two different strategies to build new methods. In both cases, the respective methods (4) and (5) will be compositions of exponentials of elements of the Lie algebra generated by $\alpha_{1},\>\alpha_{2},\>\alpha_{3}$ whose coefficients will be determined in such a way that the compositions coincide with (13) up to the desired order (six in this case).

2.2 Magnus-decomposition integrators

Firstly, we examine the structure of the Lie algebra generated by the $\alpha_{i}$. We immediately notice that, $[ij]=0$ for $i,j>1$. Furthermore,

[TABLE]

We distinguish the following types of exponentials that may appear as elements of the Lie algebra generated by $\alpha_{i},\ i=1,2,3$:

[TABLE]

Clearly, we want to avoid the computation of the full matrix exponential $E_{1}$ and instead focus on types $E_{2},\>E_{3}$. The latter is nilpotent and its exponential comes virtually for free.

In [5], it was obtained methods of order 4, 6 and 8 for the Hill equation, denoted by $\Upsilon_{k}^{[p]}$ (a $p^{th}$-order method containing $k$ exponentials of type $E_{2}$ above). In particular

[TABLE]

[TABLE]

where $C_{i}^{[p]}$ are linear combinations of $M(t)$ evaluated at a set of quadrature points of order $p$ or higher and $D_{i}^{[p]}$ are linear combinations of $M(t)$ that additionally contain one product of such linear combinations. Due to the exponentials, these schemes are specially appropriate when the solution is oscillatory or stiff.

Since the averaged matrices $D_{i}^{[p]}$ can be considered as constant matrices in each time subinterval, the following result is useful when computing the corresponding exponentials for a real-valued matrix $C=-M$ [7, sec. 11.3.3]:

[TABLE]

Although there exist efficient methods to compute matrix trigonometric functions [1, 2] as well as matrix exponentials, such as Padé approximants, Krylov/Lanczos methods, Chebyshev method and others, for this particular case a more efficient procedure is obtained by decomposing the exponential into a product of simple matrices. When $C$ is symmetric, the resulting approximations preserve the symplectic structure by construction. Specifically, if $h\rho(\sqrt{C})<\pi$, where $\rho(\sqrt{C})$ is the spectral radius of $\sqrt{C}$, then [3]

[TABLE]

where

[TABLE]

The truncated series expansions of $Q$ and $R$ up to order $q+2$ and $q$ in $h$, denoted by $Q^{[q+2]}$ and $R^{[q]}$, respectively, can be simultaneously computed with only $k=\left\lfloor\frac{q-1}{2}\right\rfloor$ products (e.g., $Q^{[6]},\>P^{[4]}$ can be computed with only one product). We denote by $\Phi^{[q]}$ a $q^{th}$-order approximation to $\Phi(h)$ obtained by replacing $Q$ and $R$ by $Q^{[q+2]}$ and $R^{[q]}$. Notice that if $C$ is a symmetric matrix then $Q^{[q+2]},\>R^{[q]}$ are also symmetric matrices and, by construction, $\Phi^{[q]}$ is a symplectic matrix $\forall q$. Taking into account all these considerations, we substitute eqs. 18 and 19 into $\Upsilon_{k}^{[p]}$ given by (15) and (16) and combine commuting matrices to get finally $\Upsilon_{1}^{[4,q]},\Upsilon_{2}^{[6,q]}$ given in Table 1. The algorithm proceeds by multiplying $\Upsilon_{k}^{[p,q]}$ by the result from the previous computational step. As the last matrix in a step commutes with the first one in the following step, some products can be saved, hence for an even $q\geq 6$ the computational cost of $\Upsilon^{[4,q]}_{1}$ and $\Upsilon^{[6,q]}_{2}$ are $(1+q/2)\mathcal{C}$ and $(7+q)\mathcal{C}$, respectively. The schemes and the relevant parameters are collected in Table 1.

2.3 Magnus–splitting integrators

For deriving the second class of schemes considered in this work, we first split the matrix $A(t)$ of eq. 3 as

[TABLE]

and denote

[TABLE]

where $\widetilde{B}(s)$ is the interpolating polynomial of $B(s)$ in the interval $[t_{n},t_{n}+h]$, and $\alpha_{1}=\delta_{1}+\beta_{1}$ $\alpha_{i}=\beta_{i},\ i>1$. It is easy to check that $[\beta_{i},\beta_{j}]=0$ and $[\delta_{1},\delta_{1},\delta_{1},\beta_{i}]=[\beta_{i},\beta_{j},\beta_{k},\delta_{1}]=0$ for any value of $i,j,k$. As a consequence, the formal solution eq. 13 simplifies to

[TABLE]

It is then clear that a composition of type (5) can be recovered by considering the following product of exponentials:

[TABLE]

$y_{m+1,j}=0,\ j=1,2,3$. Specifically, taking to account (24), we have

[TABLE]

where $M_{j}=M(t_{n}+c_{j}h),\ j=1,2,3$ and $b_{m+1,j}=0$, so that, in practice, (26) corresponds to a $m$-stage composition. All methods are symplectic when applied to Hamiltonian systems and, moreover, the coefficients are chosen so that time symmetric is preserved.

To obtain particular methods we extend the analysis carried out in [9, 10], where several 6th-order schemes were derived for the more general problem

[TABLE]

Although all of them can be applied to the present problem (3) simply by letting $M(t)=I$, new methods have also been obtained by taking into consideration the simpler structure that this system possesses. Specifically, we have taken $a_{i},\>y_{i,1}$ in (25) as the coefficients of an optimised 11-stage 6th-order method designed in [10] for (27) when $M$ and $N$ are constant. In this way the commutators involving only $\delta_{1}$ and $\beta_{1}$ (e.g. $[\delta_{1},\beta_{1},\delta_{1}]$, $[\beta_{1},\delta_{1},\beta_{1}]$, etc.) vanish up to order six, whereas the higher-order contributions are minimised by considering more stages than strictly necessary to solve the order conditions; this extra cost is compensated by a much improved accuracy and stability in an autonomous case. Next, we look for new coefficients $y_{i,2},\>y_{i,3}$, which now have to satisfy a much reduced set of order conditions. The 11-stage composition has five coefficients $y_{i,2}$ and six coefficients $y_{i,3}$ to solve seven equations. Taking into account the structure of the equations, this leaves one of the $y_{i,2}$ and three of the $y_{i,3}$ as free parameters, that are chosen in order to minimize the objective function $\sum_{i}(|y_{i,2}|^{2}+|y_{i,3}|^{2})$.

Once the coefficients $y_{i,j}$ are chosen, the matrices $\beta_{i},\ i=1,2,3$ are replaced by the corresponding linear combinations of $M(t_{n}+c_{i}h),\ i=1,2,3$ evaluated at the quadrature rule of order 6, so that one ends up with a composition of the form (26). Specifically, the following 11-stage 6th-order method is obtained:

[TABLE]

where

[TABLE]

$M_{j}=M(t_{n}+c_{j}h)$ and $c_{i}$ are the nodes of the 6th-order Gauss–Legendre quadrature rule, i.e.,

[TABLE]

The corresponding coefficients $a_{i},b_{i,j}$ are:

[TABLE]

[TABLE]

verifying the time-symmetry condition $a_{13-i}=a_{i},\quad i=1,\ldots,6$, and

[TABLE]

3 Numerical experiment

In the following section we study and demonstrate the performance of the new methods with respect to well established explicit and implicit standard Runge–Kutta (RK) and explicit symplectic Runge–Kutta–Nyström (RKN) methods from the literature. Six types of methods are considered:

$\Upsilon_{k}^{[p]}$ from [4] and [5]: symplectic $p^{th}$-order methods requiring computing of $k$ matrix exponentials.

• 2.

$\Upsilon_{k}^{[p,q]}$: new sets of methods obtained from $\Upsilon^{[p]}_{k}$ by decomposing matrix exponentials and taking Taylor series expansion up to ${\cal O}\left(h^{q}\right)$.

• 3.

$\Psi_{11}^{[6]}$: the new 11-stage 6th-order Magnus–splitting method (eqs. 28 and 29).

• 4.

$\widehat{\Psi}_{11}^{[6]}$: an 11-stage 6th-order Magnus–splitting method for non-autonomous systems from [9, 10].

• 5.

$\mathrm{RK}_{k}^{[p]}$: an explicit $k$-stage $p^{th}$-order method that uses the Radau quadrature rule of order six for the time-dependent part.

• 6.

$\mathrm{RKGL}_{s}^{[p]}$: implicit $k$-stage $p^{th}$-order symplectic Gauss–Legendre methods.

• 7.

$\mathrm{RKNb}_{k}^{[p]}$: $k$-stage $p^{th}$-order explicit symplectic methods from [12]. $\mathrm{RKNb}_{11}^{[6]}$ has been selected instead of the more effective $\mathrm{RKNa}_{14}^{[6]}$ (ibid.) because the former has the same number of stages as $\Psi_{11}^{[6]}$.

The cost of each method is estimated for two different problem types. The first one is when a numerical method acts on the fundamental matrix $\Phi$. Let $r\times r=\dim M(t)$. Then, the cost is expressed in the number of matrix – matrix products $\mathcal{C}$, required to propagate for one time step $h$. Evaluations of $M(t)$, scalar – matrix multiplications and linear combinations of matrices are not included to the cost.

The second case is when the same method is used to integrate a system whose state is represented by a vector $(v,w)^{T},\dim v=\dim w=r$. Similarly, the cost is expressed in the number of matrix – vector products $\mathcal{V}$. The methods’ costs are summarized in Table 2. The parameter $\varrho^{[p]}$ in the implicit methods refers to the number of iterations per step for a $p^{th}$-order method. Typically, $\varrho^{[p]}=4...7$ to attain convergence and preservation of symplecticity to high accuracy. For the numerical experiments in this paper, we assume them to be $\varrho^{[4]}=4$ for the 4th-order method and $\varrho^{[6]}=6$ for the 6th-order one.

However, some of these methods are not optimised for problems where matrix–matrix are exceedingly costly, and they are not used in these problems. Moreover, the schemes $\Upsilon_{k}^{[p]}$ acting on a vector could be carried out using propagators like Krylov methods to compute the action of the exponential of a matrix on a vector, but this is not considered in this work.

The reference solutions are obtained numerically using sufficiently small time steps.

3.1 Mathieu equation

The first performance test is executed employing the Mathieu equation

[TABLE]

written as a first-order system (3).

At first, we family-wise compare the performance of the new methods $\Upsilon^{[p,q]}_{k}$. Then, we integrate for $t\in[0,\pi]$ with the identity matrix as the initial condition and then measure the $L_{1}$-norm of the error of the fundamental matrix at the final time. This procedure is repeated for different time steps and different choices of $k,p,q$. In Table 3 we show the values of $q$ that provided the best performances for the choices $\varepsilon=\{0.1,\>1\}$ and $\omega=5^{j},\ j=-3,\ldots,3$. We observe that the choices $q=6,\,8$ are generally better in the 4th-order family, and $q=8,\,12$ are better among the 6th-order methods. These selected schemes will be considered in the following numerical examples.

We select now $\varepsilon\in\{0.1,\>1\}$ and $\omega\in\{1/5,\>5\}$ for plotting the $L_{1}$-norm of the error of the fundamental matrix at the final time versus the computation cost in units of $\mathcal{C}$ for all methods previously considered. Figures 2 and 3 show methods’ performance. We observe that the standard explicit and implicit RK methods perform considerably worse while the new methods show the best performances for oscillatory cases.

To illustrate how the accuracy depends on the frequency $\omega$, we take $\varepsilon=1$, the time step $h=\pi/20$ and measure the $L_{1}$-norm of the error in the fundamental matrix solution for $\omega\in{[0,10]}$. The results are shown in Figure 1. We observe that the new methods show smaller error growth as frequency of the problem increases (oscillatory problems).

3.2 Hill equation

The second benchmark to consider is the matrix Hill equation:

[TABLE]

where $A,B_{1},B_{2}\in\mathbb{R}^{r\times r}$. We assume $A=r^{2}I+D$, where $D$ is a Pascal matrix:

[TABLE]

We set $B_{1}=\varepsilon I,B_{2}=\frac{1}{10}\varepsilon I$, $\varepsilon=r$ and $\varepsilon=\frac{1}{10}r$ and compute solutions for $r=5$ and $r=7$ on the interval $t\in[0,\pi]$, and then we measure the $L_{1}$-norm of the error in the fundamental matrix solution at the final time.

In [5], it was determined that for matrix Hill-type problems $\Upsilon_{k}^{[p]}$ perform no worse than $\mathrm{RKNb}_{11}^{[6]}$. Figures 4 and 5 demonstrate that, thanks to lower computational cost, the new $\Upsilon_{1}^{[4,q]},\Upsilon_{2}^{[6,q]}$ and $\Psi_{11}^{[6]}$ methods consistently produce better results in both oscillatory (larger $r$) and nearly autonomous (small $\varepsilon$) cases.

3.3 Time-dependent wave equation

To analyse the performance of the methods that only involve matrix–vector products we consider the following trapped wave equation

[TABLE]

equipped with initial conditions $u(x,0)=\sigma\mathrm{e}^{-x^{2}/2}$, and $u_{t}(x,0)=0$. The solution for $g(x,t)=0$ can be easily be obtained by separation of variables and it is given by $u_{0}(x,t)=\sigma\cos(t)\mathrm{e}^{-x^{2}/2}$.

When an external interaction appears, $g(x,t)\neq 0$, the equation has no analytical solution in general and one has to consider, for example, a numerical scheme. We assume the solution is confined to a region $x\in[x_{0},x_{N}]$ and hence the solution and all spatial derivatives vanish at these boundaries. This allows us to treat the problem as periodic and spectral methods can be used. We divide the spatial region into $N$ intervals of length $\Delta x=(x_{N}-x_{0})/N$ and, after spatial discretization, we obtain an equation similar to eq. 1 that we write as the first order system eq. 3 where $z=(v,w)^{T}$ and $v_{i}(t)\approx u(x_{i},t),w_{i}(t)\approx u_{t}(x_{i},t)$.

For the numerical experiments we take $N=128,x_{0}=-10,x_{N}=10$ and the external interaction

[TABLE]

We take $\delta\in\{0.2,\>1\}$, $\varepsilon\in\{0.1,0.2,\>0.4,\>0.5\}$, and integrate for the period $t\in[0,20\pi/\delta]$. The reference solution is obtained numerically with a sufficiently small time step and we measure the $L_{1}$-norm of the solution versus the number of matrix – vector products for each method. The results are shown in the Figure 6 where $\Psi_{11}^{[6]}$, with only three evaluations of the time-dependent functions per step, shows the best performance for smooth time-dependencies and it is nearly the best one when such time-dependency increases.

4 Conclusions

Starting from the well know Magnus expansion for linear time-dependent differential equations, in this work we have presented two families of time-integrators especially designed for the second order linear system (1). Both share with the exact solution relevant qualitative properties (in particular, its symplectic character) and thus provide an accurate description of the system over long time intervals, due to their favourable error propagation, but they differ in significant aspects. Whereas Magnus-decomposition methods are addressed to problems with small to moderate dimensions where the numerical computation of matrix–matrix products (and thus the fundamental matrix) is feasible, Magnus–splitting methods are advantageous when the dimension of the system is exceedingly high, as is the case when (1) results from a linear time-dependent wave equation previously discretized in space. The numerical experiments reported here clearly illustrate this difference: whereas the first family performs more efficiently on the scalar Mathieu equation and a low-dimensional and oscillatory matrix Hill equation, it is the second class of integrators which shows a better behaviour on a discretised trapped wave equation, where algorithms based only in matrix-vector products are advisable.

Magnus-decomposition methods can be considered as an adaptation of the schemes presented in [5] involving appropriate truncated expansions of the otherwise computationally costly matrix exponentials initially present in them and useful combinations of matrices. Magnus–splitting methods, on the other hand, belong to the class of integrators analysed in [9, 10]: taking as a starting point an efficient splitting method for the autonomous case, the scheme is formed by averaging the matrix $M(t)$ at each step (with different weights). The new coefficients are obtained by solving the additional order conditions arising from the time-dependency. Compared with the general case, the number of order conditions is considerably reduced, which allows to get new schemes within this family in a rather straightforward way.

Acknowledgements

Bader, Blanes, Casas and Kopylov acknowledge the Ministerio de Economía y Competitividad (Spain) for financial support through the coordinated project MTM2013-46553-C3-3-P. Additionally, Kopylov has been partly supported by fellowship GRISOLIA/2015/A/137 from the Generalitat Valenciana.