pith. sign in

arxiv: 2604.11558 · v1 · submitted 2026-04-13 · 🧮 math.NA · cs.NA

A tensor-based exponential integrator for diffusion--reaction equations in common curvilinear coordinates

Marco Caliari , Fabio Cassini This is my paper

Pith reviewed 2026-05-10 14:52 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords diffusion-reaction equationsexponential integratorstensor methodscurvilinear coordinatesfinite differencesKronecker productsTuring patterns
0
0 comments X

The pith

A novel split exponential Euler method with tensor computations efficiently integrates stiff diffusion-reaction equations on curvilinear coordinate domains.

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

The authors present a tensor-based numerical method for solving diffusion-reaction equations on domains such as disks, balls, and cylinders using polar, spherical, or cylindrical coordinates. They discretize the Laplace operators with finite differences to obtain matrices expressed as sums of Kronecker products. A split variant of the exponential Euler method is introduced to manage stiffness while using mu-mode tensor-matrix products to apply the required phi1 matrix function efficiently. This setup enables stable simulations of large systems, up to a million degrees of freedom, that produce Turing patterns.

Core claim

By choosing a finite difference discretization that represents the Laplace operators as sums of Kronecker products in common curvilinear coordinates, and employing a split exponential Euler time integrator, the action of the phi1 function can be computed via tensor-matrix products in a mu-mode framework, avoiding the time step restrictions of explicit methods for stiff problems.

What carries the argument

Split exponential Euler method whose phi1 action is evaluated through mu-mode tensor-matrix products on Kronecker-sum structured operators from curvilinear finite differences.

Load-bearing premise

The finite difference discretization in the chosen curvilinear coordinates must produce Laplace operators that are exactly representable as sums of Kronecker products.

What would settle it

Running the method on a small 2D disk grid and comparing the computed solution after one step to a direct matrix exponential computation; mismatch in the phi1 application would indicate failure of the tensor approach.

Figures

Figures reproduced from arXiv: 2604.11558 by Fabio Cassini, Marco Caliari.

Figure 1
Figure 1. Figure 1: Results of the experiment in Section 5.1 for varying number of time [PITH_FULL_IMAGE:figures/full_fig_p016_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Results of the experiment in Section 5.2 using integrator (23). The [PITH_FULL_IMAGE:figures/full_fig_p017_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Results of the experiment in Section 5.3 using a suitable adaptation [PITH_FULL_IMAGE:figures/full_fig_p019_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Results of the experiment in Section 5.4 using the adaptation of [PITH_FULL_IMAGE:figures/full_fig_p020_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Results of the experiment in Section 5.5 using the adaptation of [PITH_FULL_IMAGE:figures/full_fig_p022_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Graphical representation of the spatial domain Ω for the experiment [PITH_FULL_IMAGE:figures/full_fig_p022_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Results for the r component of the experiment in Section 5.6 using the adaptation of integrators (20a) and (20d) to coupled systems. The DOF per direction are (nρ, nθ, nz) = (160, 160, 20) (total number of DOF 2 · 512000 + 2 · 25600), and the number of time steps is m = 8000. Left: r component at the final time. Right: evolution of integral mean ⟨r⟩. The total simulation wall-clock time is about 4 minutes.… view at source ↗
Figure 8
Figure 8. Figure 8: Results for the u component of the experiment in Section 5.6 using the adaptation of integrators (20a) and (20d) to coupled systems, see also [PITH_FULL_IMAGE:figures/full_fig_p024_8.png] view at source ↗
read the original abstract

In this paper, we study a tensor-based method for the numerical solution of a class of diffusion--reaction equations defined on spatial domains that admit common curvilinear coordinate representations. Typical examples in 2D include disks (polar coordinates), and in 3D balls or cylinders (spherical or cylindrical coordinates) as well as spheres for problems involving the Laplace--Beltrami operator. The proposed approach is based on a carefully chosen finite difference discretization of the Laplace operators that yields matrices with a structured representation as sums of Kronecker products. For the time integration, we introduce a novel split variant of the exponential Euler method that effectively handles the stiffness and avoids the severe time step size restriction of classical explicit methods. By exploiting the peculiar form of the obtained discretized operators and the chosen splitting strategy, we compute the needed action of the $\varphi_1$ matrix function through suitable tensor-matrix products in a $\mu$-mode framework. We demonstrate the efficiency the approach on a wide range of physically relevant 2D and 3D examples of coupled diffusion--reaction systems generating Turing patterns with up to $10^6$ degrees of freedom.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The manuscript introduces a tensor-based method for diffusion-reaction equations on domains admitting common curvilinear coordinates (e.g., polar, spherical, cylindrical). A finite-difference discretization is chosen so that the resulting Laplace operators admit a representation as sums of Kronecker products; a novel split variant of the exponential Euler method is then proposed whose phi_1 action is realized exactly by mu-mode tensor-matrix products. The approach is demonstrated on 2D and 3D Turing-pattern systems with up to 10^6 degrees of freedom.

Significance. If the Kronecker-sum structure and the splitting strategy are realized as claimed, the method supplies an efficient, structure-exploiting integrator that removes the explicit-method time-step restriction while remaining applicable to large-scale stiff reaction-diffusion problems on non-Cartesian domains. The reported experiments on Turing patterns are consistent with this efficiency gain.

major comments (2)
  1. [time integration / method statement] The central claim that the split exponential Euler method 'effectively handles the stiffness' rests on the chosen splitting and the exact phi_1 evaluation, yet no stability analysis, contractivity estimate, or error bound is supplied for the integrator (see the time-integration section and the statement of the method). Without such analysis the claim that the method avoids 'severe time step size restriction' cannot be assessed beyond the numerical demonstrations.
  2. [discretization section] The finite-difference discretization is asserted to produce Laplace operators as sums of Kronecker products that are preserved under the mu-mode framework, but the explicit construction for the curvilinear Laplace-Beltrami operator (especially on spheres) is not shown in sufficient detail to verify that the Kronecker structure survives the coordinate transformation and boundary conditions.
minor comments (2)
  1. [Abstract] Abstract, last sentence: 'We demonstrate the efficiency the approach' is missing the preposition 'of'.
  2. [numerical experiments] The numerical experiments report results up to 10^6 degrees of freedom but do not include quantitative error tables or direct runtime comparisons against a standard implicit or exponential integrator on the same meshes; adding such data would strengthen the efficiency claim.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading and constructive comments on our manuscript. We address each major comment point by point below and outline the revisions we will make.

read point-by-point responses

  1. Referee: [time integration / method statement] The central claim that the split exponential Euler method 'effectively handles the stiffness' rests on the chosen splitting and the exact phi_1 evaluation, yet no stability analysis, contractivity estimate, or error bound is supplied for the integrator (see the time-integration section and the statement of the method). Without such analysis the claim that the method avoids 'severe time step size restriction' cannot be assessed beyond the numerical demonstrations.

    Authors: We agree that a dedicated stability analysis would strengthen the presentation. The split exponential Euler method is constructed so that the linear diffusion operator (the source of stiffness) is treated exactly through the phi_1 function, while the nonlinear reaction term is handled explicitly; this is a standard strategy in exponential integrators for reaction-diffusion problems and removes the diffusive CFL restriction by design. The numerical experiments already illustrate stable behavior at time steps far larger than those permitted by explicit methods. To address the referee's point, the revised manuscript will include a brief discussion in the time-integration section that recalls the local error bound for the exponential Euler method and explains why the chosen splitting preserves the exact treatment of the linear part. A full contractivity analysis lies outside the scope of the present work, but we will note this limitation explicitly. revision: partial

  2. Referee: [discretization section] The finite-difference discretization is asserted to produce Laplace operators as sums of Kronecker products that are preserved under the mu-mode framework, but the explicit construction for the curvilinear Laplace-Beltrami operator (especially on spheres) is not shown in sufficient detail to verify that the Kronecker structure survives the coordinate transformation and boundary conditions.

    Authors: We thank the referee for highlighting this gap. In the revised manuscript we will expand the discretization section with the explicit finite-difference stencils for the Laplace-Beltrami operator in spherical coordinates, including the metric coefficients that appear after discretization. We will show, term by term, how the resulting matrix decomposes as a sum of Kronecker products and how standard boundary conditions (periodic or homogeneous Neumann) are imposed while retaining the Kronecker-sum structure. This added detail will allow direct verification that the mu-mode tensor-matrix products remain applicable. revision: yes

Circularity Check

0 steps flagged

No significant circularity; derivation exploits discretization structure without reduction to inputs

full rationale

The paper's core chain begins with a standard finite-difference discretization of Laplace operators in curvilinear coordinates, which produces the expected Kronecker-sum structure by the coordinate transformation itself. A split variant of the exponential Euler method is then defined so that the required ϕ₁ action reduces exactly to μ-mode tensor-matrix products on that structure. This is a direct algebraic consequence of the chosen splitting and the discretization, not a fitted parameter renamed as a prediction or a self-citation load-bearing premise. No uniqueness theorem, ansatz smuggling, or renaming of known results is invoked; the method is presented as a straightforward exploitation of the operator form. Experiments on Turing patterns serve as external validation rather than internal fitting. The derivation therefore remains self-contained against the stated assumptions.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

Abstract provides insufficient detail to enumerate all free parameters or axioms; the core premise is the existence of Kronecker-sum structure from the chosen finite-difference scheme in curvilinear coordinates.

axioms (1)
  • domain assumption Finite difference discretization in common curvilinear coordinates produces Laplace operators expressible as sums of Kronecker products.
    Invoked as the foundation for the tensor-based efficiency.

pith-pipeline@v0.9.0 · 5499 in / 1231 out tokens · 71493 ms · 2026-05-10T14:52:50.771798+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

38 extracted references · 38 canonical work pages

  1. [1]

    A. H. Al-Mohy and N. J. Higham. Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators.SIAM J. Sci. Comput., 33(2):488–511, 2011

    work page 2011

  2. [2]

    J. L. Arag´ on, M. Torres, D. Gil, R. A. Barrio, and P. K. Maini. Turing patterns with pentagonal symmetry.Phys. Rev. E, 65:051913, 2002

    work page 2002

  3. [3]

    U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton. Implicit-Explicit Methods for Time-Dependent Partial Differential Equations.SIAM J. Numer. Anal., 32(2):797–823, 1995

    work page 1995

  4. [4]

    Benzi and V

    M. Benzi and V. Simoncini. Approximation of functions of large matrices with Kronecker structure.Numer. Math., 135:1–26, 2017

    work page 2017

  5. [5]

    H. P. Bhatt, A. Q. M. Khaliq, and B. A. Wade. Efficient Krylov-based expo- nential time differencing method in application to 3D advection-diffusion- reaction systems.Appl. Math. Comput., 338:260–273, 2018

    work page 2018

  6. [6]

    Bialecki and L

    B. Bialecki and L. Wright. A fast direct solver for a fourth order finite dif- ference scheme for Poisson’s equation on the unit disc in polar coordinates. Numer. Algorithms, 70:727–751, 2015

    work page 2015

  7. [7]

    N. F. Britton.Reaction-Diffusion Equations and Their Applications to Biology. Academic Press, New York, NY, 1986

    work page 1986

  8. [8]

    Brugnano and D

    L. Brugnano and D. Trigiante. Tridiagonal matrices: Invertibility and conditioning.Linear Algebra Appl., 166:131–150, 1992. 26

    work page 1992

  9. [9]

    Caliari and F

    M. Caliari and F. Cassini. Direction splitting ofφ-functions in exponen- tial integrators ford-dimensional problems in Kronecker form.J. Approx. Softw., 1(1):1, 2024

    work page 2024

  10. [10]

    Caliari and F

    M. Caliari and F. Cassini. A second order directional split exponential in- tegrator for systems of advection–diffusion–reaction equations.J. Comput. Phys., 498:112640, 2024

    work page 2024

  11. [11]

    Caliari and F

    M. Caliari and F. Cassini. Matrix- and tensor-oriented numerical schemes for the evolutionary space-fractional complex Ginzburg–Landau equation. arXiv preprint arXiv:2510.21394, 2025

    work page arXiv 2025

  12. [12]

    Caliari, F

    M. Caliari, F. Cassini, L. Einkemmer, A. Ostermann, and F. Zivcovich. A µ-mode integrator for solving evolution equations in Kronecker form.J. Comput. Phys., 455:110989, 2022

    work page 2022

  13. [13]

    Caliari, F

    M. Caliari, F. Cassini, and F. Zivcovich. BAMPHI: Matrix-free and transpose-free action of linear combinations ofφ-functions from exponential integrators.J. Comput. Appl. Math., 423:114973, 2023

    work page 2023

  14. [14]

    Caliari, F

    M. Caliari, F. Cassini, and F. Zivcovich. Aµ-mode BLAS approach for mul- tidimensional tensor-structured problems.Numer. Algorithms, 92(4):2483– 2508, 2023

    work page 2023

  15. [15]

    F. Cassini. Efficient third-order tensor-oriented directional splitting for exponential integrators.Electron. Trans. Numer. Anal., 60:520–540, 2024

    work page 2024

  16. [16]

    Chen and D

    M. Chen and D. Kressner. Recursive blocked algorithms for linear systems with Kronecker product structure.Numer. Algorithms, 84(3):1199–1216, 2020

    work page 2020

  17. [17]

    M. C. D’Autilia, I. Sgura, and V. Simoncini. Matrix-oriented discretiza- tion methods for reaction–diffusion PDEs: Comparisons and applications. Comput. Math. Appl., 79(7):2067–2085, 2020

    work page 2067

  18. [18]

    Frittelli, A

    M. Frittelli, A. Madzvamuse, and I. Sgura. The bulk-surface virtual element method for reaction-diffusion PDEs: Analysis and applications.Commun. Comput. Phys., 33(3):733–763, 2023

    work page 2023

  19. [19]

    Frittelli, I

    M. Frittelli, I. Sgura, and B. Bozzini. Turing patterns in a 3D morpho- chemical bulk-surface reaction-diffusion system for battery modeling.Math. Eng., 6(2):363–393, 2024

    work page 2024

  20. [20]

    Gambino, M

    G. Gambino, M. C. Lombardo, G. Rubino, and M. Sammartino. Pattern selection in the 2D FitzHugh–Nagumo model.Ric. Mat., 68:535–549, 2019

    work page 2019

  21. [21]

    Gaudreault, G

    S. Gaudreault, G. Rainwater, and M. Tokman. KIOPS: A fast adap- tive Krylov subspace solver for exponential integrators.J. Comput. Phys., 372:236–255, 2018. 27

    work page 2018

  22. [22]

    Hairer and G

    E. Hairer and G. Wanner.Solving Ordinary Differential Equations II: Stiff and Differential Algeraic Problems, volume 14 ofSpringer Series in Com- putational Mathematics. Springer Berlin, second edition, 1996

    work page 1996

  23. [23]

    Hern´ andez, E

    D. Hern´ andez, E. C. Herrera-Hern´ andez, M. N´ u˜ nez L´ opez, and H. Hern´ andez-Coronado. Self-similar Turing patterns: An anomalous dif- fusion consequence.Phys. Rev. E, 95:022210, 2017

    work page 2017

  24. [24]

    Hochbruck and A

    M. Hochbruck and A. Ostermann. Exponential integrators.Acta Numer., 19:209–286, 2010

    work page 2010

  25. [25]

    Hundsdorfer and J

    W. Hundsdorfer and J. G. Verwer.Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, volume 33 ofSpringer Series in Computational Mathematics. Springer Berlin, Heidelberg, 2003

    work page 2003

  26. [26]

    T. G. Kolda and B. W. Bader. Tensor decompositions and applications. SIAM Rev., 51(3):455–500, 2009

    work page 2009

  27. [27]

    Kulkarni, D

    D. Kulkarni, D. Schmidt, and S.-K. Tsui. Eigenvalues of tridiagonal pseudo- Toeplitz matrices.Linear Algebra Appl., 297:63–80, 1999

    work page 1999

  28. [28]

    Lacitignola, B

    D. Lacitignola, B. Bozzini, M. Frittelli, and I. Sgura. Turing pattern for- mation on the sphere for a morphochemical reaction–diffusion model for electrodeposition.Commun. Nonlinear Sci. Numer. Simulat., 28:484–508, 2017

    work page 2017

  29. [29]

    M.-C. Lai. A note on finite difference discretizations for Poisson equation on a disk.Numer. Methods Partial Differ. Equ., 17(3):199–203, 2000

    work page 2000

  30. [30]

    A. V. L´ opez, D. Hern´ andez, C. G. Aguilar-Madera, R. C. Mart´ ınez, and E. C. Herrera-Hern´ andez. Boundary conditions influence on Turing pat- terns under anomalous diffusion: A numerical exploration.Phys. D, 470:134353, 2024

    work page 2024

  31. [31]

    V. T. Luan, J. A. Pudykiewicz, and D. R. Reynolds. Further development of efficient and accurate time integration schemes for meteorological models. J. Comput. Phys., 376:817–837, 2019

    work page 2019

  32. [32]

    Meurant.Hessenberg and Tridiagonal Matrices

    G. Meurant.Hessenberg and Tridiagonal Matrices. Other Titles in Applied Mathematics. SIAM, Philadelphia, PA, first edition, 2025

    work page 2025

  33. [33]

    J. D. Murray.Mathematical Biology II: Spatial Models and Biomedical Applications, volume 18 ofInterdisciplinary Applied Mathematics. Springer, third edition, 2003

    work page 2003

  34. [34]

    Schnakenberg

    J. Schnakenberg. Simple chemical reaction systems with limit cycle be- haviour.J. Theor. Biol., 81:389–400, 1979

    work page 1979

  35. [35]

    Singh, R

    S. Singh, R. C. Mittal, S. R. Thottoli, and M. Tamsir. High-fidelity simulations for Turing pattern formation in multi-dimensional Gray–Scott reaction-diffusion system.Appl. Math. Comput, 452:128079, 2023. 28

    work page 2023

  36. [36]

    Smoller.Shock waves and reaction–diffusion equations

    J. Smoller.Shock waves and reaction–diffusion equations. Grundlehren der mathematischen Wissenschaften. Springer, New York, NY, second edition, 1994

    work page 1994

  37. [37]

    C. F. Van Loan. The ubiquitous Kronecker product.J. Comput. Appl. Math., 123:85–100, 2000

    work page 2000

  38. [38]

    Varea, J

    C. Varea, J. L. Arag´ on, and R. A. Barrio. Turing patterns on a sphere. Phys. Rev. E, 60(4):4588–4592, 1999. 29

    work page 1999