Search

New families of composition methods with processing of order 4 and 6 are presented and analyzed. They are specifically designed to be used for the numerical integration of differential equations whose vector field is separated into three or more parts which are explicitly solvable. The new schemes are shown to be more efficient than previous state-of-the-art splitting methods.

We propose a unified approach for different exponential perturbation techniques used in the treatment of time-dependent quantum mechanical problems, namely the Magnus expansion, the Floquet--Magnus expansion for periodic systems, the quantum averaging technique and the Lie--Deprit perturbative algorithms. Even the standard perturbation theory fits in this framework. The approach is based on carrying out an appropriate change of coordinates (or picture) in each case, and can be formulated for any time-dependent linear system of ordinary differential equations. All the procedures (except the standard perturbation theory) lead to approximate solutions preserving by construction unitarity when applied to the time-dependent Schr\"odinger equation.

We propose a framework where Fer and Wilcox expansions for the solution of differential equations are derived from two particular choices for the initial transformation that seeds the product expansion. In this scheme intermediate expansions can also be envisaged. Recurrence formulas are developed. A new lower bound for the convergence of the Wilcox expansion is provided as well as some applications of the results. In particular, two examples are worked out up to high order of approximation to illustrate the behavior of the Wilcox expansion.

The present work provides a comprehensive study of symmetric-conjugate operator splitting methods in the context of linear parabolic problems and demonstrates their additional benefits compared to symmetric splitting methods. Relevant applications include nonreversible systems and ground state computations for linear Schr\"odinger equations based on the imaginary time propagation. Numerical examples confirm the favourable error behaviour of higher-order symmetric-conjugate splitting methods and illustrate the usefulness of a time stepsize control, where the local error estimation relies on the computation of the imaginary parts and thus requires negligible costs., Comment: Paper to be published in Journal of Computational Dynamics

This overview is devoted to splitting methods, a class of numerical integrators intended for differential equations that can be subdivided into different problems easier to solve than the original system. Closely connected with this class of integrators are composition methods, in which one or several low-order schemes are composed to construct higher-order numerical approximations to the exact solution. We analyze in detail the order conditions that have to be satisfied by these classes of methods to achieve a given order, and provide some insight about their qualitative properties in connection with geometric numerical integration and the treatment of highly oscillatory problems. Since splitting methods have received considerable attention in the realm of partial differential equations, we also cover this subject in the present survey, with special attention to parabolic equations and their problems. An exhaustive list of methods of different orders is collected and tested on simple examples. Finally, some applications of splitting methods in different areas, ranging from celestial mechanics to statistics, are also provided., Comment: Review paper to be published in Acta Numerica 2024

We propose new linear combinations of compositions of a basic second-order scheme with appropriately chosen coefficients to construct higher order numerical integrators for differential equations. They can be considered as a generalization of extrapolation methods and multi-product expansions. A general analysis is provided and new methods up to order 8 are built and tested. The new approach is shown to reduce the latency problem when implemented in a parallel environment and leads to schemes that are significantly more efficient than standard extrapolation when the linear combination is delayed by a number of steps., Comment: 17 figures

The present work is concerned with the extension of modified potential operator splitting methods to specific classes of nonlinear evolution equations. The considered partial differential equations of Schr{\"o}dinger and parabolic type comprise the Laplacian, a potential acting as multiplication operator, and a cubic nonlinearity. Moreover, an invariance principle is deduced that has a significant impact on the efficient realisation of the resulting modified operator splitting methods for the Schr{\"o}dinger case.} Numerical illustrations for the time-dependent Gross--Pitaevskii equation in the physically most relevant case of three space dimensions and for its parabolic counterpart related to ground state and excited state computations confirm the benefits of the proposed fourth-order modified operator splitting method in comparison with standard splitting methods. The presented results are novel and of particular interest from both, a theoretical perspective to inspire future investigations of modified operator splitting methods for other classes of nonlinear evolution equations and a practical perspective to advance the reliable and efficient simulation of Gross--Pitaevskii systems in real and imaginary time., Comment: 30 pages, 6 figures

We analyze the preservation properties of a family of reversible splitting methods when they are applied to the numerical time integration of linear differential equations defined in the unitary group. The schemes involve complex coefficients and are conjugated to unitary transformations for sufficiently small values of the time step-size. New and efficient methods up to order six are constructed and tested on the linear Schr\"odinger equation., Comment: 24 pages, 6 figures

For systems of the form $\dot q = M^{-1} p$, $\dot p = -Aq+f(q)$, common in many applications, we analyze splitting integrators based on the (linear/nonlinear) split systems $\dot q = M^{-1} p$, $\dot p = -Aq$ and $\dot q = 0$, $\dot p = f(q)$. We show that the well-known Strang splitting is optimally stable in the sense that, when applied to a relevant model problem, it has a larger stability region than alternative integrators. This generalizes a well-known property of the common St\"{o}rmer/Verlet/leapfrog algorithm, which of course arises from Strang splitting based on the (kinetic/potential) split systems $\dot q = M^{-1} p$, $\dot p = 0$ and $\dot q = 0$, $\dot p = -Aq+f(q)$., Comment: 2 figures

A procedure to numerically integrate non-autonomous linear delay differential equations is presented. It is based on the use of an spectral discretization of the delayed part to transform the original problem into a matrix linear ordinary differential equation which is subsequently solved with numerical integrators obtained from the Magnus expansion. The algorithm can be used in the periodic case to get both accurate approximations of the characteristic multipliers and the solution itself. In addition, it can be extended to deal with certain quasilinear delay equations.

We study Hamiltonian Monte Carlo (HMC) samplers based on splitting the Hamiltonian $H$ as $H_0(\theta,p)+U_1(\theta)$, where $H_0$ is quadratic and $U_1$ small. We show that, in general, such samplers suffer from stepsize stability restrictions similar to those of algorithms based on the standard leapfrog integrator. The restrictions may be circumvented by preconditioning the dynamics. Numerical experiments show that, when the $H_0(\theta,p)+U_1(\theta)$ splitting is combined with preconditioning, it is possible to construct samplers far more efficient than standard leapfrog HMC., Comment: 25 pages, 6 figures

A new family of methods involving complex coefficients for the numerical integration of differential equations is presented and analyzed. They are constructed as linear combinations of symmetric-conjugate compositions obtained from a basic time-symmetric integrator of order 2n (n $\ge$ 1). The new integrators are of order 2(n + k), k = 1, 2, ..., and preserve time-symmetry up to order 4n + 3 when applied to differential equations with real vector fields. If in addition the system is Hamiltonian and the basic scheme is symplectic, then they also preserve symplecticity up to order 4n + 3. We show that these integrators are well suited for a parallel implementation, thus improving their efficiency. Methods up to order 10 based on a 4th-order integrator are built and tested in comparison with other standard procedures to increase the order of a basic scheme., Comment: Accepted for publication in Applied Mathematics and Computation

We present a practical algorithm to approximate the exponential of skew-Hermitian matrices up to round-off error based on an efficient computation of Chebyshev polynomials of matrices and the corresponding error analysis. It is based on Chebyshev polynomials of degrees 2, 4, 8, 12 and 18 which are computed with only 1, 2, 3, 4 and 5 matrix-matrix products, respectively. For problems of the form $\exp(-iA)$, with $A$ a real and symmetric matrix, an improved version is presented that computes the sine and cosine of $A$ with a reduced computational cost. The theoretical analysis, supported by numerical experiments, indicates that the new methods are more efficient than schemes based on rational Pad\'e approximants and Taylor polynomials for all tolerances and time interval lengths. The new procedure is particularly recommended to be used in conjunction with exponential integrators for the numerical time integration of the Schr\"odinger equation.

We analyze composition methods with complex coefficients exhibiting the so-called ``symmetry-conjugate'' pattern in their distribution. In particular, we study their behavior with respect to preservation of qualitative properties when projected on the real axis and we compare them with the usual left-right palindromic compositions. New schemes within this family up to order 8 are proposed and their efficiency is tested on several examples. Our analysis shows that higher-order schemes are more efficient even when time step sizes are relatively large., Comment: 24 pages, 4 figures

A new procedure is presented for computing the matrix cosine and sine simultaneously by means of Taylor polynomial approximations. These are factorized so as to reduce the number of matrix products involved. Two versions are developed to be used in single and double precision arithmetic. The resulting algorithms are more efficient than schemes based on Pad\'e approximations for a wide range of norm matrices.

We get compact expressions for the Baker--Campbell--Hausdorff series $Z = \log(\e^X \, \e^Y)$ in terms of right-nested commutators. The reduction in the number of terms originates from two facts: (i) we use as a starting point an explicit expression directly involving independent commutators and (ii) we derive a complete set of identities arising among right-nested commutators. The procedure allows us to obtain the series with fewer terms than when expressed in the classical Hall basis at least up to terms of grade 10., Comment: 13 pages, no figures

New families of fourth-order composition methods for the numerical integration of initial value problems defined by ordinary differential equations are proposed. They are designed when the problem can be separated into three parts in such a way that each part is explicitly solvable. The methods are obtained by applying different optimization criteria and preserve geometric properties of the continuous problem by construction. Different numerical examples exhibit their improved performance with respect to previous splitting methods in the literature.

In this paper, we are concerned with the construction and analysis of a new class of methods obtained as double jump compositions with complex coefficients and projection on the real axis. It is shown in particular that the new integrators are symmetric and symplectic up to high orders if one uses a symmetric and symplectic basic method. In terms of efficiency, the aforementioned technique requires fewer stages than standard compositions of the same orders and is thus expected to lead to faster methods., Comment: 19 pages, 7 figures

The objective of this work is the introduction and investigation of favourable time integration methods for the Gross--Pitaevskii equation with rotation term. Employing a reformulation in rotating Lagrangian coordinates, the equation takes the form of a nonlinear Schr{\"o}dinger equation involving a space-time-dependent potential. A natural approach that combines commutator-free quasi-Magnus exponential integrators with operator splitting methods and Fourier spectral space discretisations is proposed. Furthermore, the special structure of the Hamilton operator permits the design of specifically tailored schemes. Numerical experiments confirm the good performance of the resulting exponential integrators., Comment: 24 pages, 13 figures

In this paper, we are concerned with a formulation of Magnus and Floquet-Magnus expansions for general nonlinear differential equations. To this aim, we introduce suitable continuous variable transformations generated by operators. As an application of the simple formulas so-obtained, we explicitly compute the first terms of the Floquet-Magnus expansion for the Van der Pol oscillator and the nonlinear Schr\"odinger equation on the torus., Comment: 19 pages, 2 figures

We propose new local error estimators for splitting and composition methods. They are based on the construction of lower order schemes obtained at each step as a linear combination of the intermediate stages of the integrator, so that the additional computational cost required for their evaluation is almost insignificant. These estimators can be subsequently used to adapt the step size along the integration. Numerical examples show the efficiency of the procedure., Comment: 23 pages, 4 figures

We propose and analyze a symmetric version of the Zassenhaus formula for disentangling the exponential of two non-commuting operators. A recursive procedure for generating the expansion up to any order is presented which also allows one to get an enlarged domain of convergence when it is formulated for matrices. It is shown that the approximations obtained by truncating the infinite expansion considerably improve those arising from the standard Zassenhaus formula., Comment: 18 pages, 3 figures

In this work, we present a new way to compute the Taylor polynomial of the matrix exponential which reduces the number of matrix multiplications in comparison with the de-facto standard Patterson-Stockmeyer method. This reduction is sufficient to make the method superior in performance to Pad\'e approximants by 10-30% over a range of values for the matrix norms and thus we propose its replacement in standard software kits. Numerical experiments show the performance of the method and illustrate its stability., Comment: 14pages, 2 figures

We present a general expression for any term of the Magnus series as an iterated integral of a linear combination of independent right-nested commutators with given coefficients. The relation with the Malvenuto--Reutenauer Hopf algebra of permutations is also discussed., Comment: 16 pages

An algorithm is presented for generating successive approximations to trigonometric functions of sums of non-commuting matrices. The resulting expressions involve nested commutators of the respective matrices. The procedure is shown to converge in the convergent domain of the Zassenhaus formula and can be useful in the perturbative treatment of quantum mechanical problems, where exponentials of sums of non-commuting skew-Hermitian matrices frequently appear., Comment: 11 pages, 1 figure

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.

A special purpose solver, based on the Magnus expansion, well suited for the integration of the linear three neutrino oscillations equations in matter is proposed. The computations are speeded up to two orders of magnitude with respect to a general numerical integrator, a fact that could smooth the way for massive numerical integration concomitant with experimental data analyses. Detailed illustrations about numerical procedure and computer time costs are provided., Comment: 10 pages, 3 figures

We consider the Vlasov-Poisson equation in a Hamiltonian framework and derive new time splitting methods based on the decomposition of the Hamiltonian functional between the kinetic and electric energy. Assuming smoothness of the solutions, we study the order conditions of such methods. It appears that these conditions are of Runge-Kutta-Nystr{\"o}m type. In the one dimensional case, the order conditions can be further simplified, and efficient methods of order 6 with a reduced number of stages can be constructed. In the general case, high-order methods can also be constructed using explicit computations of commutators. Numerical results are performed and show the benefit of using high-order splitting schemes in that context. Complete and self-contained proofs of convergence results and rigorous error estimates are also given.

KEEN waves are nonlinear, non-stationary, self-organized asymptotic states in Vlasov plasmas outside the scope or purview of linear theory constructs such as electron plasma waves or ion acoustic waves. Nonlinear stationary mode theories such as those leading to BGK modes also do not apply. The range in velocity that is strongly perturbed by KEEN waves depends on the amplitude and duration of the ponderomotive force used to drive them. Smaller amplitude drives create highly localized structures attempting to coalesce into KEEN waves. These cases have much more chaotic and intricate time histories than strongly driven ones. The narrow range in which one must maintain adequate velocity resolution in the weakly driven cases challenges xed grid numerical schemes. What is missing there is the capability of resolving locally in velocity while maintaining a coarse grid outside the highly perturbed region of phase space. We here report on a new Semi-Lagrangian Vlasov-Poisson solver based on conservative non-uniform cubic splines in velocity that tackles this problem head on. An additional feature of our approach is the use of a new high-order time-splitting scheme which allows much longer simulations per computational e ort. This is needed for low amplitude runs which take a long time to set up KEEN waves, if they are able to do so at all. The new code's performance is compared to uniform grid simulations and the advantages quanti ed. The birth pains associated with KEEN waves which are weakly driven is captured in these simulations. These techniques allow the e cient simulation of KEEN waves in multiple dimensions which will be tackled next as well as generalizations to Vlasov-Maxwell codes which are essential to understanding the impact of KEEN waves in practice.

We construct numerical integrators for Hamiltonian problems that may advantageously replace the standard Verlet time-stepper within Hybrid Monte Carlo and related simulations. Past attempts have often aimed at boosting the order of accuracy of the integrator and/or reducing the size of its error constants; order and error constant are relevant concepts in the limit of vanishing step-length. We propose an alternative methodology based on the performance of the integrator when sampling from Gaussian distributions with not necessarily small step-lengths. We construct new splitting formulae that require two, three or four force evaluations per time-step. Limited, proof-of-concept numerical experiments suggest that the new integrators may provide an improvement on the efficiency of the standard Verlet method, especially in problems with high dimensionality., Comment: 30 pages, 5 figures