Table of Contents
Fetching ...

Improving the stability and efficiency of high-order operator-splitting methods

Siqi Wei, Victoria Guenter, Raymond J. Spiteri

TL;DR

This work addresses instability and efficiency bottlenecks in high-order operator-splitting methods for diffusion–reaction systems by developing design principles that optimize both accuracy and stability. It introduces a new four-stage, third-order, 2-split method OS$_{2}(4,3)7$ with seven sub-integrations that is augmented by stability-oriented coefficient design and operator ordering, and further enhances stability by replacing backward implicit steps with explicit sub-integrators when appropriate. Through Niederer cardiac electrophysiology benchmarks, the authors demonstrate up to ~36% improvements in efficiency over classical third-order methods, and show that the chosen operator ordering and sub-integrator choices can significantly affect stability and runtime. The results provide practical strategies to realize stable, efficient high-order operator-splitting methods in diffusion-dominated, multi-physics simulations, with potential applicability beyond cardiac models and to problems with similar linear stability constraints.

Abstract

Operator-splitting methods are widely used to solve differential equations, especially those that arise from multi-scale or multi-physics models, because a monolithic (single-method) approach may be inefficient or even infeasible. The most common operator-splitting methods are the first-order Lie--Trotter (or Godunov) and the second-order Strang (Strang--Marchuk) splitting methods. High-order splitting methods with real coefficients require backward-in-time integration in each operator and hence may be adversely impacted by instability for certain operators such as diffusion. However, besides the method coefficients, there are many other ancillary aspects to an overall operator-splitting method that are important but often overlooked. For example, the operator ordering and the choice of sub-integration methods can significantly affect the stability and efficiency of an operator-splitting method. In this paper, we investigate some design principles for the construction of operator-splitting methods, including minimization of local error measure, choice of sub-integration method, maximization of linear stability, and minimization of overall computational cost. We propose a new four-stage, third-order, 2-split operator-splitting method with seven sub-integrations per step and optimized linear stability for a benchmark problem from cardiac electrophysiology. We then propose a general principle to further improve stability and efficiency of such operator-splitting methods by using low-order, explicit sub-integrators for unstable sub-integrations. We demonstrate an almost 30\% improvement in the performance of methods derived from these design principles compared to the best-known third-order methods.

Improving the stability and efficiency of high-order operator-splitting methods

TL;DR

This work addresses instability and efficiency bottlenecks in high-order operator-splitting methods for diffusion–reaction systems by developing design principles that optimize both accuracy and stability. It introduces a new four-stage, third-order, 2-split method OS with seven sub-integrations that is augmented by stability-oriented coefficient design and operator ordering, and further enhances stability by replacing backward implicit steps with explicit sub-integrators when appropriate. Through Niederer cardiac electrophysiology benchmarks, the authors demonstrate up to ~36% improvements in efficiency over classical third-order methods, and show that the chosen operator ordering and sub-integrator choices can significantly affect stability and runtime. The results provide practical strategies to realize stable, efficient high-order operator-splitting methods in diffusion-dominated, multi-physics simulations, with potential applicability beyond cardiac models and to problems with similar linear stability constraints.

Abstract

Operator-splitting methods are widely used to solve differential equations, especially those that arise from multi-scale or multi-physics models, because a monolithic (single-method) approach may be inefficient or even infeasible. The most common operator-splitting methods are the first-order Lie--Trotter (or Godunov) and the second-order Strang (Strang--Marchuk) splitting methods. High-order splitting methods with real coefficients require backward-in-time integration in each operator and hence may be adversely impacted by instability for certain operators such as diffusion. However, besides the method coefficients, there are many other ancillary aspects to an overall operator-splitting method that are important but often overlooked. For example, the operator ordering and the choice of sub-integration methods can significantly affect the stability and efficiency of an operator-splitting method. In this paper, we investigate some design principles for the construction of operator-splitting methods, including minimization of local error measure, choice of sub-integration method, maximization of linear stability, and minimization of overall computational cost. We propose a new four-stage, third-order, 2-split operator-splitting method with seven sub-integrations per step and optimized linear stability for a benchmark problem from cardiac electrophysiology. We then propose a general principle to further improve stability and efficiency of such operator-splitting methods by using low-order, explicit sub-integrators for unstable sub-integrations. We demonstrate an almost 30\% improvement in the performance of methods derived from these design principles compared to the best-known third-order methods.
Paper Structure (16 sections, 1 theorem, 22 equations, 5 figures, 7 tables)

This paper contains 16 sections, 1 theorem, 22 equations, 5 figures, 7 tables.

Key Result

Theorem 2.1

Let $\{\alpha_{k}^{[\ell]}\}_{k=1,\dots,s}^{\ell=1,2}$ be the coefficients of an operator-splitting method $\Psi$ and let $\{{\alpha}_{k}^{\ast\,[\ell]}\}_{k=1,\dots,s}^{\ell=1,2}$ be the coefficients of the adjoint method $\Psi^\ast$. Let $RK_k^{[D]}$ ($RK_k^{\ast [D]}$) and $RK_k^{[R]}$ ($RK_k^{\a

Figures (5)

  • Figure 1: Examples of stability regions of FRSK methods that use one implicit and one explicit sub-integrator. Sub-regions relevant to stability of the method in practice when solving method-of-lines ODEs are shaded. \ref{['fig:explain_relevant_region1']}: When a stability region consists of multiple disjointed parts, only the part that contains the origin is relevant to the stability of the method in practice. \ref{['fig:explain_relevant_region2']}: If a hole of instability is present in the part of the stability region containing the origin, the right-most $x$-co-ordinate of the hole determines the largest stable step size in the presence of dominant negative real eigenvalues.
  • Figure 1: Stability regions of the OS$_{2}(4,3)7$$_\text{minLEM}$ method applied to the Niederer benchmark problem with different operator orderings.
  • Figure 2: Stability regions of the Ruth, AKS3, and OS$_{2}(4,3)7$$_{DR\hat{x}}$ methods applied to the Niederer benchmark problem with different operator orderings.
  • Figure 3: Stability regions of the Ruth and AKS3 applied to the Niederer benchmark problem with different operator orderings and FE applied to both negative steps.
  • Figure 4: Stability regions of the OS$_{2}(4,3)7$$_{DR\hat{x}}$ applied to the Niederer benchmark problem with different operator orderings and FE applied to both negative sub-integrations.

Theorems & Definitions (4)

  • Remark 2.1
  • Theorem 2.1
  • Proof 1
  • Remark 4.1