RSS 1.0Numerical 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
-
Multiple-Relaxation Runge Kutta Methods for Conservative Dynamical Systems(arXiv, 2023-02-10) [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 Runge-Kutta 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 Implicit-Explicit Runge-Kutta time integration for the Korteweg-de Vries equation and investigate the feasibility and effect of conserving multiple invariants for multi-soliton solutions.
-
Computing with B-series(ACM Transactions on Mathematical Software, Association for Computing Machinery (ACM), 2022-12-02) [Article]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.
-
Very High-Order A-stable Stiffly Accurate Diagonally Implicit Runge-Kutta Methods(arXiv, 2022-11-26) [Preprint]A numerical search approach is used to design high-order diagonally implicit Runge-Kutta (DIRK) schemes suitable for stiff and oscillatory systems. We present new A-stable 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 A-stable as well as one that is stiffly accurate and therefore L-stable. The stiffly accurate schemes require more stages but can be expected to give better results for highly stiff problems and differential-algebraic equations. The development of eighth-order 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, 2022-09-25) [Preprint]The Riemann problem for first-order 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 Convection-diffusion Equations(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(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.(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(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(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.
-
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, 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:(Github, 2021-09-15) [Software]
-
A review of matrix sir arino epidemic models(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(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(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(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(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 https://github.com/ketch/RK-Opt, with documentation at http://numerics.kaust.edu.sa/RK-Opt/. 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(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"(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(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(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.
