On the Spatiotemporal Pattern Formation in Nonlinear Coupled Reaction–Diffusion Systems

Abstract

Nonlinear coupled reaction–diffusion (NCRD) systems have played a crucial role in the emergence of spatiotemporal patterns across various scientific and engineering domains. The NCRD systems considered in this study encompass various models, such as linear, Gray–Scott, Brusselator, isothermal chemical, and Schnakenberg, with the aim of capturing the spatiotemporal patterns they generate. These models cover a diverse range of intricate spatiotemporal patterns found in nature, including spots, spot replication, stripes, hexagons, and more. A mixed-type modal discontinuous Galerkin approach is employed for solving one- and two-dimensional NCRD systems. This approach introduces a mathematical formulation to handle the occurrence of second-order derivatives in diffusion terms. For spatial discretization, hierarchical modal basis functions premised on orthogonal scaled Legendre polynomials are used. Moreover, a novel reaction term treatment is proposed for the NCRD systems, demonstrating an intrinsic feature of the new DG scheme and preventing erroneous solutions due to extremely nonlinear reaction terms. The proposed approach reduces the NCRD systems into a framework of ordinary differential equations in time, which are addressed by an explicit third-order TVD Runge–Kutta algorithm. The spatiotemporal patterns generated with the present approach are comparable to those found in the literature. This approach can readily be expanded to handle large multi-dimensional problems that appear as model equations in developed biological and chemical applications.

1. Introduction

Over the last few decades, the studies of spatiotemporal pattern formation in nature have been fascinating subjects for different areas of scientific research. These research areas include biology [1], ecology [2], semiconductor physics [3], material science [4], hydrodynamics [5], astrophysics [6], computational social sciences [7], ecological networks [8], and many more. In 1952, mathematician Alan Turing [9] developed a nonlinear coupled reaction–diffusion (NCRD) system-based theory that featured an activator–inhibitor pair to describe the early development of spatiotemporal patterns in an embryo. In this hypothesis, he looked at how the activator’s and inhibitor’s diffusivities differed, leading to spontaneous symmetry breaking and spatially different stable patterns. After that, different types of Turing-type spatiotemporal pattern formations have been widely discovered in biology, chemistry, and population dynamics. Typical examples of biological applications of spatiotemporal patterns are hydra skin patterns [10], animal coat patterns [11], butterfly wing patterns [12], zebra coat patterns [13], skeletal patterning in limb evolution [14], and sea shell coloring patterns [15]. Several mathematical models based on spatiotemporal patterns utilized in biology, ecology, and biochemistry involve organisms reacting in the presence of diffusion, resulting in the NCRD system. A few examples of well-known NCRD models are Gierer and Meinhardt [10], Brusselator [16], Schnakenberg [17], isothermal [18], and Gray–Scott [19], and there are many others.

The focus of this study is restricted to two-dimensional NCRD systems containing two chemical species, which are described in a general form as

$$\frac{\partial \mathbf{U}}{\partial t}=\mathbf{D}{\nabla }^2\mathbf{U}+\mathbf{S}\left(\mathbf{U}\right),\forall (\mathbf{x},t)\in \Omega \times [0,T],$$

(1)

with initial conditions

$$\mathbf{U}(\mathbf{x},0)={\mathbf{U}}_0\left(\mathbf{x}\right),\forall \mathbf{x}\in \Omega ,$$

(2)

and Dirichlet or Neumann boundary conditions. In the above-mentioned expression,

$$\mathbf{U}=\left(\begin{array}{c}u\\ v\end{array}\right),\mathbf{D}=\left(\begin{array}{cc}{\mu }_1& 0\\ 0& {\mu }_2\end{array}\right),\mathbf{S}\left(\mathbf{U}\right)=\left(\begin{array}{c}f(u,v)\\ g(u,v)\end{array}\right),\mathrm{and}\mathbf{x}=(x,y)\in \Omega ,$$

(3)

where $\mathbf{U}$ is the unknown vector of two chemical-species components, $u,v$; $\mathbf{D}$ is the diffusion constant matrix with ${\mu }_1,{\mu }_2$ parameters; and ${\nabla }^2$ is the Laplacian operator, while $\mathbf{S}\left(\mathbf{U}\right)$ represents the reaction matrix of two nonlinear reaction kinetics $f(u,v)$ and $g(u,v)$. Notably, both an external force and a coupling between diffusion coefficients are not taken into consideration in this work. The most general form of the nonlinear reaction kinetics might be represented as

$$\begin{array}{cc}\hfill f(u,v)& =a_1{u}^{2}v+a_{2}u{v}^2+a_{3}uv+a_{4}u+a_{5}v+a_6,\hfill \\ \hfill g(u,v)& =b_1{u}^{2}v+b_{2}u{v}^2+b_{3}uv+b_{4}u+b_{5}v+b_6,\hfill \end{array}$$

(4)

where $a_i$, $b_i$, $i=1,6$ represent the genuine parameters. The above-mentioned expression (4) encompasses a number of well-known mathematical models that are taken into consideration in this study, including Gray–Scott [19], Brusselator [16], Schnakenberg [17], and isothermal [18] models. The genuine parameters for these models are chosen as

• Gray–Scott: $a_1=a_3=a_5=0,a_2=-1,a_4=-k_1<0,a_6=k_1;b_1=b_3=b_4=b_6=0,$$b_2=1,b_5=-(k_1+k_2)<0$.
• Brusselator: $a_1=1,a_2=a_3=a_5=0,a_4=-(1+k_2),a_6=k_1;b_1=-1,b_2=b_3=b_5=b_6=0,b_4=k_2$.
• Schnakenberg: $a_1=1,a_2=a_3=a_5=0,a_4=-1,a_6=k_1;b_1=-1,b_2=b_3=b_4=b_5=0,b_6=k_2$.
• Isothermal: $a_1=a_2=a_4=a_5=a_6=0,a_3=-1;b_1=b_2=b_4=b_6=0,b_3=1,$$b_5=-k_1$.

To understand and predict the complex spatiotemporal dynamical structures, including the traveling wave and self-organized patterns found in nature such as spots, spot replication, stripes, hexagons, and other dissipative solutions [19,20,21], the mathematical modeling and numerical simulation of the NCRD systems are essential. The solutions to these systems have been a significant and demanding area of research because of the challenging combination of the stiff diffusion term and the highly nonlinear reaction term. In the field of computational biology, it is necessary to have a practical and trustworthy numerical method for solving NCRD systems during parametric and numerical investigations.

In recent decades, numerous research efforts have been devoted to developing numerical methods that can effectively approximate the solutions of NCRD systems with high accuracy. Ruuth [22] used implicit–explicit schemes in conjunction with spectral methods to describe the spatiotemporal patterns in biological and chemical processes. Zegeling and Kok [23] employed an adaptive moving mesh approach for the spatiotemporal patterns of one- and two-dimensional NCRD systems with chemistry applications. Lopez and Ramos [24] used linearized methods in conjunction with operator-splitting techniques to obtain the solutions of the isothermal NCRD system. Madzvamuse [25] employed moving grid-based finite element approaches to derive the spatiotemporal patterns of the NCRD systems on fixed and continually increasing domains using time-stepping procedures. Mittal and Rohila [26] adopted the modified cubic B-spline-based differential quadrature approach for obtaining the spatiotemporal patterns of various one-dimensional NCRD systems. Subsequently, Jiwari et al. [27] captured the different types of spatiotemporal patterns in several NCRD systems through modified trigonometric cubic B-spline functions, which are based on the differential quadrature approach. Later, Yadav and Jiwari [28] illustrated a finite element scheme to investigate the computation simulations for the spatiotemporal patterns of NCRD systems. Tamsir et al. [29] proposed a hybrid numerical method of cubic trigonometric B-spline base functions and a differential quadrature algorithm for solving Fisher-type NCRD systems. Onarcan et al. [30] presented a trigonometric cubic B-spline-based collocation scheme for obtaining the spatiotemporal dynamics of NCRD systems. Recently, Manzoor [31] presented an effective and simple approximation scheme to analytically approximate the solution behavior of a multi-dimensional Brusselator NCRD system. Up to this point, the literature has documented several widely recognized computational techniques such as finite difference, positive finite volume, finite volume spectral element, and meshless local Petrov–Galerkin methods, which were employed to investigate the formation of spatiotemporal patterns in biological systems [32,33,34,35].

In the present study, the authors utilized the discontinuous Galerkin (DG) method to examine the formation of spatiotemporal patterns in NCRD systems. In 1973, Reed and Hill [36] introduced the DG method as a solution for neutron transport, and it has since gained widespread popularity as a versatile approach for approximating various types of partial differential equations. This method employs a completely discontinuous piece-wise polynomial space for both the numerical approximations and the test functions. Over the past years, Cockburn et al. devised a framework for solving PDEs-based hyperbolic conservation laws in a series of works [37,38,39,40,41,42] that contributed significantly to the development of the DG technique. This method has been effectively applied to several problems, including fluid flows, multi-phase flows, quantum physics, magneto-hydrodynamics, and many others [43,44,45,46,47,48,49,50]; it has also recently made its way into developmental biology. The motivation behind this implementation stems from progress in the approach and recent breakthroughs in generating spatiotemporal patterns. These advancements have expanded the applicability of the DG method to a wider range of biological and chemical needs.

Zhu et al. [51] merged the DG approach along with the Strang-type symmetrical operator splitting methodology to address NCRD systems. Subsequently, this technique has been applied to moving meshes [52]. Zhang et al. [53] provided a direct DG method, which is attributed to the direct weak formulation of the NCRD system and the construction of an appropriate numerical flux scheme on the elemental boundaries. Later, Zhang et al. [54] illustrated a newly nonlinear Galerkin approach subjected to finite element discretization for the numerical approximations of NCRD systems on wide intervals. In addition to these investigations, Singh [55] introduced an alternative approach known as the mixed-type modal DG method for solving general FitzHugh–Nagumo NCRD equations. Later, this work was extended to Fisher-KPP-type NRCD equations [56,57]. Recently, Singh et al. [58] presented high-fidelity numerical simulations for capturing the Turing pattern in a multi-dimensional Gray–Scott reaction–diffusion system.

The objective of this study is to establish a systematic framework utilizing the DG method to simulate two-dimensional NCRD systems. The primary focus of this framework lies in addressing the second-order derivatives present in the diffusion term and handling the nonlinear reaction term in the NCRD systems using the Galerkin framework. Our focus is on the development of an explicit mixed-type modal DG method, utilizing two-dimensional structured meshes, to study spatiotemporal pattern formation. Furthermore, we introduce a novel approach to handle the reaction terms in NCRD systems, showcasing an inherent characteristic of the new DG scheme. We employ a spatial discretization technique that utilizes hierarchical modal ansatz functions based on orthogonal scaled Legendre polynomials. This approach differs from previous research conducted on NCRD systems. Lastly, the proposed technique is innovative for one- and two-dimensional NCRD systems, specifically Gray–Scott, Brusselator, Schnakenberg, and isothermal models. A graphical abstract for modeling the NCRD systems is illustrated in Figure 1.

The outlines of this study are as follows. In Section 2, the mixed-type DG spatial discretization is reported in detail for solving the NCRD system. A novel reaction term treatment is also presented, demonstrating an intrinsic feature of the new DG scheme. In Section 3, we demonstrate the accuracy and validation of the proposed numerical approach for the NCRD system. In Section 4, the proposed scheme is employed in 1D and 2D NCRD systems utilized in developed biology to generate spatiotemporal patterns. We demonstrate the influence of changing parameters on the resulting patterns through the various simulations of these models. Finally, in Section 5, a conclusion is offered, along with recommendations for further research in light of the current study.

2. Mixed-Type Modal Discontinuous Galerkin Method

This section demonstrates the numerical procedure for solving NCRD systems (1) using a mixed-type modal discontinuous Galerkin (DG) formulation [49,58,59]. In this mixed formulation, a new additional auxiliary variable, $\Theta$, is introduced to handle the second-order derivatives appearing in the diffusion term. This auxiliary variable is defined as the derivative of the unknown variable, $\mathbf{U}$, in the NCRD systems (1). To apply this formulation, the NCRD system (1) is regarded as a connected framework of $\mathbf{U}$ and $\Theta$.

$$\begin{array}{cc}\hfill \Theta -\mathbf{D}\nabla \mathbf{U}& =0,\hfill \\ \hfill \frac{\partial \mathbf{U}}{\partial t}-\nabla \mathbf{F}(\Theta )& =\mathbf{S}\left(\mathbf{U}\right),\hfill \end{array}$$

(5)

where $\mathbf{F}(\nabla \mathbf{U})=\nabla \mathbf{U}$ is the diffusive term.

2.1. Discretization Process in Modal DG Framework

In this section, we introduce a simplified approach for the spatial discretization of two-dimensional NCRD systems based on structured meshes. Notably, one-dimensional NCRD systems can also be addressed using a similar approach. To discretize the governing Equation (1), the computational domain $D\in {\Re }^2$ is subdivided into non-overlapping square elements, $E_1,\cdots ,E_K$, and each of these elements is mapped to the reference element, $E_0={[-1,1]}^2$. The bilinear mappings, $\overrightarrow{X_k}:E_0\to E_k$ with $\overrightarrow{X_k}(\xi ,\eta )={(x,y)}^T$, are defined as

$$\begin{array}{c}\hfill \overrightarrow{X_k}(\xi ,\eta )=\frac{1}{4}[\overrightarrow{x_1}(1-\xi )(1-\eta )+\overrightarrow{x_2}(1-\xi )(1+\eta )\\ \hfill +\overrightarrow{x_3}(1+\xi )(1+\eta )+\overrightarrow{x_4}(1+\xi )(1-\eta )],\end{array}$$

(6)

for $k=1,\cdots K$, where $\{\overrightarrow{x_1},\overrightarrow{x_2},\overrightarrow{x_3},\overrightarrow{x_4}\}$ are the four corners of the element, $E_k$. Due to the uniform Cartesian mesh, the Jacobian of the adopted mappings is to be constant, i.e.,

$${\mathbf{J}}^v=\frac{1}{4}\Delta x\Delta y=\frac{\Delta x^2}{4},$$

(7)

for equal element side lengths, $\Delta x=\Delta x_k=\Delta y_k=\Delta x$, $k=1,\cdots ,K$.

Further, we define a finite element space as

$$V_h^l=\{\nu \in \left[L^2\left(D\right)\right]:\nu |{}_{E_k}\in \left[P^l\left(E_k\right)\right],k=1,\cdots ,K\},$$

(8)

where $P^l\left(E_k\right)$ denotes the set of all polynomials of degree at most l on $E_k$. Now, on every element, the solutions of $\mathbf{U}$ and $\Theta$ are approximated by a polynomial of degree, l, in each spatial direction on $E_0$, such as

$$\begin{array}{cc}\hfill {\Theta }_h{(x,y,t)|}_{E_k}& ={\Theta }_h(\xi ,\eta ,t)\approx \sum _{i=0}^{N_l-1}{\Theta }_h^i\left(t\right){\phi }_i(\xi ,\eta ),\hfill \\ \hfill {\mathbf{U}}_h{(x,y,t)|}_{E_k}& ={\mathbf{U}}_h(\xi ,\eta ,t)\approx \sum _{i=0}^{N_l-1}{\mathbf{U}}_h^i\left(t\right){\phi }_i(\xi ,\eta ),\hfill \end{array}$$

(9)

where ${\Theta }_h^i$ and ${\mathbf{U}}_h^i$ are the local time-dependent degrees of freedom for $\mathbf{U}$ and $\Theta$, respectively, in the considered element, $E_k$. The number of required basis functions for $l–$exact DG approximation is calculated by $N_l=(l+1)(l+2)/2$. For this work, the interpolating 1D and 2D orthogonal scaled Legendre basis functions are adopted, which are also summarized in Appendices Appendix A.1 and Appendix A.2. Interestingly, the advantages of scaled Legendre basis functions lie in their flexibility, orthogonality, numerical stability, efficiency, convergence properties, and compatibility with numerical techniques. These properties make them a valuable tool for approximating functions defined on non-standard intervals and solving a wide range of numerical problems.

At this stage, the discretization process can be started by establishing the weak DG form of the NCRD systems. This form is obtained by multiplying Equation (5) by a test function, ${\phi }_h$, integrating over the domain, $E_k$, and performing an integration by parts:

$$\begin{array}{cc}\hfill \underset{E_k}{\int }{\Theta }_h{\phi }_{h}dV+\mathbf{D}\underset{E_k}{\int }\nabla {\phi }_h·{\mathbf{U}}_{h}dV-\mathbf{D}\underset{\partial E_k}{\int }{\phi }_h{\mathbf{U}}_h·\mathbf{n}d\Gamma & =0,\hfill \\ \hfill \frac{\partial }{\partial t}\underset{E_k}{\int }{\mathbf{U}}_h{\phi }_{h}dV+\underset{E_k}{\int }\nabla {\phi }_h·\mathbf{F}\left({\Theta }_h\right)dV-\underset{\partial E_k}{\int }{\phi }_h\mathbf{F}\left({\Theta }_h\right)·\mathbf{n}d\Gamma & =\underset{E_k}{\int }{\phi }_h\mathbf{S}\left({\mathbf{U}}_h\right)dV.\hfill \end{array}$$

(10)

where $\mathbf{n}$ denotes the outward unit normal vector. V and $\Gamma$ are the volume and boundary of the rectangular element, $E_k$, respectively. The numerical solutions ${\mathbf{U}}_h$ and ${\Theta }_h$ are discontinuous across element interfaces; therefore, the interface fluxes are not uniquely defined. The flux functions ${\mathbf{U}}_h·\mathbf{n}$, and $\mathbf{F}\left({\Theta }_h\right)·\mathbf{n}$ emerging in Equation (10) are substituted by a numerical flux function illustrated at the cell interfaces, which are denoted by ${\mathbf{H}}^{aux}$ and ${\mathbf{H}}^{alt}$, respectively. In this work, we used an alternating scheme [40] for viscous and auxiliary flux functions, illustrated as

$$\begin{array}{cc}\hfill {\mathbf{H}}^{aux}({\mathbf{U}}_h^{int},{\mathbf{U}}_h^{ext})& =\frac{1}{2}\left[{\mathbf{U}}_h^{int}+{\mathbf{U}}_h^{ext}\right]+C_{11}\left[{\mathbf{U}}_h^{int}-{\mathbf{U}}_h^{ext}\right],\hfill \\ \hfill {\mathbf{H}}^{alt}({\Theta }_h^{int},{\Theta }_h^{ext})& =\frac{1}{2}\left[\mathbf{F}\left({\Theta }_h^{int}\right)+\mathbf{F}\left({\Theta }_h^{ext}\right)\right]-C_{11}\left[{\mathbf{U}}_h^{int}-{\mathbf{U}}_h^{ext}\right]+C_{12}\left[{\Theta }_h^{int}-{\Theta }_h^{ext}\right],\hfill \end{array}$$

(11)

where the superscripts “int” and “ext” convey the interior and exterior states of the elemental interface. When $C_{11}=0$ and $C_{12}=1/2$ or $-1/2$, the fluxes are called the alternating fluxes [40]. It can be noted that the mentioned flux scheme becomes central flux with $C_{11}=C_{12}=0$. Thus, the aforementioned weak formulation (10) can be expressed as follows:

$$\begin{array}{cc}\hfill \underset{E_k}{\int }{\Theta }_h{\phi }_{h}dV+\mathbf{D}\underset{E_k}{\int }\nabla {\phi }_h·{\mathbf{U}}_{h}dV-\mathbf{D}\underset{\partial E_k}{\int }{\phi }_h{\mathbf{H}}^{aux}d\Gamma & =0,\hfill \\ \hfill \frac{\partial }{\partial t}\underset{E_k}{\int }{\mathbf{U}}_h{\phi }_{h}dV+\underset{E_k}{\int }\nabla {\phi }_h·\mathbf{F}\left({\Theta }_h\right)dV-\underset{\partial E_k}{\int }{\phi }_h{\mathbf{H}}^{alt}d\Gamma & =\underset{E_k}{\int }{\phi }_h\mathbf{S}\left({\mathbf{U}}_h\right)dV.\hfill \end{array}$$

(12)

All volume and surface integrals appearing in the DG weak formulation (12) are calculated using the Gauss–Legendre quadrature rule, with a proper number of integration points consistent with the accuracy required. According to the analysis, if $P^l$ methods are utilized, the quadrature rules for the boundaries and interior (volume) of the elements must be precise for polynomials of degree $2l+1$ and $2l$, respectively [41]. The Gauss–Legendre quadrature of order $N_q$ for the reference rectangular element in the standard interval ${[-1,1]}^2$ is expressed as

$$\begin{array}{c}\hfill \underset{E_k}{\int }f(x,y)dxdy={\mathbf{J}}^v·\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}f(\xi ,\eta )d\xi d\eta ={\mathbf{J}}^v·\sum _{i=1}^{N_q}\sum _{j=1}^{N_q}{\omega }_i{\omega }_{j}f({\xi }_i^{★},{\eta }_j^{★}).\end{array}$$

(13)

In this expression, ${\mathbf{J}}^v$ denotes the Jacobian of the transformation for the reference element, respectively; ${\xi }_i^{★},{\eta }_j^{★}$ are the evaluation points on the reference element; ${\omega }_i,{\omega }_j$ are the weights related to the evaluation points; and $N_q$ represents the number of evaluation points for the reference element.

Finally, by assembling all the elemental contributions together, the modal DG formulation for the NCRD systems leads to a semi-discrete system of ODEs in time,

$$\mathbf{M}\frac{d{\mathbf{U}}_h}{dt}=\mathbf{L}\left({\mathbf{U}}_h\right),$$

(14)

where $\mathbf{M}$ and $\mathbf{L}\left(\mathbf{U}\right)$ are the block diagonal mass matrix and the residual vector, respectively. Due to the diagonal nature of the mass matrix, its inverse may be evaluated by hand for each element. For solving this ODE system, we employ an explicit third-order TVD Runge–Kutta scheme [41].

$$\begin{array}{cc}\hfill {\mathbf{U}}_h^{\left(1\right)}& ={\mathbf{U}}_h^n+\Delta t{\mathbf{M}}^{-1}\mathbf{L}\left({\mathbf{U}}_h^n\right),\hfill \\ \hfill {\mathbf{U}}_h^{\left(2\right)}& =\frac{3}{4}{\mathbf{U}}_h^n+\frac{1}{4}{\mathbf{U}}_h^{\left(1\right)}+\frac{1}{4}\Delta t{\mathbf{M}}^{-1}\mathbf{L}\left({\mathbf{U}}_h^{\left(1\right)}\right),\hfill \\ \hfill {\mathbf{U}}_h^{n+1}& =\frac{1}{3}{\mathbf{U}}_h^n+\frac{2}{3}{\mathbf{U}}_h^{\left(2\right)}+\frac{2}{3}\Delta t{\mathbf{M}}^{-1}\mathbf{L}\left({\mathbf{U}}_h^{\left(2\right)}\right),\hfill \end{array}$$

(15)

where $\mathbf{L}\left({\mathbf{U}}_h^n\right)$ is the residual of approximate solutions at time $t_n$ and $\Delta t$ is the suitable time-stepping length.

2.2. Reaction Term Treatment

The inclusion of reaction terms in strictly hyperbolic systems is known to result in prolonged relaxation periods, which presents significant numerical challenges [60]. Furthermore, incorrect solutions can arise when the resolution of reaction terms is inadequate. Thus, it is essential to address this matter appropriately. This subsection elucidates how the distinctive characteristics of the modal DG scheme eliminate the necessity for inefficient handling of reaction terms in conventional techniques. Let us consider the reaction terms associated with the elemental formulation (12) for a single variable, $u_h$.

$$\frac{d}{dt}\underset{E_k}{\int }{u}_h{\phi }_{h}dV+\underset{E_k}{\int }\nabla {\phi }_h·F\left(u_h\right)dV-\underset{\partial E_k}{\int }{\phi }_{h}F\left(u_h\right)·nd\Gamma =\underset{E_k}{\int }{\phi }_{h}S\left(u_h\right)dV.$$

(16)

The above-mentioned equation can also be represented in a matrix format by fascinating $\mathbf{U}$ as the global vector of modal coefficients:

$$\mathbf{M}\frac{d\mathbf{U}}{dt}+\mathbf{KU}-H^{alt}\left(u_h\right)\Phi -S\left(u_h\right){\Phi }^{{}^{\prime }}=0.$$

(17)

After some calculations, this equation may be read as

$$\begin{array}{cc}\hfill \frac{d\mathbf{U}}{dt}& =L\left(\mathbf{U}\right),\hfill \\ \hfill \mathbf{U}& ={(U^{\left(1\right)},U^{\left(2\right)},U^{\left(3\right)},\dots U^{\left(N\right)})}^T,\hfill \\ \hfill L\left(\mathbf{U}\right)& ={\mathbf{M}}^{-1}\mathbf{KU}+H^{alt}\left(u_h\right){\mathbf{M}}^{-1}\Phi -S\left(u_h\right){\mathbf{M}}^{-1}{\Phi }^{{}^{\prime }},\hfill \end{array}$$

(18)

where $\mathbf{M}$ is the mass matrix, $\mathbf{K}$ is the stiffness matrix, $\Phi$ is the vector associated with the boundary term, and ${\Phi }^{{}^{\prime }}$ is the vector associated with the source terms. Following are the definitions for these matrices:

$$\begin{array}{cc}\hfill \mathbf{M}& =\underset{E_k}{\int }{\phi }_i{\phi }_{j}dV\forall 0\le i,j\le N=N_l-1,\hfill \\ \hfill & =\left[\begin{array}{cccc}\underset{E_k}{\int }{\phi }_0{\phi }_{0}dV& \underset{E_k}{\int }{\phi }_0{\phi }_{1}dV& \cdots & \underset{E_k}{\int }{\phi }_0{\phi }_{N}dV\\ \underset{E_k}{\int }{\phi }_1{\phi }_{0}dV& \underset{E_k}{\int }{\phi }_1{\phi }_{1}dV& \cdots & \underset{E_k}{\int }{\phi }_1{\phi }_{N}dV\\ ⋮& ⋮& \ddots & ⋮\\ \underset{E_k}{\int }{\phi }_N{\phi }_{0}dV& \underset{E_k}{\int }{\phi }_N{\phi }_{1}dV& \cdots & \underset{E_k}{\int }{\phi }_N{\phi }_{N}dV\end{array}\right]\hfill \end{array}$$

(19)

Due to the orthogonal nature of the modal basic functions,

$$\mathbf{M}=\left\{\begin{array}{c}{C}_{ij}i=j\hfill \\ 0i\ne j\hfill \end{array}\right.,$$

$$\begin{array}{cc}\hfill \mathbf{K}& =\underset{E_k}{\int }\nabla {\phi }_i{\phi }_{j}dV\forall 0\le i,j\le N=N_l-1,\hfill \\ \hfill & =\left[\begin{array}{cccc}\underset{E_k}{\int }\nabla {\phi }_0{\phi }_{0}dV& \underset{E_k}{\int }\nabla {\phi }_0{\phi }_{1}dV& \cdots & \underset{E_k}{\int }\nabla {\phi }_0{\phi }_{N}dV\\ \underset{E_k}{\int }\nabla {\phi }_1{\phi }_{0}dV& \underset{E_k}{\int }\nabla {\phi }_1{\phi }_{1}dV& \cdots & \underset{E_k}{\int }\nabla {\phi }_1{\phi }_{N}dV\\ ⋮& ⋮& \ddots & ⋮\\ \underset{E_k}{\int }\nabla {\phi }_N{\phi }_{0}dV& \underset{E_k}{\int }\nabla {\phi }_N{\phi }_{1}dV& \cdots & \underset{E_k}{\int }\nabla {\phi }_N{\phi }_{N}dV\end{array}\right]\hfill \end{array}$$

(20)

$$\Phi =\left[\begin{array}{c}\underset{\partial E_k}{\int }{\phi }_0{\mathbf{J}}^{s}d\Gamma \\ \underset{\partial E_k}{\int }{\phi }_1{\mathbf{J}}^{s}d\Gamma \\ ⋮\\ \underset{\partial E_k}{\int }{\phi }_N{\mathbf{J}}^{s}d\Gamma \end{array}\right],{\Phi }^{{}^{\prime }}=\left[\begin{array}{c}\underset{E_k}{\int }{\phi }_0{\mathbf{J}}^{v}dV\\ \underset{E_k}{\int }{\phi }_1{\mathbf{J}}^{v}dV\\ ⋮\\ \underset{E_k}{\int }{\phi }_N{\mathbf{J}}^{v}dV\end{array}\right].$$

(21)

where ${\mathbf{J}}^s$ and ${\mathbf{J}}^v$ denote the Jacobian of the transformation for the reference line and rectangular elements, respectively. In the modal DG technique, the preference for orthogonal polynomial functions considerably facilitates the beneficence of the high-order moments of the approximate solutions to the reaction terms associated with vector ${\Phi }^{{}^{\prime }}$. When employing the modal polynomial functions $\left({\phi }_i\right)$, numerous integral terms appearing in ${\Phi }^{{}^{\prime }}$ become nullified within the interval $[-1,1]$. For instance, in the case of a third-order DG approximation, the computation of ${\Phi }^{{}^{\prime }}$ can be expressed as follows:

$${\Theta }^{{}^{\prime }}=\left[\begin{array}{c}\underset{E_k}{\int }{\phi }_0{\mathbf{J}}^{v}dV\\ \underset{E_k}{\int }{\phi }_1{\mathbf{J}}^{v}dV\\ \underset{E_k}{\int }{\phi }_2{\mathbf{J}}^{v}dV\\ \underset{E_k}{\int }{\phi }_3{\mathbf{J}}^{v}dV\\ \underset{E_k}{\int }{\phi }_4{\mathbf{J}}^{v}dV\\ \underset{E_k}{\int }{\phi }_5{\mathbf{J}}^{v}dV\end{array}\right]=\left[\begin{array}{c}4{\mathbf{J}}^v\\ 0\\ 0\\ -2{\mathbf{J}}^v/3\\ 0\\ -2{\mathbf{J}}^v/3\end{array}\right].$$

(22)

The reaction term handling is substantially simplified in this unique way, exactly as it is in the first-order ($P^0$) case. In short, when the left sides of Equations (12) and (16) are derived using high-order polynomial approximations, this leads to the participation of the cell average solutions dominating the reaction terms in the modal DG construction.

3. Validation Study

3.1. Accuracy Test for Solver

To verify the accuracy of the employed numerical solver, we examine a 1D linear reaction–diffusion system with the domain $\Omega =[0,\pi /2]$, as described below:

$$\left\{\begin{array}{cc}\frac{\partial u}{\partial t}={\mu }_1\frac{{\partial }^{2}u}{\partial x^2}-k_{1}u+v,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ \frac{\partial v}{\partial t}={\mu }_2\frac{{\partial }^{2}v}{\partial x^2}-k_{2}v,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ u_x(0,t)=u(\pi /2,t)=0,\hfill & t\in \times [0,T],\hfill \\ v_x(0,t)=v(\pi /2,t)=0,\hfill & t\in \times [0,T],\hfill \end{array}\right.$$

(23)

for some real parameters, $k_1$ and $k_2$. For ${\mu }_1={\mu }_2$, the exact solution of the above system is given by Chou et al. [61]:

$$\left\{\begin{array}{cc}u(x,t)=\left[e^{-(k_1+{\mu }_1)t}+e^{-(k_2+{\mu }_1)t}\right]cosx,\hfill & \\ v(x,t)=(k_1-k_2)e^{-(k_1+{\mu }_1)t}cosx.\hfill & \end{array}\right.$$

(24)

The initial conditions are obtained from the exact solutions of the system (23) by taking $t=0$. We compute the numerical solutions for this RD system with the parameters ${\mu }_1=1.0$, $k_1=0.1$, and $k_2=0.01$, which shows a diffusion-dominated case. $\Delta t=0.001$ is considered for the time discretization, and simulations are run up to $t=1$. All the computations are performed with $N=100$ and the DG scheme with third-order accuracy. In

Table 1. Comparison of $L_{\infty }$-error for u with ${\mu }_1=1.0$, $k_1=0.1$, $k_2=0.01$ in a 1D linear RD system.

Time	CN-MG [61]	IIF2 [61]	ECSCM [62]	DQM [63]	Present
0.04	$1.09\times {10}^{-4}$	$2.73\times {10}^{-5}$	1.10E-04	$2.02\times {10}^{-5}$	$1.85\times {10}^{-5}$
0.02	$2.67\times {10}^{-5}$	$6.40\times {10}^{-6}$	$2.84\times {10}^{-5}$	$2.07\times {10}^{-5}$	$2.76\times {10}^{-6}$
0.01	$6.27\times {10}^{-6}$	$1.19\times {10}^{-6}$	$7.91\times {10}^{-6}$	$2.09\times {10}^{-6}$	$1.01\times {10}^{-6}$
0.005	$1.16\times {10}^{-6}$	$1.07\times {10}^{-7}$	$2.79\times {10}^{-6}$	$2.10\times {10}^{-6}$	$1.03\times {10}^{-6}$

, $L_{\infty }$-norms are estimated for the u variable, and the computed results are compared against the existing numerical results [61,62,63]. In addition, a profile comparison between the exact and numerical solutions for u and v at $t=1$ is illustrated in Figure 2. The present numerical results are found to be quite similar to the exact solutions, which demonstrates that the proposed numerical scheme is efficient in capturing the spatiotemporal wave patterns accurately. Furthermore, the order of accuracy of the proposed numerical scheme is also estimated. For this purpose, simulations are run till $t=1.0$ and up to the third-degree of polynomials (i.e., $1\le l\le 3$) by varying the number of grid points, $N_x(=20,40,80,160)$. In

Table 2. Error estimation with ${\mu }_1=0.001$, $k_1=0.1$, $k_2=0.01$ at $t=1$ in a 1D linear RD system.

l	${\mathit{N}}_{\mathit{x}}$	${\mathit{L}}_2$-Norm	Order	${\mathit{L}}_{\mathit{\infty }}$-Norm	Order
1	20	$6.751\times {10}^{-2}$	−	$4.213\times {10}^{-2}$	−
	40	$1.804\times {10}^{-2}$	$1.91$	$1.105\times {10}^{-2}$	$1.85$
	80	$3.613\times {10}^{-3}$	$1.87$	$2.410\times {10}^{-3}$	$1.96$
	160	$2.138\times {10}^{-3}$	$1.79$	$1.036\times {10}^{-3}$	$2.01$
2	20	$1.041\times {10}^{-3}$	−	$1.248\times {10}^{-3}$	−
	40	$1.252\times {10}^{-4}$	$3.01$	$1.592\times {10}^{-4}$	$2.65$
	80	$1.732\times {10}^{-5}$	$3.10$	$2.007\times {10}^{-5}$	$2.75$
	160	$2.100\times {10}^{-6}$	$3.02$	$2.510\times {10}^{-6}$	$2.98$
3	20	$1.767\times {10}^{-5}$	−	$2.701\times {10}^{-5}$	−
	40	$1.101\times {10}^{-6}$	$4.02$	$1.701\times {10}^{-6}$	$3.69$
	80	$7.015\times {10}^{-8}$	$4.00$	$3.126\times {10}^{-7}$	$3.71$
	160	$5.252\times {10}^{-9}$	$4.01$	$5.125\times {10}^{-8}$	$3.92$

, the convergence analysis and the order of accuracy have been studied. This investigation shows that the current numerical approach reached the necessary order of accuracy $(l+1)$.

3.2. Validation for Spatiotemporal Pattern Formation

In our earlier study [55], the present numerical method was validated for the generalized 1D FitzHugh–Nagumo RD model. In the current study, to validate the spatiotemporal patterns arising in NCRD systems, the present numerical results are compared with the computational results of Zegeling and Kok [23] for the 1D Gray–Scott NCRD system. For this case, the following parameters are considered: ${\mu }_1={10}^{-4},{\mu }_2={10}^{-6},k_1=0.024,k_2=0.026$. The initial conditions are $u(x,0)=1-\frac{1}{2}{\sin }^{100}\left(\pi x\right),v(x,0)=\frac{1}{4}{\sin }^{100}\left(\pi x\right)$, along with Dirichlet boundary conditions on the domain $[0,1]$. The simulations are run with 300 grid points on the computational domain $[0,1]$, and $\Delta t=0.01$ is used as the time step. Figure 3 illustrates the comparison of solutions’ trajectories for the u variable between the computed results of Zegeling and Kok [23] and the present numerical results up to $t=2000$. It can be observed that this system forms a pulse-splitting wave pattern. The wave patterns propagate, split, and overlap endlessly until they encompass the whole computational domain. During the pulse-splitting process, a few additional pulses are produced in the back portion of the dominant pulses. To summarize, the figure indicates a satisfactory level of agreement between the present and previous computational results.

Furthermore, for the other validation case, we consider a 2D Gray–Scott NCRD system with the following parameters: ${\mu }_1=8\times {10}^{-5},{\mu }_2=4\times {10}^{-5},k_1=0.024,k_2=0.026$. The initial conditions are defined on the domain $[0,1]\times [0,1]$ as

$$\begin{array}{cc}\hfill u(x,y,0),v(x,y,0)& =\left\{\begin{array}{cc}0.5,0.25\hfill & \mathrm{if}0.3\le x,y\le 0.7,\hfill \\ 1,0\hfill & \mathrm{elsewhere},\hfill \end{array}\right.\hfill \end{array}$$

(25)

along with Dirichlet boundary conditions. The simulations are run with $61\times 61$ grid points on the computational domain $[0,1]\times [0,1]$, and $\Delta t=0.01$ is used as the time step. The comparison of the u variable solution between the current results and the results of Zegeling and Kok [23] for $t=50$, 100, and 150 is shown in Figure 4. These numerical simulations share the same initial conditions and spatiotemporal pattern formations. The present spatiotemporal patterns, including four spots, are in good agreement with the computational results of Zegeling and Kok [23]. Consequently, the current numerical solver can accurately capture the formation of spatiotemporal patterns.

4. Results of Spatiotemporal Pattern Formations in NCRD Systems and Discussion

This section illustrates an extensive investigation of the spatiotemporal pattern formations in various NCRD systems, including Gray–Scott, Brusselator, Schnakenberg, and isothermal models. Computational simulations on the spatiotemporal pattern formation are conducted on an Intel Core i7-12700 [Manufactured by Arizona (USA), Equipment sourced in Aachen (Germany) ] workstation using a single processor. These simulations are run with Visual Studio 2019 in the Fortran 95 forum. Additionally, Tecplot 360 2019R1 software is used to visualize the numerical results.

4.1. 1D Gray–Scott NCRD System

The Gray–Scott NCRD system, originally proposed by Gray and Scott [19], describes the various patterns and phenomena existing in nature, including Turing, butterfly wings, multiple spots, and embryo. In this model, the two kinematic reactions are studied as [20]

$$\begin{array}{c}\hfill U+2V\to 3V,\\ \hfill V\to P,\end{array}$$

where U and V are the chemical species and P denotes an inert outcome of these species. Both kinematic reactions are irreversible; therefore, the densities of U and V remain consistent. Both chemical species are being removed in the process and the resulting dimensionless Gray–Scott NCRD system with domain $\Omega =[-50,50]$ becomes

$$\left\{\begin{array}{cc}\frac{\partial u}{\partial t}={\mu }_1\frac{{\partial }^{2}u}{\partial x^2}-u{v}^2+k_1(1-u),\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ \frac{\partial v}{\partial t}={\mu }_2\frac{{\partial }^{2}v}{\partial x^2}+u{v}^2-(k_1+k_2)v,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ u(-50,t)=u(50,t)=1,\hfill & t\in \times [0,T],\hfill \\ v(-50,t)=v(50,t)=0,\hfill & t\in \times [0,T].\hfill \end{array}\right.$$

(26)

For this system, we consider the following initial conditions:

$$\begin{array}{cc}\hfill u(x,0)& =1-\frac{1}{2}{\mathit{sin}}^{100}\left(\pi \frac{(x-50)}{100}\right),\hfill \\ \hfill v(x,0)& =\frac{1}{4}{\mathit{sin}}^{100}\left(\pi \frac{(x-50)}{100}\right).\hfill \end{array}$$

(27)

In the aforementioned expression, u and v denote the concentrations of chemical species U and V, respectively; ${\mu }_1$, ${\mu }_2$ denote the diffusion rates; $k_1$ represents the feed rate of u into the system; and $k_2$ represents the reaction rate of the second reaction. For this system, we evaluate the solutions with the following selected parameters: ${\mu }_1=1.0$, ${\mu }_2=0.01$, $k_1=0.064$, and $k_2=0.062$. All the numerical experiments are carried out with 400 elements on the spatial domain $\Omega =[-50,50]$ with $\Delta t=0.01$, and the simulations are run until $t=2500$ so that the stationary wave patterns are obtained. The growth, separation, and replication of the solution patterns are shown in Figure 5 and Figure 6. The wave patterns are continuously growing, splitting, and replicating themselves until the entire computational domain is covered by these waves. During the pulse-splitting process, some new pulses are developed in the rear part of the prevailing pulses, as demonstrated in Figure 6.

4.2. 1D Brusselator NCRD System

The Brusselator NCRD system is crucial in real-world physical issues such as Turing pattern creation on animal skin, ozone synthesis from atomic oxygen, enzyme reactions, and laser and plasma physics due to multiple mode coupling. This system was originated by Prigogine and Lefever [16]. It uses highly nonlinear oscillations to illustrate the chemical RD process by the following tri-molecular chemical system [64]:

$$\begin{array}{cc}\hfill k_1& \to U,\hfill \\ \hfill k_2+U& \to V+D,\hfill \\ \hfill 2U+V& \to 3U,\hfill \\ \hfill U& \to E,\hfill \end{array}$$

where $k_1$ and $k_2$ are the input chemicals, D and E are the output chemicals, and U and V represent the intermediates. The second reaction represents the bi-molecular reaction, while, the third one denotes the autocatalytic tri-molecular reaction. The resulting dimensionless 1D Brusselator NCRD system with domain $\Omega =[0,1]$ is defined as [23]

$$\left\{\begin{array}{cc}\frac{\partial u}{\partial t}={\mu }_1\frac{{\partial }^{2}u}{\partial x^2}+u^{2}v-(1+k_2)u+k_1,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ \frac{\partial v}{\partial t}={\mu }_2\frac{{\partial }^{2}v}{\partial x^2}-u^{2}v+k_{2}u,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ u_x(0,t)=u_x(1,t)=0,\hfill & t\in \times [0,T],\hfill \\ v_x(0,t)=v_x(1,t)=0,\hfill & t\in \times [0,T],\hfill \end{array}\right.$$

(28)

with the following initial conditions

$$\left\{\begin{array}{c}u(x,0)=0.5,\hfill \\ v(x,0)=1+5x.\hfill \end{array}\right.$$

(29)

For the numerical experiments, we chose the following parameters, ${\mu }_1={\mu }_2={10}^{-4}$, $k_1=1.0$, $k_2=3.4$, based on previous studies [23]. The computational work has been performed with 200 elements on the spatial domain $\Omega =[0,1]$; the used time-step is $\Delta t=0.01$ and simulations are run until $t=20$. Figure 7 and Figure 8 show the spatiotemporal pattern formation of the Brusselator NCRD system. These patterns of the Brusselator system were found to be very similar to previous studies [23,26,27]. In Figure 6, the physical behavior of the u and v solution profiles during the molecular interactions are observed at different time instants. At $t=0$, these profiles were found to be high on the right side of the domain. As time increases, these spatiotemporal profiles move toward the left side of the domain.

4.3. 1D Isothermal Chemical NCRD System

The isothermal chemical NCRD system was developed with auto-catalytic reactions [24]:

$$\begin{array}{cc}\hfill U+2V& \to 3V,\hfill \\ \hfill V+V& \to P,\hfill \end{array}$$

where U and V illustrate the chemical reactants and P denotes the output reactant. The resulting dimensionless 1D isothermal chemical NCRD system with $\Omega =[0,200]$ is defined as [24]

$$\left\{\begin{array}{cc}\frac{\partial u}{\partial t}={\mu }_1\frac{{\partial }^{2}u}{\partial x^2}-uv,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ \frac{\partial v}{\partial t}={\mu }_2\frac{{\partial }^{2}v}{\partial x^2}+uv-k_{1}v,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ u_x(0,t)=u_x(200,t)=0,\hfill & t\in \times [0,T],\hfill \\ v_x(0,t)=v_x(200,t)=0,\hfill & t\in \times [0,T].\hfill \end{array}\right.$$

(30)

For this system, we adopt the following initial conditions:

$$\left\{\begin{array}{c}u(x,0)=1,\hfill \\ v(x,0)=e^{-x^2}.\hfill \end{array}\right.$$

(31)

In an isothermal chemical NCRD system, $k_1>0$ represents a dimensionless variable that expresses the intensity of the crumbling step to the autocatalytic fashioning phase. It was demonstrated using this system that almost-isothermal flames emerging from quadratic branching in the carbon sulfide oxygen reaction could be characterized using basic cubic autocatalysis [18]. Following previous studies [24,26], the numerical experiments are carried out with the parameters ${\mu }_1={\mu }_2=1$ and $k_1=0.1$, $0.5$, $0.9$. For this NCRD system, we performed the numerical simulations with 300 elements over the region $\Omega =[0,200]$ with $\Delta t=0.01$, and the simulations are run till time $t=300$. The impact of different $k_1$ values on the spatiotemporal pattern formations of u and v are illustrated in Figure 9 and Figure 10. For all $k_1$ values, the spatiotemporal profiles for the u and v solutions transport in the direction of the right side of the domain, with time increasing as shown in Figure 9, and the concentration of u remains constant towards the boundary, whereas v decreases to zero. Further, the physical behavior of these profiles can also be revealed by the spatiotemporal contours, as shown in Figure 10.

4.4. 1D Schnakenberg NCRD System

The Schnakenberg NCRD system [17] is a fundamental model that describes an autocatalytic chemical reaction with possible oscillatory behavior. The tri-molecular reactions between two chemical products, $U,V$, and two chemical source, $A,B$, are formulated as

$$\begin{array}{cc}\hfill A& \rightleftharpoons U,\hfill \\ \hfill B& \to V,\hfill \\ \hfill 2U+V& \to 3U.\hfill \end{array}$$

A system of two nonlinear RD equations for the concentrations u and v of the chemical products U and V can be obtained through the law of mass action. A 1D Schnakenberg NCRD system in dimensionless form with domain $\Omega =[0,3\pi ]$ can defined as [1,17]

$$\left\{\begin{array}{cc}\frac{\partial u}{\partial t}={\mu }_1\frac{{\partial }^{2}u}{\partial x^2}+u^{2}v-u+k_1,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ \frac{\partial v}{\partial t}={\mu }_2\frac{{\partial }^{2}v}{\partial x^2}-u^{2}v+k_2,\hfill & (x,t)\in \Omega \times [0,T],\hfill \\ u_x(0,t)=u_x(3\pi ,t)=0,\hfill & t\in \times [0,T],\hfill \\ v_x(0,t)=v_x(3\pi ,t)=0,\hfill & t\in \times [0,T],\hfill \end{array}\right.$$

(32)

with two initial conditions:

$$\mathrm{I}.\mathrm{Cs}.\left(\mathrm{i}\right)=\left\{\begin{array}{c}u(x,0)=0.8+0.1cosx,\hfill \\ v(x,0)=1.03+0.1cosx.\hfill \end{array}\right.\mathrm{I}.\mathrm{Cs}.\left(\mathrm{ii}\right)=\left\{\begin{array}{c}u(x,0)=0.3+0.001sin3x,\hfill \\ v(x,0)=1.778+0.001cos2x.\hfill \end{array}\right.$$

(33)

Following the study of Jiwari et al. [27], we have selected the following parameters for computational simulations: (a) ${\mu }_1=0.2,{\mu }_2=0.1$, $k_1=0.14$, $k_2=0.66$, with initial conditions $\left(\mathrm{I}.\mathrm{Cs}.\left(\mathrm{i}\right)\right)$ and (b) ${\mu }_1=0.01,{\mu }_2=1.0$, $k_1=0.14$, $k_2=0.16$, with initial conditions $\left(\mathrm{I}.\mathrm{Cs}.\left(\mathrm{ii}\right)\right)$. For this system, we performed the numerical simulations with 300 elements over the region $\Omega =[0,3\pi ]$ with a time-step of $\Delta t=0.001$, and DG solver is run to time $t=300$. The spatiotemporal patterns with two different initial conditions and different parameters are illustrated in Figure 11 and Figure 12. The numerical results show that the solution converges to a spatially homogeneous periodic orbit with the first initial condition (33), whereas with the second initial condition, the solutions display a spatially periodic spatiotemporal pattern with wavelength $\pi$. Further, the physical behavior of the spatiotemporal patterns can also be revealed by the space–time plots, as shown in Figure 12. The spatiotemporal patterns obtained by the present DG scheme are quite similar to the results of Jiwari et al. [27]. It can also be observed that the spatiotemporal patterns of the system strictly depend upon the initial conditions and parameters, as shown in Figure 11 and Figure 12, which are drastically different.

4.5. 2D Gray–Scott NCRD System

For capturing the spatiotemporal pattern formation, the 2D NCRD system is more obvious than the 1D system. Therefore, we consider a 2D Gray–Scott NCRD system in this test case, as follows [23]:

$$\left\{\begin{array}{cc}\frac{\partial u}{\partial t}={\mu }_1\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]-u{v}^2+k_1(1-u),\hfill & (x,y,t)\in \Omega \times [0,T],\hfill \\ \frac{\partial v}{\partial t}={\mu }_2\left[\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right]+u{v}^2-(k_1+k_2)v,\hfill & (x,y,t)\in \Omega \times [0,T].\hfill \end{array}\right.$$

(34)

For this system, we consider two block functions as the initial condition defined in the domain $\Omega =[0,1]\times [0,1]$ as

$$\begin{array}{cc}\hfill u(x,y,0)& =\left\{\begin{array}{cc}0.5\hfill & \mathrm{if}0.3\le x\le 0.7\mathrm{and}0.3\le y\le 0.7,\hfill \\ 1\hfill & \mathrm{otherwise},\hfill \end{array}\right.\hfill \\ \hfill v(x,y,0)& =\left\{\begin{array}{cc}0.25\hfill & \mathrm{if}0.3\le x\le 0.7\mathrm{and}0.3\le y\le 0.7,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

(35)

For this NCRD system, the two-dimensional domain is partitioned with the grid points $151\times 151$, and the Dirichlet conditions are used for the boundary implementation. For simulations, the time step of $\Delta t={10}^{-3}$ is used. The following parameters are considered for computations: ${\mu }_1=8\times {10}^{-5}$, ${\mu }_2=4\times {10}^{-5}$, $k_1=0.024$, $k_2=0.06$. Figure 13 illustrates the spatiotemporal contours for u and v at distinct time intervals $t=0,50,100,150,300,400,500,$$600,800,$ and 1000. The two-dimensional system splits the original block functions into four spots and then eight spots in terms of the first and second concentration components. From the results, it is found that the spatiotemporal patterns obtained by the present DG scheme are quite similar to the results of Zegeling and Kok [23].

4.6. 2D Brusselator NCRD System

Finally, a 2D Brusselator NCRD system is examined, as follows [23]:

$$\left\{\begin{array}{cc}\frac{\partial u}{\partial t}={\mu }_1\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]+u{v}^2-(1+k_1)u+k_2,\hfill & (x,y,t)\in \Omega \times [0,T],\hfill \\ \frac{\partial v}{\partial t}={\mu }_2\left[\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right]-u{v}^2+k_{1}u,\hfill & (x,y,t)\in \Omega \times [0,T].\hfill \end{array}\right.$$

(36)

For this NCRD system, the following initial conditions are considered, along with Neumann boundary conditions in the computational $\Omega =[0,1]\times [0,1]$:

$$\begin{array}{cc}\hfill u(x,y,0)& =0.5+y,\hfill \\ \hfill v(x,y,0)& =1+5x.\hfill \end{array}$$

(37)

Numerical simulations are performed with the following parameters: ${\mu }_1={\mu }_2=2.0\times {10}^{-3}$, $k_1=1$, and $k_2=3.4$. Figure 14 illustrates the spatiotemporal contours for the concentrations u and v at distinct time intervals. From the results, it is obvious that due to $k_2>1+k_1$, the solutions do not converge, and they oscillate on a regular basis. The spatiotemporal patterns travel from the y-axis to the x-axis, as shown in the contour plots. These plots, like those in Zegeling and Kok [23], provide similar behaviors to the 2D Brusselator NCRD system.

5. Concluding Remarks

In this study, we investigate the spatiotemporal pattern formations in nonlinear coupled reaction–diffusion (NCRD) systems using a mixed-type modal discontinuous Galerkin method. This method illustrates a new formulation to handle the second-order derivatives arising in the diffusion terms. Towards the spatial discretization process, hierarchical modal basis functions premised on orthogonal scaled Legendre polynomials are employed. Furthermore, a novel approach for handling the reaction term in the NCRD systems is also proposed. The employed method simplifies the NCRD systems to a set of time-dependent ordinary differential equations, which are effectively solved by an explicit third-order TVD Runge–Kutta scheme. We verify the reliability and accuracy of the used numerical solver by comparing the current results with earlier computational results, which show a good degree of agreement. In order to capture the spatiotemporal patterns, several 1D/2D NCRD systems, including Gray–Scott, Brusselator, isothermal, and Schnakenberg models, are investigated. The numerical simulations reveal that the spatiotemporal pattern formations are highly reliant on the initial condition and parameters of the model. For example, three intriguing features—growth, separation, and replication—are displayed in spatiotemporal patterns with different parameters in the 1D Gray–Scott model. Meanwhile, a spatially periodic spatiotemporal pattern with wavelength $\pi$ is observed in the 1D Schnakenberg model, with specific initial conditions. Interestingly, in the 2D Gray–Scott model, the spatiotemporal patterns split the original block functions into four spots and then eight spots in terms of the first and second concentration components.

The primary objective of the current study is to investigate the generation of spatiotemporal patterns in 1D/2D NCRD systems. Additional forces, such as electric and magnetic fields, play a prominent role in the formation of spatiotemporal patterns in multi-dimensional NCRD systems. The cross-diffusion rate is also expected to have a substantial impact on the development of the spatiotemporal pattern in NCRD systems. In this regard, future work will examine the effects of electromagnetic fields, particularly the cross-diffusion rate, on spatiotemporal pattern formation in multi-dimensional varied NCRD systems.