Abstract:Suitable discretizations through tensor product formulas of popular multidimensional operators (diffusion or diffusion--advection, for instance) lead to matrices with $d$-dimensional Kronecker sum structure. For evolutionary Partial Differential Equations containing such operators and integrated in time with exponential integrators, it is then of paramount importance to efficiently approximate the actions of $\varphi$-functions of the arising matrices. In this work, we show how to produce directional split approximations of third order with respect to the time step size. They conveniently employ tensor-matrix products (the so-called $\mu$-mode product and related Tucker operator, realized in practice with high performance level 3 BLAS), and allow for the effective usage of exponential Runge--Kutta integrators up to order three. The technique can also be efficiently implemented on modern computer hardware such as Graphic Processing Units. The approach has been successfully tested against state-of-the-art techniques on two well-known physical models that lead to Turing patterns, namely the 2D Schnakenberg and the 3D FitzHugh--Nagumo systems, on different hardware and software architectures.
What problem does this paper attempt to address?
The main problem that this paper attempts to solve is to improve the efficiency of exponential integrators when solving multi - dimensional stiff ordinary differential equations (ODEs) systems with Kronecker sum structures. Specifically, the author aims to develop an efficient directional splitting approximation method to approximate the action of the φ - function of the matrix with third - order accuracy, so that high - order exponential Runge - Kutta integrators can be more effectively applied to these systems.
### Problem Background
Many partial differential equations (PDEs) generate stiff ODEs systems with Kronecker sum structures after spatial discretization. For example:
\[ u'(t) = Ku(t) + g(t, u(t)) = f(t, u(t)), \quad u(0) = u_0, \]
where \( K \in \mathbb{C}^{N \times N} \) is a matrix with d - dimensional Kronecker sum structure:
\[ K = A_d \oplus \cdots \oplus A_1 = A \otimes I + \cdots + I \otimes A. \]
For such systems, when performing time integration using exponential integrators, it is necessary to efficiently calculate the action of the φ - function of the matrix. In particular, for third - order exponential Runge - Kutta integrators, how to efficiently implement directional splitting approximation is a key issue.
### Solution
The author proposes a directional splitting method based on tensor operations, which can achieve a third - order accurate φ - function approximation. This method utilizes tensor - matrix products (such as µ - mode products and Tucker operators) and achieves efficient computation through high - performance BLAS libraries. Specifically, the author derives the following formula:
\[ \varphi_\ell(\tau K)v=\ell! d^{-1}\left(\varphi_\ell(\tau A_d)\otimes\cdots\otimes\varphi_\ell(\tau A_1)\right)v + O(\tau^2). \]
In addition, the author also proposes two different splitting methods:
1. **Two - step two - dimensional splitting**: Applicable to two - dimensional cases and implemented through two Tucker operators.
2. **Three - step multi - dimensional splitting**: Applicable to higher - dimensional cases and implemented through three Tucker operators, and all coefficients are real numbers.
### Application Examples
To verify the effectiveness of the proposed method, the author conducted numerical experiments on different hardware architectures, including CPUs and GPUs. Specifically, two well - known physical models - the Schnakenberg model and the FitzHugh - Nagumo model, which can generate Turing patterns, were applied. The experimental results show that the new method significantly improves computational efficiency while maintaining high accuracy.
### Conclusion
By introducing an efficient directional splitting approximation method, the author successfully realizes the application of third - order exponential integrators in solving stiff ODEs systems with Kronecker sum structures. This method not only improves computational efficiency but can also be efficiently implemented on modern computer hardware (such as GPUs).