19 research outputs found
Evaluating the Evans function: Order reduction in numerical methods
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 reference
Three-dimensional solutions for the geostrophic flow in the Earth's core
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 figure
A Krylov subspace algorithm for evaluating the phi-functions appearing in exponential integrators
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 refere
Computing stability of multi-dimensional travelling waves
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 page
Convergence of the Magnus series
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 correction
On an asymptotic method for computing the modified energy for symplectic methods
We revisit an algorithm by Skeel et al. [5,16] for computing the modified, or shadow, energy associated with symplectic discretizations of Hamiltonian systems. We amend the algorithm to use Richardson extrapolation in order to obtain arbitrarily high order of accuracy. Error estimates show that the new method captures the exponentially small drift associated with such discretizations. Several numerical examples illustrate the theory
On the global error of discretization methods for ordinary differential equations
Discretization methods for ordinary differential equations are usually not exact; they commit an error at every step of the algorithm. All these errors combine to form the global error, which is the error in the final result. The global error is the subject of this thesis. In the first half of the thesis, accurate a priori estimates of the global error are derived. Three different approaches are followed: to combine the effects of the errors committed at every step, to expand the global error in an asymptotic series in the step size, and to use the theory of modified equations. The last approach, which is often the most useful one, yields an estimate which is correct up to a term of order h 2p, where h denotes the step size and p the order of the numerical method. This result is then applied to estimate the global error for the Airy equa-tion (and related oscillators that obey the Liouville–Green approximation) and the Emden–Fowler equation. The latter example has the interesting feature that it is not sufficient to consider only the leading global error term, because subsequent terms of higher order in the step size may grow faster in time. The second half of the thesis concentrates on minimizing the global error by varying the step size. It is argued that the correct objective function is the norm of the global error over the entire integration interval. Specifically, the L2 norm and the L ∞ norm are studied. In the former case, Pontryagin’s Minimum Principle converts the problem to a boundary value problem, which may be solved analyti-cally or numerically. When the L ∞ norm is used, a boundary value problem with a complementarity condition results. Alternatively, the Exterior Penalty Method may be employed to get a boundary value problem without complementarity condition, which can be solved by standard numerical software. The theory is illustrated by calculating the optimal step size for solving the Dahlquist test equation and the Kepler problem.