Efficient Numerical Schemes for a Heterogeneous Reaction–Diffusion System with Applications

Abstract

In this study, a class of nonlinear heterogeneous reaction–diffusion system (RDS) has been considered that arises in modeling epidemiological interactions, environmental sciences, and chemical and ecological systems. Numerical and analytic solutions for this kind of variable medium nonlinear RDS are challenging. This article developed a few highly accurate numerical schemes for such problems. For the spatial integration of the heterogeneous RDS, a few finite difference schemes, a Bernstein collocation scheme, and a Fourier transform scheme were employed. The stability and accuracy analysis of the spectral schemes were studied to confirm the order of convergence of the approximation. A few methods of lines were then used for the temporal integration of the resulting semidiscrete model. It was confirmed theoretically that the spectral/pseudo-spectral method is very efficient and highly accurate for such a model. A fast and efficient solver for the resulting full discrete system is highly desired. A Newton–GMRES–Multigrid solver was applied for the resulting full discrete system. It is demonstrated in tabular form that a multigrid accelerated Newton–GMRES solver outperforms most linear solvers for such a model.

1. Introduction

The spatial characteristics of the environment can be represented by various PDEs of different kinds, such as ecology, composed of chemical, biological, geological, and physical, with various types and levels of heterogeneity. Heterogeneous PDEs, in which reaction–diffusion has been used, are employed in biological and chemical systems, in population dynamics, and in epidemics. Heterogeneity is associated with uneven changes in a system with general rules and also variations in external conditions. These systems may exist in the initial values of the problem modeled (dependent variables), in space, or by time varying the values of properties (see, e.g., [1]). Furthermore, heterogeneity in reaction–diffusion systems increases the difficulty of obtaining the exact solutions. In 1952, Alan Turing first developed a seminal mechanism for self-organizing emergence, such as the development of different patterns in various biological systems. He considered a reaction–diffusion system in which two concentrations diffuse and interact with each other with reaction terms and determined that a small variation in their concentrations near the dynamic equilibrium can be heterogeneously amplified for proper reaction terms. The resulting behavior of the Turing model with fixed constant coefficients is generally highly periodic and usually fails to generate real scenarios with many real-life applications [2,3].

The reaction–diffusion system (RDS) model is a time-dependent partial differential equation that describes physical realities. A detailed mathematical initiative, especially numerical schemes, may have numerous applications in understanding and classifying structures and patterns in engineering and biological systems [3,4]. The numerical range of parameter values and the physical meaning of the boundary conditions have been explained in detail in [3,5] and some references therein. A class of highly accurate numerical schemes is available for mathematical models that contain diffusion operators in a homogeneous medium [6,7]. There are several works that explain the applications of RDS in epidemics [8,9]. In [8,9], the authors analyzed a few RDS models used for describing epidemic processes. The main attention in the studies was on an RDS that explains COVID-19 pandemic modeling. Here, they studied traveling wave solutions, exact solutions, and parameter estimations. The authors considered a class of RDSs that were capable of modeling the development and continuation of the COVID-19 pandemic in both space and time in [8]. The study confirmed that some phenomena cannot be captured using ODEs, and thus a reaction–diffusion system with a network structure within an RDS may help to better understand the progression of an epidemic. In [10], the authors considered stochastic reaction–diffusion processes using graph neural networks and used cheap Monte Carlo simulations of reaction–diffusion processes. They also analyzed the stability and accuracy analysis of the schemes. Thus, a rigorous numerical study of RDS is highly needed; hence, the ideas developed by these authors can be implemented for the problem considered in this study, which will lead to a good future research direction.

The dynamics of the progression of epidemics and chemical and ecological systems may depend on heterogeneous ordering [3,5,8,9]. Heterogeneities are an underlying feature of environmental and biological reaction–diffusion systems. There is a growing interest in mathematical models where heterogeneity plays a significant role, including heterogeneous reaction–diffusion systems, where the diffusion parameters depend on space, and the reaction terms depend explicitly on space and/or time. These types of systems can also be found in the Belousov–Zhabotinsky cardiac cell model and neuron cell dynamics.

Recently, it has become clear that spatial heterogeneity has a significant impact on ecological dynamics [3]. This work focuses on the study of the numerical approximations of two-component RDSs [4,8,11]:

$$\begin{array}{c}{\displaystyle \frac{\partial u}{\partial t}}={\delta }_1(x,y)▵u+f(u,v),\\ {\displaystyle \frac{\partial v}{\partial t}}={\delta }_2(x,y)▵v+g(u,v),\end{array}$$

(1)

where $▵={\displaystyle \frac{{\partial }^2}{\partial x^2}}+{\displaystyle \frac{{\partial }^2}{\partial y^2}}$, f and g are nonlinearities, and $(x,y)\in \Omega \subset {\mathbb{R}}^2$. One of the most used such nonlinearities is as follows:

$$f(u,v)=ru(1-{\displaystyle \frac{u}{w}})-pvh\left(ku\right),g(u,v)=qvh\left(ku\right)-sv,$$

where $u(x,t)$ and $v(x,t)$ are the concentrations of the populations of prey and predators at a particular t and spatial position $x\in \Omega$. Here, the model parameters $r,w,p,$${\delta }_i$, $i=1,2$, $k,q,s$ are in $(0,\infty )$. ${\delta }_i$, $i=1,2$ may depend on spatial variables. The first term determines how quickly the quantity u increases. The first term, ${\delta }_1▵u$, is the diffusion term. Here, the functional response $h\in C^2[a,b]$ satisfies [4]

• $h\left(x\right)=0$ at $x=0$.
• $h\left(x\right)$ approaches 1 as x approaches ∞.
• $h\left(x\right)>0$.

One such choice of h is $h\left(x\right)={\displaystyle \frac{x}{\alpha +x}}$ (forms kinetics 1), which is an example of the resulting nonlinearity $(f,g)=\left(u(1-u-{\displaystyle \frac{v}{\alpha +u}}),u({\displaystyle \frac{\beta v}{\alpha +u}}-\gamma )\right)$, whereas another option may be the function $h\left(x\right)=(1-e^{-x})$ (forms kinetics 2) [4].

In addition, the following reaction kinetics can be found in the Gray–Scott model:

$$\begin{array}{c}f(u,v)=F(1-u)-u{v}^2,\\ g(u,v)=u{v}^2-(F+k)v,\end{array}$$

(2)

where u and v are the concentrations of the two chemical components, $F>0$ denotes the dimensionless flow rate of the reaction, and $k>0$ indicates the rate constant of the activator. Thus, a few different types of kinetics are considered in this study. Stability issues/time–long dynamics of the RDS with such two nonlinearities have been well discussed in [3,4].

An extensive field of study, this kind of heterogeneous RDS replicates diverse predator–prey interactions and represents a broad spectrum of biologically, chemically, and environmentally relevant behavior arising from intrinsic variables alone. A reaction diffusion system was considered in [4], where the authors detailed a finite difference securitization and illustrated numerical solutions. The authors also studied the parameter-sensitive dynamical behavior of the solutions of the RDS. The study  [11] contained the same problem as in [4] and focused on developing an efficient linear system solver for the RDS. In this study, the author used a similar finite difference scheme to [4] and employed a few linear system solvers for the resulting system of equations. It is noticed from this study that the multigrid accelerated GMRES solver was a time–efficient technique for the homogeneous RDS considered.

Detailed mathematical and physical studies about the space–varying diffusion parameter in an RDS can be found in [3], where the author rigorously studied the effect of the heterogeneous diffusion coefficients in pattern recognition in biological systems and chemical systems. Here it was confirmed that the use of space-varying diffusion coefficients can be used to modify classical Turing patterns. The author used the Galerkin method for the space discretization. The analytic solution properties were discussed in detail. The author studied the generic instability of the solutions and the physical mechanism of real-life problems that lead to heterogeneity. The heterogeneous nonlinear problem of type (1) also arises in the progression of epidemics as well [8,10].

Studying numerical methods for heterogeneous RDS is important, as well as very challenging, due to the variability of the parameter values. A few studies have addressed the numerical solutions of this problem [3]. To the best of our knowledge, no existing solver is very fast and efficient for a heterogeneous RDS of type (1). The implementation of a spectral method for such a system of nonlinear PDEs is also challenging; a multigrid preconditioned linear solver for such a model has not been studied yet. Thus a smart algorithm will undoubtedly assist applied scientists to employ these kinds of models for research. Motivated by the studies described above, in this study, the Bernstein polynomial scheme, Fourier transform scheme, and a few finite difference schemes were employed for the spatial integration of (1) followed by a few methods of lines for time integrals. The stability and accuracy analysis of the spectral schemes were well studied to confirm the order of convergence of the approximation. Theoretically, it has been established that the spectral and pseudo-spectral methods are highly accurate and efficient for this type of model. The schemes have been illustrated through several graphical representations of solutions. Now, a combination of these schemes contains a full discrete system of nonlinear equations, which is very challenging to investigate. A preconditioned solver was employed to solve the abovementioned system. To be specific, a Newton–GMRES–Multigrid solver was applied for the resulting full non-linear system of equations. Here, our illustrations confirm that a multigrid accelerated Newton–GMRES solver outperforms most linear solvers for the full discrete heterogeneous RDS.

The rest of the article has been organized in the following manner. Section 2 deals with the Bernstein polynomial and Fourier transform schemes for the spatial domain. We also analyze the stability and accuracy of such approximations in this section. Section 3 contains a few full discrete finite difference schemes with the multigrid accelerated Newton–GMRES method for the resulting nonlinear system. A detailed numerical illustration for the schemes discussed in Section 2 and Section 3 are presented and analyzed in Section 4. A few concluding remarks are included in Section 5.

2. Spectral Schemes for the RDS

In contrast to finite difference methods, spectral methods are global, which means that the calculation at any given location is based on information from the entire domain, not just the surrounding points. Spectral approaches converge dramatically, making them more accurate than local methods. When the integral/solution changes significantly in time or space, when extremely high spatial resolution is required, or when long-term integration is required, global methods are preferred over local methods.

2.1. Bernstein Collocation Scheme

In the mathematical field of numerical analysis, a Bernstein polynomial is typically defined by a superposition of Bernstein basis functions or polynomials (BBPs) [12]. Here, our goal is to use the Bernstein collocation method for spatial integration to approximate the model. The accuracy problems with such a polynomial approximation have been well examined in [13,14,15]. Before we move onto the main discussion, a brief explanation of Bernstein polynomials over $[0,1]$ domain and $[a,b]$ domain is included here.

• Bernstein Basis Polynomials (BBPs) in $[0,1]$

BBPs of degree n over $[0,1]$ can be defined as [12,16]

$$b_{n,k}\left(x\right):=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)x^k{(1-x)}^{n-k},k=0,1,2,\cdots ,n,$$

where ${}^{n}C_r$ are the binomial coefficients. The respective Bernstein polynomials of degree n can be defined by [12,17]

$$B_n\left(x\right)=\sum _{k=0}^n{c}_k{b}_{n,k}\left(x\right),$$

where $c_k$ are known as Bernstein coefficients. These polynomials satisfy [12,13,17]

• $b_{n,k}=0,$ if $k<0$ or $k>n$,
• $b_{n,k}\ge 0,$ for $x\in [0,1]$,
• $b_{n,k}(1-x)=b_{n,n-k}\left(x\right)$,

and the derivatives of such a polynomial satisfies

$$b_{n,k}^{\prime }\left(x\right)=n\left(b_{n-1,k-1}\left(x\right)-b_{n-1,k}\left(x\right)\right),$$

$$b_{n,k}^{\left(l\right)}\left(0\right)={\displaystyle \frac{n!}{(n-l)!}}\left({\displaystyle \genfrac{}{}{0pt}{}{l}{k}}\right){(-1)}^{k+l},$$

$$b_{n,k}^{\left(l\right)}\left(1\right)={(-1)}^l{b}_{n,n-k}^{\left(l\right)}\left(0\right).$$

Let us consider $F\in C[0,1]$. Then, one of the best uses of such a polynomial is to approximate f efficiently. The efficiency is guaranteed by the famous Weierstrass approximation theorem, providing the following results [12,14,18]:

$$B_n\left(F\right)\left(x\right)=\sum _{k=0}^{n}F\left({\displaystyle \frac{k}{n}}\right)b_{n,k}\left(x\right),$$

$${∥B_n{\left(F\right)}^{\left(l\right)}∥}_{\infty }\le {\displaystyle \frac{{\left(n\right)}_l}{n^l}}{∥F^{\left(l\right)}∥}_{\infty },$$

and

$${∥F^{\left(l\right)}-B_n{\left(F\right)}^{\left(l\right)}∥}_{\infty }\to 0.$$

• Bernstein Basis Polynomials (BBPs) in $[a,b]$

Now, we aim to define such polynomials over an arbitrary domain $[a,b]$. Such polynomials can be defined over $[a,b]$ as [12]

$$b_{n,j}\left(x\right){=}^n{C}_r{\displaystyle \frac{{(x-a)}^j{(b-x)}^{n-j}}{{(b-a)}^n}},j\in \mathbb{N},\mathrm{and}0\le j\le n,$$

and it satisfies the properties stated above for the same polynomials over $[0,1]$. Now, letting $h={\displaystyle \frac{b-a}{N}}$, and $x_i=x_0+ih$, $i=0,1,2,3,\cdots ,n$, a general form of Bernstein polynomials of degree n over $[a,b]$ can be written as

$${\beta }_{i,n}\left(x\right)=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right){\displaystyle \frac{{(x-x_0)}^i{(x_n-x)}^{n-i}}{{(x_n-x_0)}^n}},$$

where $i=0,1,\cdots ,n$. These polynomials satisfy the properties stated above. These polynomials also satisfy the following recurrence relations [12]:

$${\beta }_{i,n}\left(x\right)={\displaystyle \frac{x_n-x}{x_n-x_0}}{\beta }_{i,n-1}\left(x\right)+{\displaystyle \frac{x-x_0}{x_n-x_0}}{\beta }_{i-1,n-1}\left(x\right),$$

$${\beta }_{i,n}^{\prime }\left(x\right)={\displaystyle \frac{n}{x_n-x_0}}\left({\beta }_{i-1,n-1}\left(x\right)-{\beta }_{i,n-1}\left(x\right)\right),$$

$${\beta }_{i,n}^{″}\left(x\right)={\displaystyle \frac{n(n-1)}{{(x_n-x_0)}^2}}\left({\beta }_{i-2,n-2}\left(x\right)-2{\beta }_{i-1,n-2}\left(x\right)+{\beta }_{i,n-2}\left(x\right)\right).$$

Thus, it is evident that

$$\beta \left(m\right)\left(x\right)=D^m{\beta }_{n-m}\left(x\right),\mathrm{where}\beta \left(m\right)={[{\beta }_{0,n}^m{b}_{1,n}^{*m}\cdots b_{n,n}^{*m}]}^T,$$

and $D^m$ is a derivative operator of order m [12,19]:

$$D^m{b}_{N,j}\left(x\right)={\displaystyle \frac{1}{{(b-a)}^m}}\left({\displaystyle \genfrac{}{}{0pt}{}{N}{m}}\right)\sum _{l=max(0,j+m-N)}^{min(j,m)}\left({\displaystyle \genfrac{}{}{0pt}{}{p}{l}}\right)b_{j-l,N-m}\left(x\right).$$

Here, we recall the one dimensional version of the RDS (1) by

$$\begin{array}{c}{\displaystyle \frac{\partial u}{\partial t}}={\delta }_1\left(x\right)▵u+f(u,v),\\ {\displaystyle \frac{\partial v}{\partial t}}={\delta }_2\left(x\right)▵v+g(u,v).\end{array}$$

(3)

Let the approximate solutions for (3) be

$$u_h(x,t)=\sum _{j=0}^M{C}_j^u\left(t\right){\beta }_{j,M}\left(x\right),\mathrm{and}{v}_h(x,t)=\sum _{j=0}^M{C}_j^v\left(t\right){\beta }_{j,M}\left(x\right),$$

where ${\beta }_{j,M}\left(x\right)$ are the Bernstein basis polynomials defined above and at $x_k$, $k=0\left(1\right)M$,

$$u_h(x_k,t)=\sum _{j=0}^M{C}_j^u\left(t\right){\beta }_{j,M}\left(x_k\right),v_h(x_k,t)=\sum _{j=0}^M{C}_j^v\left(t\right){\beta }_{j,M}\left(x_k\right),$$

and it yields

$${\underline{u}}_h\left(t\right)=\mathcal{B}{C}^u\left(t\right),\mathrm{where}\mathcal{B}=\left({\beta }_{j,M}\left(x_k\right)\right),\mathrm{and}{C}^u=\left(C_j^u\right),$$

and

$${\underline{v}}_h\left(t\right)=\mathcal{B}{C}^v\left(t\right),\mathrm{where}\mathcal{B}=\left({\beta }_{j,M}\left(x_k\right)\right),\mathrm{and}{C}^v=\left(C_j^v\right).$$

Then, system (3) can be approximated by the following system of ordinary differential equations through the use of product approximation of the nonlinear part  [20]

$$\begin{array}{cc}\hfill \mathcal{B}{\displaystyle \frac{d}{dt}}{C}^u\left(t\right)& =M_{\alpha }{D}^2{C}^u\left(t\right)+\mathcal{B}f(C^u,C^v),\hfill \\ \hfill \mathcal{B}{\displaystyle \frac{d}{dt}}{C}^v\left(t\right)& =M_{\beta }{D}^2{C}^v\left(t\right)+\mathcal{B}g(C^u,C^v),\hfill \end{array}$$

(4)

where $M_{\alpha }=diag\left({\delta }_1\left(x_j\right)\right)$, and $M_{\beta }=diag\left({\delta }_2\left(x_j\right)\right)$. System (4) can be written in a matrix-vector form as

$$\begin{array}{c}\hfill {\mathcal{B}}^{*}{\displaystyle \frac{d\mathcal{C}\left(t\right)}{dt}}=\mathcal{M}\mathcal{C}\left(t\right)+\mathcal{B}\mathcal{F},\end{array}$$

(5)

where ${\mathcal{B}}^{*}=\left(\begin{array}{cc}\mathcal{B}& 0\\ 0& \mathcal{B}\end{array}\right),$ $\mathcal{M}=\left(\begin{array}{cc}{M}_{\alpha }{D}^2& 0\\ 0& M_{\beta }{D}^2\end{array}\right),$ $\mathcal{C}\left(t\right)={\left(C^u\left(t\right)C^v\left(t\right)\right)}^{\prime }$, and $\mathcal{F}={\left(f(C^u,C^v)g(C^u,C^v)\right)}^{\prime }$.

• Two-Space Dimensional Implementation of the Scheme

Now, a two-space dimensional version of Bernstein interpolating polynomials can be defined through tensor product approximation by

$$u_h(x,y)=\sum _{j=0}^N\sum _{l=0}^M{\gamma }_{j,l}{b}_{N,j}\left(x\right)b_{M,l}\left(y\right),$$

(6)

where the coefficients ${\gamma }_{i,j}$ can be computed at a sufficient number of collocation points. This idea consists of substituting the partial differential operators by the derivatives of the interpolating polynomials (6). Through vectorization, the matrix $D^r$ can be implemented similar to the Tchebychev pseudospectral method discussed in [14,21]. Now, the 2D version of the discrete diffusion operator on a square domain can be defined by

$$L_h^1=I\otimes \left(M_{\alpha }{D}^2\right)+\left(M_{\alpha }{D}^2\right)\otimes I,L_h^2=I\otimes \left(M_{\beta }{D}^2\right)+\left(M_{\beta }{D}^2\right)\otimes I,$$

where I is an identity operator of same order as of D.

Thus, the resulting semi-discrete RDS can be written as

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{u}_h}{dt}}& =L_h^1{u}_h+f_h(u_h,v_h),\hfill \\ \hfill {\displaystyle \frac{d{v}_h}{dt}}& =L_h^2{v}_h+g_h(u_h,v_h),\hfill \end{array}$$

with a set of suitable initial conditions, and thus, the resulting coupled IVP can be solved using any standard implicit/explicit solvers. For more implementation details, we refer the reader to [14] (Chap. 13).

2.2. Stability and Convergence of Bernstein Collocation Scheme

The analysis of the error resulting from the previously outlined numerical method is put forward in this section. Here, we follow [22] for the theoretical analysis in this article. Consequently, the vector form of the system (3) can be expressed as follows:

$${\displaystyle \frac{\partial \mathbf{u}}{\partial t}}=\delta \left(x\right){\nabla }^2\mathbf{u}+\mathbf{f}\left(\mathbf{u}\right),x\in \Omega \subset {\mathbb{R}}^2,$$

(7)

where ${\nabla }^2={\displaystyle \frac{{\partial }^2}{\partial x^2}}+{\displaystyle \frac{{\partial }^2}{\partial y^2}}$. The vectors u and f are associated with a relevant coefficient matrix D. The following approximation features of this problem have a standard demonstration for Dirichlet boundary conditions.

Theorem 1.

Assume that u is the exact solution to the problem (7), given the necessary initial and boundary conditions. Let f be an n continuously differentiable function, and assume once more that ${\mathbf{u}}^N\in X_N$ is the approximate solution of the Bernstein collocation method. The relationship that follows therefore holds:

$${\parallel {\mathbf{u}}^N\parallel }^2+{\int }_0^t\delta {\nabla }^2{\parallel {\mathbf{u}}^N\parallel }_E^{2}ds\le {\parallel {\mathbf{u}}_0^N\parallel }^2+\mathcal{C}{\int }_0^t\parallel f^n\parallel ds,$$

with $\mathcal{C}$ functioning as a constant.

Proof. 

Let $\mathbf{u}=(u_1,u_2)$, $\mathbf{\delta }=({\mathbf{\delta }}_1\left(x\right),{\mathbf{\delta }}_2\left(x\right))$, and $u_0$ be the given initial conditions. Determine an expression of $\mathbf{k}=(k_1,k_2)$ using any pair of positive or negative integers. Hence, the Euclidean inner product ${\mathbb{R}}^2$ is represented by the expression $\mathbf{k}.\mathbf{x}=k_1{x}_1+k_2{x}_2$. Additionally, we denote the set of multiindices $\mathbf{k}=(k_1,k_2)$ by J. Assume that X is a Hilbert space and that, for a positive integer N, $X_N$ and $Y_N$ are finite-dimensional subspaces of X with the same dimension, and

$$X_N=span\left\{b_{k,N}\right|\mathbf{u}(\partial \Omega )=0,\mathbf{k}\in J\}.$$

Let $P_N$ be the orthogonal projection of $X_N$ onto the $\omega$-weighted $L^2$. The function that represents the Bernstein collocation approximation is

$${\mathbf{u}}^N(\mathbf{x},t)=\sum _{k\in J}{u}^N({\mathbf{x}}_k,t)b_{k,n},$$

where $b_{k,n}$ denotes the characteristic of Bernstein polynomials. These polynomials satisfy the following equation at the collocation points

$$\left[{\mathbf{u}}_t^N-{\mathcal{L}}_N{\mathbf{u}}^N-\mathbf{f}\right]({\mathbf{x}}_{\mathbf{k}},t)=0,forallt>0,\mathrm{and}\mathbf{k}\in J.$$

Here, ${\mathcal{L}}_N=\delta {\nabla }^2\mathbf{u}$, for all $\mathbf{u}\in X_N$; this refers to the interpolation approximation of $\mathcal{L}u$. The coercivity and continuity conditions are satisfied with $E=H_{\omega ,0}^1(\Omega )$, which is a Hilbert space in X, and the norm is [22]

$${\parallel u\parallel }_E={\left({\int }_{\Omega }{\left|{\mathbf{u}}_x\right|}^{2}d\Omega \right)}^{{\displaystyle \frac{1}{2}}}.$$

The operator ${\mathcal{L}}_N=\delta {\nabla }^2\mathbf{u}=\mathbf{\delta }({\displaystyle \frac{{\partial }^2}{\partial x^2}}+{\displaystyle \frac{{\partial }^2}{\partial y^2}})\mathbf{u}$ in this scheme satisfies the subsequent energy identity:

$$(\mathcal{L}\mathbf{u},\mathbf{v})=\alpha (\mathbf{u},\mathbf{v})={\int }_{\Omega }{\left|{\mathbf{u}}_x\right|}^{2}d\Omega .$$

Simply put, the square root on the right side of the equation is the semi–norm of the Hilbert space. Once more, we suppose that, for any $N>0$, $X_N$ is contained in E. Additionally, it is assumed that, for any $\mathbf{u}\in X_N$, $\parallel \mathbf{u}\parallel \le \mathcal{C}{\parallel \mathbf{u}\parallel }_E$, where $\mathcal{C}>0$ is not dependent on N [22]. The following is a coercivity requirement for the approximation (7). If there is a constant $\alpha >0$ (that is, independent of N), such that

$$\alpha {\parallel \mathbf{u}\parallel }_E^2\le {({\mathcal{L}}_N\mathbf{u},\mathbf{u})}_N\forall u\in X_N,$$

then it is clear that the approximation is stable. This means that the estimate

$${\parallel {\mathbf{u}}^N\parallel }_E\le {\displaystyle \frac{\mathcal{C}}{\alpha }}{\parallel f^n\parallel }_N$$

holds true. In fact, it has

$$\alpha {\parallel {\mathbf{u}}^N\parallel }_E^2\le {({\mathcal{Q}}_n{\mathcal{L}}_N{\mathbf{u}}^N,{\mathbf{u}}^N)}_N={({\mathcal{Q}}_{N}f,{\mathbf{u}}^N)}_N\le {\parallel {\mathcal{Q}}_{N}f\parallel }_N{\parallel {\mathbf{u}}^N\parallel }_N\le \mathcal{C}{\parallel f^n\parallel }_N.$$

Here, we exploit the fact that ${\mathcal{Q}}_N$ is the orthogonal projection operator upon $Y_N$ with respect to the discrete inner product ${(\mathbf{u},\mathbf{v})}_N$. From this point forward, we assume that the continuity requirement and the coercivity conditions, which are the hypotheses presented in [22], satisfy the problem. Initially, we examine an inner product using the provided equation as

$$({\displaystyle \frac{\partial }{\partial t}}{\mathbf{u}}^N,\mathbf{v})=(\delta {\nabla }^2{\mathbf{u}}^N,\mathbf{v})+(f\left({\mathbf{u}}^N\right),\mathbf{v}).$$

Consider $\mathbf{v}={\mathbf{u}}^N\left(t\right)$. For each $t>0$, we obtain [20]

$${\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}{\parallel {\mathbf{u}}^N\parallel }^2=\delta {\nabla }^2{\parallel {\mathbf{u}}^N\parallel }_E^2+\parallel f^n\parallel {\parallel {\mathbf{u}}^N\parallel }^{2}ds.$$

After integration, we have

$${\parallel {\mathbf{u}}^N\parallel }^2+{\int }_0^t\delta {\nabla }^2{\parallel {\mathbf{u}}^N\parallel }_E^{2}ds\le {\parallel {\mathbf{u}}_0^N\parallel }^2+\mathcal{C}{\int }_0^t\parallel f^n\parallel ds.$$

This establishes the Bernstein approximation’s stability.    □

For the investigation of the spectral approximations of boundary-value problems, the approximation results in the Sobolev norms are significant. In this instance, projecting onto the space of polynomials satisfying the boundary data rather than merely the space of polynomials would be more appropriate. Consider the Dirichlet conditions at both ends of the interval $[-1,1]$. The functions of $H^1(-1,1)$ are the Sobolev space that meets these requirements, which constitutes a subspace, which is typically represented by $H_0^1(-1,1)$.

Theorem 2.

Let u represent the exact solution to (7) with the Dirichlet boundary condition. Suppose that f is n times continuously differentiable function, and assume once more that ${\mathbf{u}}^N\in X_N$ is the approximate solution of the Equation (7). After that, the inequality presented below is valid:

$${\parallel \mathbf{u}-{\mathbf{u}}^N\parallel }_E\le {\parallel e\left(0\right)\parallel }_N^2+\mathcal{C}{\left(N^{1-n}{\int }_0^t{\left|{\mathbf{u}}_t\left(s\right)\right|}_{H^{n-2}}^{2}ds+N^{1-n}{\left|\mathbf{u}\left(\mathbf{s}\right)\right|}_{H^n}^2\right)}^{{\displaystyle \frac{1}{2}}}$$

$$+\mathcal{C}\left(N^{1-n}{\int }_0^t{\left|{\mathbf{u}}_t\left(s\right)\right|}_{H^{n-2}}^{2}ds\right)+\mathcal{C}{N}^{2-n}{\left|\mathbf{u}\right|}_{H^{n-2}}+\mathcal{C}{N}^{2-n}{\left|f^n\right|}_{H^{2-n}}.$$

Proof. 

To test its convergence, let us define $e_r\left(t\right)={\mathcal{R}}_N\mathbf{u}\left(t\right)-{\mathbf{u}}^N\left(t\right)$, where the projection operator is denoted by ${\mathcal{R}}_N$. The inequality is thus satisfied by the error function $e_r\left(t\right)$ as [22]

$${\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}{\parallel e_r\parallel }^2+\delta {\parallel e_r\parallel }_E^2\le |({\mathbf{u}}_t-{\mathcal{R}}_N{\mathbf{u}}_t,e_r)+(\mathbf{u}-{\mathcal{R}}_N\mathbf{u},e_r)|.$$

For any function $p\in X$, we can define the new norm as

$${\parallel p\parallel }_{E^{*}}=\underset{v\in E,v\ne 0}{sup}{\displaystyle \frac{(p,\mathbf{v})}{{\parallel \mathbf{v}\parallel }_E}}.$$

(8)

According to [22]$,{\parallel p\parallel }_{E^{*}}$ is the norm of p in the dual space $E^{*}$ of E. Keep in mind that for any ${\parallel p\parallel }_{E^{*}}\le \mathcal{C}{\parallel p\parallel }_X$, ${\parallel \mathbf{v}\parallel }_X\le \mathcal{C}{\parallel \mathbf{v}\parallel }_E$, for all $\mathbf{v}\in E$. From the formulation above and the continuity of the operator $\mathcal{L}$, it follows that

$$|({\mathbf{u}}_t-{\mathcal{R}}_N{\mathbf{u}}_t,e_r)+(\mathbf{u}-{\mathcal{R}}_N\mathbf{u},e_r)\le \mathcal{C}\left({\parallel {\mathbf{u}}_t-{\mathcal{R}}_N{\mathbf{u}}_t\parallel }_{E^{*}}+{\parallel \mathbf{u}-{\mathcal{R}}_N\mathbf{u}\parallel }_E\right){\parallel e_r\parallel }_E.$$

Thus, for each $t>0$, the following error bound can be inferred as

$${\parallel e_r\left(t\right)\parallel }_N^2+\delta {\int }_0^t{\displaystyle \frac{d}{dt}}{\parallel e_r\left(s\right)\parallel }_E^{2}ds\le {\parallel e_r\left(0\right)\parallel }_N^2$$

$$+\mathcal{C}\left({\int }_0^t{\parallel {\mathbf{u}}_t-{\mathcal{R}}_N{\mathbf{u}}_t\parallel }_{E^{*}}^{2}ds+{\int }_0^t{\parallel \mathbf{u}-{\mathcal{R}}_N\mathbf{u}\parallel }_E^{2}ds\right)+\mathcal{C}{\int }_0^t{\left({\displaystyle \frac{({\mathcal{R}}_N{\mathbf{u}}_t,e_r)-{({\mathcal{R}}_N{\mathbf{u}}_t,e_r)}_N}{{\parallel e_r\parallel }_E}}\right)}^{2}ds$$

$$+{\int }_0^t{\left({\displaystyle \frac{({\mathcal{L}}_N\mathbf{u},e_r)-{({\mathcal{L}}_N{\mathcal{R}}_N\mathbf{u},e_r)}_N}{{\parallel e_r\parallel }_E}}\right)}^{2}ds+\int {\left({\displaystyle \frac{(f,e_r)-{(f,e_r)}_N}{{\parallel e_r\parallel }_E}}\right)}^{2}ds.$$

(9)

Thus,

$${\parallel \mathbf{u}-{\mathbf{u}}^N\parallel }_E\le {\parallel e_r\left(0\right)\parallel }_N^2+\mathcal{C}{\left(N^{1-n}{\int }_0^t{\left|{\mathbf{u}}_t\left(s\right)\right|}_{H^{n-2}}^{2}ds+N^{1-n}{\left|\mathbf{u}\left(\mathbf{s}\right)\right|}_{H^n}^2\right)}^{{\displaystyle \frac{1}{2}}}$$

$$+\mathcal{C}\left(N^{1-n}{\int }_0^t{\left|{\mathbf{u}}_t\left(s\right)\right|}_{H^{n-2}}^{2}ds\right)+\mathcal{C}{N}^{2-n}{\left|\mathbf{u}\right|}_{H^{n-2}}+\mathcal{C}{N}^{2-n}{\left|f^n\right|}_{H^{2-n}}.$$

This completes the proof of convergence.    □

2.3. Fourier Spectral Method

In this section, motivated by the study [14,23,24], the Fourier spectral method was employed for the spatial integration of PDEs (1), considering a two- space dimensional case in a homogeneous medium ($\delta$ constants). Ref. [23] was followed to generate mathematical formula and computational algorithms in this study. A one-dimensional formula is similar to 2D cases, and it is basically straightforward. To that end, we consider the 2D system of PDEs in Equation (1). Since boundary conditions are imposed, Fast Fourier Transformation (FFT) was employed to approximate the spatial derivatives ($u_{xx}$, $u_{yy}$, $v_{xx}$, and $v_{yy}$) which are called spectral derivatives [14]. The problem is considered in a periodic domain to facilitate FFT and IFFTs.

Now, from the properties of the discrete Fourier transformation formulae [14], it is evident that

$$\begin{array}{cc}\hfill \mathcal{F}\left(u_{xx}(x,t)\right)& ={\left(i\xi \right)}^2\widehat{u}(\xi ,t)=-{\xi }^2\widehat{u}(\xi ,t),\hfill \\ \hfill \mathcal{F}\left(v_{xx}(x,t)\right)& ={\left(i\xi \right)}^2\widehat{v}(\xi ,t)=-{\xi }^2\widehat{v}(\xi ,t),\hfill \end{array}$$

and that

$$\begin{array}{cc}\hfill {\widehat{\mathbf{u}}}_{\xi \xi }^{\left(0\right)}& =\daleth {\widehat{\mathbf{u}}}^{\left(0\right)},{\widehat{\mathbf{v}}}_{\xi \xi }^{\left(0\right)}=\daleth {\widehat{\mathbf{v}}}^{\left(0\right)},\hfill \end{array}$$

where $\daleth =[-{\xi }_1^2-{\xi }_2^2\cdots -{\xi }_m^2]$, ${\xi }_k={\displaystyle \frac{2\pi }{L}}k$, for all $k=-{\displaystyle \frac{N}{2}}\left(1\right){\displaystyle \frac{N}{2}}-1,$ where L = length of space domain, and N = number of discretization points in space. For further detail, we refer the reader to [14]. Then, we obtain our approximate ${\tilde{\mathbf{u}}}_{xx}^{\left(0\right)}$ and ${\tilde{\mathbf{v}}}_{xx}^{\left(0\right)}$ vectors by applying the inverse FFT [24], i.e.,

$$\begin{array}{c}\hfill {\tilde{\mathbf{u}}}_{xx}^{\left(0\right)}=ifft\left({\widehat{\mathbf{u}}}_{\xi \xi }^{\left(0\right)}\right),{\tilde{\mathbf{v}}}_{xx}^{\left(0\right)}=ifft\left({\widehat{\mathbf{v}}}_{\xi \xi }^{\left(0\right)}\right).\end{array}$$

Let us recall the RDS (1) in 2D, where ${\delta }_i$, $i=1,2$ (considering a homogeneous medium) are constants

$$\begin{array}{cc}\hfill u_t& ={\delta }_1(u_{xx}+u_{yy})+f(u,v),\hfill \\ \hfill v_t& ={\delta }_2(v_{xx}+v_{yy})+g(u,v),\mathbf{x}\in \Omega ,t>0.\hfill \end{array}$$

(10)

Applying Fourier transformation on Equation (10) yields

$$\begin{array}{cc}\hfill {\widehat{u}}_t& =-{\delta }_1({\xi }_x^2+{\xi }_y^2)\widehat{u}+\widehat{f(u,v)},\hfill \\ \hfill {\widehat{v}}_t& =-{\delta }_2({\xi }_x^2+{\xi }_y^2)\widehat{v}+\widehat{g(u,v)},\hfill \end{array}$$

(11)

which is a first order time dependent semi-linear system of coupled ODE, which can be solved using integrating factors, where $\mathcal{F}\left(u(x,y,t)\right)=\widehat{u}({\xi }_x,{\xi }_y,t).$ The idea was used to integrate the problem in time, and then inverse Fourier transform (ifft) was employed to generate solutions in the spatio-temporal domain. By employing a semi-implicit scheme over time on (11), it can be approximated by

$$\begin{array}{cc}\hfill & {\widehat{u}}_{n+1}={\displaystyle \frac{{\widehat{u}}_n+\mathsf{\Delta }t{\widehat{f(u,v)}}_n}{1+\mathsf{\Delta }t{\delta }_1({\xi }_x^2+{\xi }_y^2)}},\hfill \\ \hfill & {\widehat{v}}_{n+1}={\displaystyle \frac{{\widehat{u}}_n+\mathsf{\Delta }t{\widehat{g(u,v)}}_n}{1+\mathsf{\Delta }t{\delta }_2({\xi }_x^2+{\xi }_y^2)}}.\hfill \end{array}$$

(12)

After approximating $\widehat{u}$ and $\widehat{v}$, the inverse Fourier transform (ifft) was employed to generate solutions in the spatio-temporal domain, where

$$\begin{array}{c}\hfill u(x,y,t)\approx {\mathcal{F}}^{-1}\left(\widehat{u}({\xi }_x,{\xi }_y,t)\right),\\ \hfill v(x,y,t)\approx {\mathcal{F}}^{-1}\left(\widehat{v}({\xi }_x,{\xi }_y,t)\right).\end{array}$$

It is well known that the convergence rate of the spectral method is $\mathcal{O}\left(N^{-n}\right)$ for any smooth function, and it is $\mathcal{O}\left(c^N\right)$ for an analytic function, where $c\in (0,1)$ [14]. To deduce these correlations, we apply the Fourier transform in two steps. First, a smooth function has a transform that decays swiftly. A smooth function changes slowly, and because high wave numbers correlate to rapidly oscillating waves, it has very little energy at such frequencies. Second, if a function’s Fourier transform decays quickly, the mistakes caused by discretization are minimal. These faults are generated by aliasing from high to low wave numbers [23].

2.4. Stability and Convergence of Fourier Collocation Method

We need to address a few preliminary issues that stem from [22] before we begin this discussion. Note that, we follow [22] and some references therein for the theoretical analysis in this section. A polynomial of degree N that is trigonometric due to Fourier collocation is represented by the space $S_N$:

$$S_N=span\left\{e^{(j+k)x}\right|-N\le k<N\}.$$

It is important to remember that the ${\mathbf{L}}_p$-norm of a function u over $(0,2\pi )$ has the following definition:

$${\parallel \mathbf{u}\parallel }_{{\mathbf{L}}^p(0,2\pi )}={\left({\int }_0^{2\pi }{\left|\mathbf{u}\left(x\right)\right|}^{p}dx\right)}^{{\displaystyle \frac{1}{p}}},1\le p<\infty ,$$

and

$${\parallel \mathbf{u}\parallel }_{{\mathbf{L}}^{\infty }(0,2\pi )}=\underset{0\le x\le 2\pi }{sup}\left|\mathbf{u}\left(x\right)\right|,p=\infty .$$

Again, we take the Sobolev norm of order $n\ge 0$ ($n\in \mathbb{Z}$), given by

$${\parallel \mathbf{u}\parallel }_{H^n(0,2\pi )}={\left(\sum _{k=0}^n{\int }_0^{2\pi }{\left|{\mathbf{u}}^k\left(x\right)\right|}^{2}dx\right)}^{{\displaystyle \frac{1}{2}}},$$

with the norm

$${\parallel \mathbf{u}\parallel }_n={\left(\sum _{k=-\infty }^{\infty }(1+{\left|k\right|}^{2n}){\left|\widehat{{\mathbf{u}}_k}\left(x\right)\right|}^2\right)}^{{\displaystyle \frac{1}{2}}}.$$

Here, the interpolant of the function u at the nodes $x_j={\displaystyle \frac{\pi j}{N}},j=0,\cdots ,2N-1$ is denoted by $I_N\mathbf{u}\in S_N$. We compute the approximation results for $\mathbf{u}-I_N\mathbf{u}$, the interpolation error. We have

$${\parallel \mathbf{u}-I_N\mathbf{u}\parallel }_{H^l(0,2\pi )}\le C{N}^{l-n}{\parallel {\mathbf{u}}^n\parallel }_{L^2(0,2\pi )}\mathrm{for}0\le l\le n,$$

$$and\mathbf{u}\in H_P^n(0,2\pi )\mathrm{with} n\ge 1.$$

(13)

Theorem 3.

Let u represent the exact solution to the problem with periodic boundary conditions that is suggested. Assume once more that the approximate solution of the Fourier collocation method is ${\mathbf{u}}^N\in S_N$, and that the function f is n times continuously differentiable. Next, the problem’s stability is described as

$${\parallel {\mathbf{u}}^N\parallel }^2+{\int }_0^t\delta {\nabla }^2{\parallel {\mathbf{u}}^N\parallel }_E^{2}ds\le {\parallel {\mathbf{u}}_0^N\parallel }^2+\mathcal{C}{\int }_0^t\parallel f^n\parallel ds,$$

where $\mathcal{C}$ is considered as a constant.

Proof. 

The collocation points are as follows: ${\mathbf{x}}_{rk}=(x_r,x_k),0\le r,k\le 2N-1$; $x_r={\displaystyle \frac{\pi r}{N}}$; the interpolant of u at these points is indicated by $I_N\mathbf{u}\in S_N$. The function ${\mathbf{u}}_N(\mathbf{x},t)={\sum }_{rk}{\mathbf{u}}_N(x_{rk},t){\varphi }_{rk}\left(x\right)$ represents the Fourier collocation approximation to u. ${\varphi }_{rk}$ represents the standard Lagrange trigonometric polynomials at the collocation points that fulfill the equations

$$[{\left({\mathbf{u}}_N\right)}_t-\mathcal{L}{\mathbf{u}}_N-f]\left({\mathbf{x}}_{rk,t}\right)=0\forall t>0\mathrm{and}0\le r,k\le 2N-1$$

$${\mathbf{u}}_N({\mathbf{x}}_{rk},0)={\mathbf{u}}_0\left({\mathbf{x}}_{rk}\right),\forall t>0\mathrm{and}0\le r,k\le 2N-1,$$

where ${\mathcal{L}}_N\mathbf{u}=\mathbf{\delta }{\nabla }^2\mathbf{u}$ for all $\mathbf{u}\in S_N$; this is an approximation for the interpolation of $\mathcal{L}\mathbf{u}$. The general form of this approach can be written by choosing

$${(\mathbf{u},\mathbf{v})}_N={\left({\displaystyle \frac{\pi }{N}}\right)}^2\sum _{0\le r,k\le 2N-1}\mathbf{u}\left({\mathbf{x}}_{rk}\right)\overline{\mathbf{v}\left({\mathbf{x}}_{rk}\right)}.$$

The stability of collocation is established by following the relations stated in [22]:

$${(Q_N{\mathcal{L}}_N\mathbf{u},\mathbf{u})}_N={({\mathcal{L}}_N\mathbf{u},\mathbf{u})}_N\forall \mathbf{u}\in X_N,$$

where $Q_N$ represents the orthogonal projection onto $S_N$ for the discrete inner product ${(\mathbf{u},\mathbf{v})}_N$. Additionally, $\mathbf{u}\in S_N$,

$${({\mathcal{L}}_N\mathbf{u},\mathbf{u})}_N={\int }_{\Omega }{I}_N(\delta \left(x\right)\nabla \mathbf{u})\mathbf{u}dxdy$$

$$={\int }_{\Omega }\delta \left(x\right){|\nabla \mathbf{u}|}^{2}dxdy$$

$$\ge {\delta }_0{\int }_{\Omega }{|\nabla \mathbf{u}|}^{2}dxdy.$$

In the Hilbert space $H_P^1(0,2\pi )$, as given in ([22]), the integral on the right-hand side is precisely the square of the norm ${\parallel \mathbf{u}\parallel }_m$.

Thus, using $E=H_P^n(0,2\pi )$, the stability requirement is confirmed, and the approximation is stable in accordance with

$${\parallel {\mathbf{u}}_N\parallel }_E\le {\displaystyle \frac{C}{{\delta }_0}}{\parallel f^n\parallel }_N$$

$$\alpha {\parallel {\mathbf{u}}^N\parallel }_E^2\le {({\mathcal{Q}}_n{\mathcal{L}}_N{\mathbf{u}}^N,{\mathbf{u}}^N)}_N={({\mathcal{Q}}_{N}f,{\mathbf{u}}^N)}_N\le {\parallel {\mathcal{Q}}_{N}f\parallel }_N{\parallel {\mathbf{u}}^N\parallel }_N\le \mathcal{C}{\parallel f^n\parallel }_N.$$

Thus,

$${\parallel {\mathbf{u}}^N\parallel }^2+{\int }_0^t\delta {\nabla }^2{\parallel {\mathbf{u}}^N\parallel }_E^{2}ds\le {\parallel {\mathbf{u}}_0^N\parallel }^2++\mathcal{C}{\int }_0^t\parallel f^n\parallel ds.$$

   □

Theorem 4.

Let u be the exact solution, according to the periodic boundary conditions, to the given problem. Let f be an n times continuously differentiable function. Assume again that the approximate solution of Equation (10) is ${\mathbf{u}}^N\in S_N$. According to the Fourier collocation method, the solution converges as follows:

$${\parallel \mathbf{u}-{\mathbf{u}}^N\parallel }_{H^1(0,2\pi )}\le {\parallel e\left(0\right)\parallel }_{L^2(0,2\pi )}+\mathcal{C}{\left(N^{1-n}{\int }_0^t{\left|{\mathbf{u}}^n\right|}_{L^2(0,2\pi )}ds+N^{1-n}{\left|{\mathbf{u}}^n\right|}_{L^2(0,2\pi )}\right)}^{{\displaystyle \frac{1}{2}}}$$

$$+\mathcal{C}\left(N^{1-n}{\int }_0^t{\left|{\mathbf{u}}^n\right|}_{L^2(0,2\pi )}ds\right)+\mathcal{C}{N}^{2-n}{\left|{\mathbf{u}}^n\right|}_{L^2(0,2\pi )}+\mathcal{C}{N}^{2-n}{\left|f^n\right|}_{L^2(0,2\pi )},$$

which holds.

Proof. 

For the Fourier collocation solution, the convergence estimate can be demonstrated by using Equation (9) and assuming $R_N\mathbf{u}=I_N\mathbf{u}$ as well. Then, applying the following approximation features of this operator, one can estimate the aliasing error using result (13).

$$\begin{array}{cc}\hfill {\parallel \mathbf{u}-I_N\mathbf{u}\parallel }_{H^1(0,2\pi )}& \le {\parallel \mathbf{u}-P_N\mathbf{u}\parallel }_{H^1(0,2\pi )}+{\parallel R_N\mathbf{u}\parallel }_{H^1(0,2\pi )}\hfill \\ \hfill & \le C{N}^{l-m}{\parallel {\mathbf{u}}^m\parallel }_{{\mathbf{L}}^2(0,2\pi )}+C{N}^l{\parallel R_N\mathbf{u}\parallel }_{{\mathbf{L}}^2(0,2\pi )}\hfill \\ \hfill & \le C{N}^{l-m}{\parallel {\mathbf{u}}^m\parallel }_{{\mathbf{L}}^2(0,2\pi )}.\hfill \end{array}$$

Consequently, the following error bound can be deduced for every $t>0$ as

$$\begin{array}{cc}\hfill & {\parallel e_r\left(t\right)\parallel }_{H^1(0,2\pi )}+\delta {\int }_0^t{\displaystyle \frac{d}{dt}}{\parallel e_r\left(s\right)\parallel }_{H^1(0,2\pi )}ds\hfill \\ \hfill & \le {\parallel e_r\left(0\right)\parallel }_{H^1(0,2\pi )}+\mathcal{C}\left({\int }_0^t{\parallel {\mathbf{u}}_t-R_N{\mathbf{u}}_t\parallel }_{H^1(0,2\pi )}ds+{\int }_0^t{\parallel \mathbf{u}-R_N\mathbf{u}\parallel }_{H^1(0,2\pi )}ds\right)\hfill \\ \hfill & +\mathcal{C}{\int }_0^t{\left({\displaystyle \frac{{(R_N{\mathbf{u}}_t,e)}_r-{(R_N{\mathbf{u}}_t,e_r)}_{H^1(0,2\pi )}}{{\parallel e_r\parallel }_{H^1(0,2\pi )}}}\right)}^{2}ds\hfill \\ \hfill & +{\int }_0^t{\left({\displaystyle \frac{({\mathcal{L}}_N\mathbf{u},e_r)-{({\mathcal{L}}_N{R}_N\mathbf{u},e_r)}_{H^1(0,2\pi )}}{{\parallel e_r\parallel }_{H^1(0,2\pi )}}}\right)}^{2}ds+\int {\left({\displaystyle \frac{(f,e_r)-{(f,e_r)}_{H^1(0,2\pi )}}{{\parallel e_r\parallel }_{H^1(0,2\pi )}}}\right)}^{2}ds.\hfill \end{array}$$

Thus,

$${\parallel \mathbf{u}-{\mathbf{u}}^N\parallel }_{H^1(0,2\pi )}\le {\parallel e\left(0\right)\parallel }_{{\mathbf{L}}^2(0,2\pi )}+\mathcal{C}{\left(N^{1-n}{\int }_0^t{\left|{\mathbf{u}}^n\right|}_{{\mathbf{L}}^2(0,2\pi )}ds+N^{1-n}{\left|{\mathbf{u}}^n\right|}_{{\mathbf{L}}^2(0,2\pi )}\right)}^{{\displaystyle \frac{1}{2}}}$$

$$+\mathcal{C}\left(N^{1-n}{\int }_0^t{\left(\left|{\mathbf{u}}^n\right|\right)}_{{\mathbf{L}}^2(0,2\pi )}ds\right)+\mathcal{C}{N}^{2-n}{\left|{\mathbf{u}}^n\right|}_{{\mathbf{L}}^2(0,2\pi )}+\mathcal{C}{N}^{2-n}{\left|f^n\right|}_{{\mathbf{L}}^2(0,2\pi )}.$$

This completes the proof of convergence.    □

3. A Multigrid Preconditioned Finite Difference Scheme

In this section, a few numerical schemes are employed to solve the class of system of heterogeneous PDEs in Equation  (1) with Dirichlet/Neumann boundary conditions. We start by employing finite difference schemes for spatial integration followed by a method of lines for the temporal approximation. Then, we employ an efficient solver for the resulting discrete model. A preconditioned Newton solver is studied here for the nonlinear system of equations.

We study a two-space dimensional version of (1). Let us consider a spatial domain $\Omega =[a,b]\times [c,d]$. Let $h_x={\displaystyle \frac{b-a}{n_x}},n_x\in \mathbb{N}$, $h_y={\displaystyle \frac{d-c}{n_y}},n_y\in \mathbb{N}$, $x_i=a+i{h}_x$, $i=0,1,2,\cdots ,n_x$, and $y_j=a+j{h}_y$, $j=0,1,2,\cdots ,n_y$. Let $u(x,y,t)\in C^{2,1}[\Omega ;[0,T]],T>0$. Then, a simple Taylor expansion for $u(x,y,t)\in C^{2,1}[\Omega ;[0,T]],T>0$ yields [1]

$$u_{xx}(x_i,y_j,t^n)\approx [u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n]/h_x^2,$$

$$u_{yy}(x_i,y_j,t^n)\approx [u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n]/h_y^2,$$

and

$$u_t(x_i,y_j,t^n)\approx {\displaystyle \frac{1}{\mathsf{\Delta }t}}[u_{i,j}^{n+1}-u_{i,j}^n],\mathrm{or}{u}_t(x_i,y_j,t^n)\approx {\displaystyle \frac{1}{\mathsf{\Delta }t}}[u_{i,j}^n-u_{i,j}^{n-1}],$$

where $\mathsf{\Delta }t=T/m,m\in \mathbb{N}.$ Let us recall the RDS (1) with heterogeneous diffusion coefficients (${\delta }_i(x,y),i=1,2$) of the form

$$\begin{array}{c}{\displaystyle \frac{\partial u}{\partial t}}={\delta }_1(x,y)▵u+f(u,v),\\ {\displaystyle \frac{\partial v}{\partial t}}={\delta }_2(x,y)▵v+g(u,v),\end{array}$$

(14)

where $▵u={\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}$, along with the Dirichlet/Neumann boundary conditions required by the problem modeled. Then at $(x_i,y_j)$, the system (14) can be approximated by the following implicit Euler scheme:

$$\begin{array}{cc}\hfill {\displaystyle \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\mathsf{\Delta }t}}& ={\displaystyle \frac{{\delta }_1(x_i,y_j)}{h^2}}[u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}-4u_{i,j}^{n+1}]+f(u_{i,j}^{n+1},v_{i,j}^{n+1}),\hfill \\ \hfill {\displaystyle \frac{v_{i,j}^{n+1}-v_{i,j}^n}{\mathsf{\Delta }t}}& ={\displaystyle \frac{{\delta }_2(x_i,y_j)}{h^2}}[v_{i+1,j}^{n+1}+v_{i-1,j}^{n+1}+v_{i,j+1}^{n+1}+v_{i,j+1}^{n+1}-4v_{i,j}^{n+1}]+g(u_i^{n+1},v_i^{n+1}),\hfill \end{array}$$

(15)

which is a full discrete system of nonlinear equations. Now, to facilitate the computations, two different semi-implicit schemes were employed, so that the equations can be solved separately. For example, two semi-implicit schemes are as follows:

• Semi-implicit scheme 1

  $$\begin{array}{cc}\hfill {\displaystyle \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\mathsf{\Delta }t}}& ={\displaystyle \frac{{\delta }_1(x_i,y_j)}{h^2}}[u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}-4u_{i,j}^{n+1}]+f(u_{i,j}^n,v_{i,j}^n),\hfill \\ \hfill {\displaystyle \frac{v_{i,j}^{n+1}-v_{i,j}^n}{\mathsf{\Delta }t}}& ={\displaystyle \frac{{\delta }_2(x_i,y_j)}{h^2}}[v_{i+1,j}^{n+1}+v_{i-1,j}^{n+1}+v_{i,j+1}^{n+1}+v_{i,j+1}^{n+1}-4v_{i,j}^{n+1}]+g(u_i^{n+1},v_i^n),\hfill \end{array}$$

  (16)

  and
• Semi-implicit scheme 2

  $$\begin{array}{cc}\hfill {\displaystyle \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\mathsf{\Delta }t}}& ={\displaystyle \frac{{\delta }_1(x_i,y_j)}{h^2}}[u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}-4u_{i,j}^{n+1}]+f(u_{i,j}^{n+1},v_{i,j}^n),\hfill \\ \hfill {\displaystyle \frac{v_{i,j}^{n+1}-v_{i,j}^n}{\mathsf{\Delta }t}}& ={\displaystyle \frac{{\delta }_2(x_i,y_j)}{h^2}}[v_{i+1,j}^{n+1}+v_{i-1,j}^{n+1}+v_{i,j+1}^{n+1}+v_{i,j+1}^{n+1}-4v_{i,j}^{n+1}]+g(u_i^{n+1},v_i^{n+1}).\hfill \end{array}$$

  (17)

The Euler formula in time and a five–point formula in space for one- and two–dimensional models give an computational accuracy of $\mathcal{O}(h^2+\mathsf{\Delta }t)$. Instead, the use of the Crank–Nicolson scheme is accurate for $\mathcal{O}(h^2+\mathsf{\Delta }{t}^2)$. For the explicit Euler approximation, the stability of the scheme needs to follow the Courant–Friedrichs–Lewy (CFL) conditions for choosing h and $\mathsf{\Delta }t$, whereas the implicit Euler scheme is unconditionally stable [1,25,26].

We are interested in approximating the solutions of the full discrete model efficiently by using a fast Newton-based algorithm. In this case, $u_{i,j}^{n+1}$ is approximated first, followed by $v_{i,j}^{n+1}$, which is easier to execute than utilizing Newton’s technique for (15). The primary challenge in implementing Newton’s method for (15), (16), or (17) is efficiently solving the resultant system of linear equations at each step. To be more explicit, the inverse of the Jacobian matrix must be effectively approximated in each step of Newton’s algorithm. Iterative solvers are commonly used to estimate the inverse of the Jacobian of the nonlinear equations. The most important step in completing such an algorithm for our issue is efficiently solving the resultant system of linear equations at each stage. We compare two iterative strategies for solving the system of linear equations that develops during each iteration of the Newton algorithm for nonlinear systems. In short, the main part in each step of Newton’s algorithm for a nonlinear system $F\left(x\right)=0$ works in the following way:

Step 1:

Solve $DF\left(x_i\right)\mathsf{\Delta }{x}_i=-F\left(x_i\right)$, where $DF\left(x_i\right)$, $i=0,1,2,\cdots$ stands for the Jacobian of F at $x_i\in {\mathbb{R}}^n,n\ge 1$.

Step 2:

Update $x_{i+1}=x_i+\mathsf{\Delta }{x}_i.$

In this case, any iterative solvers inside the Newton algorithm oscillate or converge very slowly. This is because, while the minimal eigenvalue is bounded, the maximal eigenvalue of the discrete operator acting on each equation of (17) or (16) grows exponentially. The condition number $cond\left(\mathcal{A}\right)=\mathcal{O}\left(N^2\right)=\mathcal{O}\left(2^{2l}\right)$ grows exponentially, in this case, for some $l>1$. Here, $\mathcal{A}$ represents the Jacobian of the equation parts of the nonlinear system (16) or (17). Thus, preconditioning is necessary to tackle such problems [27]. A multigrid is a very good candidate for such a situation [28]. So here, we study and employ a multigrid preconditioned GMRES solver on each iteration of the Newton algorithm [1] for the solution of a nonlinear system of equations.

The multigrid scheme is a numerical method used for solving partial differential equations (PDEs), particularly in the context of computational fluid dynamics and other areas requiring efficient and accurate solutions to complex problems. The scheme is designed to address challenges related to large-scale computations, including issues of stability, accuracy, and computational efficiency. This scheme can be used within the framework of the generalized minimal residual method (GMRES), which is an iterative method, as a preconditioning to speed up the computational process [27,28]. Here, we approximate the inverse of the Jacobian of each nonlinear system of equations we encountered, while calculating the solutions using a multigrid preconditioned GMRES algorithm. We refer to it as the Newton multigrid accelerated GMRES solver [27,28].

In this instance, the plan is to carry out a few smoothing iterations in a fine grid, then move to a coarser level and carry out a few more iterations, and so on. We refer to these as coarse grid adjustments. Following adjustments, one returns to the fine grid and carries out a few post-smoothing operations. Using the Newton technique for the nonlinear system given above, we use the two-grid iterative linear system solution technique as a preconditioner inside the GMRES solver within the MATLAB R2024b environment using the following algorithm:

Multigrid Algorithm

In short, a multigrid scheme works in the following way:

• Grid Hierarchy: In the multigrid approach, the problem is solved on a hierarchy of grids. The coarse grids capture the large-scale features of the solution, while the finer grids handle the small-scale details. By transferring information between these grids, the method can efficiently reduce errors and improve the convergence rates.
• Smoothing and Correction: The multigrid scheme involves two main components:
  • Smoothing: This step reduces the high-frequency errors on a given grid. It usually involves iterative solvers like the Gauss–Seidel or Jacobi relaxation methods.
  • Correction: After smoothing, the residual (error) is corrected using information from coarser grids. This process helps address low-frequency errors that are not effectively reduced by smoothing alone.

It can be implemented using the following Algorithm 1.

Algorithm 1 Multigrid V-Cycle method for $\mathcal{A}\underline{d}=\underline{f}$

Here is a more formal outline of the multigrid algorithm for solving $\mathcal{A}\underline{d}=\underline{f}$. Initialization  Start with an initial guess $d^{\left(0\right)}$ for the solution. OUTPUT  The approximate solutions ${\underline{d}}^h.$ Step 1  Pre-Smoothing: Apply a smoothing or relaxation technique on the fine grid to reduce high-frequency errors (apply k iterations of a smoothing method to $d^{\left(i\right)}$). Step 2  Residual Computation: Compute the residual vector ($r^{\left(i\right)}=b-\mathcal{A}{d}^{\left(i\right)}$). This represents the error in the current approximation. Step 3  Restriction: Restrict $r^{\left(i\right)}$ to a coarser grid to obtain $r_{coarse}$. Step 4  Coarse Grid Correction: Solve the linear system ${\mathcal{A}}_{coarse}{d}_{coarse}=r_{coarse}$ approximately on the coarse grid to obtain a correction vector. This solution addresses the low-frequency errors that are not well-captured by the fine grid. Step 5  Prolongation: Transfer (or prolong) the correction vector $d_{coarse}$ from the coarse grid back to the fine grid and add the result to $d^{\left(i\right)}$. This step interpolates the coarse grid solution to the finer grid and adds it to the current approximation. Step 6  Post-Smoothing: Apply k iterations of the smoothing technique again on the fine grid to the updated $d^{\left(i\right)}$ to further reduce any remaining high-frequency errors after the correction. Then, check whether the residual norm $\parallel r^{\left(i\right)}\parallel$ is below a predetermined tolerance. If not, iterate. STOP

4. Numerical Results and Discussions

Here, we present a few experimented results for the schemes studied earlier in this article. Graphical and tabular illustration gives us a clear view of the problem and the acceptability of the solutions generated. It helps us to analyze the physical reality neatly. In the previous sections, a few class of schemes were studied. Here, we are motivated to demonstrate approximate solutions of the RDS by the numerical schemes discussed in the previous sections. The results in one- and two-space dimensions are presented here in the preceding segments of the section. To obtain the solution to the problems, a few linear system solvers like the multigrid preconditioned Newton–GMRES (generalized minimal residual) and Newton–ILU factorization algorithm are employed.

Here, from the graphical representations in Figure 1 and Figure 2 for the BBP schemes and overall, it is observed that the outputs are numerically convergent. Here, for all the computations the same set of parameter values were used. All the schemes can be good alternatives for further investigations for the problems in population dynamics, following [4]. Note that, for the experiments presented in Figure 1 and Figure 2, the parameter values $\alpha =0.3$, $\beta =2.0$, $\gamma =0.8$, $\delta =1.0$, $a=0$, $b=100$ $h=1/2$, $u_0\left(x\right)=exp(-(x-0.5b){.}^2)$, and $v_0\left(x\right)=2/5$ were considered. As shown in Figure 1, the computation was performed using the Bernstein polynomials, and both the schemes are compatible for the model problem with a long time interval. In Figure 2, the numerical results obtained by the full implicit finite difference method, the Bernstein spectral method, and the Fourier spectral methods are compared. All the schemes are compatible, and the Bernstein spectral method outperforms the other two in terms of storage consideration, as it requires fewer numbers of space grid points compared to the other two schemes. For all the one-space dimensional cases, Kinetics 1 was considered. Here, in Figure 3 and Figure 4, the results for heterogeneity are presented, and the Figures show the effect of inhomogeneity in the RDS. Then, we move to the two space-dimensional cases in Figure 5. Figure 5 and Figure 6 confirm that for two space dimensions, the schemes are shown off apparent convergence to each other, which is a positive approach to use for further investigations in population dynamics [4]. In Figure 5, the solution dynamics illustrate that these schemes in two-dimensional space need some significant analysis to make the system user-friendly. In Figure 6 and Figure 7, we show the same experiment with different choices of the diffusion coefficient, $\alpha$, and initial data to further demonstrate the solution’s behavior. In Figure 8, we use a system size $[0,2]\times [0,2]$. The self-repeated replication of the initial spots is observed in the solution of the Gary–Scott model with the parameters. This type of self-replicating pattern is formed due to the complex interplay between the activator and inhibitor chemicals aided by the reaction and diffusion components. This type of development is found in proto-organisms. Here, the entire space is occupied by making replicas of cell-like localized structures.

Then, we move on to implement a multigrid preconditioned Newton–GMRES scheme for the full nonlinear system of equations (semi-implicit finite difference scheme 1 and semi-implicit finite difference scheme 2). The same experiments can be performed for the other nonlinear systems developed in this study. Here, the schemes were implemented within the framework of MATLAB. It is well observed from the earlier study [28] that the multigrid accelerated scheme performs better for the homogeneous RDS. So, we are interested in implementing and computing the results for GMRES and the multigrid accelerated GMRES solver only for the heterogeneous counterpart, to compare our results with that of the unaccelerated GMRES. Here, from the Tables (

Table 1. Here, we demonstrate the total iteration numbers to reach T. We consider the semi-implicit scheme 1 with nonlinearity type 1 and ${\delta }_1(x,y)={\delta }_2(x,y)=0.0001(1+exp(-{(x-0.5)}^2-{(y-0.5)}^2)$.

$\mathcal{N}$$\downarrow |\mathsf{\Delta }\mathit{t}\to$	0.01	0.1	0.2	0.5
$u_h(x,y,T)$, GMRES	6375	1156	860	630
$u_h(x,y,T)$, MGGMRES	3800	610	480	262
$v_h(x,y,T)$, GMRES	4220	640	295	262
$v_h(x,y,T)$, MGGMRES	2250	422	176	242

,

Table 2. Here, we demonstrate the total iteration numbers to reach T. We consider the semi-implicit scheme 2, nonlinearity type 1, and ${\delta }_1(x,y)={\delta }_2(x,y)=0.0001(1+exp(-{(x-0.5)}^2-{(y-0.5)}^2)$.

$\mathcal{N}$$\downarrow |\mathsf{\Delta }\mathit{t}\to$	0.001	0.01	0.05	0.1	0.2	0.5	1
$u_h(x,y,T)$, GMRES	132,672	NC	NC	NC	NC	NC	NC
$u_h(x,y,T)$, MGGMRES	12,320	4476	2015	780	450	115	67
$v_h(x,y,T)$, GMRES	121,482	NC	NC	NC	NC	NC	NC
$v_h(x,y,T)$, MGGMRES	12,320	4432	1286	760	540	115	68

and

Table 3. Number of iterations needed to reach T. We consider the semi-implicit scheme 2, nonlinearity type 2, and ${\delta }_1(x,y)={\delta }_2(x,y)=0.0001(1+exp(-{(x-25)}^2-{(y-25)}^2)$.

$\mathcal{N}$$\downarrow |\mathsf{\Delta }\mathit{t}\to$	0.01	0.1	0.2	0.5	1
$u_h(x,y,T)$, GMRES	4575	2286	475	NCon	NCon
$u_h(x,y,T)$, MGGMRES	750	430	242	162	112
$v_h(x,y,T)$, GMRES	4424	2254	1054	NCon	NCon
$v_h(x,y,T)$, MGGMRES	650	424	255	168	108

) presented below, it is observed that multigrid preconditioned Newton–GMRES outperforms the Newton–GMRES scheme for both the choices of nonlinearity and the semi-implicit schemes. In fact, the multigrid preconditioned Newton–GMRES works better than the solver without preconditioning. Note that, for Table 1 and Table 2, we consider $\Omega ={[0,1]}^2$, $T=1$, $h=0.01$, and $tol=1\times {10}^{-10}$. In the tables,

• MGGMRES means the Newton multigrid preconditioned GMRES scheme was used to solve the resulting full discrete system,
• GMRES means the unaccelerated Newton GMRES scheme was used to solve the resulting full discrete system,
• “NC” means the Newton GMRES did not converge as desired [11].

In Table 3, we consider $\omega =[0,50]\times [0,50]$, $h=0.1$, and $tol={10}^{-10}$. Here, $\mathcal{N}$ means the number of iterations.

5. Conclusions

In this study, a system of heterogeneous PDEs/RDS (1) was considered and studied numerically to observe the solution dynamics. Three different classes of schemes, namely the Bernstein polynomial schemes, Fourier pseudo-spectral scheme, and a few finite difference-based schemes were employed for spatial approximations. The resulting semi-discrete model was investigated using a few method of lines for time integration. By studying the numerical solutions of (1) by three different classes of schemes, it has been demonstrated that the Bernstein polynomial scheme works very well for the heterogeneous medium RDS. The Fourier transform schemes (FTS) for such models are one of the smartest schemes, if it can be applied, as it is very simple, easy to use, and a fast scheme. The main drawback of FTS is the use of convolution in the Fourier domain while working with the spatially variable medium models. It was confirmed theoretically that the spectral/pseudo-spectral methods are very efficient and highly accurate for such a model. The schemes were illustrated through several graphical representations of solutions. From our study, in one dimension, through experiments, all the schemes we considered in this study were compatible. The semi-explicit schemes require smaller $\mathsf{\Delta }t$ compared to the full implicit scheme, as expected, and the spectral methods are well converged for a reasonably small number of spacial points compared to the finite difference schemes. It is noted that the multigrid preconditioned Newton–GMRES solver for a full discrete nonlinear system outperformed the Newton–GMRES solver in terms of iteration numbers. So, to study such kinds of heterogeneous problems, one may use this idea to generate the solutions quickly to analyze their real-life problems. Here, the Fourier scheme for the heterogeneous problem remains as a future research direction. The three-space dimensional problems and effects of the nonlocal behavior of the problems were not considered in this study and, thus, remain a future research direction.