Second-order time integrators with the Fourier spectral method in application to multidimensional space-fractional FitzHugh-Nagumo model

1.   Introduction

In the field of excitable media, FitzHugh ^[1] introduced the FitzHugh model as a mathematical model to understand the generation and propagation of electrical signals in neuronal dynamics. Initially, the FitzHugh model was formulated as a nonlinear ordinary differential equation. Later, Nagumo in ^[2] extended this model to create the FitzHugh-Nagumo (FHN) model. The FitzHugh-Nagumo model is a well known two components classical reaction-diffusion model that describes the dynamics of excitable systems such as neurons and cardiac tissues ^[3]. However, when heterogeneity and complex connectivity properties characterize the medium in which the transport phenomenon is observed on a very large number of scales, it cannot be adequately described using the classical FHN model. Hence, the use of space-fractional FHN (SFFHN) model becomes necessary. The SFFHN model is obtained by replacing the standard Laplacian operator in the FHN model by a fractional one. In this study, we consider the following SFFHN model.

with the periodic, homogeneous Neumann, or Dirichlet boundary conditions along with appropriate initial conditions

where $u$ is a normalized transmembrane potential, $v$ is a dimensionless time-dependent recovery variable, $D_u$ is the diffusion coefficient, $d$ is the dimension of the problem, $F_1(u, v) = u(u-\mu)(1-u) - v, F_2(u, v) = \epsilon(\beta u - \gamma v - \delta)$ are nonlinear functions that satisfy the local Lipschitz condition, $\epsilon, \mu, \beta, \gamma$, and $\delta$ are constants that describe the model behavior and $(-\Delta)^{\frac{\alpha}{2}}, 1 < \alpha\leq 2$ is thel fractional Laplacian. As this model often exhibits complex and nonlinear behavior, analytical solutions to this model are often not available or they are very challenging to obtain. This has led to the development of a variety of efficient and accurate numerical methods to simulate the propagation of electrical signals, with varying degrees of accuracy and computational efficiency. In recent years, there has been an increasing interest in the development of various-order numerical methods for the simulation of partial differential equation models such as the SFFHN model. To name a few of them, Liu et al. ^[4] developed an implicit operator splitting method with a shifted Grünwald-Letnikov approximation to solve the two-dimensional (2D) SFFHN model with zero Dirichlet boundary conditions. Bueno-Orovio et al. ^[5] proposed a backward Euler method along with a Fourier spectral method as accurate and efficient numerical stencils for the SFFHN model, with the advantage of fully diagonal discretization that can scale to any spatial dimension. Bu et al. in ^[6] introduced a Crank–Nicolson (CN) alternating direction implicit Galerkin finite element method for a 2D fractional FitzHugh–Nagumo monodomain model with homogeneous Dirichlet boundary conditions. In another study, Li and Song ^[7] introduced a second-order operator splitting spectral element method for solving the space fractional reaction-diffusion (SFRD) equations including the SFFHN. Liu and Zhang ^[8] developed a semi-alternating direction Galerkin finite element method to solve a 2D fractional FHN monodomain model on an irregular domain. Li et al. ^[9] used a fourth-order Runge-Kutta method along with a Fourier spectral method to solve the FitzHugh-Nagumo model with Riesz fractional-derivatives. Zhang et al. in ^[10] introduced a numerical scheme based on the CN method in time and a Legendre spectral method in space for the SFFHN model. In ^[11,12], authors proposed a second-order stabilized semi-implicit method and a second-order operator splitting Fourier spectral method respectively to solve SFRD equations, including the 2D SFFHN model. Wang et al. ^[13] proposed a ghost structure finite difference method for a fractional FHN monodomain model on a moving irregular domain. Che et al. in ^[14] introduced the Fourier transform for spatial discretization and the Runge-Kutta method for time discretization for the time dependent 2D SFFHN mode. Authors in ^[15] proposed a fast high-order method for the multidimensional SFRD equation with general boundary conditions, including the 2D SFFHN model.

Second-order numerical methods in particular have shown promise in providing increased accuracy and efficiency over first-order methods. However, the application of these methods to multidimensional (especially 3D) SFFHN models with general boundary conditions (periodic, Dirichlet and Neumann) remains relatively unexplored, despite the potential benefits that they could provide. Therefore, the primary objective of this research project is to investigate the effectiveness of second-order time integrators in combination with the Fourier spectral method for solving the multidimensional SFFHN model (1.1), with given initial and boundary conditions. We aim to explore the accuracy, stability, and computational efficiency of second-order numerical methods in capturing the dynamics of the SFFHN model for varying $\alpha$, comparing the results based on two testing parameters: Accuracy and CPU time consumption.

The paper is organized as follows: Section 2 describes the spatial discretization approach using the Fourier spectral method for different boundary conditions. Section 3 introduces four second-order time integrators for solving the SFFHN model. The linear stability analysis of the proposed methods is presented in Section 4. Section 5 provides numerical results and discussions to evaluate the accuracy and computational efficiency of the methods. Finally, Section 6 concludes the work and outlines potential avenues for future research.

2.   Fourier spectral method

There are various definitions available for the fractional Laplacian operator $(-\Delta)^{\frac{\alpha}{2}}$, and for more detailed information, readers are encouraged to refer to the references ^[16,17]. The selection of an appropriate numerical approximation method depends on the specific definition being employed. In this study, we define the fractional Laplacian $(-\Delta)^{\alpha/2}$ operator in terms of the spectral decomposition of the Laplacian $(-\Delta)$ operator.

Suppose the $d$ dimensional Laplacian $(-\Delta)$ has a complete set of orthonormal eigenfunctions $\phi_{\eta_1, \cdots, \eta_d}$ corresponding to eigenvalues $\lambda_{\eta_1, \cdots, \eta_d}$ on the bounded region $\Omega = [a, b]^d$ with periodic, homogeneous Dirichlet or homogeneous Neumann boundary conditions, that is $(-\Delta)\phi_{\eta_1, \cdots, \eta_d} = \lambda_{\eta_1, \cdots, \eta_d}\phi_{\eta_1, \cdots, \eta_d}$, then $(-\Delta)^{\alpha/2}$ is defined by

in which $u$ can be expressed as follows

where $\langle u, \phi_{\eta_1} \rangle = \int_{\Omega} u \phi_{\eta_1} dx$ for $d = 1$ and eigenvalues $\lambda_{\eta_1, \cdots, \eta_d}$ and eigenvectors $\phi_{\eta_1, \cdots, \eta_d}$ will depend on the provided boundary conditions:

Case 1: For periodic boundary conditions:

Case 2: For homogeneous Dirichlet boundary conditions:

Case 3 For homogeneous Neumann boundary conditions:

where $L = b-a.$ If we consider a finite number of the orthonormal trigonometric eigenfunctions $\lbrace\phi_{i}\rbrace$(equal to the number of discretization points $N$), then $(-\Delta)^{\frac{\alpha}{2}}u$ is approximated by a truncated series:

Efficient computation of the coefficients $\hat{u}_{\eta_1, \cdots, \eta_d}$ in Eq (2.1) and the inverse reconstruction of $u$ in physical space can be accomplished using fast and robust algorithms such as direct and inverse Discrete Sine Transforms for the homogeneous Dirichlet boundary, direct and inverse Discrete Cosine Transforms for the homogeneous Neumann boundary and direct and inverse Discrete Fourier Transforms for the periodic boundary data, as detailed in ^[18].

To illustrate the idea of the Fourier spectral method applied to the fractional Laplacian, let's consider the following equation on one-dimensional (1D) space for simplicity:

The following is the representation of Eq (2.2) in Fourier space by using the definition of the fractional Laplacian and applying the Fourier transform:

where $A = D(\lambda_{\eta_1})^\frac{\alpha}{2}, {\bf U}: = [u_1(t), u_2(t)\cdots, u_N(t)]^\intercal, {\bf F({\bf U}}) = [f(x_1, t, u), f(x_2, t, u), \cdots, f(x_N, t, u)]^\intercal$.

In the actual calculation, the choice of the mesh will depend on the specified boundary data, and we generate the grids in MATLAB as follows. For simplicity, we consider a 1D space. The 2D/three-dimensional (3D) case is defined analogously.

For periodic boundary conditions,

For homogeneous Neumann boundary conditions,

For homogeneous Dirichlet boundary conditions,

3.   Second-order time integrators

In this section, we briefly discuss four second-order time integration methods for solving the SFFHN model. To describe the methods, let us consider the following system of nonlinear ordinary differential equations obtained by approximating the fractional Laplacian in (2.2) with the Fourier spectral technique discussed above:

where $\hat{{\bf U}} = \mathcal{F}({\bf U})$ and $\hat{\mathbf{\mathcal{F}}} = \mathcal{F}({\bf F({\bf U})}).$ Let $k = t_{n+1}-t_n$ be the time step size, then using a variation of constant formula, we come up with the following formula:

The Eq (3.2) serves as the foundation for various time integrators, including the exponential time differential (ETD) and the integrating factor (IF), which vary based on how we approximate the matrix exponential function and nonlinear function.

3.1.   Extrapolated exponential time differencing backward Euler (EETD-BE) method

To get the first order ETD method, the simplest approximation to the integral is to impose that $\hat{\mathbf{\mathcal{F}}}$ is constant for $t\in [t_n, t_{n+1}]$ i.e., $\hat{\mathbf{\mathcal{F}}}\approx \hat{\mathbf{\mathcal{F}}}_n$. Denoting the semi-discrete approximation to $\hat{{\bf U}}(t_{n})$ by $\hat{{\bf U}}_{n}$ yields:

However, the effectiveness of this semi-discrete method hinges on discretizing the matrix exponential. By incorporating the $(0, 1)$-Pad$\acute{\mathrm{e}}$ approximation, the scheme becomes more practical and computationally efficient. This approximation facilitates a direct discretization of the matrix exponential and eliminates the need for explicit computation of the inverse matrix A.

Definition 3.1. The $(r+s)^{th}$ order rational Pad$\acute{\mathrm{e}}$ approximation to $e^{-z}$ is defined as:

where

Implementing $(0, 1)$-Pad$\acute{\mathrm{e}}$ approximation $R_{0, 1}(k{ A}) = (1+k{ A})^{-1}$ into Eq (3.3), yields the ETD backward Euler method (ETD-BE):

The ETD-BE method is only first order accurate in time. So, in order to enhance the temporal accuracy to second-order, we implement a local extrapolation approach proposed by Lawson and Morris ^[19]:

Considering Eq (3.4) written over a double time step $2k$ and denoting the result with $\hat{{\bf U}}^{(2)},$ we have:

Alternatively, if Eq (3.4) is computed over two single time steps and if we denote the result with $\hat{{\bf U}}^{(1)}$ we have:

As a result, Equations (3.5) and (3.6) represent two distinct ETD-BE schemes for computing the solution at double time steps. Neither (3.5) nor (3.6) achieves an accuracy of $O(k^2)$. Nonetheless, by expanding $\hat{{\bf U}}^{(1)}$ and $\hat{{\bf U}}^{(2)}$ to $O(k^2)$ and combining these expansions through the operation of subtracting $\hat{{\bf U}}^{(2)}$ from $2\hat{{\bf U}}^{(1)}$, the local extrapolation of ETD-BE is obtained ^[20]:

For the efficient implementation of EETD-BE, in the following Algorithm the fast Fourier transform (FFT) computations of the EETD-BE method are provided for the case of 3D with Neumann boundary conditions:

3.2.   Exponential time differencing Crank-Nicolson (ETD-CN) method

Using ETD-BE as a prediction of $\hat{{\bf U}}_{n+1}$ and approximating the nonlinear term $\hat{\mathbf{\mathcal{F}}}$ ^[21] in (3.2) with:

where $\hat{{\bf a}}_{n} = e^{-k{\bf A}}\hat{{\bf U}}_n-{ A}^{-1}\left(e^{-k{\bf A}}-1\right)\hat{\mathbf{\mathcal{F}}}_n$ yields

Replacing $e^{-kA}$ in (3.8) with $(1, 1)$-Pad$\acute{\mathrm{e}}$ approximation $R_{1, 1}(kA) = (2-k{ A})(2+kA)^{-1}$ results the ETD-CN method:

where

For the efficient implementation of ETD-CN, the FFT computations of the method are provided for the case of 3D with periodic boundary conditions:

3.3.   ETD-Pad$\acute{\mathrm{e}}\; (0, 2)$

Replacing $e^{-kA}$ in (3.8) with $(0, 2)$-Pad$\acute{\mathrm{e}}$ approximation $R_{0, 2}(kA) = 2\left(2+2kA+(kA)^2\right)^{-1}$ yields the ETD-Pad$\acute{\mathrm{e}}\; (0, 2)$ method:

where

For the efficient implementation of ETD-Pad$\acute{\mathrm{e}}\; (0, 2)$, the FFT computations of the method are provided for the case of 3D with Dirichlet boundary conditions:

3.4.   IIF2 method

To construct the second-order implicit integration factor (IIF2) method ^[22], we approximate $e^{kA\tau} \hat{\mathbf{\mathcal{F}}}({\bf U}(t_n+k\tau), t_n+k\tau)$ using second-order Lagrange interpolation polynomial involving points $t_{n+1}, t_n,$ and $t_{n-1}$:

With the above approximation, the direct evaluation of the integral in (3.2) leads to the IIF2 method given as below:

For the efficient implementation of the IIF2 method, the FFT computations of the method are provided for the case of 3D with Dirichlet boundary conditions:

4.   Stability regions

We consider the following nonlinear autonomous ordinary differential equation (ODE):

In order to get the stability regions of an EETD-BE, we linearize Eq (4.1) about a fixed point $u_0$ such that $-cu_0+F(u_0) = 0$ to obtain:

where $u$ is perturbation of $u_0$ and $\lambda = F^{'}(u_0)$. In order to keep a fixed point $u_0$ stable we need to have Re$(\lambda-c) < 0$ (note that the fixed points of the ETD schemes are the same as those of the ODE (4.1), in contrast to the IF methods, which do not preserve the fixed points for the ODE that they discretize) ^[23].

The boundaries of the stability regions (a family of curves for different values of $ck$, where $k$ is the time step) based on Eq (4.2) are presented below for EETD-BE, ETD-CN and ETD-Pad$\acute{\mathrm{e}}\; (0, 2)$.

To obtain the stability regions, we applied EETD-BE, ETD-CN and ETD-Pad$\acute{\mathrm{e}}\; (0, 2)$ to (4.2), which leads to the following amplification factors:

4.1.   Amplification factor of EETD-BE

where

4.2.   Amplification factor of ETD-CN

where

4.3.   Amplification factor of ETD-Pad$\acute{\mathrm{e}}\; (0, 2)$

where

In general, the parameters $c$ and $\lambda$ may both be complex-valued. The stability region is four dimensional (4D) and, therefore, difficult to represent ^[23]. However, if both $c$ and $\lambda$ are pure imaginary ^[24] or real ^[23] or if $\lambda$ is complex and $c$ is fixed and real ^[25], then the stability region is 2D. The boundaries of the stability region are obtained by substituting $r = e^{i\theta}, \ \theta\in [0, 2\pi]$ into Eq (4.2) and solving for $x$. In order to plot a 2D stability region in the complex $x-$plane, where the solution is bounded as $n\to \infty$, we start our analysis by assuming $c$ to be fixed and real and $\lambda$ as complex.

Figure 1 illustrates the stability regions of the methods for different values of $y = 0, -3, -5, -10, -15,$ and $-20$ in the complex $x$-plane. According to Beylkin et al. ^[25], it is crucial for a method to have stability regions that grow as $ck$ becomes larger to be useful. From the stability regions shown in Figure 1, we observed that as $y\rightarrow -\infty$, the methods maintain their shape and exhibit larger stability regions. This characteristic enables the methods to utilize larger time-step sizes as $|c|\to \infty)$ while solving the problem. This observation indicates the stability of the methods EETD-BE, ETD-CN and ETD-Padé$(0, 2)$. Comparing the stability regions in Figure 1, we can note that ETD-CN exhibits slightly better stability than EETD and ETD-Padé$(0, 2)$, as its stability region encompasses larger values along the real axis.

4.4.   Amplification factor of IIF2

For the linear stability analysis and amplification factor of the IIF2 method, readers are referred to the work of ^[22]. In their study, the authors stated that the IIF2 method is unconditionally stable and exhibits excellent stability properties when used to solve stiff reaction systems.

5.   Numerical experiments and discussions

In this section, five test problems are considered to demonstrate the effectiveness of the methods based on the two key testing parameters: Accuracy and CPU time consumption for the integrating SFFHN model with different fractional powers $\alpha$, different initial conditions along with periodic, homogeneous Dirichlet, or Neumann boundary conditions. The accuracy of the method is measured in terms of maximum error norm $\left\lVert{\cdot}\right\rVert_\infty$ and the empirical rates of convergence (in time direction) of the methods is computed utilizing the following formula:

where $E_k = \left\lVert{u_k-u_{2k}}\right\rVert_\infty$ and $E_{\frac{k}{2}} = \left\lVert{u_{\frac{k}{2}}-u_{k}}\right\rVert_\infty$ are maximum norm errors at $k$ and $\frac{k}{2}$. Throughout this study, we fixed the parameters $D_u = 10^{-4}, \epsilon = 0.01, \mu = 0.1, \beta = 0.1, \gamma = 1,$ and $\delta = 0,$ which are known to generate stable patterns in the system. All the computations are conducted on the MATLAB 2020b platform based on a MacBook Pro with 2.6 GHz 6-Core Intel Core i7 CPU and 16 GB RAM.

Example 1 (2D SFFHN model with Neumann Boundary): In this example, we consider the 2D SFFHN model with homogeneous Neumann boundary conditions, along with the following initial conditions on domain $[0, 2.5]^2$:

To test the order of accuracy of the methods mentioned in Section 3 for varying $\alpha$ along with Neumann boundary conditions, we simulated the solution profile of the component $u$ up to the final time $t = 10.$ In the experiment we set the spectral resolution $N = 256$ large enough to ensure that the errors primarily arise from the time integration process rather than spatial discretization and repeatedly halving the temporal step size in each simulation, starting with the initial value $k = 0.1$. The results of our experiments are illustrated in Figure 2 for various values of $\alpha = 2, 1.7$ and $1.5$. Figure 2 confirms second-order convergence for the ETD-CN, ETD-Pad$\acute{\mathrm{e}}\; (0, 2)$, IIF2, and EETD methods. These second-order methods effectively integrate the system for varying $\alpha$ values.

In Figure 3, we extended the comparison of the methods by analyzing their accuracy in relation to the CPU time. The results depicted in Figure 3 highlight a noticeable variation in CPU time among the methods, even at a comparable level of accuracy. Consequently, it can be concluded that the computational time required by each method is not similar, indicating the presence of significant differences in efficiency. The performance of the ETD-CN method closely resembles that of the ETD-Pad$\acute{\mathrm{e}}\; (0, 2)$, as evident from the identical errors that scale with the expected second-order as shown in Figure 2. However, it should be noted that the IIF2 method is the most computationally expensive due to the involvement of a nonlinear solver, while the EETD method is the least demanding in terms of computational resources.

Considering both accuracy and computational cost when selecting a second-order method for time stepping for Example 1, it becomes evident that the ETD-CN method stands out as the most efficient choice in this regard.

In this experiment, we utilized the ETD-CN method to investigate the effect of the fractional order in the SFFHN model. Figure 4 displays the stable rotation solutions of the component $u$ at $t = 2000$ with $k = 1$ and $N = 256$ for four different fractional orders $\alpha = 2, 1.8, 1.7$ and $1.5$. As we observed in Figure 4, decreasing the fractional power has a notable effect on the width of the excitation wavefront and the wavelength of the system. A smaller value of $\alpha$ leads to a significantly reduced width of the wavefront and a shorter wavelength, allowing the domain to accommodate a larger number of wavefronts. These visualizations align with the findings presented by Bueno-Orovio et al. in ^[5]. In addition, these results highlight the application of fractional diffusion as a modeling tool to characterize intermediate dynamic states that cannot be solely described by pure diffusion mechanisms.

Example 2 (2D SFFHN model with Periodic Boundary): In this example, we consider the 2D SFFHN model with periodic boundary conditions, along with following initial conditions on domain $[-2.5, 2.5]^2$ to analyze the effectiveness and accuracy of the methods for periodic boundary conditions and varying $\alpha$ values.

In this experiment, we conducted a comparative analysis of the methods to assess their accuracy, computational cost, and effectiveness for solving SFFHN in the case of periodic boundary conditions and varying fractional order $\alpha$. We simulated the solution profile of the component $u$ up to the final time $t = 10$ by setting the fixed spectral resolution $N = 512$, which was chosen to be large enough to ensure that the errors primarily originated from the time integration process rather than spatial discretization.

In each simulation, we systematically reduced the temporal step size by halving it, starting with an initial value of $k = 0.1$. This step allowed us to investigate the impact of different temporal resolutions on the accuracy of the results and assess the convergence behavior of the numerical method. In Figures 5 and 6, we plotted the maximum error as a function of temporal step and the CPU time, respectively. From these figures, it is evident that all the methods achieved expected second-order convergence in the temporal direction for varying $\alpha$. Additionally, there is a noticeable variation in the CPU time among the methods, even when they achieve a comparable level of accuracy.

When considering both accuracy and computational cost in the selection of a second-order method for time stepping in Example 2, the ETD-CN method emerges as the most efficient choice, making it a favorable option when balancing accuracy and computational cost.

After ensuring that the proposed methods provide second-order accuracy for Example 2, in this sets of experiments we simulated the spatio-temporal evolution profile of the component $v$ and studied the effects of the super-diffusion on the SFFHN model. All the plots shown in Figure 7 are the aerial views of the numerical solutions $v$ and are obtained via the ETD-CN method with $N = 512, k = 1$ and $\alpha = 2, 1.8, 1.7, 1.5$ for various times. Since the evolution profile of the component $u$ is similar to the $v$, we will not discuss it in this paper. From Figure 7, we observe that as the fractional order decreases, both the width of the excitation wavefront and the wavelength of the system become significantly smaller. Consequently, the domain is capable of accommodating a greater number of wavefronts for lower orders. In the case of super-diffusion, which is characterized by the presence of a fractional Laplacian operator, the wave propagation exhibits long-tailed characteristics. This means that the influence of the wave extends over a larger spatial range compared to traditional diffusion. As a result, for wavefronts with similar widths, the distance between consecutive wavelengths becomes larger in the case of super-diffusion.These observations are in accordance with previous investigations conducted by researchers in ^[9].

Example 3 (2D SFFHN model with Dirichlet Boundary): In this example, we consider the 2D SFFHN model with homogeneous Dirichlet boundary conditions, along with the following initial conditions on domain $[-1.25, 1.25]^2$ to compare accuracy, efficiency and effectiveness of the methods for homogeneous Dirichlet boundary conditions and varying $\alpha$.

In order to analyze the empirical rates of convergence, computational efficiency, and effectiveness of the methods for the case of homogenous Dirichlet boundary conditions and super diffusion, we simulated the solution profile of the component $u$ with the methods up to time $t = 10.$ We fixed $N = 256$ and ran the simulations for various values of temporal step size starting with $k = 0.1$ and illustrated the maximum errors as a function of temporal step and CPU time in Figures 8 and 9, respectively. From these figures, it is evident that all the methods were able to capture second-order accuracy in the time direction for varying fractional orders $\alpha$. The convergence rates were consistent with the expected behavior. Moreover, there was noticeable variation in the CPU time among the methods, similar to what was observed in Examples 1 and 2, even when achieving a comparable level of accuracy.

Considering the case of Dirichlet boundary conditions with the specified initial conditions, again the ETD-CN method emerged as the most efficient choice, taking into account both accuracy and computational efficiency. The results depicted in Figures 8 and 9 further support the superiority of the ETD-CN method in this particular scenario.

Example 4 (3D SFFHN model with periodic boundary): In this example, we consider the 3D SFFHN model with periodic boundary conditions, along with the following initial conditions on domain $[-2.5, 2.5]^3$:

To further demonstrate the effectiveness and applicability of the methods, we applied them to investigate pattern formation in the 3D SFFHN model. In our analysis, we aimed to compare the accuracy and efficiency of the methods by conducting simulations on Example 4 until time $t = 5$. We kept the spectral resolution fixed at $N = 128$ and systematically reduced the temporal step size by halving it in each simulation, starting with an initial value of $k = 0.0625$.

The results regarding the order of convergence and computational time are depicted in Figures 10 and 11, respectively. From these figures, it is evident that all the methods achieved second-order convergence in the time direction for varying $\alpha$ values, as expected. Similar to the 2D case, the ETD-CN method stood out as the most efficient choice, considering both accuracy and computational efficiency. It demonstrated excellent performance in solving the 3D SFFHN model with periodic boundary conditions and varying fractional orders. These findings further emphasize the superiority of the ETD-CN method in capturing the dynamics and patterns of the 3D SFFHN model accurately while maintaining computational efficiency.

After confirming the accuracy and efficiency of the ETD-CN method for solving the 3D SFFHN model, we proceeded with a series of experiments to simulate the spatio-temporal evolution profiles of the component $u$ and $v$ for varying $\alpha = 2, 1.9, 1.8, 1.7$ to investigate the effects of super-diffusion on the SFFHN model. In these simulations, we used $N = 128, k = 1$ and ran the simulations until different times $t = 400,500, 1000,$ and $2000$. The complex patterns obtained from these simulations are illustrated in Figures 12 and 13. Notably, some of the patterns observed in these figures differ from those obtained in previous numerical work ^[9]. These unique patterns highlight the ability of the methods to capture novel and distinct dynamical behavior in the context of the SFFHN model. The results presented in Figures 12 and 13 provide valuable insights into the intricate spatio-temporal dynamics and pattern formation in the 3D SFFHN model under the influence of super-diffusion.

Example 5 (3D SFFHN model with Dirichlet boundary): In this example, we consider the 3D SFFHN model with homogeneous Dirichlet boundary conditions, along with the following initial conditions on domain $[-1.25, 1.25]^3$.

In this experiment we conducted a series of simulations using the ETD-Padé$(0, 2)$ method with a fixed spectral resolution of $N = 128$, temporal step size of $k = 1$ and varying fractional orders $\alpha = 2, 1.8, 1.7, 1.5$. The simulations were performed until time $t = 2000$ to investigate the effect of the fractional operator on the wave propagation process in the 3D SFFHN model. Figure 14 provides visual representation of the system's evolution, highlighting the generation of a spiral wave from the initial data. The process begins with the emergence of a traveling pulse, which subsequently starts rotating in a clockwise direction. Notably, as the fractional power decreases, the wave propagation extends toward the edge of the region. This phenomenon is accompanied by a noticeable reduction in both the width of the excitation wavefront and the wavelength of the system. The observed generation of spiral waves and the subsequent propagation dynamics shed light on the complex relationship between the fractional order and the spatio-temporal patterns in the system.

6.   Conclusions

In this manuscript, we have introduced four second-order time integrators along with the Fourier spectral method to investigate the influence of the fractional operator $\alpha$ in the 2D and 3D SFFHN model with periodic, Dirichlet, or Neumann boundary conditions. The effectiveness of these methods was compared based on two key parameters: Accuracy and computational efficiency. The experimental results have demonstrated that all the considered methods accurately captured the long-time spatio-temporal solution profiles of the system for varying $\alpha$, although with noticeable variations in CPU time. Considering both accuracy and computational cost when selecting a second-order method for time stepping the SFFHN model suggests that certain methods are more preferable than others, although all the methods yielded satisfactory results. Notably, the ETD-CN method emerged as the most efficient choice in this regard.

However, it is crucial to note that these conclusions are specific to the studies conducted on the SFFHN model with the specified initial and boundary conditions. The performance of the methods may vary depending on the specific case, such as the choice of initial conditions or other problem scenarios. It is important to emphasize that these results are not universally applicable and may differ for different sets of initial conditions, boundary conditions, and problem scenarios.

The methods introduced in this paper have the potential to be extended for solving fractional nonlinear reaction-diffusion models with variable coefficients or models with a fractional Laplacian of complex order ^[26] posed in higher dimensions, which will be the focus of our future research project.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The author is grateful to the scholarly activities committee of Utah Valley University for the summer research award. In addition, the author is grateful to the editor and anonymous referees for their valuable comments and suggestions, which helped us to improve the quality of this manuscript.

Conflict of interest

There is no conflict of interest to declare.