Pith. sign in

REVIEW 3 major objections 5 minor 34 references

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

T0 review · 3 major / 5 minor · reviewed 2026-08-10 · deepseek-v4-flash

Pith's one-line read This paper constructs a four-stage, third-order operator-splitting method with optimized linear stability that runs the Niederer cardiac benchmark about 36% more efficiently than Ruth, and about 29% faster when backward steps use forward…

desk verdict A useful, honest numerical study that delivers a new high-order splitting method and a practical stability trick, but its central efficiency claim rests on one benchmark the method was optimized for. read the letter →

arxiv 2501.02310 v1 pith:VWJAN7YJ submitted 2025-01-04 math.NA cs.NA

classification math.NAcs.NA MSC 65L0565L0665L20
keywords operatorsplittingfractional-stepmethodsRunge-Kuttalinearstabilityanalysisbackward-in-timeintegrationcardiacelectrophysiologyNiedererbenchmarkmethodoflines
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper is trying to establish that high-order operator-splitting methods, long avoided for problems like diffusion because third-order real-coefficient splittings require backward-in-time integration, can be made practical by treating stability and implementation choices as design variables. It shows that two design levers, choosing the operator ordering and replacing unstable implicit backward sub-integrations with forward Euler, together produce a method with a larger practical stability region and lower cost per step. On the Niederer cardiac benchmark, the proposed OS2(4,3)7DR xhat method achieves about a 36% efficiency gain over the Ruth method, and with forward Euler substitution it is about 29% faster than the best Ruth implementation while still meeting the error tolerance. If correct, high-order splitting methods become viable for stability-constrained reaction-diffusion problems.

What carries the argument

The argument is carried by the scalar linear stability function of a fractional-step Runge-Kutta method, R(z[1],z[2]) = product over stages of each sub-integrator's stability function evaluated at the splitting-scaled arguments. For reaction-diffusion systems with simultaneously diagonalizable Jacobians, the two variables collapse to one: with z = lambda_R $\Delta$ t and the fixed eigenvalue ratio lambda_D/lambda_R, the stability function becomes a one-variable product whose practical stability region is the intersection of |R(z)| <= 1 containing the origin with the negative real axis, and the largest usable step is estimated by the right-most negative x-intercept xhat. The design target is to minimize xhat, which yields OS2(4,3)7DR xhat when combined with a favorable operator ordering. The second mechanism is replacing an implicit backward sub-integrator with explicit forward Euler, whose polynomial stability function has no pole in the left half-plane, thereby removing holes of instability and reducing the cost per step.

What would settle it

Run the OS2(4,3)7DR xhat coefficients on a reaction-diffusion benchmark with a different diffusion-to-reaction eigenvalue ratio or with anisotropic diffusion; if the measured largest stable step falls below that of the Ruth method, or if the [MRMS]_v error exceeds 0.05 at the predicted step, the scalar-ratio design assumption is refuted.

Watch

Extended reading notes

Core claim

The central claim is that a four-stage, third-order, two-split fractional-step Runge-Kutta method with seven sub-integrations per step, denoted OS2(4,3)7DR xhat, has a larger practical linear stability interval on the negative real axis than the Ruth and AKS3 methods for the Niederer benchmark, whose Jacobian eigenvalue ratio lambda_D/lambda_R is about 1.92/1260 and makes the problem stability-constrained. The method's coefficients are optimized to push the right-most negative x-intercept xhat of |R(z)|=1 as far from the origin as possible, and the operator ordering is chosen so that the stiff reaction operator is integrated backward for as short an interval as possible. When the backward-in-time sub-integrations, which normally create poles and instability, are replaced by forward Euler, the poles are removed entirely, improving stability and cutting computational cost per step while errors rise but stay below the accepted threshold [MRMS]_v = 0.05. The measured consequence is a 36% efficiency gain over the Ruth method in standard form and a 29% speedup over the best Ruth implementation when forward Euler is used for the negative sub-steps.

Load-bearing premise

The load-bearing premise is that the scalar stability function with the fixed eigenvalue ratio lambda_D/lambda_R from the Niederer benchmark predicts the stable step for the full method-of-lines system, and that coefficients tuned to that ratio keep their advantage when the ratio or diagonalizability changes.

Editorial extensions

If this is right

  • Third-order operator splitting becomes a practical option for stability-constrained reaction-diffusion problems, where the conventional advice to avoid high-order methods because of backward-in-time integration is shown to be too pessimistic.
  • Operator ordering is not an ancillary detail: the paper shows that the ordering placing the backward integration of the stiff reaction operator over the shortest interval substantially increases the stable step size.
  • Replacing backward implicit sub-integrations with forward Euler can simultaneously improve stability and reduce cost, provided the problem is stability-constrained so the resulting accuracy loss remains within tolerance.
  • The OS2(4,3)7DR xhat method is about 36% more efficient than the Ruth method on the Niederer benchmark in standard form, and about 29% faster than the optimal Ruth implementation when forward Euler handles the negative steps.
  • AKS3's smaller local error measure does not translate into a practical advantage on this benchmark because the step size is limited by stability, not accuracy, supporting the paper's emphasis on stability optimization.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the scalar-ratio design mechanism generalizes, one could tune splitting coefficients to the stiffness ratio of any given problem, treating the eigenvalue ratio as a problem-specific design input rather than a fixed constant.
  • The forward Euler substitution formally lowers the order of the overall method, so for accuracy-constrained problems, rather than stability-constrained ones, the error increase could outweigh the speed gain; testing on a problem with a much stricter error tolerance would reveal the boundary of the strategy.
  • The stability analysis assumes simultaneously diagonalizable Jacobians; for non-normal or non-diagonalizable operators, the right-most negative x-intercept may not predict the practical stable step, so applying the optimized coefficients to such systems would test the limits of the design principle.
  • Because the optimized stability regions can contain holes, an adaptive-step implementation would need to guard against stepping into an unstable island; the paper's fixed-step results suggest but do not establish how adaptive controllers would behave.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper studies design choices for high-order fractional-step Runge–Kutta (FSRK) operator-splitting methods applied to reaction–diffusion systems. It reviews order conditions, local error measures, and stability functions, then proposes a new four-stage, third-order, two-split method with seven sub-integrations per step, OS2(4,3)7DRxhat, obtained by minimizing the right-most negative real-axis intercept of |R(z)|=1 for the eigenvalue ratio lambda_D/lambda_R approximately 1.92/1260 taken from the Niederer cardiac benchmark. It also proposes replacing backward-in-time implicit sub-integrations with forward Euler, and reports that on the Niederer benchmark the new method takes dt=0.011 versus dt=0.0062 for the best Ruth implementation, giving roughly 36% efficiency gain without forward Euler and 29% with forward Euler while satisfying the stated MRMS error threshold. A low-LEM method, OS2(4,3)7minLEM, is also constructed but is shown to underperform on this stability-constrained problem.

Significance. If validated, the paper gives a practical recipe for improving high-order splitting methods: choose operator ordering and sub-integrators to enhance linear stability, and use forward Euler for backward substeps to remove poles and cut cost. The order-condition and stability-function derivations are standard, the method coefficients are reported to 15 decimals, and the numerical experiments are clearly described, with CPU times as minima of three runs. The principal caveat is that the headline speedup is demonstrated for a single benchmark on which the coefficients were optimized, so the generality of the claimed design principle is asserted rather than independently tested. With an out-of-sample test or a clearly scoped claim, the paper would be a useful and citable contribution to practical operator-splitting methods.

major comments (3)
  1. [§3.2, §4.2–4.4] The coefficients of OS2(4,3)7DRxhat in Table 5 are obtained by minimizing the right-most negative x-intercept xhat for the fixed eigenvalue ratio lambda_D/lambda_R approximately 1.92/1260, and the same ratio and the same Niederer benchmark are then used to report the efficiency gains in Tables 6 and 7. The observed 36% and 29% gains therefore confirm the optimization objective rather than providing an independent validation of the proposed design principle. I recommend either adding a clearly out-of-sample test (for example, a different mesh, cell model, or eigenvalue ratio) or explicitly restricting the claim to the benchmark-specific optimized method.
  2. [§2.4, §4.4] The paper itself cautions that the scalar stability analysis based on equation (2.6) cannot be directly applied to the method-of-lines system and that stability plots must be interpreted in the sense of extremes. For the ten Tusscher–Panfilov model in the Niederer benchmark, the reaction Jacobian is state-dependent, so lambda_R = -1260 is not a constant during the action potential, and the diffusion eigenvalue lambda_D changes if the spatial mesh is refined. The stability advantage reported in Table 6 (dt = 0.011 for OS2 versus dt = 0.0062 for Ruth) is thus demonstrated for one mesh and one design linearization point. The manuscript should either test the robustness of the advantage over the range of eigenvalue ratios encountered during the simulation and under mesh refinement, or state the benchmark-specific scope of the stability claim more prominently.
  3. [§3.2] The label optimized linear stability rests on a numerical search using MATLAB GlobalSearch with 100 random initial guesses and a coefficient bound [-1,1]. No proof of global optimality is given, and no seed or search script is supplied. This is not fatal to the design recipe, but the paper should qualify optimal as best among the candidates found in a restricted search and report the search details needed for reproducibility.
minor comments (5)
  1. [§2.4] The caption of Figure 1 contains the typo FRSK; it should read FSRK.
  2. [§2.4] In Theorem 2.1, the notation RK∗[R]] contains a stray bracket; the theorem would also be easier to read if the factorization of the stability function were written out explicitly.
  3. [§4.2] The statement that the best OS2(3,3)6 method has an xhat only marginally more negative than Ruth is not quantified; reporting the xhat values for Ruth, AKS3, and the new methods would let the reader verify the predicted stable step sizes in Tables 6 and 7.
  4. [§4.5] The effective order of the overall method after replacing backward sub-steps with forward Euler is not analyzed; the paper states that the order is negatively impacted, but a short quantitative statement of the resulting order would be helpful.
  5. [§4.4] The speedup is tolerance-specific because both FE variants in Table 7 have errors near the threshold of 0.05 (0.0414 and 0.041); a small sweep of tolerances would strengthen the efficiency comparison.

Circularity Check

1 steps flagged · score 5.0 of 10

Reported Niederer efficiency gain restates the linear-stability optimization objective, because the benchmark's eigenvalue ratio fixes the design target and the computation is stability-constrained.

  1. fitted input called prediction [Sections 3.2, 4, and 4.4 (Table 6)]
    "by designing a method that minimizes x̂, the right-most negative x-intercept of |R(z)| = 1. When the Jacobians of the operators are simultaneously diagonalizable, this value is a good predictor of the largest stable step size in practice. ... For the Niederer benchmark problem, the most negative eigenvalues of the Jacobians of the diffusion and reaction operators are λ[D] ≈ −1.92 and λ[R] = −1260 [24], respectively. ... We note that this computation is stability constrained because the [MRMS] v errors for the allowed step sizes are well below 0.05."

    OS2(4,3)7DRx is selected by minimizing x̂ for the stability function (2.6), and that stability function is constructed with the fixed ratio λD/λR ≈ 1.92/1260 taken from the Niederer benchmark. The same benchmark is then used in Table 6 to measure the largest stable step size, and the paper explicitly states the computation is stability-constrained. Hence the measured 36% efficiency gain over Ruth is the optimization objective converted into a step size for the very eigenvalue ratio used in the design, rather than an independent out-of-sample validation of a general stability principle.

full rationale

The paper contains one significant in-sample design/evaluation loop: the splitting coefficients are optimized for the linear stability bound x̂ using the Niederer benchmark's eigenvalue ratio, and the headline efficiency comparison is then performed on that same stability-constrained benchmark. This fits the 'fitted input called prediction' pattern and explains the moderate circularity score. The remaining derivation is largely self-contained: the order conditions, LEM minimization, and FSRK stability-function product are standard mathematical developments, and the self-citations to prior work by the same authors (e.g., [26] for the stability-function formula, [24] for the eigenvalue −1260) either state reproducible model properties or standard formulas rather than importing an unverified uniqueness theorem. The paper also reports an honest negative result for LEM-optimized methods, which reduces the concern that the framework is rigged. Overall, the design principles have independent content, but the 36% benchmark gain should be read as an in-sample demonstration of the optimization target rather than as an independent prediction.

Assumptions & free parameters 4 free parameters · 5 assumptions · 0 invented entities

The central claim rests on the optimized coefficients (fitted to the benchmark's stiffness ratio), the scalar test-equation stability model, and the empirical assumption that forward Euler on backward steps retains acceptable accuracy. No new physical or mathematical entities are postulated.

free parameters (4)
  • OS2(4,3)7DR\hat{x} coefficients = Table 5 (e.g., alpha[1]_2=0.5115, alpha[2]_2=0.6687)
    Seven nonzero coefficients chosen by GlobalSearch to minimize the right-most negative x-intercept of the stability function for the Niederer benchmark's eigenvalue ratio and sub-integrators (SDIRK(2,3), RK3).
  • OS2(4,3)7minLEM coefficients = Table 3
    Seven nonzero coefficients chosen by GlobalSearch to minimize the local error measure (LEM) for third order. Reported but not used in final performance comparison.
  • Coefficient design bound [-1,1] = -1 to 1
    Hand-chosen restriction on sub-step sizes to avoid excessively large backward steps; affects the optimized coefficients.
  • Error tolerance [MRMS]_v = 0.05
    User-chosen acceptable error level used to determine max step sizes in efficiency comparison.
assumptions (5)
  • standard math Baker-Campbell-Hausdorff formula yields the order conditions for operator-splitting methods
    Used in Section 2.1 to derive order conditions (2.1).
  • domain assumption Jacobians of reaction and diffusion operators are simultaneously diagonalizable near the solution
    Section 2.4: this allows reducing the stability function to a single-variable function of the eigenvalue ratio.
  • domain assumption For method-of-lines ODEs, the largest stable step size is determined by the right-most negative x-intercept of |R(z)|=1 in the stability region containing the origin
    Section 2.4: used to convert stability regions into step-size predictions.
  • domain assumption The eigenvalue ratio λD/λR ≈ 1.92/1260 is representative for the Niederer benchmark
    Section 4: taken from Spiteri and Dean [24]; used to construct RDR and RRD.
  • ad hoc to paper Replacing backward-in-time sub-integrators with forward Euler does not destroy the practical accuracy for the tested regime
    Section 4.5: empirically observed (errors remain below tolerance), but no order or error analysis is provided.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Improving the stability and efficiency of high-order operator-splitting methods." pith.science (2026). https://pith.science/paper/VWJAN7YJ

@misc{pith2026250102310,
  author       = {Pith},
  title        = {Pith review of: Improving the stability and efficiency of high-order operator-splitting methods},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/VWJAN7YJ}},
  note         = {Machine review of arXiv:2501.02310}
}
read the original 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.

Figures

Figures reproduced from arXiv: 2501.02310 by the authors.

Figure 1
Figure 1. Examples of stability regions of FRSK methods that use one implicit and one [PITH_FULL_IMAGE:figures/full_fig_p010_1.png] view at source ↗
Figure 2
Figure 2. Stability regions of the OS2(4, 3)7minLEM method applied to the Niederer benchmark problem with different operator orderings. 4.2. Optimized linear stability. For the case of OS2(3, 3)6 methods, we find the best result is a method whose ˆx is marginally more negative than that of the Ruth method, which implies that it will allow a slightly larger step-size than that allowed by the Ruth method. However, because the c… view at source ↗
Figure 3
Figure 3. Stability regions of the Ruth, AKS3, and OS [PITH_FULL_IMAGE:figures/full_fig_p016_3.png] view at source ↗
Figures from the paper (2 more)
Figure 4
Figure 4. Figure 4: Stability regions of the Ruth and AKS3 applied to the Niederer benchmark [PITH_FULL_IMAGE:figures/full_fig_p018_4.png]
Figure 5
Figure 5. Figure 5: Stability regions of the OS2(4, 3)7DRxˆ applied to the Niederer benchmark problem with different operator orderings and FE applied to both negative sub￾integrations. [7] J. Cervi and R. J. Spiteri, A comparison of fourth-order operator splitting methods for cardiac sim…

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

34 extracted references · 33 canonical work pages

  1. [1]

    U. M. Ascher, S. J. Ruuth, and R. J. Spiteri , Implicit-explicit Runge–Kutta Methods for Time-dependent Partial Differential Equations , Applied Numerical Mathematics, 25 (1997), pp. 151–167

  2. [2]

    Auzinger and W

    W. Auzinger and W. Herfort , Local error structures and order conditions in terms of Lie elements for exponential splitting schemes, Opuscula Mathematica, 34 (2014), pp. 243–255

  3. [3]

    Auzinger, H

    W. Auzinger, H. Hofst¨atter, D. Ketcheson, and O. Koch, Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes , BIT Numerical Mathematics, 57 (2017), pp. 55–74

  4. [4]

    Blanes, F

    S. Blanes, F. Casas, and A. Murua , Splitting methods for differential equations , arXiv.org, (2024)

  5. [5]

    Cervi and R

    J. Cervi and R. J. Spiteri , High-order operator splitting for the bidomain and monodomain models, SIAM Journal on Scientific Computing, 40 (2018), pp. A769–A786

  6. [6]

    Cervi and R

    J. Cervi and R. J. Spiteri, High-order operator-splitting methods for the bidomain and mon- odomain models, in Mathematical and Numerical Modeling of the Cardiovascular System and Applications, Springer, 2018, pp. 23–40. 19 Fig. 5: Stability regions of the OS 2(4, 3)7DR ˆx applied to the Niederer benchmark problem with different operator orderings and FE ...

  7. [7]

    Cervi and R

    J. Cervi and R. J. Spiteri, A comparison of fourth-order operator splitting methods for cardiac simulations, Applied Numerical Mathematics, 145 (2019), pp. 227–235

  8. [8]

    Cooper, R

    J. Cooper, R. J. Spiteri, and G. R. Mirams , Cellular cardiac electrophysiology modeling with chaste and cellml , Frontiers in physiology, 5 (2015), pp. 511–511

Show all 34 references
  1. [9]

    S. K. Godunov , A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics , Matematicheskii Sbornik, 89 (1959), pp. 271–306

  2. [10]

    Goldman and T

    G. Goldman and T. J. Kaper , Nth-order operator splitting schemes and nonreversible sys- tems, SIAM Journal on Numerical Analysis, 33 (1996), pp. 349–367

  3. [11]

    Guenter, S

    V. Guenter, S. Wei, and R. J. Spiteri, pythos: A Python library for solving IVPs by operator splitting, https://arxiv.org/abs/2407.05475, https://arxiv.org/abs/https://arxiv.org/abs/ 2407.05475

  4. [12]

    Hairer, C

    E. Hairer, C. Lubich, and G. W anner , Geometric numerical integration: structure- preserving algorithms for ordinary differential equations , vol. 31, Springer Science & Busi- ness Media, 2006

  5. [13]

    Hairer, S

    E. Hairer, S. P. Nørsett, and G. W anner, Solving ordinary differential equations. I, vol. 8 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, second ed., 1993. Nonstiff problems

  6. [14]

    Hansen and A

    E. Hansen and A. Ostermann , High order splitting methods for analytic semigroups exist , BIT, 49 (2009), pp. 527–542

  7. [15]

    C. A. Kennedy and M. H. Carpenter , Additive Runge–Kutta schemes for convection– diffusion–reaction equations, Applied Numerical Mathematics, 44 (2003), pp. 139–181

  8. [16]

    D. I. Ketcheson and H. Ranocha , Computing with b-series , ACM transactions on mathe- matical software, 49 (2023), pp. 1–23

  9. [17]

    Kutta, Beitrag zur naherungsweisen integration totaler differentialgleichungen , Z

    W. Kutta, Beitrag zur naherungsweisen integration totaler differentialgleichungen , Z. Math. Phys., 46 (1901), pp. 435–453

  10. [18]

    G. I. Marchuk , On the theory of the splitting-up method. , in Numerical Solution of Partial Differential Equations-II, Academic Press, 1971, pp. 469 – 500

  11. [19]

    M. E. Marsh, S. Ziaratgahi, and R. J. Spiteri , The Secrets to the Success of the Rush– Larsen Method and its Generalizations , IEEE Transactions on Biomedical Engineering, 59 (2012), pp. 2506–2515

  12. [20]

    S. A. Niederer, E. Kerfoot, A. P. Benson, M. O. Bernabeu, O. Bernus, C. Bradley, E. M. Cherry, R. Clayton, F. H. Fenton, A. Gary, et al., Verification of cardiac tissue electrophysiology simulators using an N-version benchmark , Philosophical Transactions of the Royal Society ...

  13. [21]

    Plank, A

    G. Plank, A. Loewe, A. Neic, C. Augustin, Y.-L. Huang, M. A. Gsell, E. Kara- 20 belas, M. Nothstein, A. J. Prassl, J. S ´anchez, G. Seemann, and E. J. Vig- mond, The opencarp simulation environment for cardiac electrophysiology , Computer Methods and Programs in Biomedicine, 2...

  14. [22]

    R. D. Ruth , A canonical integration technique , IEEE Transactions on Nuclear Science, 30 (1983), pp. 2669–2671

  15. [23]

    Sandu and M

    A. Sandu and M. G ¨unther, A generalized-structure approach to additive Runge–Kutta meth- ods, SIAM Journal on Numerical Analysis, 53 (2015), pp. 17–42

  16. [24]

    R. J. Spiteri and R. C. Dean, Stiffness analysis of cardiac electrophysiological models, Annals of Biomedical Engineering, 38 (2010), pp. 3592–3604

  17. [25]

    R. J. Spiteri, A. Tavassoli, S. Wei, and A. Smolyakov, Practical 3-splitting beyond Strang, 2023, https://arxiv.org/abs/2302.08034

  18. [26]

    R. J. Spiteri and S. Wei , Fractional-step Runge–Kutta methods: Representation and linear stability analysis , Journal of computational physics, 476 (2023), p. 111900

  19. [27]

    Strang , On the construction and comparison of difference schemes , SIAM Journal on Numerical Analysis, 5 (1968), pp

    G. Strang , On the construction and comparison of difference schemes , SIAM Journal on Numerical Analysis, 5 (1968), pp. 506–517

  20. [28]

    Sundnes, G

    J. Sundnes, G. T. Lines, and X. Cai, Computing the Electrical Activity in the Heart, Springer, 2006

  21. [29]

    K. H. ten Tusscher and A. V. Panfilov, Alternans and spiral breakup in a human ventricular tissue model , American Journal of Physiology - Heart and Circulatory Physiology, 291 (2006), pp. H1088–H1100

  22. [30]

    H. F. Trotter , Approximation of semi-groups of operators. , Pacific Journal of Mathematics, 8 (1958), pp. 887–919

  23. [31]

    Tyson, L

    R. Tyson, L. G. Stern, and R. J. LeVeque , Fractional step methods applied to a chemotaxis model, J. Math. Biol., 41 (2000), pp. 455–475, https://doi.org/10.1007/s002850000038

  24. [32]

    Wei and R

    S. Wei and R. J. Spiteri , Qualitative property preservation of high-order operator splitting for the sir model , Applied Numerical Mathematics, 172 (2022), pp. 332–350

  25. [33]

    N. N. Yanenko , The method of fractional steps. The solution of problems of mathematical physics in several variables , Springer-Verlag, New York-Heidelberg, 1971. Translated from the Russian by T. Cheron. English translation edited by M. Holt

  26. [34]

    Yoshida, Construction of higher order symplectic integrators, Physics Letters A, 150 (1990), pp

    H. Yoshida, Construction of higher order symplectic integrators, Physics Letters A, 150 (1990), pp. 262–268

Pith tools

Reviewed August 10, 2026 · model on record in the stance chip above.