Search

In this paper, we consider the application of exponential integrators to problems that are advection dominated, either on the entire or on a subset of the domain. In this context, we compare Leja and Krylov based methods to compute the action of exponential and related matrix functions. We set up a performance model by counting the different operations needed to implement the considered algorithms. This model assumes that the evaluation of the right-hand side is memory bound and allows us to evaluate performance in a hardware independent way. We find that exponential integrators perform comparably to explicit Runge-Kutta schemes for problems that are advection dominated in the entire domain. Moreover, they are able to outperform explicit methods in situations where small parts of the domain are diffusion dominated. We generally observe that Leja based methods outperform Krylov iterations in the problems considered. This is in particular true if computing inner products is expensive.

Nonlinear diffusion-reaction systems model a multitude of physical phenomena. A common situation is biological development modeling where such systems have been widely used to study spatiotemporal phenomena in cell biology. Systems of coupled diffusion-reaction equations are usually subject to some complicated features directly related to their multiphysics nature. Moreover, the presence of advection is source of numerical instabilities, in general, and adds another challenge to these systems. In this study, we propose a NURBS-based isogeometric analysis (IgA) combined with a second-order Strang operator splitting to deal with the multiphysics nature of the problem. The advection part is treated in a semi-Lagrangian framework and the resulting diffusion-reaction equations are then solved using an efficient time-stepping algorithm based on operator splitting. The accuracy of the method is studied by means of a advection-diffusion-reaction system with analytical solution. To further examine the performance of the new method on complex geometries, the well-known Schnakenberg-Turing problem is considered with and without advection. Finally, a Gray-Scott system on a circular domain is also presented. The results obtained demonstrate the efficiency of our new algorithm to accurately reproduce the solution in the presence of complex patterns on complex geometries. Moreover, the new method clarifies the effect of geometry on Turing patterns.

The exponential trapezoidal rule is proposed and analyzed for the numerical integration of semilinear integro-differential equations. Although the method is implicit, the numerical solution is easily obtained by standard fixed-point iteration, making its implementation straightforward. Second-order convergence in time is shown in an abstract Hilbert space framework under reasonable assumptions on the problem. Numerical experiments illustrate the proven order of convergence.

We investigate a filtered Lie-Trotter splitting scheme for the ``good" Boussinesq equation and derive an error estimate for initial data with very low regularity. Through the use of discrete Bourgain spaces, our analysis extends to initial data in $H^{s}$ for $01/2$ imposed by the bilinear estimate in smooth Sobolev spaces. We establish convergence rates of order $\tau^{s/2}$ in $L^2$ for such levels of regularity. Our analytical findings are supported by numerical experiments., Comment: 3 figures

The filtered Lie splitting scheme is an established method for the numerical integration of the periodic nonlinear Schr\"{o}dinger equation at low regularity. Its temporal convergence was recently analyzed in a framework of discrete Bourgain spaces in one and two space dimensions for initial data in $H^s$ with $0

For the numerical solution of the cubic nonlinear Schr\"{o}dinger equation with periodic boundary conditions, a pseudospectral method in space combined with a filtered Lie splitting scheme in time is considered. This scheme is shown to converge even for initial data with very low regularity. In particular, for data in $H^s(\mathbb T^2)$, where $s>0$, convergence of order $\mathcal O(\tau^{s/2}+N^{-s})$ is proved in $L^2$. Here $\tau$ denotes the time step size and $N$ the number of Fourier modes considered. The proof of this result is carried out in an abstract framework of discrete Bourgain spaces, the final convergence result, however, is given in $L^2$. The stated convergence behavior is illustrated by several numerical examples., Comment: arXiv admin note: substantial text overlap with arXiv:2301.10639

In this paper we consider an approach to improve the performance of exponential Runge--Kutta integrators and Lawson schemes} in cases where the solution of a related, but usually much simpler, problem can be computed efficiently. While for implicit methods such an approach is common (e.g. by using preconditioners), for exponential integrators this has proven more challenging. Here we propose to extract a constant coefficient differential operator from the semilinear advection-diffusion-reaction equation for which, in many situations, efficient methods are known to compute the required matrix functions. Both a linear stability analysis and {\color{black} extensive} numerical experiments show that the resulting schemes can be unconditionally stable. In fact, we find that exponential integrators of Runge--Kutta type and Lawson schemes can have better stability properties than similarly constructed implicit-explicit schemes. We also derive two new Lawson type integrators that further improve on these stability properties. The overall effectiveness of the approach is highlighted by a number of performance comparisons on examples in two and three space dimensions.

A filtered Lie splitting scheme is proposed for the time integration of the cubic nonlinear Schr\"odinger equation on the two-dimensional torus $\mathbb{T}^2$. The scheme is analyzed in a framework of discrete Bourgain spaces, which allows us to consider initial data with low regularity; more precisely initial data in $H^s(\mathbb{T}^2)$ with $s>0$. In this way, the usual stability restriction to smooth Sobolev spaces with index $s>1$ is overcome. Rates of convergence of order $\tau^{s/2}$ in $L^2(\mathbb{T}^2)$ at this regularity level are proved. Numerical examples illustrate that these convergence results are sharp.

Dynamical low-rank approximation, as has been demonstrated recently, can be extremely efficient in solving kinetic equations. However, a major deficiency is that they do not preserve the structure of the underlying physical problem. For example, the classic dynamical low-rank methods violate mass, momentum, and energy conservation. In [L. Einkemmer, I. Joseph, J. Comput. Phys. 443:110495, 2021] a conservative dynamical low-rank approach has been proposed. However, directly integrating the resulting equations of motion, similar to the classic dynamical low-rank approach, results in an ill-posed scheme. In this work we propose a robust, i.e. well-posed, integrator for the conservative dynamical low-rank approach that conserves mass and momentum (up to machine precision) and significantly improves energy conservation. We also report improved qualitative results for some problems and show how the approach can be combined with a rank adaptive scheme.

The aim of this paper is to construct and analyze explicit exponential Runge-Kutta methods for the temporal discretization of linear and semilinear integro-differential equations. By expanding the errors of the numerical method in terms of the solution, we derive order conditions that form the basis of our error bounds for integro-differential equations. The order conditions are further used for constructing numerical methods. The convergence analysis is performed in a Hilbert space setting, where the smoothing effect of the resolvent family is heavily used. For the linear case, we derive the order conditions for general order $p$ and prove convergence of order $p$, whenever these conditions are satisfied. In the semilinear case, we consider in addition spatial discretization by a spectral Galerkin method, and we require locally Lipschitz continuous nonlinearities. We derive the order conditions for orders one and two, construct methods satisfying these conditions and prove their convergence. Finally, some numerical experiments illustrating our theoretical results are given.

Two types of second-order in time partial differential equations (PDEs), namely semilinear wave equations and semilinear beam equations are considered. To solve these equations with exponential integrators, we present an approach to compute efficiently the action of the matrix exponential as well as those of related matrix functions. Various numerical simulations are presented that illustrate this approach., Comment: 19 pages, 10 figures

In this paper, we analyse a new exponential-type integrator for the nonlinear cubic Schr\"odinger equation on the $d$ dimensional torus $\mathbb T^d$. The scheme has recently also been derived in a wider context of decorated trees in [Y. Bruned and K. Schratz, arXiv:2005.01649]. It is explicit and efficient to implement. Here, we present an alternative derivation, and we give a rigorous error analysis. In particular, we prove second-order convergence in $H^\gamma(\mathbb T^d)$ for initial data in $H^{\gamma+2}(\mathbb T^d)$ for any $\gamma > d/2$. This improves the previous work in [Kn\"oller, A. Ostermann, and K. Schratz, SIAM J. Numer. Anal. 57 (2019), 1967-1986]. The design of the scheme is based on a new method to approximate the nonlinear frequency interaction. This allows us to deal with the complex resonance structure in arbitrary dimensions. Numerical experiments that are in line with the theoretical result complement this work.

For the solution of the cubic nonlinear Schr\"odinger equation in one space dimension, we propose and analyse a fully discrete low-regularity integrator. The scheme is explicit and can easily be implemented using the fast Fourier transform with a complexity of $\mathcal{O}(N\log N)$ operations per time step, where $N$ denotes the degrees of freedom in the spatial discretisation. We prove that the new scheme provides an $\mathcal{O}(\tau^{\frac32\gamma-\frac12-\varepsilon}+N^{-\gamma})$ error bound in $L^2$ for any initial data belonging to $H^\gamma$, $\frac12<\gamma\leq 1$, where $\tau$ denotes the temporal step size. Numerical examples illustrate this convergence behavior.

Recently a splitting approach has been presented for the simulation of sonic-boom propagation. Splitting methods allow one to divide complicated partial differential equations into simpler parts that are solved by specifically tailored numerical schemes. The present work proposes a second order exponential integrator for the numerical solution of sonic-boom propagation modelled through a dispersive equation with Burgers' nonlinearity. The linear terms are efficiently solved in frequency space through FFT, while the nonlinear terms are efficiently solved by a WENO scheme. The numerical method is designed to be highly parallelisable and therefore takes full advantage of modern computer hardware. The new approach also improves the accuracy compared to the splitting method and it reduces oscillations. The enclosed numerical results illustrate that parallelisation on a CPU results in a speedup of 22 times faster than the straightforward sequential version. The GPU implementation further accelerates the runtime by a factor 3, which improves to 5 when single precision is used instead of double precision.

In this paper, we propose a $\mu$-mode integrator for computing the solution of stiff evolution equations. The integrator is based on a $d$-dimensional splitting approach and uses exact (usually precomputed) one-dimensional matrix exponentials. We show that the action of the exponentials, i.e. the corresponding batched matrix-vector products, can be implemented efficiently on modern computer systems. We further explain how $\mu$-mode products can be used to compute spectral transforms efficiently even if no fast transform is available. We illustrate the performance of the new integrator by solving, among the others, three-dimensional linear and nonlinear Schr\"odinger equations, and we show that the $\mu$-mode integrator can significantly outperform numerical methods well established in the field. We also discuss how to efficiently implement this integrator on both multi-core CPUs and GPUs. Finally, the numerical experiments show that using GPUs results in performance improvements between a factor of $10$ and $20$, depending on the problem.

We study a filtered Lie splitting scheme for the cubic nonlinear Schr\"{o}dinger equation. We establish error estimates at low regularity by using discrete Bourgain spaces. This allows us to handle data in $H^s$ with $01/2$ . More precisely, we prove convergence rates of order $\tau^{s/2}$ in $L^2$ at this level of regularity.

The space fractional Cahn-Hilliard phase-field model is more adequate and accurate in the description of the formation and phase change mechanism than the classical Cahn-Hilliard model. In this article, we propose a temporal second-order energy stable scheme for the space fractional Cahn-Hilliard model. The scheme is based on the second-order backward differentiation formula in time and a finite difference method in space. Energy stability and convergence of the scheme are analyzed, and the optimal convergence orders in time and space are illustrated numerically. Note that the coefficient matrix of the scheme is a $2 \times 2$ block matrix with a Toeplitz-like structure in each block. Combining the advantages of this special structure with a Krylov subspace method, a preconditioning technique is designed to solve the system efficiently. Numerical examples are reported to illustrate the performance of the preconditioned iteration., Comment: 26 pages, 5 tables, 4 figures

Fractional Ginzburg-Landau equations as the generalization of the classical one have been used to describe various physical phenomena. In this paper, we propose a numerical integration method for solving space fractional Ginzburg-Landau equations based on a dynamical low-rank approximation. We first approximate the space fractional derivatives by using a fractional centered difference method. Then, the resulting matrix differential equation is split into a stiff linear part and a nonstiff (nonlinear) one. For solving these two subproblems, a dynamical low-rank approach is used. The convergence of our method is proved rigorously. Numerical examples are reported which show that the proposed method is robust and accurate., Comment: 17 pages; 3 figures; 4 tables

Rational exponential integrators (REXI) are a class of numerical methods that are well suited for the time integration of linear partial differential equations with imaginary eigenvalues. Since these methods can be parallelized in time (in addition to the spatial parallelization that is commonly performed) they are well suited to exploit modern high performance computing systems. In this paper, we propose a novel REXI scheme that drastically improves accuracy and efficiency. The chosen approach will also allow us to easily determine how many terms are required in the approximation in order to obtain accurate results. We provide comparative numerical simulations for a shallow water equation that highlight the efficiency of our approach and demonstrate that REXI schemes can be efficiently implemented on graphic processing units., Comment: 25 pages, 5 figures, submitted to Journal of Computational Physics

Volterra subdiffusion problems with weakly singular kernel describe the dynamics of subdiffusion processes well.The graded $L1$ scheme is often chosen to discretize such problems since it can handle the singularity of the solution near $t = 0$. In this paper, we propose a modification. We first split the time interval $[0, T]$ into $[0, T_0]$ and $[T_0, T]$, where $T_0$ ($0 < T_0 < T$) is reasonably small. Then, the graded $L1$ scheme is applied in $[0, T_0]$, while the uniform one is used in $[T_0, T]$. Our all-at-once system is derived based on this strategy. In order to solve the arising system efficiently, we split it into two subproblems and design two preconditioners. Some properties of these two preconditioners are also investigated. Moreover, we extend our method to solve semilinear subdiffusion problems. Numerical results are reported to show the efficiency of our method., Comment: 3 figures; 3 tables. This preprint has been accepted by Journal of Scientific Computing

In this paper, we propose a new scheme for the integration of the periodic nonlinear Schr\"{o}dinger equation and rigorously prove convergence rates at low regularity. The new integrator has decisive advantages over standard schemes at low regularity. In particular, it is able to handle initial data in $H^s$ for $0 < s\le 1$. The key feature of the integrator is its ability to distinguish between low and medium frequencies in the solution and to treat them differently in the discretization. This new approach requires a well-balanced filtering procedure which is carried out in Fourier space. The convergence analysis of the proposed scheme is based on discrete (in time) Bourgain space estimates which we introduce in this paper. A numerical experiment illustrates the superiority of the new integrator over standard schemes for rough initial data., Comment: Revised version with improved main result

The present work proposes a second-order time splitting scheme for a linear dispersive equation with a variable advection coefficient subject to transparent boundary conditions. For its spatial discretization, a dual Petrov--Galerkin method is considered which gives spectral accuracy. The main difficulty in constructing a second-order splitting scheme in such a situation lies in the compatibility condition at the boundaries of the sub-problems. In particular, the presence of an inflow boundary condition in the advection part results in order reduction. To overcome this issue a modified Strang splitting scheme is introduced that retains second-order accuracy. For this numerical scheme a stability analysis is conducted. In addition, numerical results are shown to support the theoretical derivations.

The goal of the present work is to solve a linear dispersive equation with variable coefficient advection on an unbounded domain. In this setting, transparent boundary conditions are vital to allow waves to leave (or even re-enter) the, necessarily finite, computational domain. To obtain an efficient numerical scheme we discretize space using a spectral method. This allows us to drastically reduce the number of grid points required for a given accuracy. Applying a fully implicit time integrator, however, would require us to invert full matrices. This is addressed by performing an operator splitting scheme and only treating the third order differential operator, stemming from the dispersive part, implicitly; this approach can also be interpreted as an implicit-explicit scheme. However, the fact that the transparent boundary conditions are non-homogeneous and depend implicitly on the numerical solution presents a significant obstacle for the splitting/pseudo-spectral approach investigated here. We show how to overcome these difficulties and demonstrate the proposed numerical scheme by performing a number of numerical simulations.

In this paper, we combine deterministic splitting methods with a polynomial chaos expansion method for solving stochastic parabolic evolution problems. The stochastic differential equation is reduced to a system of deterministic equations that we solve explicitly by splitting methods. The method can be applied to a wide class of problems where the related stochastic processes are given uniquely in terms of stochastic polynomials. A comprehensive convergence analysis is provided and numerical experiments validate our approach., Comment: 28 pages, 3 figures, 1 table

We introduce two exponential-type integrators for the "good" Bousinessq equation. They are of orders one and two, respectively, and they require lower regularity of the solution compared to the classical exponential integrators. More precisely, we will prove first-order convergence in Hrfor solutions in H^{r+1} with r > 1/2 for the derived first-order scheme. For the second integrator, we prove second-order convergence in Hrfor solutions in H^{r+3} with r > 1/2 and convergence in L2for solutions in H^3. Numerical results are reported to illustrate the established error estimates. The experiments clearly demonstrate that the new exponential-type integrators are favorable over classical exponential integrators for initial data with low regularity., Comment: 24 pages, 12 figures

We present a new filtered low-regularity Fourier integrator for the cubic nonlinear Schr\"odinger equation based on recent time discretization and filtering techniques. For this new scheme, we perform a rigorous error analysis and establish better convergence rates at low regularity than known for classical schemes in the literature so far. In our error estimates, we combine the better local error properties of the new scheme with a stability analysis based on general discrete Strichartz-type estimates. The latter allow us to handle a much rougher class of solutions as the error analysis can be carried out directly at the level of $L^2$ compared to classical results \black in dimension $d$, \black which are limited to higher-order (sufficiently smooth) Sobolev spaces $H^s$ with $s>d/2$. In particular, we are able to establish a global error estimate in $L^2$ for $H^1$ solutions which is roughly of order $\tau^{ {1\over 2} + { 5-d \over 12} }$ in dimension $d \leq 3$ ($\tau$ denoting the time discretization parameter). This breaks the "natural order barrier" of $\tau^{1/2}$ for $H^1$ solutions which holds for classical numerical schemes (even in combination with suitable filter functions).

The Vlasov--Maxwell equations are used for the kinetic description of magnetized plasmas. As they are posed in an up to 3+3 dimensional phase space, solving this problem is extremely expensive from a computational point of view. In this paper, we exploit the low-rank structure in the solution of the Vlasov equation. More specifically, we consider the Vlasov--Maxwell system and propose a dynamic low-rank integrator. The key idea is to approximate the dynamics of the system by constraining it to a low-rank manifold. This is accomplished by a projection onto the tangent space. There, the dynamics is represented by the low-rank factors, which are determined by solving lower-dimensional partial differential equations. The proposed scheme performs well in numerical experiments and succeeds in capturing the main features of the plasma dynamics. We demonstrate this good behavior for a range of test problems. The coupling of the Vlasov equation with the Maxwell system, however, introduces additional challenges. In particular, the divergence of the electric field resulting from Maxwell's equations is not consistent with the charge density computed from the Vlasov equation. We propose a correction based on Lagrange multipliers which enforces Gauss' law up to machine precision.

We propose an explicit numerical method for the periodic Korteweg-de Vries equation. Our method is based on a Lawson-type exponential integrator for time integration and the Rusanov scheme for Burgers' nonlinearity. We prove first-order convergence in both space and time under a mild Courant-Friedrichs-Lewy condition $\tau=O(h)$, where $\tau$ and $h$ represent the time step and mesh size, respectively. Numerical examples illustrating our convergence result are given., Comment: 17 pages, 4 figures

Standard numerical integrators suffer from an order reduction when applied to nonlinear Schr\"{o}dinger equations with low-regularity initial data. For example, standard Strang splitting requires the boundedness of the solution in $H^{r+4}$ in order to be second-order convergent in $H^r$, i.e., it requires the boundedness of four additional derivatives of the solution. We present a new type of integrator that is based on the variation-of-constants formula and makes use of certain resonance based approximations in Fourier space. The latter can be efficiently evaluated by fast Fourier methods. For second-order convergence, the new integrator requires two additional derivatives of the solution in one space dimension, and three derivatives in higher space dimensions. Numerical examples illustrating our convergence results are included. These examples demonstrate the clear advantage of the Fourier integrator over standard Strang splitting for initial data with low regularity.

We propose a numerical integrator for determining low-rank approximations to solutions of large-scale matrix differential equations. The considered differential equations are semilinear and stiff. Our method consists of first splitting the differential equation into a stiff and a non-stiff part, respectively, and then following a dynamical low-rank approach. We conduct an error analysis of the proposed procedure, which is independent of the stiffness and robust with respect to possibly small singular values in the approximation matrix. Following the proposed method, we show how to obtain low-rank approximations for differential Lyapunov and for differential Riccati equations. Our theory is illustrated by numerical experiments.

Strang splitting is a well established tool for the numerical integration of evolution equations. It allows the application of tailored integrators for different parts of the vector field. However, it is also prone to order reduction in the case of non-trivial boundary conditions. This order reduction can be remedied by correcting the boundary values of the intermediate splitting step. In this paper, three different approaches for constructing such a correction in the case of inhomogeneous Dirichlet, Neumann, and mixed boundary conditions are presented. Numerical examples that illustrate the effectivity and benefits of these corrections are included.

In the present paper we consider numerical methods to solve the discrete Schr\"odinger equation with a time dependent Hamiltonian (motivated by problems encountered in the study of spin systems). We will consider both short-range interactions, which lead to evolution equations involving sparse matrices, and long-range interactions, which lead to dense matrices. Both of these settings show very different computational characteristics. We use Magnus integrators for time integration and employ a framework based on Leja interpolation to compute the resulting action of the matrix exponential. We consider both traditional Magnus integrators (which are extensively used for these types of problems in the literature) as well as the recently developed commutator-free Magnus integrators and implement them on modern CPU and GPU (graphics processing unit) based systems. We find that GPUs can yield a significant speed-up (up to a factor of $10$ in the dense case) for these types of problems. In the sparse case GPUs are only advantageous for large problem sizes and the achieved speed-ups are more modest. In most cases the commutator-free variant is superior but especially on the GPU this advantage is rather small. In fact, none of the advantage of commutator-free methods on GPUs (and on multi-core CPUs) is due to the elimination of commutators. This has important consequences for the design of more efficient numerical methods.

Since their introduction in 1967, Lawson methods have achieved constant interest in the time discretization of evolution equations. The methods were originally devised for the numerical solution of stiff differential equations. Meanwhile, they constitute a well-established class of exponential integrators. The popularity of Lawson methods is in some contrast to the fact that they may have a bad convergence behaviour, since they do not satisfy any of the stiff order conditions. The aim of this paper is to explain this discrepancy. It is shown that non-stiff order conditions together with appropriate regularity assumptions imply high-order convergence of Lawson methods. Note, however, that the term regularity here includes the behaviour of the solution at the boundary. For instance, Lawson methods will behave well in the case of periodic boundary conditions, but they will show a dramatic order reduction for, e.g., Dirichlet boundary conditions. The precise regularity assumptions required for high-order convergence are worked out in this paper and related to the corresponding assumptions for splitting schemes. In contrast to previous work, the analysis is based on expansions of the exact and the numerical solution along the flow of the homogeneous problem. Numerical examples for the Schr\"odinger equation are included.

The efficient numerical integration of large-scale matrix differential equations is a topical problem in numerical analysis and of great importance in many applications. Standard numerical methods applied to such problems require an unduly amount of computing time and memory, in general. Based on a dynamical low-rank approximation of the solution, a new splitting integrator is proposed for a quite general class of stiff matrix differential equations. This class comprises differential Lyapunov and differential Riccati equations that arise from spatial discretizations of partial differential equations. The proposed integrator handles stiffness in an efficient way, and it preserves the symmetry and positive semidefiniteness of solutions of differential Lyapunov equations. Numerical examples that illustrate the benefits of this new method are given. In particular, numerical results for the efficient simulation of the weather phenomenon El Ni\~no are presented.

In this paper we propose a method to solve the Kadomtsev--Petviashvili equation based on splitting the linear part of the equation from the nonlinear part. The linear part is treated using FFTs, while the nonlinear part is approximated using a semi-Lagrangian discontinuous Galerkin approach of arbitrary order. We demonstrate the efficiency and accuracy of the numerical method by providing a range of numerical simulations. In particular, we find that our approach can outperform the numerical methods considered in the literature by up to a factor of five. Although we focus on the Kadomtsev--Petviashvili equation in this paper, the proposed numerical scheme can be extended to a range of related models as well.

In this paper we consider splitting methods in the presence of non-homogeneous boundary conditions. In particular, we consider the corrections that have been described and analyzed in Einkemmer, Ostermann 2015 and Alonso-Mallo, Cano, Reguera 2016. The latter method is extended to the non-linear case, and a rigorous convergence analysis is provided. We perform numerical simulations for diffusion-reaction, advection-reaction, and dispersion-reaction equations in order to evaluate the relative performance of these two corrections. Furthermore, we introduce an extension of both methods to obtain order three locally and evaluate under what circumstances this is beneficial.

We consider Lie and Strang splitting for the time integration of constrained partial differential equations with a nonlinear reaction term. Since such systems are known to be sensitive with respect to perturbations, the splitting procedure seems promising as we can treat the nonlinearity separately. This has some computational advantages, since we only have to solve a linear constrained system and a nonlinear ODE. However, Strang splitting suffers from order reduction which limits its efficiency. This is caused by the fact that the nonlinear subsystem produces inconsistent initial values for the constrained subsystem. The incorporation of an additional correction term resolves this problem without increasing the computational cost. Numerical examples including a coupled mechanical system illustrate the proven convergence results.

The Schr\"odinger equation in the presence of an external electromagnetic field is an important problem in computational quantum mechanics. It also provides a nice example of a differential equation whose flow can be split with benefit into three parts. After presenting a splitting approach for three operators with two of them being unbounded, we exemplarily prove first-order convergence of Lie splitting in this framework. The result is then applied to the magnetic Schr\"odinger equation, which is split into its potential, kinetic and advective parts. The latter requires special treatment in order not to lose the conservation properties of the scheme. We discuss several options. Numerical examples in one, two and three space dimensions show that the method of characteristics coupled with a nonequispaced fast Fourier transform (NFFT) provides a fast and reliable technique for achieving mass conservation at the discrete level.

We introduce low regularity exponential-type integrators for nonlinear Schr\"odinger equations for which first-order convergence only requires the boundedness of one additional derivative of the solution. More precisely, we will prove first-order convergence in $H^r$ for solutions in $H^{r+1}$ ($r>d/2$) of the derived schemes. This allows us lower regularity assumptions on the data in the energy space than for instance required for classical splitting or exponential integration schemes. For one dimensional quadratic Schr\"odinger equations we can even prove first-order convergence without any loss of regularity. Numerical experiments underline the favorable error behavior of the newly introduced exponential-type integrators for low regularity solutions compared to classical splitting and exponential integration schemes.

Splitting methods constitute a well-established class of numerical schemes for the time integration of partial differential equations. Their main advantages over more traditional schemes are computational efficiency and superior geometric properties. In the presence of non-trivial boundary conditions, however, splitting methods usually suffer from order reduction and some additional loss of accuracy. For diffusion-reaction equations with inhomogeneous oblique boundary conditions, a modification of the classic second-order Strang splitting is proposed that successfully resolves the problem of order reduction. The same correction also improves the accuracy of the classic first-order Lie splitting. The proposed modification only depends on the available boundary data and, in the case of non Dirichlet boundary conditions, on the computed numerical solution. Consequently, this modification can be implemented in an efficient way, which makes the modified splitting schemes superior to their classic versions. The framework employed in our error analysis also allows us to explain the fractional orders of convergence that are often encountered for classic Strang splitting. Numerical experiments that illustrate the theory are provided.

The equations of motion of a single particle subject to an arbitrary electric and a static magnetic field form a Poisson system. We present a second-order time integration method which preserves well the Poisson structure and compare it to commonly used algorithms, such as the Boris scheme. All the methods are represented in a general framework of splitting methods. We use the so-called $\phi$ functions, which give efficient ways for both analyzing and implementing the algorithms. Numerical experiments show an excellent long term stability for the new method considered., Comment: 14 pages, 8 figures

The Leja method is a polynomial interpolation procedure that can be used to compute matrix functions. In particular, computing the action of the matrix exponential on a given vector is a typical application. This quantity is required, e.g., in exponential integrators. The Leja method essentially depends on three parameters: the scaling parameter, the location of the interpolation points, and the degree of interpolation. We present here a backward error analysis that allows us to determine these three parameters as a function of the prescribed accuracy. Additional aspects that are required for an efficient and reliable implementation are discussed. Numerical examples that illustrate the performance of our Matlab code are included.

For diffusion-reaction equations employing a splitting procedure is attractive as it reduces the computational demand and facilitates a parallel implementation. Moreover, it opens up the possibility to construct second-order integrators that preserve positivity. However, for boundary conditions that are neither periodic nor of homogeneous Dirichlet type order reduction limits its usefulness. In the situation described the Strang splitting procedure is not more accurate than Lie splitting. In this paper, we propose a splitting procedure that, while retaining all the favorable properties of the original method, does not suffer from order reduction. We demonstrate our results by conducting numerical simulations in one and two space dimensions with inhomogeneous and time dependent Dirichlet boundary conditions. In addition, a mathematical rigorous convergence analysis is conducted that confirms the results observed in the numerical simulations.

We consider a splitting approach for the Kadomtsev--Petviashvili equation with periodic boundary conditions and show that the necessary interpolation procedure can be efficiently implemented. The error made by this numerical scheme is compared to exponential integrators which have been shown in Klein and Roidot (SIAM J. Sci. Comput., 2011) to perform best for stiff solutions of the Kadomtsev--Petviashvili equation. Since many classic high order splitting methods do not perform well, we propose a stable extrapolation method in order to construct an efficient numerical scheme of order four. In addition, the conservation properties and the possibility of order reduction for certain initial values for the numerical schemes under consideration is investigated.

In this paper we study the error propagation of numerical schemes for the advection equation in the case where high precision is desired. The numerical methods considered are based on the fast Fourier transform, polynomial interpolation (semi-Lagrangian methods using a Lagrange or spline interpolation), and a discontinuous Galerkin semi-Lagrangian approach (which is conservative and has to store more than a single value per cell). We demonstrate, by carrying out numerical experiments, that the worst case error estimates given in the literature provide a good explanation for the error propagation of the interpolation-based semi-Lagrangian methods. For the discontinuous Galerkin semi-Lagrangian method, however, we find that the characteristic property of semi-Lagrangian error estimates (namely the fact that the error increases proportionally to the number of time steps) is not observed. We provide an explanation for this behavior and conduct numerical simulations that corroborate the different qualitative features of the error in the two respective types of semi-Lagrangian methods. The method based on the fast Fourier transform is exact but, due to round-off errors, susceptible to a linear increase of the error in the number of time steps. We show how to modify the Cooley--Tukey algorithm in order to obtain an error growth that is proportional to the square root of the number of time steps. Finally, we show, for a simple model, that our conclusions hold true if the advection solver is used as part of a splitting scheme., Comment: submitted to Computers & Mathematics with Applications

In this paper we propose a strategy to suppress the recurrence effect present in grid-based Vlasov solvers. This method is formulated by introducing a cutoff frequency in Fourier space. Since this cutoff only has to be performed after a number of time steps, the scheme can be implemented efficiently and can relatively easily be incorporated into existing Vlasov solvers. Furthermore, the scheme proposed retains the advantage of grid-based methods in that high accuracy can be achieved. This is due to the fact that in contrast to the scheme proposed by Abbasi et al. no statistical noise is introduced into the simulation. We will illustrate the utility of the method proposed by performing a number of numerical simulations, including the plasma echo phenomenon, using a discontinuous Galerkin approximation in space and a Strang splitting based time integration.

In this paper we consider splitting methods for the time integration of parabolic and certain classes of hyperbolic partial differential equations, where one partial flow can not be computed exactly. Instead, we use a numerical approximation based on the linearization of the vector field. This is of interest in applications as it allows us to apply splitting methods to a wider class of problems from the sciences. However, in the situation described the classic Strang splitting scheme, while still a method of second order, is not longer symmetric. This, in turn, implies that the construction of higher order methods by composition is limited to order three only. To remedy this situation, based on previous work in the context of ordinary differential equations, we construct a class of Strang splitting schemes that are symmetric up to a desired order. We show rigorously that, under suitable assumptions on the nonlinearity, these methods are of second order and can then be used to construct higher order methods by composition. In addition, we illustrate the theoretical results by conducting numerical experiments for the Brusselator system and the KdV equation.

Exponential Runge-Kutta methods constitute efficient integrators for semilinear stiff problems. So far, however, explicit exponential Runge-Kutta methods are available in the literature up to order 4 only. The aim of this paper is to construct a fifth-order method. For this purpose, we make use of a novel approach to derive the stiff order conditions for high-order exponential methods. This allows us to obtain the conditions for a method of order 5 in an elegant way. After stating the conditions, we first show that there does not exist an explicit exponential Runge-Kutta method of order 5 with less than or equal to 6 stages. Then, we construct a fifth-order method with 8 stages and prove its convergence for semilinear parabolic problems. Finally, a numerical example is given that illustrates our convergence bound., Comment: Revised version 01.07.2013

In this paper we consider splitting methods for nonlinear ordinary differential equations in which one of the (partial) flows that results from the splitting procedure can not be computed exactly. Instead, we insert a well-chosen state $y_{\star}$ into the corresponding nonlinearity $b(y)y$, which results in a linear term $b(y_{\star})y$ whose exact flow can be determined efficiently. Therefore, in the spirit of splitting methods, it is still possible for the numerical simulation to satisfy certain properties of the exact flow. However, Strang splitting is no longer symmetric (even though it is still a second order method) and thus high order composition methods are not easily attainable. We will show that an iterated Strang splitting scheme can be constructed which yields a method that is symmetric up to a given order. This method can then be used to attain high order composition schemes. We will illustrate our theoretical results, up to order six, by conducting numerical experiments for a charged particle in an inhomogeneous electric field, a post-Newtonian computation in celestial mechanics, and a nonlinear population model and show that the methods constructed yield superior efficiency as compared to Strang splitting. For the first example we also perform a comparison with the standard fourth order Runge--Kutta methods and find significant gains in efficiency as well better conservation properties.

We analyze the convergence of the exponential Lie and exponential Strang splitting applied to inhomogeneous second-order parabolic equations with Dirichlet boundary conditions. A recent result on the smoothing properties of these methods allows us to prove sharp convergence results in the case of homogeneous Dirichlet boundary conditions. When no source term is present and natural regularity assumptions are imposed on the initial value, we show full-order convergence of both methods. For inhomogeneous equations, we prove full-order convergence for the exponential Lie splitting, whereas order reduction to 1.25 for the exponential Strang splitting. Furthermore, we give sufficient conditions on the inhomogeneity for full-order convergence of both methods. Moreover our theoretical convergence results explain the severe order reduction to 0.25 of splitting methods applied to problems involving inhomogeneous Dirichlet boundary conditions. We include numerical experiments to underline the sharpness of our theoretical convergence results.