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

MultipleRelaxation Runge Kutta Methods for Conservative Dynamical Systems(arXiv, 20230210) [Preprint]We generalize the idea of relaxation time stepping methods in order to preserve multiple nonlinear conserved quantities of a dynamical system by projecting along directions defined by multiple time stepping algorithms. Similar to the directional projection method of Calvo et. al., we use embedded RungeKutta methods to facilitate this in a computationally efficient manner. Proof of the accuracy of the modified RK methods and the existence of valid relaxation parameters are given, under some restrictions. Among other examples, we apply this technique to ImplicitExplicit RungeKutta time integration for the Kortewegde Vries equation and investigate the feasibility and effect of conserving multiple invariants for multisoliton solutions.

Computing with Bseries(ACM Transactions on Mathematical Software, Association for Computing Machinery (ACM), 20221202) [Article]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.

Very HighOrder Astable Stiffly Accurate Diagonally Implicit RungeKutta Methods(arXiv, 20221126) [Preprint]A numerical search approach is used to design highorder diagonally implicit RungeKutta (DIRK) schemes suitable for stiff and oscillatory systems. We present new Astable schemes of orders six (the highest order of previously designed DIRK schemes) up to eight. For each order, we include one scheme that is only Astable as well as one that is stiffly accurate and therefore Lstable. The stiffly accurate schemes require more stages but can be expected to give better results for highly stiff problems and differentialalgebraic equations. The development of eighthorder schemes requires the highly accurate numerical solution of a system of 200 equations in over 100 variables, which is accomplished via a combination of global and local optimization. The accuracy and stability of the schemes is analyzed and tested on diverse problems.

A Comparative Study of Iterative Riemann Solvers for the Shallow Water and Euler Equations(arXiv, 20220925) [Preprint]The Riemann problem for firstorder hyperbolic systems of partial differential equations is of fundamental importance for both theoretical and numerical purposes. Many approximate solvers have been developed for such systems; exact solution algorithms have received less attention because computation of the exact solution typically requires iterative solution of algebraic equations. Iterative algorithms may be less computationally efficient or might fail to converge in some cases. We investigate the achievable efficiency of robust iterative Riemann solvers for relatively simple systems, focusing on the shallow water and Euler equations. We consider a range of initial guesses and iterative schemes applied to an ensemble of test Riemann problems. For the shallow water equations, we find that Newton's method with a simple modification converges quickly and reliably. For the Euler equations we obtain similar results; however, when the required precision is high, a combination of Ostrowski and Newton iterations converges faster. These solvers are slower than standard approximate solvers like Roe and HLLE, but come within a factor of two in speed. We also provide a preliminary comparison of the accuracy of a finite volume discretization using an exact solver versus standard approximate solvers.

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.

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.