Jitse Niesen

Algorithm 919: A Krylov Subspace Algorithm for Evaluating the ϕ-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 ϕ-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.

Convergence of the Magnus Series

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

On a new procedure for finding nonclassical symmetries

A new technique for deriving the determining equations of nonclassical symmetries associated with a partial differential equation system is introduced. The problem is reduced to computing the determining equations of the classical symmetries associated with a related equation with coefficients which depend on the nonclassical symmetry operator. As a consequence, all the symbolic manipulation programs designed for the latter task can also be used to find the determining equations of the nonclassical symmetries, without any adaptation of the program. The algorithm was implemented as the maple routine gendefnc and uses the maple package desolv (authors Carminati and Vu). As an example, we consider the Huxley partial differential equation.

Computing stability of multidimensional traveling waves

We present a numerical method for computing the pure-point spectrum associated with the linear stability of multidimensional traveling 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 ...

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 travel- ling 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 Schrodinger 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 conver- gence in the end result. Finally we discuss the relative merits of both methods as spectral tools.

Localized Patterns in Periodically Forced Systems

Spatially localized, time-periodic structures are common in pattern-forming systems, appearing in fluid mechanics, chemical reactions, and granular media. We examine the existence of oscillatory localized states in a PDE model with single frequency time dependent forcing, introduced in [A. M. Rucklidge and M. Silber, SIAM J. Appl. Math., 8 (2009), pp. 298--347] as a phenomenological model of the Faraday wave experiment. In this study, we reduce the PDE model to the forced complex Ginzburg--Landau equation in the limit of weak forcing and weak damping. This allows us to use the known localized solutions found in [J. Burke, A. Yochelis, and E. Knobloch, SIAM J. Appl. Dyn. Syst., 7 (2008), pp. 651--711]. We reduce the forced complex Ginzburg--Landau equation to the Allen--Cahn equation near onset, obtaining an asymptotically exact expression for localized solutions. We also extend this analysis to the strong forcing case, recovering the Allen--Cahn equation directly without the intermediate step. We find exc...

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.

A NEW CLASS OF SYMMETRY REDUCTIONS FOR PARAMETER IDENTIFICATION PROBLEMS

This paper introduces a new type of symmetry reductions called extended nonclassical symmetries that can be studied for parameter identification problems described by partial differential equations. Including the data function in the parameter space, we show that specific data and parameter classes that lead to a reduced dimension model can be found. More exactly, since the extended nonclassical symmetries relate the forward and inverse problems, the dimension of the studied equation may be reduced by expressing the data and parameter in terms of the group invariants. The main advantage of these new symmetries is that they may be incorporated into the boundary conditions as well, and, consequently, the dimension reduction problem can be analyzed on new types of domains. Special group-invariant solutions or additional information on the parameter can be obtained. Besides, in the case of the first-order partial differential equations, this symmetry reduction method might be an effective alternative tool for finding particular analytical solutions to the studied model, especially when the Maple subroutine pdsolve does not output satisfactory results. As an example, we consider the nonlinear stationary heat conduction equation. 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 extended nonclassical symmetries.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.

Determination of the instantaneous geostrophic flow within the three-dimensional magnetostrophic regime

In his seminal work, Taylor (1963 Proc. R. Soc. Lond. A 274, 274–283. (doi:10.1098/rspa.1963.0130).) 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 Taylors 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 Taylors 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 generalization of Taylors method that enables correct evaluation of the instantaneous geostrophic flow for any three-dimensional 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 three-dimensional case we show that this is absent and indeed the geostrophic flow appears to be everywhere regular.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.

Closed-form modified Hamiltonians for integrable numerical integration schemes

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 twodegrees 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 light of these examples, the potential importance of integrable mappings to the field of geometric numerical integration is discussed.