Numerical Mathematics Group
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 Convectiondiffusion Equations(Journal of Scientific Computing, Springer Science and Business Media LLC, 20220801) [Article]We provide a framework for highorder discretizations of nonlinear scalar convectiondiffusion 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 maximumprinciplepreserving (MPP) with no step size restriction. The schemes are based on a twotiered limiting strategy, starting with a highorder limiterbased method that may have small oscillations or maximumprinciple 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(arXiv, 20220424) [Preprint]RungeKutta (RK) methods may exhibit order reduction when applied to certain stiff problems. While fully implicit RK schemes exist that avoid order reduction via highstage order, DIRK (diagonally implicit RungeKutta) 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 highorder DIRK schemes with WSO 4 and above. The resulting DIRK schemes are stiffly accurate, Lstable, have optimized error coefficients, and are demonstrated to perform well on a portfolio of relevant ODE and PDE test problems.

Rapid evolution of SARSCoV2 challenges human defenses.(Scientific reports, Springer Science and Business Media LLC, 20220419) [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 SARSCoV2, which can be monitored through unprecedented realtime tracking of SARSCoV2 population genomics at the global scale. Here we quantify the accelerating evolution of SARSCoV2 by tracking the SARSCoV2 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 RungeKutta Methods(arXiv, 20220407) [Preprint]RungeKutta (RK) methods may exhibit order reduction when applied to stiff problems. For linear problems with timeindependent 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 fullyimplicit 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.

Boundpreserving Flux Limiting for HighOrder Explicit RungeKutta Time Discretizations of Hyperbolic Conservation Laws(JOURNAL OF SCIENTIFIC COMPUTING, Springer Science and Business Media LLC, 20220304) [Article]We introduce a general framework for enforcing local or global maximum principles in highorder spacetime discretizations of a scalar hyperbolic conservation law. We begin with sufficient conditions for a space discretization to be bound preserving (BP) and satisfy a semidiscrete maximum principle. Next, we propose a global monolithic convex (GMC) flux limiter which has the structure of a fluxcorrected transport (FCT) algorithm but is applicable to spatial semidiscretizations 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 highorder RK method using GMCtype 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 Bseries(arXiv, 20211123) [Preprint]We present BSeries.jl, a Julia package for the computation and manipulation of Bseries, which are a versatile theoretical tool for understanding and designing discretizations of differential equations. We give a short introduction to the theory of Bseries 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 Bseries of high order.

On the rate of error growth in time for numerical solutions of nonlinear dispersive wave equations(Partial Differential Equations and Applications, Springer Science and Business Media LLC, 20211018) [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.

manuelquezada/BP_Lim_for_imp_RK_Methods:(Github, 20210915) [Software]

A review of matrix sir arino epidemic models(Mathematics, MDPI AG, 20210628) [Article]Many of the models used nowadays in mathematical epidemiology, in particular in COVID19 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 SIRtype epidemics”. However, many papers continue to reprove them in particular instances. This motivated us to redraw attention to these basic laws and provide a selfcontained reference of related formulas for (x, y, z) models. For the case of one susceptible class, we propose to use the name SIRPH, due to a simple probabilistic interpretation as SIR models where the exponential infection time has been replaced by a PHtype distribution. Note that to each SIRPH 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 finitetime nonpharmaceutical intervention(Journal of Mathematical Biology, Springer Science and Business Media LLC, 20210626) [Article]We consider the problem of controlling an SIRmodel 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 longtime limit, subject to some cost function for the control. We first consider the nocost 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 HamiltonJacobiBellman equation. Finally, we provide some examples related directly to the current pandemic.

Positivity preservation of implicit discretizations of the advection equation(arXiv, 20210516) [Preprint]We analyze, from the viewpoint of positivity preservation, certain discretizations of a fundamental partial differential equation, the onedimensional advection equation with periodic boundary condition. The full discretization is obtained by coupling a finite difference spatial semidiscretization (the second and some higherorder 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 secondorder method by choosing $\theta\in [0,1]$ suitably). The full discretization generates a twoparameter 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 nonnegativity 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 loworder 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 semidiscretizations studied are \emph{not} positivity preserving.

Rapid Evolution of SARSCoV2 Challenges Human Defenses(Cold Spring Harbor Laboratory, 20210328) [Preprint]Evolutionary ecology theory provides an avenue to anticipate the future behavior of SARSCoV2. Here we quantify the accelerating evolution of SARSCoV2 by tracking the SARSCoV2 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. OneSentence Summary Accelerating evolution of SARSCoV2 demands formulating universal vaccines and treatments based on bigdata simulations of possible new variants.

RKOpt: A package for the design of numerical ODE solvers(Journal of Open Source Software, The Open Journal, 20201030) [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. RKOpt is a software package for designing numerical ODE solvers with coefficients optimally chosen to provide desired properties. It is available from https://github.com/ketch/RKOpt, with documentation at http://numerics.kaust.edu.sa/RKOpt/. The primary focus of the package is on the design of RungeKutta methods, but some routines for designing other classes of methods such as multistep RungeKutta and general linear methods are also included.

General relaxation methods for initialvalue problems with application to multistep schemes(Numerische Mathematik, Springer Nature, 20201028) [Article]Recently, an approach known as relaxation has been developed for preserving the correct evolution of a functional in the numerical solution of initialvalue 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 highorder 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.

manuelquezada/BP_Lim_for_RK_Methods: This repository contains the code to reproduce the results in "Boundpreserving convex limiting for highorder RungeKutta time discretizations of hyperbolic conservation laws"(Github, 20200828) [Software]This repository contains the code to reproduce the results in "Boundpreserving convex limiting for highorder RungeKutta time discretizations of hyperbolic conservation laws"

Dirk schemes with high weak stage order(Springer Nature, 20200812) [Conference Paper]RungeKutta timestepping 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, diagonallyimplicit RungeKutta (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 rankinehugoniot conditions for shock waves in periodic media(Communications in Mathematical Sciences, International Press of Boston, 20200728) [Article]Solutions of firstorder 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 twodimensional 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 RankineHugoniot conditions applied to a leadingorder 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(Journal of Scientific Computing, Springer Nature, 20200709) [Article]The recentlyintroduced 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(Society for Industrial & Applied Mathematics (SIAM), 20200626) [Book]This book addresses an important class of mathematical problems—the Riemann problem—for firstorder 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 genderdependence of the severity and case fatality rates of COVID19 disease in Spain(Wellcome Open Research, F1000 Research Ltd, 20200602) [Article]Background: The assessment of the severity and case fatality rates of coronavirus disease 2019 (COVID19) 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 demographyadjusted underascertainment 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: COVID19 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 decisionmaking.