David I. Ketcheson

David I. Ketcheson, Matteo Parsani, Randall J. LeVeque

We present a finite volume method that is applicable to hyperbolic PDEs including spatially varying and semilinear nonconservative systems. The spatial discretization, like that of the well-known Clawpack software, is based on solving Riemann problems and calculating fluctuations (not fluxes). The implementation employs weighted essentially non-oscillatory reconstruction in space and strong stability preserving Runge-Kutta integration in time. The method can be extended to arbitrarily high order of accuracy and allows a well-balanced implementation for capturing solutions of balance laws near steady state. This well-balancing is achieved through the $f$-wave Riemann solver and a novel wave-slope WENO reconstruction procedure. The wide applicability and advantageous properties of the method are demonstrated through numerical examples, including problems in nonconservative form, problems with spatially varying fluxes, and problems involving near-equilibrium solutions of balance laws.

David I. Ketcheson, Kyle T. Mandli, Aron Ahmadia et al.

Development of scientific software involves tradeoffs between ease of use, generality, and performance. We describe the design of a general hyperbolic PDE solver that can be operated with the convenience of MATLAB yet achieves efficiency near that of hand-coded Fortran and scales to the largest supercomputers. This is achieved by using Python for most of the code while employing automatically-wrapped Fortran kernels for computationally intensive routines, and using Python bindings to interface with a parallel computing library and other numerical packages. The software described here is PyClaw, a Python-based structured grid solver for general systems of hyperbolic PDEs \cite{pyclaw}. PyClaw provides a powerful and intuitive interface to the algorithms of the existing Fortran codes Clawpack and SharpClaw, simplifying code development and use while providing massive parallelism and scalable solvers via the PETSc library. The package is further augmented by use of PyWENO for generation of efficient high-order weighted essentially non-oscillatory reconstruction code. The simplicity, capability, and performance of this approach are demonstrated through application to example problems in shallow water flow, compressible flow and elasticity.

David I. Ketcheson

We further develop a simple modification of Runge--Kutta methods that guarantees conservation or stability with respect to any inner-product norm. The modified methods can be explicit and retain the accuracy and stability properties of the unmodified Runge--Kutta method. We study the properties of the modified methods and show their effectiveness through numerical examples, including application to entropy-stability for first-order hyperbolic PDEs.

David I. Ketcheson, Sigal Gottlieb, Colin B. Macdonald

We investigate the strong stability preserving (SSP) property of two-step Runge-Kutta (TSRK) methods. We prove that all SSP TSRK methods belong to a particularly simple subclass of TSRK methods, in which stages from the previous step are not used. We derive simple order conditions for this subclass. Whereas explicit SSP Runge-Kutta methods have order at most four, we prove that explicit SSP TSRK methods have order at most eight. We present TSRK methods of up to eighth order that were found by numerical search. These methods have larger SSP coefficients than any known methods of the same order of accuracy, and may be implemented in a form with relatively modest storage requirements. The usefulness of the TSRK methods is demonstrated through numerical examples, including integration of very high order WENO discretizations.

Hendrik Ranocha, David I. Ketcheson

We propose and study a class of arbitrarily high-order numerical discretizations that preserve multiple invariants and are essentially explicit (they do not require the solution of any large systems of algebraic equations). In space, we use Fourier Galerkin methods, while in time we use a combination of orthogonal projection and relaxation. We prove and numerically demonstrate the conservation properties of the method by applying it to the Benjamin-Bona-Mahony, Korteweg-de Vries, and nonlinear Schrödinger (NLS) PDEs as well as a hyperbolic approximation of NLS. For each of these equations, the proposed schemes conserve mass, momentum, and energy up to numerical precision. We show that this conservation leads to reduced growth of numerical errors for long-term simulations.

David I. Ketcheson, Aron J. Ahmadia

We consider the problem of finding optimally stable polynomial approximations to the exponential for application to one-step integration of initial value ordinary and partial differential equations. The objective is to find the largest stable step size and corresponding method for a given problem when the spectrum of the initial value problem is known. The problem is expressed in terms of a general least deviation feasibility problem. Its solution is obtained by a new fast, accurate, and robust algorithm based on convex optimization techniques. Global convergence of the algorithm is proven in the case that the order of approximation is one and in the case that the spectrum encloses a starlike region. Examples demonstrate the effectiveness of the proposed algorithm even when these conditions are not satisfied.

David I. Ketcheson, Colin B. Macdonald, Steven J. Ruuth

We study spatially partitioned embedded Runge--Kutta (SPERK) schemes for partial differential equations (PDEs), in which each of the component schemes is applied over a different part of the spatial domain. Such methods may be convenient for problems in which the smoothness of the solution or the magnitudes of the PDE coefficients vary strongly in space. We focus on embedded partitioned methods as they offer greater efficiency and avoid the order reduction that may occur in non-embedded schemes. We demonstrate that the lack of conservation in partitioned schemes can lead to non-physical effects and propose conservative additive schemes based on partitioning the fluxes rather than the ordinary differential equations. A variety of SPERK schemes are presented, including an embedded pair suitable for the time evolution of fifth-order weighted non-oscillatory (WENO) spatial discretizations. Numerical experiments are provided to support the theory.

David I. Ketcheson

Strong stability preserving (SSP) integrators for initial value ODEs preserve temporal monotonicity solution properties in arbitrary norms. All existing SSP methods, including implicit methods, either require small step sizes or achieve only first order accuracy. It is possible to achieve more relaxed step size restrictions in the discretization of hyperbolic PDEs through the use of both upwind- and downwind-biased semi-discretizations. We investigate bounds on the maximum SSP step size for methods that include negative coefficients and downwind-biased semi-discretizations. We prove that the downwind SSP coefficient for linear multistep methods of order greater than one is at most equal to two, while the downwind SSP coefficient for explicit Runge--Kutta methods is at most equal to the number of stages of the method. In contrast, the maximal downwind SSP coefficient for second order Runge--Kutta methods is shown to be unbounded. We present a class of such methods with arbitrarily large SSP coefficient and demonstrate that they achieve second order accuracy for large CFL number.

Manuel Quezada de Luna, David I. Ketcheson

We study two-dimensional wave propagation in materials whose properties vary periodically in one direction only. High order homogenization is carried out to derive a dispersive effective medium approximation. One-dimensional materials with constant impedance exhibit no effective dispersion. We show that a new kind of effective dispersion may arise in two dimensions, even in materials with constant impedance. This dispersion is a macroscopic effect of microscopic diffraction caused by spatial variation in the sound speed. We analyze this dispersive effect by using high-order homogenization to derive an anisotropic, dispersive effective medium. We generalize to two dimensions a homogenization approach that has been used previously for one-dimensional problems. Pseudospectral solutions of the effective medium equations agree to high accuracy with finite volume direct numerical simulations of the variable-coefficient equations.

Inmaculada Higueras, David I. Ketcheson, Tihamér A. Kocsis

Perturbed Runge--Kutta methods (also referred to as downwind Runge--Kutta methods) can guarantee monotonicity preservation under larger step sizes relative to their traditional Runge--Kutta counterparts. In this paper we study, the question of how to optimally perturb a given method in order to increase the radius of absolute monotonicity (a.m.). We prove that for methods with zero radius of a.m., it is always possible to give a perturbation with positive radius. We first study methods for linear problems and then methods for nonlinear problems. In each case, we prove upper bounds on the radius of a.m., and provide algorithms to compute optimal perturbations. We also provide optimal perturbations for many known methods.

Yiannis Hadjimichael, Colin B. Macdonald, David I. Ketcheson et al.

We apply the concept of effective order to strong stability preserving (SSP) explicit Runge-Kutta methods. Relative to classical Runge-Kutta methods, methods with an effective order of accuracy are designed to satisfy a relaxed set of order conditions, but yield higher order accuracy when composed with special starting and stopping methods. We show that this allows the construction of four-stage SSP methods with effective order four (such methods cannot have classical order four). However, we also prove that effective order five methods - like classical order five methods - require the use of non-positive weights and so cannot be SSP. By numerical optimization, we construct explicit SSP Runge-Kutta methods up to effective order four and establish the optimality of many of them. Numerical experiments demonstrate the validity of these methods in practice.

David I. Ketcheson, Lajos Lóczi, Aliya Jangabylova et al.

We investigate dense output formulae (also known as continuous extensions) for strong stability preserving (SSP) Runge-Kutta methods. We require that the dense output formula also possess the SSP property, ideally under the same step-size restriction as the method itself. A general recipe for first-order SSP dense output formulae for SSP methods is given, and second-order dense output formulae for several optimal SSP methods are developed. It is shown that SSP dense output formulae of order 3 and higher do not exist, and that in any method possessing a second-order SSP dense output, the coefficient matrix A has a zero row.

Yiannis Hadjimichael, David I. Ketcheson

The analysis of strong-stability-preserving (SSP) linear multistep methods is extended to semi-discretized problems for which different terms on the right-hand side satisfy different forward Euler (or circle) conditions. Optimal additive and perturbed monotonicity-preserving linear multistep methods are studied in the context of such problems. Optimal perturbed methods attain larger monotonicity-preserving step sizes when the different forward Euler conditions are taken into account. On the other hand, we show that optimal SSP additive methods achieve a monotonicity-preserving step-size restriction no better than that of the corresponding non-additive SSP linear multistep methods.

Lajos Loczi, David I. Ketcheson

We study the radius of absolute monotonicity R of rational functions with numerator and denominator of degree s that approximate the exponential function to order p. Such functions arise in the application of implicit s-stage, order p Runge-Kutta methods for initial value problems and the radius of absolute monotonicity governs the numerical preservation of properties like positivity and maximum-norm contractivity. We construct a function with p=2 and R>2s, disproving a conjecture of van de Griend and Kraaijevanger. We determine the maximum attainable radius for functions in several one-parameter families of rational functions. Moreover, we prove earlier conjectured optimal radii in some families with 2 or 3 parameters via uniqueness arguments for systems of polynomial inequalities. Our results also prove the optimality of some strong stability preserving implicit and singly diagonally implicit Runge-Kutta methods. Whereas previous results in this area were primarily numerical, we give all constants as exact algebraic numbers.

Imre Fekete, David I. Ketcheson, Lajos Lóczi

We propose a technique for investigating stability properties like positivity and forward invariance of an interval for method-of-lines discretizations, and apply the technique to study positivity preservation for a class of TVD semi-discretizations of 1D scalar hyperbolic conservation laws. This technique is a generalization of the approach suggested in ref. 12. We give more relaxed conditions on the time-step for positivity preservation for slope-limited semi-discretizations integrated in time with explicit Runge-Kutta methods. We show that the step-size restrictions derived are sharp in a certain sense, and that many higher-order explicit Runge-Kutta methods, including the classical 4th-order method and all non-confluent methods with a negative Butcher coefficient, cannot generally maintain positivity for these semi-discretizations under any positive step size. We also apply the proposed technique to centered finite difference discretizations of scalar hyperbolic and parabolic problems.

Manvendra Pratap Rajvanshi, David I. Ketcheson

We study relaxation-based approaches for conserving mass and energy in the numerical solution of Schrödinger-Poisson (SP) type systems. Relaxation-based methods offer a general approach that can be applied as post-time step processing to achieve conservation with any time-stepping scheme. Here we study two types of relaxation techniques applied to implicit-explicit Runge-Kutta schemes, with Fourier collocation in space. We also study SP equations with time-varying coefficients (which appear naturally in cosmology) where energy is not conserved but satisfies a balance equation. We show that the fully-discrete system conserves both mass and energy (or satisfies the balance equation in case of time-varying coefficients), up to rounding errors. The effectiveness of these methods is demonstrated via numerical examples, including a three-dimensional cosmological simulation.