1.
Introduction
Reaction-diffusion phenomena is ubiquitous occurring mostly in different fields of science and engineering in which the system components interact [1,2,3,4,5]. They are used to describe population growth, propagation of travelling waves, pattern formation, other intriguing phenomena arising in combustion theory, tumor invasion, and neural networks spanning across various temporal and spatial scales [6,7,8,9]. Mathematically, the reaction-diffusion model is a parabolic partial differential equation (PDE) involving terms which denote diffusion and local reaction kinetics as given below:
where the first term on the left side represents the temporal part, the second term represents the diffusion component, and the term on the right side of equation denotes the local reaction kinetics. For the choice of $ Q(\mathrm{Y}) = \mathrm{Y} (1-\mathrm{Y}) $, the resultant model is Fisher's equation that can be used to describe population evolution and how the wave propagates in certain medium [10]. If $ Q(\mathrm{Y}) = \mathrm{Y} (1-\mathrm{Y}^{2}) $, the reaction-diffusion is called Newell-Whitehead-Segel (NWS) equation which is primarily used to describe convection phenomena in fluid thermodynamics [11]. Zeldovich-Frank-Kamenetskii equation (ZFK) is obtained for the choice of $ Q(\mathrm{Y}) = \mathrm{Y} (1-\mathrm{Y})e^{-\beta (1-\mathrm{Y})} $. In combustion theory, the ZFK equation is used to describe how flames propagate [12].
In theoretical neuroscience, the most well-studied reaction diffusion equation is the FitzHugh-Nagumo (FHN) model [13,14]. The FHN model is the simplification of the Hodgkin-Huxley model for the action potentials in squid giant axon [15,16,17]. Although FHN is not phenomenological (the involved parameters are not biophysical), it is a good candidate for studying how an action potential is generated and propagated. The main advantage of the FHN is that the solution space is two-dimensional and hence geometrical phase plane tools can be utilized to show how the trajectory evolves, which gives rise to the excitability and spike mechanism of spike generation. In the phase plane diagram of the FHN neuronal model, one of the nullclines has cubic nonlinearity. Over the last few decades, due to its model simplicity, FHN gained a lot of attention in the scientific community. Besides the generation and propagation of a train of spikes, FHN is also used to explained the dynamics of various physical systems in several fields of science, typically in the propagation of flame, Brownian motion process, autocatalytic chemical reaction, growth of logistic population, and neurophysiology [18,19].
Consider the following scaler reaction diffusion model with constant coefficients:
where $ \eta $ is a known real number and $ \mathrm{Y}(\xi, t) $ is an unknown function to be determined. In this equation, the local reaction kinetics is denoted by the cubic function $ \mathrm{Y}(1-\mathrm{Y})(\mathrm{Y}-\eta) $. For $ \eta = -1 $, Eq (1.2) reduces to the NWS equation.
In literature, various analytical and numerical strategies have been used to solve FHN-type equation. Shih et al. [20] used approximate conditional Lie point symmetry method for the numerical solution of the perturbed FHN-type equations. Mehta et al. [21] solved such nonlinear time dependent problems by using a novel block method coupled with compact finite difference schemes. In the perturbed model, the term $ \epsilon \mathrm{Y} $ has been added to the cubic term and for $ \epsilon = 0, $ the problem reduces to an unperturbed FHN-type model. Li and Guo implemented an integral approach for the exact solution of FHN-type equations [22]. Abbasbandy proposed a homotopy analysis method (HAM) to calculate a solitary-type solution of the nonlinear FHN-type model [23]. Kawahara and Tanaka derived the exact solution for the interaction of traveling fronts [24]. Gorder and Vajravelua obtained exact solutions for FHN-type equations and Nagumo telegraph equation using a variational method for fixed initial conditions in order to control the error [25]. Ali et al. proposed a Galerkin finite element method for the numerical solutions of FHN-type equation[26]. Mehran and Sadegh used a non-conventional finite difference scheme for the numerical solutions of FHN-type equations [27]. Gorder implemented HAM for the approximate solutions of FHN-type reaction diffusion equations [28].
Reaction-diffusion of FHN-type equations discusses previously involved constant coefficients. However, in general, time-dependent coefficients and dispersion-reaction terms involved in such models are more realistic physically [29,30,31]. Due to variation in problem geometry and taking into account the factor of heterogeneity in the propagating medium, FHN-type models with time-dependent coefficients are quite sophisticated to study.
The main goal of this work is to demonstrate the numerical solutions of the generalized FitzHugh Nagumo (GFHN) equation with a linear dispersion term and time-dependent coefficients, which are given below:
where $ \alpha(t) $, $ \vartheta(t) $, and $ \delta(t) $ are time-dependent coefficients. The associated initial condition is:
and the boundary conditions are:
where $ \psi_{0}(\xi), \psi_{a}(t), $ and $ \psi_{b}(t) $ are the known functions.
In recent literature, various authors computed the numerical solutions of GFHN equation with time-dependent coefficients. For instance, tanh and the Jacobi Gauss Lobatto collocation technique [18,29], the polynomial differential quadrature method [32], and the mixed-types discontinuous Galerkin method [33,34] were used.
Recently, wavelet numerical procedures attracted different researchers in various scientific communities because of their good properties and easy implementation. Some interesting aspects of wavelet methods include, compact support, localization in time and space, and multi-resistant analysis. Among different families of wavelet, Haar wavelet (HW) is simple and popular. The importance of the HW has been explored in various articles: Lepik proposed HW-based methods for different problems [35], Jiwari [36] utilized the HW method for the solution of Burger's equation, Kumar and Pandit [37] solved coupled Burger's equations using the HW method, Shiralashetti et al. [38] used the HW scheme for the numerical solutions of the initial valued problems. For further analysis of the HW, interested readers are referred to [39,40,41]. The remaining parts of the paper are presented in the following pattern:
● The underlying motivation and preliminaries are given in Sections 2 and 3, respectively.
● The proposed methodology is presented in Section 4.
● Numerical results and stability are reported in Sections 5 and 6.
● Finally, conclusions are drawn in Section 7.
2.
Main motive
The analytical solution of the time-dependent coefficients partial differential equations is quite complicated to compute. Therefore, numerical treatment is an alternative and an efficient way to cope with this issue. According to our analysis, the method of lines using the HW is not proposed for FHN-type models. In this work, we propose an HW-based method of lines for the solutions of FHN-type equations. The theoretical stability and its numerical verification will also be a part of this work.
3.
Haar wavelet and its integrals
Here, we address some basic results. Consider an arbitrary interval $ [a, b) $ and divide it into $ 2L $ equal subintervals of length $ \partial \xi = \frac{b-a}{2L} $, where $ L = 2^\mathcal{J} $ represents the maximum level of resolution. Define the dilation and translational parameters as $ j = 0, 1, ..., \mathcal{J} $, $ \mathcal{K} = 0, 1, ..., l-1 $, where $ l = 2^{j} $, respectively. Using the dilation and translational parameters, define the wavelet number $ \iota = l+\mathcal{K}+1 $. Now the first and $ \iota $-th HW are given below [35]:
where $ \zeta_{1}(\iota) = a+2\mathcal{K}\omega\partial \xi $, $ \zeta_{1}(\iota) = a+(2\mathcal{K}+1)\omega\delta \xi $, $ \zeta_{3}(\iota) = a+2(\mathcal{K}+1)\omega\partial \xi, $ and $ \omega = \frac{\mathrm{L}}{l} $.
In our analysis, we will use the following integrals:
where $ \alpha = 1, 2, 3, ..., n $ and $ \iota = 1, 2, 3, ..., 2L $. Analytically, the evaluation of these integrals for Eqs (3.1) and (3.2) are given below:
4.
Proposed scheme
In this section, the proposed method is described in detail. Here, an integral approach is utilized, therefore, the greatest order spatial derivative in Eq (1.3) can be approximated via the HW series as:
where $ \lambda_\iota (t) $ denotes the unknown HW coefficients and $ \mathrm{h}_\iota(\xi) $ are HW basis. From Eq (4.1), we deduce the following equations via twice integration:
The evaluation of Eq (4.3) at $ \xi = b $ gives:
Plugging $ \partial _{\xi}\mathrm{Y}(a, t) $ in Eqs (4.2) and (4.3), we get:
Now, the corresponding matrix form of Eqs (4.1)–(4.5), and (4.6) by using $ \xi \rightarrow \xi_{m} $ are:
where
From Eqs (4.7)–(4.9) and (1.3), the following equation can be obtained:
where "*" represents an element wise product. In a more compact form, Eq (4.10) can be written as:
where $ \mathbf{F}(t, \mathrm{Y}) = \vartheta(t)\big(\mathrm{h}\lambda(t) \big)- \alpha(t)\Big(\mathrm{M}_{1}\lambda(t)+\mathrm{M}_{2}(t)\Big)$-$\delta(t) \Big[(\mathrm{N}_{1}\lambda(t)+\mathrm{N}_{2}(t))*\Big(1-(\mathrm{N}_{1}\lambda(t) +\mathrm{N}_{2}(t))\Big)*\Big((\mathrm{N}_{1}\lambda(t)+\mathrm{N}_{2}(t))-\eta \Big)\Big]. $
Here, Eq (4.11) represents the system of first-order ordinary differential equations, which can be solved via the RK-4 scheme discussed in a later section.
4.1. Temporal discretization
To obtain the solutions of Eq (4.11), we use the following $ \mbox{RK-4} $ scheme:
where $ \partial t $ is the time step size.
5.
Simulations and comparison
In this section, numerical solutions of GFHN with constant and time-dependent coefficient are listed. The efficiency of the proposed scheme is checked by applying different error measures, namely: $ L_2 $, $ L_\infty $, and $ L_{rms} $ as described below:
5.1. Problem 1
Here, we consider Eq (1.2) with constant coefficients and various choices of the exact solutions:
The associated initial and boundary conditions for all cases are used from the exact solutions. Simulations are done in different spatial domains to test the technique, and a comparison is made with the previously published work. The spatial domains used for case $ (i) $ are $ [0, 1] $ and $ [-10, 10] $, for case $ (ii) $ is $ [-22, 22] $, for case $ (iii) $ are $ [0, 1] $ and $ [-10, 10] $, while for case $ (iv) $ is the interval $ [-10, 10] $. The obtained results of case $ (i) $ in $ [-10, 10] $ with different values of $ \eta $ are given in Tables 1–3, while the outcomes in the interval $ [0, 1] $ with various values of $ \eta $ are given in Table 4.
In these tables, the obtained solutions are compared with the existing results in the literature [32,42,43,44]. Through comparison, it is obvious that the present technique shows better performance than the cited work [32,42,44]. The present method is based on RK-4 which gives good results using a small step size.
In Table 5, the computed values of the constants $ \mathcal{K}_{1}, \mathcal{K}_{2}, \mathcal{K}_{3}, \mbox{and}\; \mathcal{K}_{4}, $ for $ \eta $ = 0.75, $ \partial $t = 0.001, and $ L $ = 2 in the spatial domain $ [-10, 10] $ at distinct time levels are reported.
In Figure 1, numerical and exact solutions are plotted for various values of $ \eta $ for case $ (i) $ while in Figure 2, the absolute error at distinct points are plotted with various time levels. Similar simulations are obtained for case ($ ii $) and its absolute error at distinct points are compared with the existing results [33,45,46] in Table 6, which show the superiority of the method.
In Figures 3 and 4, the solutions profiles are shown for different values of $ \eta $ in the form two- and three- dimensions plots, respectively, for the same case. Both figures reveal the mutual agreement of exact and numerical solutions.
The obtained results of case $ (iii) $ and those available in [44] are given in Tables 7 and 8, where the small error norm shows the high quality of the scheme. The computed norms of case ($ iv $) are listed in Table 9.
From this table, an obvious relation between the resolution level and accuracy is observed, which predicts when resolution increases the accuracy also increases. Similarly, one can see that when time increases, the accuracy reduces due to round-off errors.
Graphically, the solutions for case $ (iii) $ are displayed in Figures 5 and 6, respectively. Likewise, solutions for case $ (iv) $ are presented in Figures 7 and 8. In both cases, the numerical and exact solutions show good agreement. From all tabulated and graphical solutions, we conclude that the presented scheme is quite suitable for solving constant coefficient FHN models.
5.2. Problem 2
Here, we considered Eq (1.3) with $ \alpha(t) = \cos(t) = \vartheta(t) $ and $ \delta(t) = 2\cos(t). $ The associated exact solution is given by:
The corresponding initial and boundary conditions are used from the given solution. This problem is solved in the spatial domain $ [-10, 10] $ for comparison purposes. The numerical experiments are conducted with different values of $ \eta $ at time. The extracted findings are reported and matched with the existing results [32,33] in Table 10. The table demonstrates noticeable accuracy versus the cited work.
Moreover, the exact and numerical solutions at different time levels are presented graphically in Figure 9 for mutual comparison. The solution comparison is also illustrated in the form of surfaces in Figure 10 for different values of $ \eta $. In Figure 11, solutions are shown for different values of $ \eta $ which discloses the travelling front as $ \eta $ increases. The corresponding error at different time levels are pictured in Figure 12. Again, it is evident that the scheme produces pretty good results for time-dependent variable coefficient problems.
6.
Stability analysis
Here, we discuss the stability of the present scheme computationally. Consider the aforementioned model as:
with the associated initial and boundary conditions. Discretization of Eq (6.1) in space by the truncated Haar wavelet series gives rise to the ordinary differential equations system in a time grid as:
where $ [\mathrm{Y}] $ describes the vector of unknowns, $ \Theta $ is the coefficient matrix, and $ \Xi $ is a vector comprised of a nonhomogeneous part and boundary conditions. The stability of Eq (6.2) based on the coefficient matrix $ \Theta $ which is as defined as follows:
(a) For a constant coefficient: $ \Theta = H + (diag(\mathrm{Y}(t)))*(1-\beta_{2})*(\beta_2-\eta) $,
(b) For a time-dependent coefficient: $ \Theta = \vartheta(t)H-\alpha(t)\beta_1-\delta(t)(diag(\mathrm{Y}(t)))*(1-\beta_{2})*(\beta_2-\eta) $, where $ * $ denotes the element-wise product. To discuss the stability of the proposed method, we debate on the eigenvalues of $ \Theta $. Let $ \lambda_i $ be the eigenvalues of $ \Theta $ and $ \partial t $ the time step size. As $ t \rightarrow \infty $, the stable solution $ \mathrm{Y} $ needs to satisfy the following conditions:
● For real eigenvalues: $ -2.78 < \partial t \lambda_\iota < 0, $
● For imaginary eigenvalues: $ -2\sqrt{2} < \partial t \lambda_\iota < 2\sqrt{2}, $
● For complex eigenvalues: $ \partial t \lambda_{\iota} $ lies in the stability region addressed in [47,48].
As stated in [48], if the eigenvalues are complex, the $ Re(\partial t \lambda_\iota) $ may be a small positive number. For different problems the eigenvalues are calculated and plotted in Figures 13–15, which show that the eigenvalues lie in the stable region. So, the method is stable and condition 3 is fulfilled. In Figures, the values of $ Re(\partial t \lambda_\iota) $ are small, which is shown by the value of $ 10^{-3} $.
7.
Concluding remarks
In this work, the Haar wavelet method of lines is implemented for the numerical solutions of the FHN reaction diffusion model with constant and time-dependent coefficients. The collocation procedure has been adopted in Haar wavelet basis for the estimation of the derivatives and solution. In this way, the nonlinear FHN models has been transformed to the initial value problems. Thereafter, the initial value problems have been solved with RK-4 scheme. The resultant outcomes have been matched with some existing literature work. Moreover, the stability of the scheme has been verified computationally. It has been noticed that the proposed numerical strategy is a good tool to estimate the solutions of the FHN models with constant and time-dependent coefficients.
Author contributions
Aslam Khan: Software, writing original draft preparation; Abdul Ghafoor: Supervision, conceptualization, methodology; Emek Khan: reviewing original draft; Kamal Shah: formal analysis, Funding acquisition; Thabet Abdeljawad: validation, project administration.
Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
The authors would like to thank Prince Sultan University for paying the APC and for their support through the TAS research lab.
Conflict of interest
The authors declare there is no conflict of interest.
DownLoad:
