Search

The paper investigates two new use cases for the Boris Spectral Deferred Corrections (Boris-SDC) time integrator for plasma simulations. First, we show that using Boris-SDC as a particle pusher in an electrostatic particle-in-cell (PIC) code can, at least in the linear regime, improve simulation accuracy compared with the standard second order Boris method. In some instances, the higher order of Boris-SDC even allows a much larger time step, leading to modest computational gains. Second, we propose a modification of Boris-SDC for the relativistic regime. Based on an implementation of Boris-SDC in the \textsc{runko} PIC code, we demonstrate for a relativistic Penning trap that Boris-SDC retains its high order of convergence for velocities ranging from $0.5c$ to $>0.99c$. We also show that for the force-free case where acceleration from electric and magnetic field cancel, Boris-SDC produces less numerical drift than Boris.

Mounting evidence from both seismology and experiments on core composition suggests the existence of a layer of stably stratified fluid at the top of Earth's outer core. In this work we examine the structure of the geomagnetic field within such a layer, building on the important but little known work of Malkus (1979). We assume (i) an idealised magnetostrophic spherical model of the geodynamo neglecting inertia, viscosity and the solid inner core, and (ii) a strongly stratified layer of constant depth immediately below the outer boundary within which there is no spherically radial flow. Due to the restricted dynamics, Malkus showed that the geomagnetic field must obey a certain condition which is a more restrictive version of the condition of Taylor (1963). The nonlinear nature of these constraints makes finding a magnetic field that obeys them, here termed a Malkus state, a challenging task. Nevertheless, such Malkus states when constrained further by geomagnetic observations have the potential to probe the interior of the core. By focusing on a particular class of magnetic fields for which the Malkus constraints are linear, we describe a constructive method that turns any purely-poloidal field into an exact Malkus state by adding a suitable toroidal field. We consider poloidal fields following a prescribed smooth profile within the core that match observation-derived models of the magnetic field in either epoch 2015 or the 10000-yr time averaged field. Multiple possible solutions for the toroidal field exist, hence we determine the Malkus state of miumum toroidal energy and we find that it has a strong azimuthal toroidal field, larger than the observed poloidal component at the core-mantle boundary. For the 2015 field for a layer of depth 300 km, we estimate a root mean squared azimuthal toroidal field of 3 mT with a pointwise maximum of 8 mT occurring at a depth of about 70 km., Comment: 27 pages, 10 figures

In his seminal work, Taylor (1963) argued that the geophysically relevant limit for dynamo action within the outer core is one of negligibly small inertia and viscosity in the magnetohydrodynamic equations. Within this limit, he showed the existence of a necessary condition, now well known as Taylor's constraint, which requires that the cylindrically-averaged Lorentz torque must everywhere vanish; magnetic fields that satisfy this condition are termed Taylor states. Taylor further showed that the requirement of this constraint being continuously satisfied through time prescribes the evolution of the geostrophic flow, the cylindrically-averaged azimuthal flow. We show that Taylor's original prescription for the geostrophic flow, as satisfying a given second order ordinary differential equation, is only valid for a small subset of Taylor states. An incomplete treatment of the boundary conditions renders his equation generally incorrect. Here, by taking proper account of the boundaries, we describe a generalisation of Taylor's method that enables correct evaluation of the instantaneous geostrophic flow for any 3D Taylor state. We present the first full-sphere examples of geostrophic flows driven by non-axisymmetric Taylor states. Although in axisymmetry the geostrophic flow admits a mild logarithmic singularity on the rotation axis, in the fully 3D case we show that this is absent and indeed the geostrophic flow appears to be everywhere regular., Comment: 29 Pages, 8 figures

Modified Hamiltonians are used in the field of geometric numerical integration to show that symplectic schemes for Hamiltonian systems are accurate over long times. For nonlinear systems the series defining the modified Hamiltonian usually diverges. In contrast, this paper constructs and analyzes explicit examples of nonlinear systems where the modified Hamiltonian has a closed-form expression and hence converges. These systems arise from the theory of discrete integrable systems. We present cases of one- and two-degrees symplectic mappings arising as reductions of nonlinear integrable lattice equations, for which the modified Hamiltonians can be computed in closed form. These modified Hamiltonians are also given as power series in the time step by Yoshida's method based on the Baker-Campbell-Hausdorff series. Another example displays an implicit dependence on the time step which could be of relevance to certain implicit schemes in numerical analysis. In the light of these examples, the potential importance of integrable mappings to the field of geometric numerical integration is discussed.

We revisit an algorithm by Skeel et al. for computing the modified, or shadow, energy associated with the symplectic discretization of Hamiltonian systems. By rephrasing the algorithm as a Richardson extrapolation scheme arbitrary high order of accuracy is obtained, and provided error estimates show that it does capture the theoretical exponentially small drift associated with such discretizations. Several numerical examples illustrate the theory., Comment: 16 pages, 7 figures

We develop an algorithm for computing the solution of a large system of linear ordinary differential equations (ODEs) with polynomial inhomogeneity. This is equivalent to computing the action of a certain matrix function on the vector representing the initial condition. The matrix function is a linear combination of the matrix exponential and other functions related to the exponential (the so-called phi-functions). Such computations are the major computational burden in the implementation of exponential integrators, which can solve general ODEs. Our approach is to compute the action of the matrix function by constructing a Krylov subspace using Arnoldi or Lanczos iteration and projecting the function on this subspace. This is combined with time-stepping to prevent the Krylov subspace from growing too large. The algorithm is fully adaptive: it varies both the size of the time steps and the dimension of the Krylov subspace to reach the required accuracy. We implement this algorithm in the Matlab function phipm and we give instructions on how to obtain and use this function. Various numerical experiments show that the phipm function is often significantly more efficient than the state-of-the-art., Comment: 20 pages, 3 colour figures, code available from http://www.maths.leeds.ac.uk/~jitse/software.html . v2: Various changes to improve presentation as suggested by the referee

A special class of symmetry reductions called nonclassical equivalence transformations is discussed in connection to a class of parameter identification problems represented by partial differential equations. These symmetry reductions relate the forward and inverse problems, reduce the dimension of the equation, yield special types of solutions, and may be incorporated into the boundary conditions as well. As an example, we discuss the nonlinear stationary heat conduction equation and show that this approach permits the study of the model on new types of domains. Our MAPLE routine GENDEFNC which uses the package DESOLV (authors Carminati and Vu) has been updated for this propose and its output is the nonlinear partial differential equation system of the determining equations of the nonclassical equivalence transformations., Comment: 18 pages

We present a numerical method for computing the pure-point spectrum associated with the linear stability of multi-dimensional travelling fronts to parabolic nonlinear systems. Our method is based on the Evans function shooting approach. Transverse to the direction of propagation we project the spectral equations onto a finite Fourier basis. This generates a large, linear, one-dimensional system of equations for the longitudinal Fourier coefficients. We construct the stable and unstable solution subspaces associated with the longitudinal far-field zero boundary conditions, retaining only the information required for matching, by integrating the Riccati equations associated with the underlying Grassmannian manifolds. The Evans function is then the matching condition measuring the linear dependence of the stable and unstable subspaces and thus determines eigenvalues. As a model application, we study the stability of two-dimensional wrinkled front solutions to a cubic autocatalysis model system. We compare our shooting approach with the continuous orthogonalization method of Humpherys and Zumbrun. We then also compare these with standard projection methods that directly project the spectral problem onto a finite multi-dimensional basis satisfying the boundary conditions., Comment: 23 pages, 9 figures (some in colour). v2: added details and other changes to presentation after referees' comments, now 26 pages

The Evans function is a tool for assessing the stability of travelling wave solutions for partial differential equations. A recent paper (math.NA/0605581) analyzes the order reduction experienced when evaluating the Evans function numerically. The details of some lengthy calculations were excluded from that paper for clarity. The purpose of this technical report is to make these details publicly available. The report is not intended to be read on its own; the reader is referred to math.NA/0605581 for background and references., Comment: 22 pages

The Magnus series is an infinite series which arises in the study of linear ordinary differential equations. If the series converges, then the matrix exponential of the sum equals the fundamental solution of the differential equation. The question considered in this paper is: When does the series converge? The main result establishes a sufficient condition for convergence, which improves on several earlier results., Comment: 11 pages; v2: added justification for conjecture, minor clarifications and corrections

We consider the numerical evaluation of the Evans function, a Wronskian-like determinant that arises in the study of the stability of travelling waves. Constructing the Evans function involves matching the solutions of a linear ordinary differential equation depending on the spectral parameter. The problem becomes stiff as the spectral parameter grows. Consequently, the Gauss--Legendre method has previously been used for such problems; however more recently, methods based on the Magnus expansion have been proposed. Here we extensively examine the stiff regime for a general scalar Schr\"odinger operator. We show that although the fourth-order Magnus method suffers from order reduction, a fortunate cancellation when computing the Evans matching function means that fourth-order convergence in the end result is preserved. The Gauss--Legendre method does not suffer from order reduction, but it does not experience the cancellation either, and thus it has the same order of convergence in the end result. Finally we discuss the relative merits of both methods as spectral tools., Comment: 21 pages, 3 figures; removed superfluous material (+/- 1 page), added paragraph to conclusion and two references

The Alekseev-Gr{\"o}bner lemma is combined with the theory of modified equations to obtain an \emph{a priori} estimate for the global error of numerical integrators. This estimate is correct up to a remainder term of order $h^{2p}$, where $h$ denotes the step size and $p$ the order of the method. It is applied to a class of nonautonomous linear oscillatory equations, which includes the Airy equation, thereby improving prior work which only gave the $h^p$ term. Next, nonlinear oscillators whose behaviour is described by the Emden-Fowler equation $y'' + t^\nu y^n = 0$ are considered, and global errors committed by Runge-Kutta methods are calculated. Numerical experiments show that the resulting estimates are generally accurate. The main conclusion is that we need to do a full calculation to obtain good estimates: the behaviour is different from the linear case, it is not sufficient to look only at the leading term, and merely considering the local error does not provide an accurate picture either., Comment: 23 pages, 2 figures. v2: generalized Section 4 and 5, and minor changes