Splitting methods for constrained diffusion-reaction systems

We consider Lie and Strang splitting for the time integration of constrained partial differential equations with a nonlinear reaction term. Since such systems are known to be sensitive with respect to perturbations, the splitting procedure seems promising as we can treat the nonlinearity separately. This has some computational advantages, since we only have to solve a linear constrained system and a nonlinear ODE. However, Strang splitting suffers from order reduction which limits its efficiency. This is caused by the fact that the nonlinear subsystem produces inconsistent initial values for the constrained subsystem. The incorporation of an additional correction term resolves this problem without increasing the computational cost. Numerical examples including a coupled mechanical system illustrate the proven convergence results

Hydrodynamic Instability Simulations Using Front-Tracking with Higher-Order Splitting Methods

The Rayleigh-Taylor Instability (RTI) is an instability that occurs at the interface of a lighter density fluid pushing onto a higher density fluid in constant or time-dependent accelerations. The Richtmyer-Meshkov Instability (RMI) occurs when two fluids of different densities are separated by a perturbed interface that is accelerated impulsively, usually by a shock wave. When the shock wave is applied, the less dense fluid will penetrate the denser fluid, forming a characteristic bubble feature in the displacement of the fluid. The displacement will initially obey a linear growth model, but as time progresses, a nonlinear model is required. Numerical studies have been performed in the past to accurately approximate this nonlinear model. A techniques called front tracking has provided an enhanced resolution and zero numerical diffusion that is helpful with the sharp discontinuities of the fluid properties in simulations involving RTI and RMI. Weighted essentially non-oscillatory (WENO) finite difference schemes are used for accurate and precise results in both early and late time of fluid mixing simulations. In more traditional projects, WENO schemes utilized Lax-Friedrichs flux splitting. However, an alternative type of splitting developed by Gilbert Strang splits a two-dimensional problem into two one-dimensional problems that are easier and faster to solve. His splitting method was shown to achieve up to second-order accuracy. For this research, such a splitting method was derived for higher-order accuracy in three-dimensional problems. RTI simulations utilizing this newly derived model were used to incorporate front tracking technique, WENO, and operator splitting in a way that has not been done for a three-dimensional problem

Overcoming the order barrier two in splitting methods when applied to semilinear parabolic problems with non-periodic boundary conditions

In general, high order splitting methods suffer from an order reduction phenomena when applied to the time integration of partial differential equations with non-periodic boundary conditions. In the last decade, there were introduced several modifications to prevent the second order Strang Splitting method from such a phenomena. In this article, inspired by these recent corrector techniques, we introduce a splitting method of order three for a class of semilinear parabolic problems that avoids order reduction in the context of non-periodic boundary conditions. We give a proof for the third order convergence of the method in a simplified linear setting and confirm the result by numerical experiments. Moreover, we show numerically that the high order convergence persists for an order four variant of a splitting method, and also for a nonlinear source term

A First-Order Explicit-Implicit Splitting Method for a Convection-Diffusion Problem

We analyze a second-order in space, first-order in time accurate finite difference method for a spatially periodic convection-diffusion problem. This method is a time stepping method based on the first-order Lie splitting of the spatially semidiscrete solution. In each time step, on an interval of length k, of this solution, the method uses the backward Euler method for the diffusion part, and then applies a stabilized explicit forward Euler approximation on m >= 1 intervals of length k/m for the convection part. With h the mesh width in space, this results in an error bound of the form C(0)h(2) + C(m)k for appropriately smooth solutions, where C-m <= C\u27 + C-\u27\u27/m. This work complements the earlier study [V. Thomee and A. S. Vasudeva Murthy, An explicit- implicit splitting method for a convection-diffusion problem, Comput. Methods Appl. Math. 19 (2019), no. 2, 283-293] based on the second-order Strang splitting

Efficient time splitting schemes for the monodomain equation in cardiac electrophysiology

Approximating the fast dynamics of depolarization waves in the human heart described by the monodomain model is numerically challenging. Splitting methods for the PDE-ODE coupling enable the computation with very fine space and time discretizations. Here, we compare different splitting approaches regarding convergence, accuracy, and efficiency. Simulations were performed for a benchmark problem with the Beeler–Reuter cell model on a truncated ellipsoid approximating the left ventricle including a localized stimulation. For this configuration, we provide a reference solution for the transmembrane potential. We found a semi-implicit approach with state variable interpolation to be the most efficient scheme. The results are transferred to a more physiological setup using a bi-ventricular domain with a complex external stimulation pattern to evaluate the accuracy of the activation time for different resolutions in space and time

Towards New High-Order Operator Splitting Time-Integration Methods

Operator splitting (OS) methods represent a powerful strategy to solve an extensive range of mathematical models in the form of differential equations. They have a long and celebrated history, having been successfully used for well over half a century to provide efficient numerical solutions to challenging problems. In fact, OS methods are often the only viable way to solve many problems in practice. The simplest, and perhaps, most well-known OS methods are Lie--Trotter--Godunov and the Strang--Marchuk methods. They compute a numerical solution that is first-, and second-order accurate in time, respectively. OS methods can be derived by imposing order conditions using the Campbell--Baker--Hausdorff formula. It follows that, by setting the appropriate order conditions, it is possible to derive OS methods of any desired order. An important observation regarding OS methods with order higher than two is that, according to the Sheng--Suzuki theorem, at least one of their defining coefficients must be negative. Therefore, the time integration with OS methods of order higher than two has not been considered suitable to solve deterministic parabolic problems, because the necessary backward time integration would cause instabilities. Throughout this thesis, we focus our attention on high-order (i.e., order higher than two) OS methods. We successfully assess the convergence properties of some higher-order OS methods on diffusion-reaction problems describing cardiac electrophysiology and on an advection-diffusion-reaction problem describing chemical combustion. Furthermore, we compare the efficiency performance of higher-order methods to second-order methods. For all the cases considered, we confirm an improved efficiency performance compared to methods of lower order. Next, we observe how, when using OS and Runge--Kutta type methods to advance the time integration, we can construct a unique extended Butcher tableau with a similar structure to the ones describing Generalized Additive Runge--Kutta (GARK) methods. We define a combination of methods to be OS-GARK methods. We apply linear stability analysis to OS-GARK methods; this allows us to conveniently analyze the stability properties of any combination of OS and Runge--Kutta methods. Doing so, we are able to perform an eigenvalue analysis to understand and improve numerically unstable solutions

An assessment of numerical methods for cardiac simulation

Solving the mathematical models of the electrical activity in the heart is difficult mainly due to the complexity of the models required to capture the electrochemical details of the organ. A variety of mathematical models has been developed to describe the electrical activity of individual heart cells. Cardiac cells respond to an electrical stimulus, causing ions to flow across the cell membrane, changing the electrical potential difference between the interior and the exterior of the cell. Cardiac cell models describe the potential difference across the cell membrane, and depending on the complexity of the model, the ion concentrations and the movement of ions through the cell membrane. This thesis studies 37 cardiac cell models and compares the efficiency of the Forward Euler and Rush-Larsen methods, as well as two generalized Rush-Larsen methods (of order one and order two). The Backward Euler method is compared for six of the 37 models and type-insensitive methods are compared for four of the 37 models. From the results, it is determined that the Rush-Larsen method is the most efficient for moderately stiff models and that the generalized Rush-Larsen methods perform well on the stiff models. The type-insensitive methods are more efficient than the single methods for each of the four models considered. The bidomain model combines a cardiac cell model with two partial differential equations that model the propagation of the electrical activity throughout the entire heart. Simulations of the bidomain model are computationally expensive, necessitating improvements to the numerical methods used to solve the model. In order to determine whether the cell model results can be applied to the bidomain model, a one-dimensional simulation of the bidomain model is considered and solved using the Chaste software package developed by the Computational Biology Group at Oxford University Computing Laboratory, together with additions written for this thesis. The cellular electrical activity within the bidomain model is modelled by eight of the 37 cell models previously studied, chosen to represent a range of the available models. From these results, it is shown that the cell model results are directly applicable to the one-dimensional bidomain problem. The first-order generalized Rush-Larsen method drastically reduces computation time for the two stiffest models, and the Rush-Larsen method performs optimally for the moderately stiff models. One of the de facto standards, the Forward Euler method, is shown to perform poorly for seven of the eight one-dimensional bidomain simulations