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 analyse 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.

This article surveys research on the application of compatible finite element methods to large-scale atmosphere and ocean simulation. Compatible finite element methods extend Arakawa’s C-grid finite difference scheme to the finite element world. They are constructed from a discrete de Rham complex, which is a sequence of finite element spaces linked by the operators of differential calculus. The use of discrete de Rham complexes to solve partial differential equations is well established, but in this article we focus on the specifics of dynamical cores for simulating weather, oceans and climate. The most important consequence of the discrete de Rham complex is the Hodge–Helmholtz decomposition, which has been used to exclude the possibility of several types of spurious oscillations from linear equations of geophysical flow. This means that compatible finite element spaces provide a useful framework for building dynamical cores. In this article we introduce the main concepts of compatible finite element spaces, and discuss their wave propagation properties. We survey some methods for discretizing the transport terms that arise in dynamical core equation systems, and provide some example discretizations, briefly discussing their iterative solution. Then we focus on the recent use of compatible finite element spaces in designing structure preserving methods, surveying variational discretizations, Poisson bracket discretizations and consistent vorticity transport.

The ongoing energy transition challenges the stability of the electrical power system. Stable operation of the electrical power grid requires both the voltage (amplitude) and the frequency to stay within operational bounds. While much research has focused on frequency dynamics and stability, the voltage dynamics has been neglected. Here, we study frequency and voltage stability in the case of simple networks via linear stability and bulk analysis. In particular, our linear stability analysis of the network shows that the frequency secondary control guarantees the stability of a particular electric network. Even more interesting, while we only consider secondary frequency control, we observe a stabilising effect on the voltage dynamics, especially in our numerical bulk analysis.

We investigate, both analytically and numerically, dispersive fractalisation and quantisation of solutions to periodic linear and nonlinear Fermi–Pasta–Ulam–Tsingou systems. When subject to periodic boundary conditions and discontinuous initial conditions, e.g., a step function, both the linearised and nonlinear continuum models for FPUT exhibit fractal solution profiles at irrational times (as determined by the coefficients and the length of the interval) and quantised profiles (piecewise constant or perturbations thereof) at rational times. We observe a similar effect in the linearised FPUT chain at times t where these models have validity, namely t = O(h−2), where h is proportional to the intermass spacing or, equivalently, the reciprocal of the number of masses. For nonlinear periodic FPUT systems, our numerical results suggest a somewhat similar behaviour in the presence of small nonlinearities, which disappears as the nonlinear force increases in magnitude. However, these phenomena are manifested on very long time intervals, posing a severe challenge for numerical integration as the number of masses increases. Even with the high-order splitting methods used here, our numerical investigations are limited to nonlinear FPUT chains with a smaller number of masses than would be needed to resolve this question unambiguously.

We develop a deep autoencoder architecture that can be used to find a coordinate transformation which turns a non-linear partial differential equation (PDE) into a linear PDE. Our architecture is motivated by the linearising transformations provided by the Cole–Hopf transform for Burgers’ equation and the inverse scattering transform for completely integrable PDEs. By leveraging a residual network architecture, a near-identity transformation can be exploited to encode intrinsic coordinates in which the dynamics are linear. The resulting dynamics are given by a Koopman operator matrix K. The decoder allows us to transform back to the original coordinates as well. Multiple time step prediction can be performed by repeated multiplication by the matrix K in the intrinsic coordinates. We demonstrate our method on a number of examples, including the heat equation and Burgers’ equation, as well as the substantially more challenging Kuramoto–Sivashinsky equation, showing that our method provides a robust architecture for discovering linearising transforms for non-linear PDEs.

In the theory of spontaneous combustion, identifying the critical value of the Frank-Kamenetskii parameter corresponds to solving a bifurcation point problem. There are two different numerical methods used to solve this problem—the direct and indirect numerical methods. The latter finds the bifurcation point by solving a partial differential equation (PDE) problem. This is a better method to find the bifurcation point for complex geometries. This paper improves the indirect numerical method by combining the grid-domain extension method with the matrix equation computation method. We calculate the critical parameters of the Frank-Kamenetskii equation for some complex geometries using the indirect numerical method. Our results show that both the curve of the outer boundary and the height of the geometries have an effect on the values of the critical Frank-Kamenetskii parameter, however, they have little effect on the critical dimensionless temperature.

Macroscale “continuum” level predictions are made by a new way to construct computationally efficient “wrappers” around fine-scale, microscopic, detailed descriptions of dynamical systems, such as molecular dynamics. It is often significantly easier to code a microscale simulator with periodicity: so the challenge addressed here is to develop a scheme that uses only a given periodic microscale simulator; specifically, one for atomistic dynamics. Numerical simulations show that applying a suitable proportional controller within “action regions” of a patch of atomistic simulation effectively predicts the macroscale transport of heat. Theoretical analysis establishes that such an approach will generally be effective and efficient, and also determines good values for the strength of the proportional controller. This work has the potential to empower systematic analysis and understanding at a macroscopic system level when only a given microscale simulator is available.

In this paper we investigate gradient estimation for a class of contracting stochastic systems on a continuous state space. We find conditions on the one-step transitions, namely differentiability and contraction in a Wasserstein distance, that guarantee differentiability of stationary costs. Then we show how to estimate the derivatives, deriving an estimator that can be seen as a generalization of the forward sensitivity analysis method used in deterministic systems. We apply the results to examples, including a neural network model.

This paper describes an adaptive preconditioner for numerical continuation of incompressible Navier–Stokes flows based on Stokes preconditioning [42] which has been used successfully in studies of pattern formation in convection. The preconditioner takes the form of the Helmholtz operator I–ΔtL which maps the identity (no preconditioner) for Δt≪1 to Laplacian preconditioning for Δt≫1. It is built on a first order Euler time-discretization scheme and is part of the family of matrix-free methods. The preconditioner is tested on two fluid configurations: three-dimensional doubly diffusive convection and a two-dimensional projection of a shear flow. In the former case, it is found that Stokes preconditioning is more efficient for , away from the values used in the literature. In the latter case, the simple use of the preconditioner is not sufficient and it is necessary to split the system of equations into two subsystems which are solved simultaneously using two different preconditioners, one of which is parameter dependent. Due to the nature of these applications and the flexibility of the approach described, this preconditioner is expected to help in a wide range of applications.

We study a novel class of numerical integrators, the adapted nested force-gradient schemes, used within the molecular dynamics step of the Hybrid Monte Carlo (HMC) algorithm. We test these methods in the Schwinger model on the lattice, a well known benchmark problem. We derive the analytical basis of nested force-gradient type methods and demonstrate the advantage of the proposed approach, namely reduced computational costs compared with other numerical integration schemes in HMC.

In this paper, we present some efficient numerical schemes to solve a two-phase hydrodynamics coupled phase field model with moving contact line boundary conditions. The model is a nonlinear coupling system, which consists the Navier-Stokes equations with the general Navier Boundary conditions or degenerated Navier Boundary conditions, and the Allen-Cahn type phase field equations with dynamical contact line boundary condition or static contact line boundary condition. The proposed schemes are linear and unconditionally energy stable, where the energy stabilities are proved rigorously. Various numerical tests are performed to show the accuracy and efficiency thereafter.

The 2D Maxwell eigenproblems are studied from a new point of view. An electromagnetic problem is cast from the Lagrangian system with single variable into the Hamiltonian system with dual variables. The electric and magnetic components transverse to the wave propagation direction are treated as dual variables to each other. Higher order curl-conforming and divergence-conforming vector basis functions are used to construct dual vector spectral elements. Numerical examples demonstrate some unique advantages of the proposed method.

We study the construction of symplectic Runge-Kutta methods for stochastic Hamiltonian systems (SHS). Three types of systems, SHS with multiplicative noise, special separable Hamiltonians and multiple additive noise, respectively, are considered in this paper. Stochastic Runge-Kutta (SRK) methods for these systems are investigated, and the corresponding conditions for SRK methods to preserve the symplectic property are given. Based on the weak/strong order and symplectic conditions, some effective schemes are derived. In particular, using the algebraic computation, we obtained two classes of high weak order symplectic Runge-Kutta methods for SHS with a single multiplicative noise, and two classes of high strong order symplectic Runge-Kutta methods for SHS with multiple multiplicative and additive noise, respectively. The numerical case studies confirm that the symplectic methods are efficient computational tools for long-term simulations.

We prove that a class of A-stable symplectic Runge–Kutta time semi-discretizations (including the Gauss–Legendre methods) applied to a class of semilinear Hamiltonian partial differential equations (PDEs) that are well posed on spaces of analytic functions with analytic initial data can be embedded into a modified Hamiltonian flow up to an exponentially small error. Consequently, such time semi-discretizations conserve the modified Hamiltonian up to an exponentially small error. The modified Hamiltonian is O(hp )-close to the original energy, where p is the order of the method and h is the time-step size. Examples of such systems are the semilinear wave equation, and the nonlinear Schrödinger equation with analytic nonlinearity and periodic boundary conditions. Standard Hamiltonian interpolation results do not apply here because of the occurrence of unbounded operators in the construction of the modified vector field. This loss of regularity in the construction can be taken care of by projecting the PDE to a subspace in which the operators occurring in the evolution equation are bounded, and by coupling the number of excited modes and the number of terms in the expansion of the modified vector field with the step size. This way we obtain exponential estimates of the form O(exp(–c/h 1/(1+q))) with c > 0 and q ⩾ 0; for the semilinear wave equation, q = 1, and for the nonlinear Schrödinger equation, q = 2. We give an example which shows that analyticity of the initial data is necessary to obtain exponential estimates.

This paper explores the discrete singular convolution method for Hamiltonian PDEs. The differential matrices corresponding to two delta type kernels of the discrete singular convolution are presented analytically, which have the properties of high-order accuracy, bandlimited structure and thus can be excellent candidates for the spatial discretizations for Hamiltonian PDEs. Taking the nonlinear Schrödinger equation and the coupled Schrödinger equations for example, we construct two symplectic integrators combining this kind of differential matrices and appropriate symplectic time integrations, which both have been proved to satisfy the square conservation laws. Comprehensive numerical experiments including comparisons with the central finite difference method, the Fourier pseudospectral method, the wavelet collocation method are given to show the advantages of the new type of symplectic integrators.

In this paper, we study the Camassa-Holm equation and the Degasperis-Procesi equation. The two equations are in the family of integrable peakon equations, and both have very rich geometric properties. Based on these geometric structures, we construct the geometric numerical integrators for simulating their soliton solutions. The Camassa-Holm equation and the Degasperis-Procesi equation have many common properties, however they also have the significant difference, for example there exist the shock wave solutions for the Degasperis-Procesi equation. By using the symplectic Fourier pseudo-spectral integrator, we simulate the peakon solutions of the two equations. To illustrate the smooth solitons and shock wave solutions of the DP equation, we use the splitting technique and combine the composition methods. In the numerical experiments, comparisons of these two kinds of methods are presented in terms of accuracy, computational cost and invariants preservation.

Many physical processes are described by elliptic or parabolic partial differential equations. For linear stability problems associated with such equations, the inverse Laplacian provides a very effective preconditioner. In addition, it is also readily available in most scientific calculations in the form of a Poisson solver or an implicit diffusive timestep. We incorporate Laplacian preconditioning into the inverse Arnoldi method, using BiCGSTAB to solve the large linear systems. Two successful implementations are described: spherical Couette flow described by the Navier-Stokes equations and Bose-Einstein condensation described by the nonlinear Schrödinger equation.

Conservation laws provide important constraints on the solutions of partial differential equations (PDEs), therefore it is important to preserve them when discretizing such equations. In this paper, a new systematic method for discretizing a PDE, so as to preserve the local form of multiple conservation laws, is presented. The technique, which uses symbolic computation, is applied to the Korteweg–de Vries (KdV) equation to find novel explicit and implicit schemes that have finite difference analogues of its first and second conservation laws and its first and third conservation laws. The resulting schemes are numerically compared with a multisymplectic scheme.