Project A7 • Numerical methods for highly oscillatory problems

Simulating phenomena like vibration and light is an omnipresent task in modern technology and research. From physics we have for example the semilinear wave equation \begin{equation} \label{eq:wave.nonlinear} q''(t) = -L q(t) + g(q(t)), \qquad q(0) = q_0, \quad q'(0) = q'_0, \end{equation} on a spatial domain $\Omega \subset \mathbb{R}^d$. The operator $L$ is positive and self-adjoint, and $g$ is a sufficiently regular nonlinearity. Such solutions can be approximated by many well established numerical algorithms. For the analysis and construction the behaviour of the wave, so the solution u, needs to be somewhat "nice", which means for example that higher derivatives of the solution have to exist and are suitably bounded. From an analytic point of view we often cannot guarantee such regularity. Also numerical experiments show that for non smooth solutions resonances in the error can be observed if the step size times a high frequency of the system happens to hit a multiple of $2\pi$.

In [GH06], [HLW06] several trigonometric integrators are proposed to solve such problems. They also note that those can be understood as a splitting method for a filtered equation, where the nonlinearity $g$ is modified by filters depending on the operator $L$ and the step size $\tau$. To gain a better insight, in [BGGHJ18] the analysis was done from this point of view. We considered a spatially discrete and linear version of \eqref{eq:wave.nonlinear} and integrated it using the Strang splitting applied to the modified equation. In [BDH21] and [Doe21], we generalized this approach to the abstract equation \eqref{eq:wave.nonlinear} and considered a large class of exponential integrators (including the Strang splitting) as well as different scales of regularity and boundary conditions. We assume the abstract finite-energy condition \begin{equation*} \|L^{1/2} q(t)\|^2 + \|q'(t)\|^2 \leq K^2, \qquad t\geq 0, \end{equation*} and establish error bounds of the order $\mathcal{O}(\tau^2)$. In a first step, one estimates the error obtained by the perturbation through the filter. In a second step, the error of the exponential integrator is analyzed. Because of the filters some critical terms can be handled and one can partly use some classical techniques. Arguments like the Lady Windermere's fan, though, need to be carefully adapted to this special situation using for example summation by parts. We point out that in contrast to previous work, we are now able to equip almost any second-order exponential method with appropriate filter functions in such a way that in spite of the low regularity we can prove second-order error bounds.

Our contributions: [BGGHJ18], [BDH21], [Doe21].

Under certain configurations of laser wavelength and plasma density, the laser-pulse becomes more focused during and after its propagation through the plasma barrier, cf. [KSHLS06, KSHS08]

It can be shown that the exact solution to the problem is highly oscillatory, which imposes difficulties in the numerical approximation. Standard schemes like the leapfrog method require a fine resolution in space in order to resolve the oscillations. This results in a matrix differential equation (MDE) of large dimension, \begin{equation*} A''(t) = D_y A(t) + A(t) D_x + f\big(A(t)\big), \qquad D_y \in \mathbb{C}^{m \times m}, \quad D_x \in \mathbb{C}^{n \times n}, \quad A(t) \in \mathbb{C}^{m \times n}, \end{equation*} hence the memory footprint of the simulation is often huge. Gautschi-type schemes, which are tailored to highly oscillatory second-order equations, often also have a large memory footprint and suffer from long runtimes.

A key observation in the laser-plasma interaction is that for all times the solution is localized in space, i.e., only a small part of the computational domain is not covered by zero. Thus, the approximation of the solution by a time-dependent matrix of low-rank seams feasible. In this situation, the numerical approximation can be computed based on a decomposition into low-rank matrices $U$, $S$, and $V$ like $$A = U S V^H, \quad U \in \mathbb{C}^{m \times r}, \quad S \in \mathbb{C}^{r \times r}, \quad V \in \mathbb{C}^{n \times r},$$ with $U,V$ being orthonormal. Particularly in the case, where the rank $r$ is very small compared to the dimension $m,n$ of the problem, the numerical integration based on the factors is significantly faster and also requires only a fraction of the original memory footprint.

For first-order MDEs, dynamical low-rank integrators (DLRIs) have been introduced in [KL07] and [LO14]. They have been successfully used to efficiently compute approximations to the solutions of a large number of differential equations, see also project B9. However, for second-order MDEs no DLRI has been proposed yet.

In [HNS23a] we constructed and analyzed a novel DLRI for such equations based on the ideas presented in [LO14] and the leapfrog scheme. The resulting integrator was named the LRLF (low-rank leapfrog) scheme. The close relation between the LRLF and leapfrog methods are also visible in the analysis, as we could show that the error of the LRLF scheme is a combination of a low-rank error and the time-discretization error of the leapfrog integrator. In absence of low-rank errors, both schemes compute the same approximations.

For the laser-plasma interaction, already with rank 4 the new DLRI yields results of the same quality as the leapfrog scheme, with significantly reduced computational effort and storage requirements. Moreover, we derived a variant of the LRLF scheme which is particularly suited for semilinear equations. It is based on a splitting ansatz, where the linear parts of the right-hand side are separated from the nonlinear part. The nonlinear part is approximately solved with the LRLF method. For the linear parts, we compute low-rank approximations by their exact solutions and the integrator proposed in [LO14]. The new scheme yields smaller error constants than the LRLF scheme, with reasonable additional effort.

During the project, it turned out that rank-adaptivity is crucial for an efficient low-rank integration. We therefore designed an adaptive strategy based on the singular values of the numerical solution and estimators for the time-discretization error. By balancing the low-rank error contributions suitably with the time-discretization error we managed to control the rank adaptively while maintaining the desirable temporal convergence. The strategy can also be adapted to first-order problems and was introduced in detail in [HNS22b].

Further examples and more new DLRIs for second-order MDEs can be found in the doctoral thesis by S. Schrammer, cf. [Sch22].

Our contributions: [HNS23a], [HNS22b], [Sch22].

To overcome these difficulties, we considered the variational approximation by Gaussian wave packets \begin{equation*} u(\cdot,x) = \exp{\Bigl(\frac{\mathrm{i}}{\varepsilon}\bigl(\frac12(x-q)^T C (x-q)+(x-q)^T p +\zeta\bigr)\Bigr) } , \end{equation*} with time-dependent parameters $q,p \in \mathbb{R}^{d}$, $C \in \mathbb{C}^{d\times d}$ symmetric with positive definite imaginary part, and $\zeta \in\mathbb{C}$, as proposed in [LL20] for the non-magnetic case.

We were able to transfer $L^2$-error bounds and error bounds for physical quantities from the non-magnetic to the magnetic case, cf. [BDHL23]. Moreover, for the numerical approximation, we derived a set of ordinary differential equations - equations of motion - for the parameters of the approximating Gaussian wave packet.

Our contribution: [BDHL23].

This is because we then have a multiscale problem with different behavior on different time scales: fast gyrorotation, guiding centre motion, slow perpendicular drift, near-conservation of the magnetic moment over very long times and conservation of energy for all times. Using modulated Fourier expansions of the exact and numerical solutions obtained with standard and newly developed numerical methods, it is analysed to which extent this behavior is reproduced by numerical integrators used with large step sizes that do not resolve the fast gyrorotations.

Our contributions: [HL20, HLS22, HLW20].

The matrices $A_1, \ldots, A_d \in \mathbb{R}^{n\times n}$ are symmetric, whereas $E\in\mathbb{R}^{n\times n} $ is skew-symmetric ($d,n\in\mathbb{N}$). Typical initial data are highly oscillatory functions of the form \begin{align} \label{A7:PDE.uu.initial.data} u(0,x) &= \left(\mathrm{e}^{\mathrm{i} (\kappa \cdot x)/\varepsilon} + \mathrm{e}^{-\mathrm{i} (\kappa \cdot x)/\varepsilon}\right) p(x) = 2\cos\left(\tfrac{\kappa \cdot x}{\varepsilon}\right)p(x), \end{align} where $\kappa \in \mathbb{R}^d\setminus\{0\}$ is a given wave vector and $p \colon \mathbb{R}^d \rightarrow \mathbb{R}^n$ is assumed to be sufficiently smooth. See Figure 3 for an illustrative example.

Solving \eqref{A7:PDE.uu} numerically is a challenging task. Physically relevant solutions of \eqref{A7:PDE.uu} oscillate rapidly in time and space due to the small parameter $\varepsilon$ which occurs both in the PDE \eqref{A7:PDE.uu} and in the initial data \eqref{A7:PDE.uu.initial.data}. Moreover, the problem is scaled in such a way that the physically relevant nonlinear and diffractive effects appear only on a long time interval $[0,t_{\text{end}}/\varepsilon]$. Hence, a standard space discretization would require an extremely fine resolution with a tremendous number of unknowns, and standard time-stepping methods could only be used with a tiny step size $\tau < \varepsilon$ and a huge number of time steps. Clearly, the total cost of such a method would be way too high to be practicable. The goal was to approximate solutions of such problems analytically.

For every $\alpha\in\mathbb{R}$ and $\beta\in\mathbb{R}^d$ the matrix \begin{align*} \mathcal{L}(\alpha,\beta) = -\alpha I + A(\beta)-\mathrm{i} E \in \mathbb{C}^{n\times n} \qquad \text{with} \qquad A(\beta)&=\sum_{\ell=1}^d \beta_\ell A_\ell \end{align*} is Hermitian. We choose $\omega=\omega(\kappa)$ to be an eigenvalue of $A(\kappa)-\mathrm{i} E$, such that $\mathcal{L}(\omega,\kappa)$ has a nontrivial kernel. The orthogonal projection from $\mathbb{C}^{n\times n}$ onto the kernel of $\mathcal{L}(\omega,\kappa)$ is denoted by $\mathcal{P}$, and we set $\mathcal{P}^\perp=I-\mathcal{P}$. We impose the polarization condition $ \mathcal{P} p = p $, which is a typical assumption in applications. Our goal is to derive an approximation of the form \begin{align} \label{A7:Ansatz} u(t,x) \; \approx \; \widetilde{u}^{(m)}(t,x) = \sum_{j\in \mathcal{J}(m)} \mathrm{e}^{\mathrm{i} j(\kappa \cdot x - \omega t)/\varepsilon} u_j(t,x), \qquad u_{-j} = \overline{u_j} \end{align} where $ \mathcal{J}(m) = \{\pm 1, \pm 3, \ldots, \pm m\} $ for some odd integer $m\in\mathbb{N}$. Substituting this ansatz into \eqref{A7:PDE.uu} and discarding certain terms yields a system of PDEs for the coefficient functions $u_j$. This ansatz is an extension of the celebrated slowly varying envelope approximation (see, e.g., [CL09]) which is obtained for $m=1$. The main advantage of this approach is that the coefficient functions $u_j$ are free from $\varepsilon$-induced oscillations in space, because the spatial oscillations are already incorporated in the ansatz \eqref{A7:Ansatz}.

For $m=3$ we proved in [BJ22] that under certain assumptions the error bound \begin{align} \label{A7:error.bound.m3a} \sup_{t\in[0,t_{\text{end}}/\varepsilon]} \| u(t) - \widetilde{u}^{(3)}(t) \|_{L^\infty} &\leq C\varepsilon^2, \end{align} holds with a constant which does not depend on $ \varepsilon $. This result is an extension of the estimate \begin{align} \label{A7:error.bound.m1a} \sup_{t\in[0,t_{\text{end}}/\varepsilon]} \| u(t) - \widetilde{u}^{(1)}(t) \|_{L^\infty} &\leq C\varepsilon \end{align} for $m=1$ shown by Colin and Lannes in [CL09]. In [Bau22] a bootstrapping argument is used to improve \eqref{A7:error.bound.m1a} and \eqref{A7:error.bound.m3a} to \begin{align*} \sup_{t\in[0,t_{\text{end}}/\varepsilon]} \| u(t) - \widetilde{u}^{(1)}(t) \|_{L^\infty} \leq C\varepsilon^2 \qquad \text{and} \qquad \sup_{t\in[0,t_{\text{end}}/\varepsilon]} \| u(t) - \widetilde{u}^{(3)}(t) \|_{L^\infty} \leq C\varepsilon^4 \end{align*} under slightly stronger regularity assumptions on the initial data. In fact, these bounds correspond to the behavior of approximations observed in our numerical examples (see Figure 4).

As a related but complementary line of research, we constructed solutions to \eqref{A7:PDE.uu} with polarized initial values by modulated Fourier expansions. In [BJL22] we showed that under suitable assumptions there is a solution $u$ on the interval $[0,t_{\text{end}}/\varepsilon]$ that admits a modulated Fourier expansion \begin{equation*} u(t,x) = \sum_{j\in \mathcal{J}(m)} \sum_{\ell=j-1}^{m} \varepsilon^{\ell} u_j^{\ell}(t,x) \, \mathrm{e}^{\mathrm{i} j (\kappa\cdot x - \omega t)/\varepsilon} + \mathcal{O}(\varepsilon^{m+1}). \end{equation*} The coefficient functions $ u_j^{\ell} $ can be computed recursively: for all $\ell$ the parts $\mathcal{P} u_1^{\ell} \in \text{ker}(\mathcal{L}(\omega,\kappa))$ are obtained from solving Schrödinger-type equations, whereas the parts $\mathcal{P}^\perp u_1^{\ell} $ and $ u_j^{\ell} $ with $ 3 \leq j \in \mathcal{J}(m) $ are determined by algebraic equations. All equations are non-oscillatory, which is a crucial advantage for numerical approximations. The initial value $ u(0,x) = 2 \cos(\tfrac{\kappa \cdot x}{\varepsilon}) p(x) + O(\varepsilon) $ is determined by $p$ uniquely up to $O(\varepsilon^{m+1})$ via the algebraic equations. Hence, the approximation result holds for special initial values $u(0,\cdot)$ of \eqref{A7:PDE.uu}, which we call polarized initial values. This is in contrast to \eqref{A7:Ansatz} where no such polarization is required, but with the consequence that the coefficient functions $ u_j $ still oscillate in time.

Our contributions: [BJ22], [Bau22], [BJL22].

It was shown in [CW20], however, that on bounded time-intervals the vector-valued solution $\psi^\varepsilon(t,x)\in\mathbb{C}^4$ has the representation \begin{align} \label{SNRLS} \psi^\varepsilon(t,x) = e^{-it/\varepsilon^2} \varphi_+^\varepsilon(t,x) + e^{it/\varepsilon^2} \varphi_-^\varepsilon(t,x) + \mathcal{O}(\varepsilon^2), \end{align} where $ \varphi_{\pm}^\varepsilon $ are the solutions of two coupled semilinear PDEs called the semi-nonrelativistic limit system. Although $ \varphi_{\pm}^\varepsilon $ do still depend on $\varepsilon$, these functions are smooth (non-oscillatory). Hence, \eqref{SNRLS} offers a way to approximate the solution of the nonlinear Dirac equation up to an error of $\mathcal{O}(\varepsilon^2)$ without having to solve a highly oscillatory problem. In [JK22] we have constructed and analyzed a non-standard exponential two-step method for the semi-nonrelativistic limit system called the exponential explicit midpoint rule (EEMR). It is second-order accurate in the time step size and time-symmetric. Compared to standard second-order exponential integrators, it is easier to implement and considerably more efficient. Numerical experiments show that the new integrator reduces the computational costs per time step to 40% and within a given runtime improves the accuracy by a factor of six. These factors can not only be observed in numerical experiments, but also justified from the construction of the methods.

Our contribution: [JK22].

The coefficient function \(\gamma\) is given by \[\begin{equation}\label{gamma} \gamma(t) = \chi(t)+\varepsilon\alpha\,,\end{equation}\] where \(\alpha\gt0\) is the mean dispersion, and where \[\begin{equation}\label{chi} \chi(t) = \left\{\begin{array}{ll}-\delta\,, \text{ if } t \in [m, m+1) \text{ for even } m \in \mathbb{N}, \\[2mm] \,\,\,\,\delta\,, \text{ if } t \in [m, m+1) \text{ for odd }\, m \in \mathbb{N} \end{array}\right.\end{equation}\] is a periodic, piecewise constant function which alternates between \(-\delta\) and \(+\delta\). A short introduction to dispersion management from a mathematical perspective can be found in the introduction of [JM18] and in the dissertation [Mik17].

Approximating the solution of \eqref{DMNLS} numerically is a challenging task. The factor \(\mathrm{i}/\varepsilon\) in the right-hand side causes a highly oscillatory behavior of typical solutions. As a consequencs, traditional numerical time-integrators (e.g. Runge-Kutta or multistep methods) yield a very poor accuracy unless a huge number of time-steps with a tiny step-size \(\tau\ll\varepsilon\) is made, which causes excessive numerical costs. Additional difficulties are caused by the fact that the time-derivative \(t \mapsto \partial_t u(t,\cdot)\) discontinuous due to the piecewise constant coefficient function \(\gamma\). This is in contradiction to standard assumptions which are typically made to prove higher order convergence of many time integrators. Last but not least, the nonlinearity \(\mathrm{i} |u(t,x)|^2 u(t,x)\) makes implicit methods prohibitively costly and complicates the construction of new integrators.

In [JM18] we have devised and analyzed a novel time-integrator for the dispersion-managed nonlinear Schrödinger equation \eqref{DMNLS} on the torus \(\mathbb{T}=\mathbb{R}/2\pi\mathbb{Z}\). This new method, called the adiabatic midpoint rule, is based on a transformation of \eqref{DMNLS} to an evolution equation which has several advantageous properties, and on the fact that certain highly oscillatory exponential functions which occur in the transformed equation can be integrated analytically. This approach is related to the ideas developed in [Jah04, JL03]. It was shown that under suitable regularity assumptions the adiabatic midpoint rule converges with order 1 with an error constant which does not depend on \(\varepsilon\) and without any \(\varepsilon\)-induced step-size restriction. Moreover, we proved that the accuracy increases if the step-size \(\tau\) is chosen in a special way: The global error reduces to \(\mathcal{O}(\varepsilon^2 + \tau^2)\) if \(\tau = \varepsilon k \) with \(k\in \mathbb{N}\), and to \(\mathcal{O}(\varepsilon\tau)\) if \(\tau = \varepsilon/k\), respectively. In both cases the error constant remains independent on \(\varepsilon\), which explains the excellent accuracy of the adiabatic midpoint rule in numerical experiments; cf. Figure 1.

Unfortunately, the \(L_2\) norm of numerical solutions provided by the adiabatic midpoint rule is in general not constant, although it can easily be shown that the norm \(\|u(t,\cdot)\|_{L_2}\) of every solution of \eqref{DMNLS} is preserved. This was our motivation to construct an exponential version of the adiabatic midpoint rule which conserves the norm of the numerical solution exactly and, in addition, has the same favourable convergence behaviour as its non-exponential counterpart. The exponential adiabatic midpoint rule presented in [JM19] is closely related to our approach in [JM18], but in the error analysis we had to cope with new challenges arising from the exponential structure of the new method.

This part of the project was successfully completed in the first funding period (2015-2019).

Our contributions: [JM19], [JM18], [Mik17].

Continue reading. Collapse content.