Research within the Numerical Mathematics Group at KAUST focuses on the design, analysis, and implementation of numerical methods for ordinary and partial differential equations, as well as application of numerical methods to problems in nonlinear wave propagation. Their work spans the numerical method development pipeline from theoretical numerical analysis, to computational experimentation, to software development and high performance computational science, in collaboration with application domain researchers. Particular emphasis is placed on robustness, wide applicability, and efficiency.

Recent Submissions

  • Maximum Principle Preserving Space and Time Flux Limiting for Diagonally Implicit Runge–Kutta Discretizations of Scalar Convection-diffusion Equations

    Quezada de Luna, Manuel; Ketcheson, David I. (Journal of Scientific Computing, Springer Science and Business Media LLC, 2022-08-01) [Article]
    We provide a framework for high-order discretizations of nonlinear scalar convection-diffusion equations that satisfy a discrete maximum principle. The resulting schemes can have arbitrarily high order accuracy in time and space, and can be stable and maximum-principle-preserving (MPP) with no step size restriction. The schemes are based on a two-tiered limiting strategy, starting with a high-order limiter-based method that may have small oscillations or maximum-principle violations, followed by an additional limiting step that removes these violations while preserving high order accuracy. The desirable properties of the resulting schemes are demonstrated through several numerical examples.
  • Design of DIRK Schemes with High Weak Stage Order

    Biswas, Abhijit; Ketcheson, David I.; Seibold, Benjamin; Shirokoff, David (arXiv, 2022-04-24) [Preprint]
    Runge-Kutta (RK) methods may exhibit order reduction when applied to certain stiff problems. While fully implicit RK schemes exist that avoid order reduction via high-stage order, DIRK (diagonally implicit Runge-Kutta) schemes are practically important due to their structural simplicity; however, these cannot possess high stage order. The concept of weak stage order (WSO) can also overcome order reduction, and it is compatible with the DIRK structure. DIRK schemes of WSO up to 3 have been proposed in the past, however, based on a simplified framework that cannot be extended beyond WSO 3. In this work a general theory of WSO is employed to overcome the prior WSO barrier and to construct practically useful high-order DIRK schemes with WSO 4 and above. The resulting DIRK schemes are stiffly accurate, L-stable, have optimized error coefficients, and are demonstrated to perform well on a portfolio of relevant ODE and PDE test problems.
  • Rapid evolution of SARS-CoV-2 challenges human defenses.

    Duarte, Carlos M.; Ketcheson, David I.; Eguíluz, Víctor M; Agusti, Susana; Fernández-Gracia, Juan; Jamil, Tahira; Laiolo, Elisa; Gojobori, Takashi; Alam, Intikhab (Scientific reports, Springer Science and Business Media LLC, 2022-04-19) [Article]
    The race between pathogens and their hosts is a major evolutionary driver, where both reshuffle their genomes to overcome and reorganize the defenses for infection, respectively. Evolutionary theory helps formulate predictions on the future evolutionary dynamics of SARS-CoV-2, which can be monitored through unprecedented real-time tracking of SARS-CoV-2 population genomics at the global scale. Here we quantify the accelerating evolution of SARS-CoV-2 by tracking the SARS-CoV-2 mutation globally, with a focus on the Receptor Binding Domain (RBD) of the spike protein determining infection success. We estimate that the > 820 million people that had been infected by October 5, 2021, produced up to 1021 copies of the virus, with 12 new effective RBD variants appearing, on average, daily. Doubling of the number of RBD variants every 89 days, followed by selection of the most infective variants challenges our defenses and calls for a shift to anticipatory, rather than reactive tactics involving collaborative global sequencing and vaccination
  • Algebraic Structure of the Weak Stage Order Conditions for Runge-Kutta Methods

    Biswas, Abhijit; Ketcheson, David I.; Seibold, Benjamin; Shirokoff, David (arXiv, 2022-04-07) [Preprint]
    Runge-Kutta (RK) methods may exhibit order reduction when applied to stiff problems. For linear problems with time-independent operators, order reduction can be avoided if the method satisfies certain weak stage order (WSO) conditions, which are less restrictive than traditional stage order conditions. This paper outlines the first algebraic theory of WSO, and establishes general order barriers that relate the WSO of a RK scheme to its order and number of stages for both fully-implicit and DIRK schemes. It is shown in several scenarios that the constructed bounds are sharp. The theory characterizes WSO in terms of orthogonal invariant subspaces and associated minimal polynomials. The resulting necessary conditions on the structure of RK methods with WSO are then shown to be of practical use for the construction of such schemes.
  • Bound-preserving Flux Limiting for High-Order Explicit Runge-Kutta Time Discretizations of Hyperbolic Conservation Laws

    Kuzmin, Dmitri; Quezada de Luna, Manuel; Ketcheson, David I.; Grull, Johanna (JOURNAL OF SCIENTIFIC COMPUTING, Springer Science and Business Media LLC, 2022-03-04) [Article]
    We introduce a general framework for enforcing local or global maximum principles in high-order space-time discretizations of a scalar hyperbolic conservation law. We begin with sufficient conditions for a space discretization to be bound preserving (BP) and satisfy a semi-discrete maximum principle. Next, we propose a global monolithic convex (GMC) flux limiter which has the structure of a flux-corrected transport (FCT) algorithm but is applicable to spatial semi-discretizations and ensures the BP property of the fully discrete scheme for strong stability preserving (SSP) Runge–Kutta time discretizations. To circumvent the order barrier for SSP time integrators, we constrain the intermediate stages and/or the final stage of a general high-order RK method using GMC-type limiters. In this work, our theoretical and numerical studies are restricted to explicit schemes which are provably BP for sufficiently small time steps. The new GMC limiting framework offers the possibility of relaxing the bounds of inequality constraints to achieve higher accuracy at the cost of more stringent time step restrictions. The ability of the presented limiters to recognize undershoots/overshoots, as well as smooth solutions, is verified numerically for three representative RK methods combined with weighted essentially nonoscillatory (WENO) finite volume space discretizations of linear and nonlinear test problems in 1D. In this context, we enforce global bounds and prove preservation of accuracy for the linear advection equation.
  • Computing with B-series

    Ketcheson, David I.; Ranocha, Hendrik (arXiv, 2021-11-23) [Preprint]
    We present BSeries.jl, a Julia package for the computation and manipulation of B-series, which are a versatile theoretical tool for understanding and designing discretizations of differential equations. We give a short introduction to the theory of B-series and associated concepts and provide examples of their use, including method composition and backward error analysis. The associated software is highly performant and makes it possible to work with B-series of high order.
  • On the rate of error growth in time for numerical solutions of nonlinear dispersive wave equations

    Ranocha, Hendrik; Quezada de Luna, Manuel; Ketcheson, David I. (Partial Differential Equations and Applications, Springer Science and Business Media LLC, 2021-10-18) [Article]
    AbstractWe study the numerical error in solitary wave solutions of nonlinear dispersive wave equations. A number of existing results for discretizations of solitary wave solutions of particular equations indicate that the error grows quadratically in time for numerical methods that do not conserve energy, but grows only linearly for conservative methods. We provide numerical experiments suggesting that this result extends to a very broad class of equations and numerical methods.
  • manuel-quezada/BP_Lim_for_imp_RK_Methods:

    Quezada de Luna, Manuel; Ketcheson, David I. (Github, 2021-09-15) [Software]
  • A review of matrix sir arino epidemic models

    Avram, Florin; Adenane, Rim; Ketcheson, David I. (Mathematics, MDPI AG, 2021-06-28) [Article]
    Many of the models used nowadays in mathematical epidemiology, in particular in COVID-19 research, belong to a certain subclass of compartmental models whose classes may be divided into three “(x, y, z)” groups, which we will call respectively “susceptible/entrance, diseased, and output” (in the classic SIR case, there is only one class of each type). Roughly, the ODE dynamics of these models contains only linear terms, with the exception of products between x and y terms. It has long been noticed that the reproduction number R has a very simple Formula in terms of the matrices which define the model, and an explicit first integral Formula is also available. These results can be traced back at least to Arino, Brauer, van den Driessche, Watmough, and Wu (2007) and to Feng (2007), respectively, and may be viewed as the “basic laws of SIR-type epidemics”. However, many papers continue to reprove them in particular instances. This motivated us to redraw attention to these basic laws and provide a self-contained reference of related formulas for (x, y, z) models. For the case of one susceptible class, we propose to use the name SIR-PH, due to a simple probabilistic interpretation as SIR models where the exponential infection time has been replaced by a PH-type distribution. Note that to each SIR-PH model, one may associate a scalar quantity Y(t) which satisfies “classic SIR relations”, which may be useful to obtain approximate control policies.
  • Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention

    Ketcheson, David I. (Journal of Mathematical Biology, Springer Science and Business Media LLC, 2021-06-26) [Article]
    We consider the problem of controlling an SIR-model epidemic by temporarily reducing the rate of contact within a population. The control takes the form of a multiplicative reduction in the contact rate of infectious individuals. The control is allowed to be applied only over a finite time interval, while the objective is to minimize the total number of individuals infected in the long-time limit, subject to some cost function for the control. We first consider the no-cost scenario and analytically determine the optimal control and solution. We then study solutions when a cost of intervention is included, as well as a cost associated with overwhelming the available medical resources. Examples are studied through the numerical solution of the associated Hamilton-Jacobi-Bellman equation. Finally, we provide some examples related directly to the current pandemic.
  • Positivity preservation of implicit discretizations of the advection equation

    Hadjimichael, Yiannis; Ketcheson, David I.; Lóczi, Lajos (arXiv, 2021-05-16) [Preprint]
    We analyze, from the viewpoint of positivity preservation, certain discretizations of a fundamental partial differential equation, the one-dimensional advection equation with periodic boundary condition. The full discretization is obtained by coupling a finite difference spatial semi-discretization (the second- and some higher-order centered difference schemes, or the Fourier spectral collocation method) with an arbitrary $\theta$-method in time (including the forward and backward Euler methods, and a second-order method by choosing $\theta\in [0,1]$ suitably). The full discretization generates a two-parameter family of circulant matrices $M\in\mathbb{R}^{m\times m}$, where each matrix entry is a rational function in $\theta$ and $\nu$. Here, $\nu$ denotes the CFL number, being proportional to the ratio between the temporal and spatial discretization step sizes. The entrywise non-negativity of the matrix $M$ -- which is equivalent to the positivity preservation of the fully discrete scheme -- is investigated via discrete Fourier analysis and also by solving some low-order parametric linear recursions. We find that positivity preservation of the fully discrete system is impossible if the number of spatial grid points $m$ is even. However, it turns out that positivity preservation of the fully discrete system is recovered for \emph{odd} values of $m$ provided that $\theta\ge 1/2$ and $\nu$ are chosen suitably. These results are interesting since the systems of ordinary differential equations obtained via the spatial semi-discretizations studied are \emph{not} positivity preserving.
  • Rapid Evolution of SARS-CoV-2 Challenges Human Defenses

    Duarte, Carlos M.; Ketcheson, David I.; Eguíluz, Víctor; Agusti, Susana; Fernández-Gracia, Juan; Jamil, Tahira; Laiolo, Elisa; Gojobori, Takashi; Alam, Intikhab (Cold Spring Harbor Laboratory, 2021-03-28) [Preprint]
    Evolutionary ecology theory provides an avenue to anticipate the future behavior of SARS-CoV-2. Here we quantify the accelerating evolution of SARS-CoV-2 by tracking the SARS-CoV-2 mutation globally, with a focus on the Receptor Binding Domain (RBD) of the spike protein believed to determine infectivity. We estimate that 384 million people were infected by March 1st, 2021, producing up to 10 21 copies of the virus, with one new RBD variant appearing for every 600,000 human infections, resulting in approximately three new effective RBD variants produced daily. Doubling the number of RBD variants every 71.67 days followed by selection of the most infective variants challenges our defenses and calls for a shift to anticipatory, rather than reactive tactics. One-Sentence Summary Accelerating evolution of SARS-CoV-2 demands formulating universal vaccines and treatments based on big-data simulations of possible new variants.
  • RK-Opt: A package for the design of numerical ODE solvers

    Ketcheson, David I.; Parsani, Matteo; Grant, Zachary; Ahmadia, Aron; Ranocha, Hendrik (Journal of Open Source Software, The Open Journal, 2020-10-30) [Article]
    Ordinary and partial differential equations (ODEs and PDEs) are used to model many important phenomena. In most cases, solutions of these models must be approximated by numerical methods. Most of the relevant algorithms fall within a few classes of methods, with the properties of individual methods determined by their coefficients. The choice of appropriate coefficients in the design of methods for specific applications is an important area of research. RK-Opt is a software package for designing numerical ODE solvers with coefficients optimally chosen to provide desired properties. It is available from, with documentation at The primary focus of the package is on the design of Runge-Kutta methods, but some routines for designing other classes of methods such as multistep Runge-Kutta and general linear methods are also included.
  • General relaxation methods for initial-value problems with application to multistep schemes

    Ranocha, Hendrik; Loczi, Lajos; Ketcheson, David I. (Numerische Mathematik, Springer Nature, 2020-10-28) [Article]
    Recently, an approach known as relaxation has been developed for preserving the correct evolution of a functional in the numerical solution of initial-value problems, using Runge–Kutta methods. We generalize this approach to multistep methods, including all general linear methods of order two or higher, and many other classes of schemes. We prove the existence of a valid relaxation parameter and high-order accuracy of the resulting method, in the context of general equations, including but not limited to conservative or dissipative systems. The theory is illustrated with several numerical examples.
  • manuel-quezada/BP_Lim_for_RK_Methods: This repository contains the code to reproduce the results in "Bound-preserving convex limiting for high-order Runge-Kutta time discretizations of hyperbolic conservation laws"

    Kuzmin, Dmitri; Quezada de Luna, Manuel; Ketcheson, David I.; Grull, Johanna (Github, 2020-08-28) [Software]
    This repository contains the code to reproduce the results in "Bound-preserving convex limiting for high-order Runge-Kutta time discretizations of hyperbolic conservation laws"
  • Dirk schemes with high weak stage order

    Ketcheson, David I.; Seibold, Benjamin; Shirokoff, David; Zhou, Dong (Springer Nature, 2020-08-12) [Conference Paper]
    Runge-Kutta time-stepping methods in general suffer from order reduction: the observed order of convergence may be less than the formal order when applied to certain stiff problems. Order reduction can be avoided by using methods with high stage order. However, diagonally-implicit Runge-Kutta (DIRK) schemes are limited to low stage order. In this paper we explore a weak stage order criterion, which for initial boundary value problems also serves to avoid order reduction, and which is compatible with a DIRK structure. We provide specific DIRK schemes of weak stage order up to 3, and demonstrate their performance in various examples.
  • Effective rankine-hugoniot conditions for shock waves in periodic media

    Ketcheson, David I.; Quezada de Luna, Manuel (Communications in Mathematical Sciences, International Press of Boston, 2020-07-28) [Article]
    Solutions of first-order nonlinear hyperbolic conservation laws typically develop shocks infinite time even from smooth initial conditions. However, in heterogeneous media with rapid spatial variation, shock formation may be delayed or avoided. When shocks do form in such media, their speed of propagation depends on the material structure. We investigate conditions for shock formation and propagation in heterogeneous media. We focus on the propagation of plane waves in two-dimensional media with a periodic structure that changes in only one direction. We propose an estimate for the speed of the shocks that is based on the Rankine-Hugoniot conditions applied to a leading-order homogenized (constant coefficient) system. We verify this estimate via numerical simulations using different nonlinear constitutive relations and layered and smoothly varying media with a periodic structure. In addition, we discuss conditions and regimes under which shocks form in this type of media.
  • Relaxation Runge–Kutta Methods for Hamiltonian Problems

    Ranocha, Hendrik; Ketcheson, David I. (Journal of Scientific Computing, Springer Nature, 2020-07-09) [Article]
    The recently-introduced relaxation approach for Runge–Kutta methods can be used to enforce conservation of energy in the integration of Hamiltonian systems. We study the behavior of implicit and explicit relaxation Runge–Kutta methods in this context. We find that, in addition to their useful conservation property, the relaxation methods yield other improvements. Experiments show that their solutions bear stronger qualitative similarity to the true solution and that the error grows more slowly in time. We also prove that these methods are superconvergent for a certain class of Hamiltonian systems.
  • Riemann Problems and Jupyter Solutions

    Ketcheson, David I.; LeVeque, Randall J.; del Razo, Mauricio J. (Society for Industrial & Applied Mathematics (SIAM), 2020-06-26) [Book]
    This book addresses an important class of mathematical problems—the Riemann problem—for first-order hyperbolic partial differential equations (PDEs), which arise when modeling wave propagation in applications such as fluid dynamics, traffic flow, acoustics, and elasticity. The solution of the Riemann problem captures essential information about these models and is the key ingredient in modern numerical methods for their solution. This book covers the fundamental ideas related to classical Riemann solutions, including their special structure and the types of waves that arise, as well as the ideas behind fast approximate solvers for the Riemann problem. The emphasis is on the general ideas, but each chapter delves into a particular application. Riemann Problems and Jupyter Solutions - is available in electronic form as a collection of Jupyter notebooks that contain executable computer code and interactive figures and animations, allowing readers to grasp how the concepts presented are affected by important parameters and to experiment by varying those parameters themselves; - is the only interactive book focused entirely on the Riemann problem; and - develops each concept in the context of a specific physical application, helping readers apply physical intuition in learning mathematical concepts. Graduate students and researchers working in the analysis and/or numerical solution of hyperbolic PDEs will find this book of interest. This includes mathematicians, as well as scientists and engineers, working with applications like fluid dynamics, water waves, traffic modeling, or electromagnetism. Educators interested in developing instructional materials using Jupyter notebooks will also find this book useful. The book is appropriate for courses in Numerical Methods for Hyperbolic PDEs and Analysis of Hyperbolic PDEs, and it can be a great supplement for courses in computational fluid dynamics, acoustics, and gas dynamics.
  • Assessing the age- and gender-dependence of the severity and case fatality rates of COVID-19 disease in Spain

    Moraga, Paula; Ketcheson, David I.; Ombao, Hernando; Duarte, Carlos M. (Wellcome Open Research, F1000 Research Ltd, 2020-06-02) [Article]
    Background: The assessment of the severity and case fatality rates of coronavirus disease 2019 (COVID-19) and the determinants of its variation is essential for planning health resources and responding to the pandemic. The interpretation of case fatality rates (CFRs) remains a challenge due to different biases associated with surveillance and reporting. For example, rates may be affected by preferential ascertainment of severe cases and time delay from disease onset to death. Using data from Spain, we demonstrate how some of these biases may be corrected when estimating severity and case fatality rates by age group and gender, and identify issues that may affect the correct interpretation of the results. Methods: Crude CFRs are estimated by dividing the total number of deaths by the total number of confirmed cases. CFRs adjusted for preferential ascertainment of severe cases are obtained by assuming a uniform attack rate in all population groups, and using demography-adjusted under-ascertainment rates. CFRs adjusted for the delay between disease onset and death are estimated by using as denominator the number of cases that could have a clinical outcome by the time rates are calculated. A sensitivity analysis is carried out to compare CFRs obtained using different levels of ascertainment and different distributions for the time from disease onset to death. Results: COVID-19 outcomes are highly influenced by age and gender. Different assumptions yield different CFR values but in all scenarios CFRs are higher in old ages and males. Conclusions: The procedures used to obtain the CFR estimates require strong assumptions and although the interpretation of their magnitude should be treated with caution, the differences observed by age and gender are fundamental underpinnings to inform decision-making.

View more