Spatiotemporal Dynamics of a Reaction Diffusive Predator-Prey Model: A Weak Nonlinear Analysis

Abstract

In the realm of ecology, species naturally strive to enhance their own survival odds. This study introduces and investigates a predator-prey model incorporating reaction-diffusion through a system of differential equations. We scrutinize how diffusion impacts the model’s stability. By analysing the stability of the model’s uniform equilibrium state, we identify a condition leading to Turing instability. The study delves into how diffusion influences pattern formation within a predator-prey system. Our findings reveal that various spatiotemporal patterns, such as patches, spots, and even chaos, emerge based on species diffusion rates. We derive the amplitude equation by employing the weak nonlinear multiple scales analysis technique and the Taylor series expansion. A novel sinc interpolation approach is introduced. Numerical simulations elucidate the interplay between diffusion and Turing parameters. In a two-dimensional domain, spatial pattern analysis illustrates population density dynamics resulting in isolated groups, spots, stripes, or labyrinthine patterns. Simulation results underscore the method’s effectiveness. The article concludes by discussing the biological implications of these outcomes.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Mathematical modeling is an effective tool for describing and comprehending the complex systems’ spatiotemporal dynamics. The behaviours of various species in the natural ecosystem reveal a wealth of dynamical traits. The habits of different species in the natural ecosystem exhibit a wide variety of dynamic characteristics [1]. Many species have threatened extinction in recent decades as a result of inadequate resources, overexploitation, pollution, and predation. Most species extinctions are caused by an imbalance in the environment or an ecological system. External support trends to prevent species extinction including refuge and population limitation to a single place [2–4]. Considering the issue of extinction, many authors [5–10] investigated the dynamical behaviour of the same system. However, because we live in a spatial environment, spatial aspects in predator-prey systems should be considered. As a result, to characterize these systems, reaction-diffusion equations [11] should be used. The interaction between prey and predators receives the most attention in the ecosystem because it is so important and happens all over the world [4, 12]. In the realm of ecology, the dynamics of predator-prey interactions have long been the focus of interest and research. For forecasting and controlling ecosystems as well as deciphering the complex interspecies balance, it is essential to comprehend the spatiotemporal dynamics of these interactions. Mathematical models may be used as a useful tool for analyzing and forecasting the behaviour of predator and prey populations, which is one way to investigate these dynamics [13].

The spatiotemporal dynamics of a response diffusive predator-prey model is an important study area for a number of reasons [14–16]. First of all, it enables us to investigate the intricate interplay of species in a vast ecosystem. We may investigate how the movement and dispersal of predators and prey affect their population dynamics, geographical distribution, and overall stability by introducing diffusion factors. Furthermore, this study issue is extremely important in the light of global environmental changes and habitat fragmentation [17–19]. Landscape changes and the loss of habitat connectivity can have a significant impact on the dynamics of predator-prey interactions. We may acquire insights into the resilience and susceptibility of ecosystems under changing conditions and improve conservation policies aiming at maintaining biodiversity by examining the spatiotemporal dynamics of predator-prey models. Finally, the study of the spatiotemporal dynamics of a response diffusive predator-prey model is an extremely interesting research area. We can acquire a better understanding of how spatial dynamics and mobility impact the stability, coexistence, and persistence of predator and prey populations by including diffusion variables into mathematical models. This study’s findings have implications for biodiversity conservation, ecosystem management, and our general knowledge of complex spatiotemporal systems in nature and beyond [20, 21].

The creation of spatiotemporal patterns is significant in describing the dynamics of interacting populations as part of a larger ecological community [22, 23]. Such populations are scattered heterogeneously over their habitats, resulting in patterns that might be stationary or nonstationary in terms of time. Turing [24] demonstrated that diffusion-driven instability might bring out the spatial patterns in a system of coupled diffusion equations. Segel and Jackson [25] proved that Turing’s theories can also be applied to species diversity. Gierer and Meinhardt [26] also demonstrated a biologically justified formulation of the Turing-type model. Levin and Segel [27] proposed spatial pattern formation as one of the possible causes of planktonic flakiness.

Nevertheless, when compared to two-species food chain models, pattern development for three-species systems has received far less attention. In the study in [28], White and Gilligan give adequate conditions for a diffusive model to induce diffusion-driven instability. Chen and Peng [29] established the existence of nonconstant positive steady conditions to prove the endurance of static patterns for a cross-diffusion model. They also showed that without cross diffusion, the corresponding model fails to deliver any stationary pattern.

The time-dependent spatial patterns can be studied using the amplitude equation. It offers us a diplomatic classification of weak nonlinear studies’ pattern formation. It is feasible to derive the given model’s amplitude equation using multiple-scale analysis.

A wide range of computational approaches has been employed for the nonlinear prey-predator model, including the finite difference method [30], B-spline method [31], finite element method [32], spectral method [33], perturbation method, variational iteration method (VIM) [34], and so on. Researchers, on the other hand, have always sought to identify the most effective numerical method [35]. A sinc function interpolation method for the prey-predator system has been developed.

In this article, our objective is to explore the intricate spatiotemporal dynamics exhibited by models involving three interacting species (specifically, two prey species and their shared predator). Our primary focus lies in analyzing the consequences of the stochastic dispersal of all three species on their respective population distributions. In addition, we aim to delve into fundamental ecological inquiries, such as the phenomenon of extended transience and the influence of habitat size on the resulting patterns. To address these questions, we thoroughly examine a temporal model based on reaction-diffusion principles. These systematic investigations hold values in discerning and contrasting the ensuing dynamics and their underlying causal factors, given that these three models form a sequence of nested systems. The following is a breakdown of the structure of our work. In Section 2, we introduce a temporal model with three species—two prey species and one predator, and discuss the dynamical behaviours briefly as well as provide stability of the model then expand by including spatial components, where we also present a broad range of spatiotemporal dynamics that the model is capable of. The necessary condition of Turing instability as well as bifurcation is identified in Section 3. In Section 4, a sinc interpolation approach has been developed using our model. We obtain the amplitude equation by applying a weak nonlinear analysis close to the parameter threshold value for the patterns in Section 5. In Section 6, we perform numerical simulations to support our analytical results. Then, we show the Turing patterns for some fixed parameter values in Section 7. We discussed and concluded our work in Section 8.

2. Temporal Model and Stability

In this section, we framed a model of prey-predator interactions involving three species associated with a system of three differential equations and discussed the stability of this model.$$\begin{array}{c}\text{(1)}& \frac{\mathrm{d}{\mathrm{a}}_1}{\text{dt}}=g_1{a}_{1}1-\frac{a_1}{{\mathbb{K}}_1}-\sigma a_1{a}_2-\nu a_{1}b,\frac{\mathrm{d}{\mathrm{a}}_2}{\text{dt}}=g_2{a}_{2}1-\frac{a_2}{{\mathbb{K}}_2}-\phi a_1{a}_2-\eta a_{2}b,\\ \frac{\text{db}}{\text{dt}}={\mu }_1\nu a_{1}b+{\mu }_2\eta a_{2}b-\mathrm{\varrho }b.\end{array}$$

The parameters used in the model are described in Table 1.

Table 1

Description of parameters.

Parameters	Description
$a_{1}t$	Prey 1 population
$a_{2}t$	Prey 2 population
$bt$	Predator population
$g_1,g_2$	Intrinsic growth rate
${\mathbb{K}}_1,K_2$	Carrying capacity of preys
$\sigma ,\phi$	Competition rates between two preys
$\nu$, $\eta$	The rate of predation on preys by predator
${\mu }_1,{\mu }_2$	Conversion efficiencies of preys
$\mathrm{\varrho }$	Predator’s natural death rate

Stability in ecology is frequently described as the capacity to resume the initial structure or functioning following external disturbances. Now, we find the equilibrium points, confirm their existence, and examine their stability in order to better understand the dynamical behaviour of the system (1). Solving the following simultaneous equations yields the equilibrium points as follows:$$\begin{array}{c}\text{(2)}& g_1{a}_{1}1-\frac{a_1}{{\mathbb{K}}_1}-\sigma a_1{a}_2-\nu a_{1}b=0,g_2{a}_{2}1-\frac{a_2}{{\mathbb{K}}_2}-\phi a_1{a}_2-\eta a_{2}b=0,\\ {\mu }_1\nu a_{1}b+{\mu }_2\eta a_{2}b-\mathrm{\varrho }b=0.\end{array}$$

System (1) admits five equilibrium points by solving the above equations.$$\begin{array}{c}\text{(3)}& E_0=\mathrm{0,0,0},E_1={\mathbb{K}}_1,\mathrm{0,0},E_2=0,{\mathbb{K}}_2,0,E_3=\frac{g_1{\mathbb{K}}_1{g}_2-\sigma {\mathbb{K}}_1{\mathbb{K}}_2{g}_2}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2},\frac{g_1{\mathbb{K}}_2{g}_2-\phi {\mathbb{K}}_1{\mathbb{K}}_2{g}_1}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2},0.\end{array}$$

System (1) yields a unique positive equilibrium $E^{\ast }=a_1^{\ast },a_2^{\ast },b^{\ast }$, where$$\begin{array}{c}\text{(4)}& a_1^{\ast }=\frac{{\mathbb{K}}_1\nu g_2\mathrm{\varrho }+{\mathbb{K}}_2\eta g_1{\mu }_2\eta -\nu g_2{\mu }_2-\sigma \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},a_2^{\ast }=\frac{{\mathbb{K}}_2\eta g_1\mathrm{\varrho }+{\mathbb{K}}_1\nu g_2{\mu }_1\nu -\eta g_2{\mu }_1-\phi \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ b^{\ast }=\frac{g_1{g}_2\nu {\mu }_1{\mathbb{K}}_1+\eta {\mu }_2{\mathbb{K}}_2-\mathrm{\varrho }+{\mathbb{K}}_1{\mathbb{K}}_2\sigma \phi \mathrm{\varrho }-\sigma \nu g_2{\mu }_1-\phi \eta g_1{\mu }_2}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2}.\end{array}$$

The Jacobian matrix is as follows:$$\begin{array}{}\text{(5)}& J{a}_1,a_2,b=\begin{array}{ccc}{g}_1-\frac{2g_1{a}_1}{{\mathbb{K}}_1}-\sigma a_2-\nu b& -\sigma a_1& -\nu a_1\\ -\phi a_2& g_2-\frac{2g_2{a}_2}{{\mathbb{K}}_2}-\phi a_1-\eta b& -\eta a_2\\ {\mu }_1\nu b& {\mu }_2\eta b& {\mu }_1\nu a_1+{\mu }_2\eta a_2-\rho \end{array}.\end{array}$$

Theorem 1.

$E_0=\mathrm{0,0,0}$ of system (1) is always a saddle point.

Proof.

The Jacobian matrix of (1) about $E_0=\mathrm{0,0,0}$ is as follows:$$\begin{array}{}\text{(6)}& J{E}_0=\begin{array}{ccc}{g}_1& 0& 0\\ 0& g_2& 0\\ 0& 0& -\rho \end{array}.\end{array}$$

The eigenvalues of the matrix are ${\mathrm{\Lambda }}_1=g_1,{\mathrm{\Lambda }}_2=g_2,{\text{and\hspace{0.17em}}\mathrm{\Lambda }}_3=-\rho$. Here, ${\mathrm{\Lambda }}_1>0\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Lambda }}_3<0$. As a result, system (1)’s equilibrium point is a saddle point.

Theorem 2.

$E_1={\mathbb{K}}_1,\mathrm{0,0}$ of system (1) is asymptotically stable if $g_2<\phi {\mathbb{K}}_1$ and ${\mathbb{K}}_1{\mu }_1\nu <\rho$. Otherwise $E_1$ is a saddle point.

Proof.

The Jacobian matrix for $E_1$ is given by the following equation:$$\begin{array}{}\text{(7)}& J{E}_1=\begin{array}{ccc}-g_1& -\sigma {\mathbb{K}}_1& -{\mathbb{K}}_1\\ 0& g_2-\phi {\mathbb{K}}_1& 0\\ 0& 0& {\mathbb{K}}_1{\mu }_1\nu -\rho \end{array}.\end{array}$$

The eigenvalues of the matrix are ${\mathrm{\Lambda }}_1=-g_1,{\mathrm{\Lambda }}_2=g_2-\phi {\mathbb{K}}_1,{\text{and\hspace{0.17em}}\mathrm{\Lambda }}_3={\mathbb{K}}_1{\mu }_1\nu -\rho$.

Hence, $E_1$ is asymptotically stable when the restrictions $g_2<\phi {\mathbb{K}}_1$ and $K{\mu }_1\nu <\rho$ hold and unstable when $g_2>\phi {\mathbb{K}}_1$ and $K{\mu }_1\nu >\rho$.

Theorem 3.

$E_2=0,{\mathbb{K}}_2,0$ of system (1) is asymptotically stable if ${\mathbb{K}}_2<2g_2$ and ${\mu }_2\nu {\mathbb{K}}_2<\rho$. Otherwise $E_2$ is a saddle point.

Proof.

The Jacobian matrix for $E_2$ is given by the following equation:$$\begin{array}{}\text{(8)}& J{E}_2=\begin{array}{ccc}-\sigma {\mathbb{K}}_2& 0& 0\\ -\phi {\mathbb{K}}_2& {\mathbb{K}}_2-2g_2& -\eta {\mathbb{K}}_2\\ 0& 0& {\mu }_2\eta {\mathbb{K}}_2-\rho \end{array}.\end{array}$$

The eigenvalues of the matrix are ${\mathrm{\Lambda }}_1=-\sigma {\mathbb{K}}_2,{\mathrm{\Lambda }}_2={\mathbb{K}}_2-2g_2,\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Lambda }}_3={\mu }_2\nu {\mathbb{K}}_2-\rho$.

Here, ${\mathrm{\Lambda }}_1<0$. If ${\mathbb{K}}_2<2g_2$, then ${\mathrm{\Lambda }}_2<0$ and if ${\mu }_2\nu {\mathbb{K}}_2<\rho$, then ${\mathrm{\Lambda }}_3<0$.

Therefore, $E_2$ is asymptotically stable if the above condition holds, otherwise $E_2$ is a saddle point.

Theorem 4.

$E_3=g_1{\mathbb{K}}_1{g}_2-\sigma {\mathbb{K}}_1{\mathbb{K}}_2{g}_2/g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2,g_1{\mathbb{K}}_2{g}_2-\phi {\mathbb{K}}_1{\mathbb{K}}_2{g}_1/g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2,0$ of system (1) is asymptotically stable if ${\upsilon }_1>0$, ${\upsilon }_2>0$, ${\upsilon }_3>0$, and ${\upsilon }_1{\upsilon }_2-{\upsilon }_3>0$.

Proof.

At $E_3$, the Jacobian matrix is given by the following equation:$$\begin{array}{}\text{(9)}& \begin{array}{ccc}{A}_{\ast }& B_{\ast }& C_{\ast }\\ D_{\ast }& E_{\ast }& F_{\ast }\\ 0& 0& G_{\ast }\end{array},\end{array}$$where$$\begin{array}{c}\text{(10)}& A_{\ast }=g_1-\frac{2g_1}{{\mathbb{K}}_1}\frac{{\mathbb{K}}_1{g}_1{g}_2-\sigma {\mathbb{K}}_1{\mathbb{K}}_2{g}_2}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2}-\sigma \frac{{\mathbb{K}}_2{g}_1{g}_2-\phi {\mathbb{K}}_1{\mathbb{K}}_2{g}_1}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2},B_{\ast }=-\sigma \frac{{\mathbb{K}}_1{g}_1{g}_2-\sigma {\mathbb{K}}_1{\mathbb{K}}_2{g}_2}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2},\\ C_{\ast }=-\nu \frac{{\mathbb{K}}_1{g}_1{g}_2-\sigma {\mathbb{K}}_1{\mathbb{K}}_2{g}_2}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2},\\ D_{\ast }=-\phi \frac{{\mathbb{K}}_2{g}_1{g}_2-\phi {\mathbb{K}}_1{\mathbb{K}}_2{g}_1}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2},\\ E_{\ast }=g_2-\frac{2g_2}{{\mathbb{K}}_1}\frac{{\mathbb{K}}_1{g}_1{g}_2-\phi {\mathbb{K}}_1{\mathbb{K}}_2{g}_1}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2}-\phi \frac{{\mathbb{K}}_1{g}_1{g}_2-\sigma {\mathbb{K}}_1{\mathbb{K}}_2{g}_2}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2},\\ F_{\ast }=-\eta \frac{{\mathbb{K}}_1{g}_1{g}_2-\phi {\mathbb{K}}_1{\mathbb{K}}_2{g}_1}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2},\\ G_{\ast }={\mu }_1\nu \frac{{\mathbb{K}}_1{g}_1{g}_2-\sigma {\mathbb{K}}_1{\mathbb{K}}_2{g}_2}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2}+{\mu }_2\eta \frac{{\mathbb{K}}_2{g}_1{g}_2-\phi {\mathbb{K}}_1{\mathbb{K}}_2{g}_1}{g_1{g}_2-\sigma \phi {\mathbb{K}}_1{\mathbb{K}}_2}-\rho .\end{array}$$

The characteristic equation is as follows:$$\begin{array}{}\text{(11)}& {\mathrm{\Lambda }}^3-A_{\ast }+E_{\ast }+G_{\ast }{\mathrm{\Lambda }}^2+E_{\ast }{G}_{\ast }+A_{\ast }{G}_{\ast }+A_{\ast }{E}_{\ast }-D_{\ast }{B}_{\ast }\mathrm{\Lambda }-A_{\ast }{E}_{\ast }{G}_{\ast }-B_{\ast }{D}_{\ast }{G}_{\ast }=0.\end{array}$$

This can be written as follows:$$\begin{array}{}\text{(12)}& {\mathrm{\Lambda }}^3+{\upsilon }_1{\mathrm{\Lambda }}^2+{\upsilon }_2\mathrm{\Lambda }+{\upsilon }_3=0,\end{array}$$where ${\upsilon }_1=-A_{\ast }+E_{\ast }+G_{\ast }$, ${\upsilon }_2=E_{\ast }{G}_{\ast }+A_{\ast }{G}_{\ast }+A_{\ast }{E}_{\ast }-D_{\ast }{B}_{\ast }$, and ${\upsilon }_3=A_{\ast }{E}_{\ast }{G}_{\ast }-B_{\ast }{D}_{\ast }{G}_{\ast }.$

The equilibrium point $E_3$ is stable by the Routh–Hurwitz criterion if ${\upsilon }_1>0$, ${\upsilon }_2>0$, ${\upsilon }_3>0$, and ${\upsilon }_1{\upsilon }_2-{\upsilon }_3>0$.

Theorem 5.

$$\begin{array}{c}\text{(13)}& E^{\ast }=a_1^{\ast }=\frac{{\mathbb{K}}_1\nu g_2\mathrm{\varrho }+{\mathbb{K}}_2\eta g_1{\mu }_2\eta -\nu g_2{\mu }_2-\sigma \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},a_2^{\ast }=\frac{{\mathbb{K}}_2\eta g_1\mathrm{\varrho }+{\mathbb{K}}_1\nu g_2{\mu }_1\nu -\eta g_2{\mu }_1-\phi \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ b^{\ast }=\frac{g_1{g}_2\nu {\mu }_1{\mathbb{K}}_1+\eta {\mu }_2{\mathbb{K}}_2-\mathrm{\varrho }+{\mathbb{K}}_1{\mathbb{K}}_2\sigma \phi \mathrm{\varrho }-\sigma \nu g_2{\mu }_1-\phi \eta g_1{\mu }_2}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\end{array}$$of system (1) is locally stable if $b_1>0$, $b_2>0$, $b_3>0$, and $b_1{b}_2-b_3>0$.

Proof.

The Jacobian matrix at $E^{\ast }$ of (1) is as follows:$$\begin{array}{}\text{(14)}& J{E}^{\ast }=\begin{array}{ccc}{M}_{\ast }& N_{\ast }& P_{\ast }\\ Q_{\ast }& R_{\ast }& S_{\ast }\\ T_{\ast }& U_{\ast }& V_{\ast }\end{array},\end{array}$$where$$\begin{array}{c}\text{(15)}& M_{\ast }=g_1-\frac{2g_1}{{\mathbb{K}}_1}\frac{{\mathbb{K}}_1\nu g_2\mathrm{\varrho }+{\mathbb{K}}_2\eta g_1{\mu }_2\eta -\nu g_2{\mu }_2-\sigma \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2}-\sigma \frac{{\mathbb{K}}_2\eta g_1\mathrm{\varrho }+{\mathbb{K}}_1\nu g_2{\mu }_1\nu -\eta g_2{\mu }_1-\phi \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2}\\ -\nu \frac{g_1{g}_2\nu {\mu }_1{\mathbb{K}}_1+\eta {\mu }_2{\mathbb{K}}_2-\mathrm{\varrho }+{\mathbb{K}}_1{\mathbb{K}}_2\sigma \phi \mathrm{\varrho }-\sigma \nu g_2{\mu }_1-\phi \eta g_1{\mu }_2}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ N_{\ast }=-\sigma \frac{{\mathbb{K}}_1\nu g_2\mathrm{\varrho }+{\mathbb{K}}_2\eta g_1{\mu }_2\eta -\nu g_2{\mu }_2-\sigma \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ P_{\ast }=-\nu \frac{{\mathbb{K}}_1\nu g_2\mathrm{\varrho }+{\mathbb{K}}_2\eta g_1{\mu }_2\eta -\nu g_2{\mu }_2-\sigma \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ Q_{\ast }=-\phi \frac{{\mathbb{K}}_2\eta g_1\mathrm{\varrho }+{\mathbb{K}}_1\nu g_2{\mu }_1\nu -\eta g_2{\mu }_1-\phi \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ R_{\ast }=g_2-\frac{2g_2}{K_2}\frac{{\mathbb{K}}_2\eta g_1\mathrm{\varrho }+{\mathbb{K}}_1\nu g_2{\mu }_1\nu -\eta g_2{\mu }_1-\phi \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2}\\ -\phi \frac{{\mathbb{K}}_1\nu g_2\mathrm{\varrho }+{\mathbb{K}}_2\eta g_1{\mu }_2\eta -\nu g_2{\mu }_2-\sigma \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2}\\ -\eta \frac{g_1{g}_2\nu {\mu }_1{\mathbb{K}}_1+\eta {\mu }_2{\mathbb{K}}_2-\mathrm{\varrho }+{\mathbb{K}}_1{\mathbb{K}}_2\sigma \phi \mathrm{\varrho }-\sigma \nu g_2{\mu }_1-\phi \eta g_1{\mu }_2}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ S_{\ast }=-\eta \frac{{\mathbb{K}}_2\eta g_1\mathrm{\varrho }+{\mathbb{K}}_1\nu g_2{\mu }_1\nu -\eta g_2{\mu }_1-\phi \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ T_{\ast }={\mu }_1\nu \frac{g_1{g}_2\nu {\mu }_1{\mathbb{K}}_1+\eta {\mu }_2{\mathbb{K}}_2-\mathrm{\varrho }+{\mathbb{K}}_1{\mathbb{K}}_2\sigma \phi \mathrm{\varrho }-\sigma \nu g_2{\mu }_1-\phi \eta g_1{\mu }_2}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ U_{\ast }={\mu }_2\eta \frac{g_1{g}_2\nu {\mu }_1{\mathbb{K}}_1+\eta {\mu }_2{\mathbb{K}}_2-\mathrm{\varrho }+{\mathbb{K}}_1{\mathbb{K}}_2\sigma \phi \mathrm{\varrho }-\sigma \nu g_2{\mu }_1-\phi \eta g_1{\mu }_2}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2},\\ V_{\ast }={\mu }_1\nu \frac{{\mathbb{K}}_1\nu g_2\mathrm{\varrho }+{\mathbb{K}}_2\eta g_1{\mu }_2\eta -\nu g_2{\mu }_2-\sigma \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2}\\ +{\mu }_2\eta \frac{{\mathbb{K}}_2\eta g_1\mathrm{\varrho }+{\mathbb{K}}_1\nu g_2{\mu }_1\nu -\eta g_2{\mu }_1-\phi \mathrm{\varrho }}{\sigma \nu {\mu }_1{\mathbb{K}}_1{\mathbb{K}}_2\eta +\phi \nu {\mu }_2{\mathbb{K}}_1{\mathbb{K}}_2\eta -{\nu }^2{g}_2{\mu }_1{\mathbb{K}}_1-{\eta }^2{g}_1{\mu }_2{\mathbb{K}}_2}-\rho .\end{array}$$

The characteristic equation is as follows:$$\begin{array}{c}\text{(16)}& {\mathrm{\Lambda }}^3-M_{\ast }+R_{\ast }+V_{\ast }{\mathrm{\Lambda }}^2+R_{\ast }{V}_{\ast }-U_{\ast }{S}_{\ast }+M_{\ast }{V}_{\ast }-T_{\ast }{P}_{\ast }+M_{\ast }{R}_{\ast }-Q_{\ast }{N}_{\ast }\mathrm{\Lambda }-M_{\ast }{R}_{\ast }{V}_{\ast }-S_{\ast }{U}_{\ast }-N_{\ast }{Q}_{\ast }{V}_{\ast }-S_{\ast }{T}_{\ast }+P_{\ast }{Q}_{\ast }{U}_{\ast }-T_{\ast }{R}_{\ast }=0.\end{array}$$

This can be written as follows:$$\begin{array}{}\text{(17)}& {\mathrm{\Lambda }}^3+b_1{\mathrm{\Lambda }}^2+b_2\mathrm{\Lambda }+b_3=0.\end{array}$$where$$\begin{array}{c}\text{(18)}& b_1=-M_{\ast }+R_{\ast }+V_{\ast },b_2=R_{\ast }{V}_{\ast }-U_{\ast }{S}_{\ast }+M_{\ast }{V}_{\ast }-T_{\ast }{P}_{\ast }+M_{\ast }{R}_{\ast }-Q_{\ast }{N}_{\ast },\\ b_3=-M_{\ast }{R}_{\ast }{V}_{\ast }-S_{\ast }{U}_{\ast }-N_{\ast }{Q}_{\ast }{V}_{\ast }-S_{\ast }{T}_{\ast }+P_{\ast }{Q}_{\ast }{U}_{\ast }-T_{\ast }{R}_{\ast }.\end{array}$$Here, $E^{\ast }$ is stable if $b_1>0$, $b_2>0$, $b_3>0$, and $b_1{b}_2-b_3>0$.

The stability of (1) in the nonspatial case is described in details above. Then, the Turing instability of system (1) will be investigated by incorporating diffusive terms.

3. Spatiotemporal Model and Bifurcation Analysis

In this section, we extend (1) by introducing the random dispersal of all three species. The associated spatiotemporal system is given by the following equation:$$\begin{array}{c}\text{(19)}& \frac{\partial a_1}{\partial t}=D_1\frac{{\partial }^2{a}_1}{\partial x^2}+g_1{a}_{1}1-\frac{a_1}{{\mathbb{K}}_1}-\sigma a_1{a}_2-\nu a_{1}b,\frac{\partial a_2}{\partial t}=D_2\frac{{\partial }^2{a}_2}{\partial x^2}+g_2{a}_{2}1-\frac{a_2}{{\mathbb{K}}_2}-\phi a_1{a}_2-\eta a_{2}b,\\ \frac{\partial b}{\partial t}=D_3\frac{{\partial }^{2}b}{\partial x^2}+{\mu }_1\nu a_{1}b+{\mu }_2\eta a_{2}b-\mathrm{\varrho }b.\end{array}$$with the following initial conditions:$$\begin{array}{c}\text{(20)}& a_{1}X,0=a_{1_0}X\ge 0,a_{2}X,0=a_{2_0}X\ge 0,\\ bX,0=b_{0}X\ge 0,X\in {\mathrm{\Lambda }}^{\ast },\end{array}$$and homogeneous boundary conditions as follows:$$\begin{array}{}\text{(21)}& \frac{\partial a_1}{\partial n}=\frac{\partial a_2}{\partial n}=\frac{\partial b}{\partial n}=0,X\in \partial {\mathrm{\Lambda }}^{\ast },t>0.\end{array}$$

Here, $\mathrm{\Lambda }$ is a bounded domain in $R^2$ with the smooth boundary $\partial {\mathrm{\Lambda }}^{\ast }$; $\partial /\partial n$ represents the outward drawn normal derivative on $\partial {\mathrm{\Lambda }}^{\ast }$ while $n$ is the outward unit normal vector of $\partial {\mathrm{\Lambda }}^{\ast }$. The densities of the three species are denoted by $a_{1}x,t$, $a_{2}x,t$, and $bx,t$, respectively, at spatial position $x$ and time $t$. The positive parameters $d_i$ (i = 1, 2, 3) denote the diffusion coefficient of the species.

The Jacobian matrix $J_0$ for the model in (2) at the equilibrium point $E^{\ast }=a_1^{\ast },a_2^{\ast },b^{\ast }$ is as follows:$$\begin{array}{}\text{(22)}& \begin{array}{ccc}{p}^{\ast }& -\sigma a_1& -\nu a_1\\ -\phi a_2& q^{\ast }& -\eta a_2\\ {\mu }_1\nu b& {\mu }_2\eta b& r^{\ast }\end{array},\end{array}$$where $p^{\ast }=g_1-2g_1{a}_1/{\mathbb{K}}_1-\sigma a_2-\nu b$, $q^{\ast }=g_2-2g_2{a}_2/{\mathbb{K}}_2-\phi a_1-\eta b$, and $r^{\ast }={\mu }_1\nu a_1+{\mu }_2\eta a_2-\mathrm{\varrho }$.

By implementing a small amplitude spatiotemporal perturbation around the coexistence steady state $E^{\ast }=a_1^{\ast },a_2^{\ast },b^{\ast }$ and by linearizing, we get the following:$$\begin{array}{}\text{(23)}& \begin{array}{ccc}{p}^{\ast }& -\sigma a_1& -\nu a_1\\ -\phi a_2& q^{\ast }& -\eta a_2\\ {\mu }_1\nu b& {\mu }_2\eta b& r^{\ast }\end{array}-k_{\ast }^2\begin{array}{ccc}{d}_1^{\ast }& 0& 0\\ 0& d_2^{\ast }& 0\\ 0& 0& d_3^{\ast }\end{array}=\begin{array}{ccc}{p}^{\ast }-d_1^{\ast }{k}_{\ast }^2& -\sigma a_1& -\nu a_1\\ -\phi a_2& q^{\ast }-d_2^{\ast }{k}_{\ast }^2& -\eta a_2\\ {\mu }_1\nu b& {\mu }_2\eta b& r^{\ast }-d_3^{\ast }{k}_{\ast }^2\end{array}.\end{array}$$

Here, $k^{\ast }=k^{\ast }x,k^{\ast }y,{\mathrm{a}\mathrm{n}\mathrm{d}k}^{\ast }z$ is the wave number vector and $k^{\ast }=k^{\ast }$ represents the wave number. The characteristic equation of (23) is as follows:$$\begin{array}{}\text{(24)}& {\mathrm{\Lambda }}^3+P{k}_{\ast }^2{\mathrm{\Lambda }}^2+Q{k}_{\ast }^2\mathrm{\Lambda }+R{k}_{\ast }^2=0.\end{array}$$

Here,$$\begin{array}{c}\text{(25)}& P{k}_{\ast }^2=p^{\ast }+q^{\ast }+r^{\ast }-d_1^{\ast }+d_2^{\ast }+d_3^{\ast }{k}_{\ast }^2,Q{k}_{\ast }^2=k_{\ast }^4{d}_1^{\ast }{d}_2^{\ast }+d_2^{\ast }{d}_3^{\ast }-d_1^{\ast }{d}_3^{\ast }\\ +k_{\ast }^2-d_3^{\ast }{q}^{\ast }-d_2^{\ast }{r}^{\ast }+d_3^{\ast }{p}^{\ast }+d_1^{\ast }{r}^{\ast }-d_2^{\ast }{p}^{\ast }-d_1^{\ast }{q}^{\ast }\\ +q^{\ast }{r}^{\ast }+{\eta }^2{a}_2{\mu }_{2}b-p^{\ast }{r}^{\ast }-{\mu }_1{\nu }^{2}b{a}_1+p^{\ast }{q}^{\ast }-\sigma \phi a_1{a}_2,\\ R{k}_{\ast }^2=-k_{\ast }^6{d}_1^{\ast }{d}_2^{\ast }{d}_3^{\ast }+k_{\ast }^4{d}_2^{\ast }{d}_3^{\ast }{p}^{\ast }+d_1^{\ast }{d}_3^{\ast }{q}^{\ast }+d_1^{\ast }{d}_2^{\ast }{r}^{\ast }\\ +k_{\ast }^2{d}_3^{\ast }\sigma \phi a_1{a}_2-d_3^{\ast }{p}^{\ast }{q}^{\ast }-d_2^{\ast }{p}^{\ast }{r}^{\ast }-d_1^{\ast }{q}^{\ast }{r}^{\ast }\\ -k_{\ast }^2{d}_1^{\ast }{\eta }^2{a}_2{\mu }_{2}b-k_{\ast }^2{d}_2^{\ast }{\nu }^2{a}_1{\mu }_{1}b+p^{\ast }{q}^{\ast }{r}^{\ast }+{\eta }^2{a}_2{\mu }_{2}b{p}^{\ast }\\ -\sigma \phi a_1{a}_2{r}^{\ast }+\sigma \eta a_1{a}_2{\mu }_1\nu b+\nu \phi a_1{a}_2{\mu }_2\eta b+{\mu }_1{\nu }^2{a}_{1}b{q}^{\ast },\end{array}$$where $P{k}_{\ast }^2>0$ if $p^{\ast }+q^{\ast }+r^{\ast }>d_1^{\ast }+d_2^{\ast }+d_3^{\ast }{k}_{\ast }^2$.

System (2) is stable if $\text{tr}{J}_{\lambda }<0$ and $\det J_{\lambda }>0$. If it fails to meet the abovementioned conditions, the system is considered unstable. As a result, $R{k}_{\ast }^2$ can be reframed as follows:$$\begin{array}{}\text{(26)}& R{k}_{\ast }^2=R_3{k_{\ast }^2}^3+R_2{k_{\ast }^2}^2+R_1{k}_{\ast }^2+R_0,\end{array}$$where$$\begin{array}{c}\text{(27)}& R_3=-d_1{d}_2{d}_3<0\text{\hspace{0.17em}as\hspace{0.17em}}{d}_1=d_2=d_3>0,R_2=d_1^{\ast }{d}_2^{\ast }{r}^{\ast }+d_1^{\ast }{d}_3^{\ast }{q}^{\ast }+d_2^{\ast }{d}_3^{\ast }{p}^{\ast }>0,\\ R_1=d_3^{\ast }\sigma \phi a_1{a}_2-d_1^{\ast }{q}^{\ast }{r}^{\ast }-d_2^{\ast }{p}^{\ast }{r}^{\ast }-d_3^{\ast }{p}^{\ast }{q}^{\ast }-d_1^{\ast }{\eta }^2{a}_2{\mu }_{2}b-d_2^{\ast }{\nu }^2{a}_1{\mu }_{1}b>0,\\ R_0=p^{\ast }{q}^{\ast }{r}^{\ast }+{\eta }^2{a}_2{\mu }_{2}b{p}^{\ast }-\sigma \phi a_1{a}_2{r}^{\ast }+\sigma \eta a_1{a}_2{\mu }_1\nu b+\nu \phi a_1{a}_2{\mu }_2\eta b+{\mu }_1{\nu }^2{a}_{1}b{q}^{\ast }.\end{array}$$

Here, $R_0>0$ if $\sigma \phi a_1{a}_2{r}^{\ast }<p^{\ast }{q}^{\ast }{r}^{\ast }+{\eta }^2{a}_2{\mu }_{2}b{p}^{\ast }-\sigma \phi a_1{a}_2{r}^{\ast }+\sigma \eta a_1{a}_2{\mu }_1\nu b+\nu \phi a_1{a}_2{\mu }_2\eta b+{\mu }_1{\nu }^2{a}_{1}b{q}^{\ast }$ and $R_1>0$ if $d_3^{\ast }\sigma \phi a_1{a}_2>d_1^{\ast }{q}^{\ast }{r}^{\ast }-d_2^{\ast }{p}^{\ast }{r}^{\ast }-d_3^{\ast }{p}^{\ast }{q}^{\ast }-d_1^{\ast }{\eta }^2{a}_2{\mu }_{2}b-d_2^{\ast }{\nu }^2{a}_1{\mu }_{1}b$.

Therefore,$$\begin{array}{}\text{(28)}& R{k}_{\ast }^2=-R_3{k_{\ast }^2}^3+R_2{k_{\ast }^2}^2+R_1{k}_{\ast }^2+R_0.\end{array}$$

Hence, $R{k}_{\ast }^2<0$ if $R_3>R_0+R_1+R_2$.

$R{k}_{\ast }^2$ is a criterion for Turing bifurcation spatial patterns if all other mentioned conditions have been met.

In the next section, system (2) can be solved by using the sinc interpolation method as follows.

4. Interpolation Method

The sinc method is a very effective numerical technique [36]. It is frequently used in many different areas of numerical analysis, including quadrature, integral and differential equation solution, and interpolation [37]. When singularities or unbounded domains are present, the method may effectively handle them. As the number of bases rises, the method’s error converges exponentially to zero.

System (2) can be solved by dividing the interval $\mathrm{0,2}\pi$ into $\mathbb{N}$ equally spaced nodes, each of which corresponds to a regular region $\xi =\mathrm{0,2}\pi \times \mathrm{0,2}\pi$. A periodic sinc function ${\mathbb{S}}_{\mathbb{N}}$ [38] is defined as follows:$$\begin{array}{}\text{(29)}& {\mathbb{S}}_{\mathbb{N}}=\frac{\sin \pi x/h}{2\pi /h\tan x/2},\end{array}$$where $h=2\pi /n$. On the other hand, ${\mathbb{S}}_{\mathbb{N}}{x}_{\alpha }-x_{\beta }$ is an $\mathbb{N}$ order unit matrix.

For $h>0$,$$\begin{array}{}\text{(30)}& \text{span}{\mathbb{S}}_{\mathbb{N}}x-\beta h,\beta =\mathrm{1,2},\dots \mathbb{N}.\end{array}$$

For functions $a_{1}x,a_{2}x$ and $bx$ defined on $\mathrm{0,2}\pi$, let ${\mathbb{I}}_{\mathbb{N}}$ be the interpolation operator. We have [39] the following equations:$$\begin{array}{c}\text{(31)}& a_{1}x,y,t\sim {\mathbb{I}}_{\mathbb{N}}{a}_{1}x,y,t={\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}x-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}y-y_{\alpha }{a}_1{x}_{\alpha },y_{\beta },t,a_{2}x,y,t\sim {\mathbb{I}}_{\mathbb{N}}{a}_{2}x,y,t\hspace{1em}={\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}x-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}y-y_{\alpha }{a}_2{x}_{\alpha },y_{\beta },t,\\ bx,y,t\sim {\mathbb{I}}_{\mathbb{N}}bx,y,t\hspace{1em}={\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}x-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}y-y_{\alpha }b{x}_{\alpha },y_{\beta },t.\end{array}$$

The following relation holds at the collocation nodes $x_u,y_v$:$$\begin{array}{c}\text{(32)}& a_1{x}_u,y_v,t={\mathbb{I}}_{\mathbb{N}}{a}_1{x}_u,y_v,t={\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}{x}_u-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}{y}_v-y_{\alpha }{a}_1{x}_{\alpha },y_{\beta },t,\\ a_2{x}_u,y_v,t={\mathbb{I}}_{\mathbb{N}}{a}_2{x}_u,y_v,t\\ ={\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}{x}_u-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}{y}_v-y_{\alpha }{a}_2{x}_{\alpha },y_{\beta },t,\\ b{x}_u,y_v,t={\mathbb{I}}_{\mathbb{N}}b{x}_u,y_v,t\\ ={\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}{x}_u-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}{y}_v-y_{\alpha }b{x}_{\alpha },y_{\beta },t.\\ \text{(33)}& {a_1}^{c,d}{x}_u,y_v,t\sim {\mathbb{I}}_{\mathbb{N}}{a}_1^{c,d}{x}_u,y_v,t=\frac{{\partial }^{c+d}{a}_1{x}_u,y_v,t}{\partial x^c\partial y^d}{\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}^c{x}_u-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}^d{y}_q-y_{\alpha }{a}_1{x}_{\alpha },y_{\beta },t,\\ {a_2}^{c,d}{x}_u,y_v,t\sim {\mathbb{I}}_{\mathbb{N}}{a}_2^{c,d}{x}_u,y_v,t\\ =\frac{{\partial }^{c+d}{a}_2{x}_u,y_v,t}{\partial x^c\partial y^d}{\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}^c{x}_u-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}^d{y}_v-y_{\alpha }{a}_2{x}_{\alpha },y_{\beta },t,\\ b^{c,d}{x}_u,y_v,t\sim {\mathbb{I}}_{\mathbb{N}}{n}^{c,d}{x}_u,y_v,t\\ =\frac{{\partial }^{c+d}n{x}_u,y_v,t}{\partial x^c\partial y^d}{\displaystyle \sum }_{\alpha =1}^{\mathbb{N}}{\displaystyle \sum }_{\beta =1}^{\mathbb{N}}{\mathbb{S}}_{\mathbb{N}}^c{x}_u-x_{\alpha }{\mathbb{S}}_{\mathbb{N}}^d{y}_v-y_{\alpha }b{x}_{\alpha },y_{\beta },t.\end{array}$$

Here, $a_1$, $a_2$, and $b$ are as follows:$$\begin{array}{c}\text{(34)}& {{a_1}_{11},\dots ,{a_1}_{\mathbb{N}1},{a_1}_{12},\dots {a_1}_{\mathbb{N}2},{a_1}_{1\mathbb{N}},{a_1}_{2\mathbb{N}},\dots {a_1}_{\mathbb{N}\mathbb{N}}}^{\mathbb{T}},{{a_2}_{11},\dots ,{a_2}_{\mathbb{N}1},{a_2}_{12},\dots {a_2}_{\mathbb{N}2},{a_2}_{1\mathbb{N}},{a_2}_{2\mathbb{N}},\dots {a_2}_{\mathbb{N}\mathbb{N}}}^{\mathbb{T}},\\ {b_{11},\dots ,b_{\mathbb{N}1},b_{12},\dots b_{\mathbb{N}2},b_{1\mathbb{N}},b_{2\mathbb{N}},\dots b_{\mathbb{N}\mathbb{N}}}^{\mathbb{T}},\end{array}$$respectively.

As a result, (33) can be expressed as follows:$$\begin{array}{c}\text{(35)}& a_1={\mathbb{D}}^{\mathrm{0,0}}{a}_1,a_2={\mathbb{D}}^{\mathrm{0,0}}{a}_2,b={\mathbb{D}}^{\mathrm{0,0}}b,{a_1}^{c,d}={\mathbb{D}}_{\mathbb{N}}^{c,d}{a}_1,{a_2}^{c,d}={\mathbb{D}}_{\mathbb{N}}^{c,d}{a}_2,b^{c,d}={\mathbb{D}}_{\mathbb{N}}^{c,d}b,\end{array}$$where ${\mathbb{D}}_{\mathbb{N}}^{c,d}={\mathbb{D}}_{\mathbb{N}}^c\otimes {\mathbb{D}}_{\mathbb{N}}^d$ is the matrix Kronecker product. ${\mathbb{D}}_{\mathbb{N}}^c$ and ${\mathbb{D}}_{\mathbb{N}}^d$, and ${\mathbb{D}}^{\mathrm{0,0}}={\mathbb{I}}_{\mathbb{N}}\otimes {\mathbb{I}}_{\mathbb{N}}$. ${\mathbb{I}}_{\mathbb{N}}$ is an $\mathbb{N}$ order unit matrix. The discrete form of (19) can be written as follows using equations (32), (34), and (35),$$\begin{array}{c}\text{(36)}& \frac{\partial }{\partial t}\begin{array}{c}{a}_1\\ a_2\\ b\end{array}=\begin{array}{ccc}{d}_1^{\ast }\mathbb{D}+g_1{\mathbb{D}}^{\mathrm{0,0}}& 0& 0\\ 0& d_2^{\ast }\mathbb{D}+g_2{\mathbb{D}}^{\mathrm{0,0}}& 0\\ 0& 0& d_3^{\ast }\mathbb{D}-\mathrm{\varrho }{\mathbb{D}}^{\mathrm{0,0}}\end{array}\begin{array}{c}{a}_1\\ a_2\\ b\end{array}=\begin{array}{c}{m}_1{a}_1,a_2,b\\ m_2{a}_1,a_2,b\\ m_3{a}_1,a_2,b\end{array}.\end{array}$$

Here,$$\begin{array}{c}\text{(37)}& a_1,a_2,b={a_1}_{11},\dots ,{a_1}_{\mathbb{N}1},{a_1}_{12},\dots {a_1}_{\mathbb{N}\mathbb{N}}{a_2}_{11},\dots ,{a_2}_{\mathbb{N}1},{a_2}_{12},\dots {a_2}_{\mathbb{N}\mathbb{N}}\\ b_{11},\dots ,b_{\mathbb{N}1},b_{12},\dots b_{\mathbb{N}\mathbb{N}},\\ \mathbb{D}={\mathbb{D}}_{\mathbb{N}}^{\mathrm{2,0}}+{\mathbb{D}}_{\mathbb{N}}^{\mathrm{0,2}},E={\mathbb{D}}^{\mathrm{0,0}},\\ m_1{a}_1,a_2,b,m_2{a}_1,a_2,b,m_3{a}_1,a_2,b\\ =m_1{a_1}_{11},{a_2}_{11},b_{11},\dots m_1{a_1}_{\mathbb{N}\mathbb{N}},{a_2}_{\mathbb{N}\mathbb{N}},b_{\mathbb{N}\mathbb{N}}\\ m_2{a_1}_{11},{a_2}_{11},b_{11},\dots m_2{a_1}_{\mathbb{N}\mathbb{N}},{a_2}_{\mathbb{N}\mathbb{N}},b_{\mathbb{N}\mathbb{N}}\\ m_3{a_1}_{11},{a_2}_{11},b_{11},\dots m_3{a_1}_{\mathbb{N}\mathbb{N}},{a_2}_{\mathbb{N}\mathbb{N}},b_{\mathbb{N}\mathbb{N}}.\end{array}$$

Therefore,$$\begin{array}{c}\text{(38)}& m_1{a}_1,a_2,b=g_1{a}_1-\frac{a_1}{{\mathbb{K}}_1}-\alpha a_1{a}_2-\nu a_{1}b,m_2{a}_1,a_2,b=g_2{a}_2-\frac{a_2}{{\mathbb{K}}_2}-\phi a_1{a}_2-\delta a_{2}b,\\ m_3{a}_1,a_2,b={\mu }_1\nu a_{1}b+{\mu }_2\eta a_{2}b.\end{array}$$

To get the numerical solution of (19), we use the numerical scheme NDSolve in Mathematica software for solving various initial conditions.

The Turing patterns near the bifurcative value were described using amplitude equations. For the aforementioned models, the dynamics obtained by varying the bifurcative values were interpreted using weak nonlinear analysis in the following section.

5. Amplitude Equation and Stability Analysis

In this section, to obtain the amplitude equations that make up the various Turing patterns, we will use the multiple-scale analysis method. We can then determine the necessary conditions for the appearance of various Turing patterns by examining the stability of the amplitude equations, allowing us to construct various Turing patterns through numerical simulation.

We use Taylor’s expansion to expand (2) at equilibrium $E^{\ast }$, since we know that the amplitude of the equation cannot be determined directly. The system solution could be expanded as follows [40]: $\mathfrak{c}={\mathfrak{c}}_0+{{\displaystyle \sum }}_{i=1}^3{Z}_i{e}^{\mathrm{i}{\mathrm{k}}_{i}r}$ + Complex conjugate. Furthermore, system (2) can be written as follows:$$\begin{array}{}\text{(39)}& \frac{\partial \mathfrak{c}}{\partial t}=L\mathfrak{c}+N\mathfrak{c},\mathfrak{c},\end{array}$$where $\mathfrak{c}=\begin{array}{c}{a}_1\\ a_2\\ b\end{array}$ is the variable and $L=B+D{\nabla }^2$ is the linear operator such that $L=\begin{array}{ccc}{\omega }_{11}& {\omega }_{12}& {\omega }_{13}\\ {\omega }_{21}& {\omega }_{22}& {\omega }_{23}\\ {\omega }_{31}& {\omega }_{32}& {\omega }_{33}\end{array}$ provided ${\omega }_{11}=g_1-\sigma a_2-\gamma b+d_1{\nabla }^2$, ${\omega }_{12}=-\sigma a_1$, ${\omega }_{13}=-\gamma a_1$, ${\omega }_{21}=-\phi a_2$, ${\omega }_{22}=g_2-\phi a_1-\eta b+d_2{\nabla }^2$, ${\omega }_{23}=-\eta a_2$, ${\omega }_{31}={\mu }_1\gamma b$, ${\omega }_{32}={\mu }_2\eta b$, and ${\omega }_{33}={\mu }_1\gamma a_1+{\mu }_2\eta a_2-\rho +d_3{\nabla }^2$, and $N=\begin{array}{c}-g_1{a}_1^2/K\\ 0\\ 0\end{array}$ is the nonlinear term.

When $\Omega$ approaches ${\Omega }_{\mathfrak{c}}$, we must examine the dynamical behaviour and then expand as follows:$$\begin{array}{}\text{(40)}& {\Omega }_{\mathfrak{c}}-\Omega =\epsilon {\Omega }_1+{ϵ}^2{\Omega }_2+\dots ,\end{array}$$where $ϵ$ is a sufficiently small parameter and $\Omega$ is the critical value.

As the series form of $ϵ$, we expand $\mathfrak{c}$ and $N$ as follows:$$\begin{array}{c}\text{(41)}& \mathfrak{c}=\begin{array}{c}{a}_1\\ a_2\\ b\end{array}=\begin{array}{c}{a_1}_1\\ {a_2}_1\\ b_1\end{array}ϵ+\begin{array}{c}{a_1}_2\\ {a_2}_2\\ b_2\end{array}{ϵ}^2+\dots ,N=\begin{array}{c}\frac{-g_1{a_1^2}_1}{K}{\epsilon }^2+\frac{-g_1{a_1}_1{a_2}_1}{K}+{a_1^3}_1{\epsilon }^3+o{\epsilon }^4\\ 0\\ 0\end{array}.\end{array}$$Here, $L$ can be expressed as $L=L_{\mathfrak{c}}+{\Omega }_{\mathfrak{c}}-\Omega$.

Let$$\begin{array}{c}\text{(42)}& T_0=t,T_1=ϵt,\\ T_2={ϵ}^{2}t,\dots ,\end{array}$$such that $\partial /\partial t=\partial /\partial T_0+ϵ\partial /\partial T_1+{ϵ}^2\partial /\partial T_2+\dots$

$T_i$ is a dependent variable that we use as a base for our time calculation, whereas amplitude is a slow variable.$$\begin{array}{}\text{(43)}& \frac{\partial W}{\partial t}=ϵ\frac{\partial W}{\partial T_1}+{ϵ}^2\frac{\partial W}{\partial T_2}+\dots ,\end{array}$$

By incorporating the prior equations into (19) and extending to different orders of $\epsilon$, we get the following:$$\begin{array}{c}\text{(44)}& \epsilon :L_c\begin{array}{c}{a_1}_1\\ {a_2}_1\\ b_1\end{array}=0,{\epsilon }^2:L_c\begin{array}{c}{a_1}_2\\ {a_2}_2\\ b_2\end{array}=\frac{\partial }{\partial T_1}\begin{array}{c}{a_1}_1\\ {a_2}_1\\ b_1\end{array}-{\Omega }_{1}M\begin{array}{c}{a_1}_1\\ {a_2}_1\\ b_1\end{array}-\begin{array}{c}\frac{-g_1{a}_1^2}{K}\\ 0\\ 0\end{array},\\ {\epsilon }^3:L_c\begin{array}{c}{a_1}_3\\ {a_2}_3\\ b_3\end{array}=\frac{\partial }{\partial T_1}\begin{array}{c}{a_1}_2\\ {a_2}_2\\ b_2\end{array}+\frac{\partial }{\partial T_2}\begin{array}{c}{a_1}_1\\ {a_2}_1\\ b_1\end{array}-{\Omega }_{1}M\begin{array}{c}{a_1}_2\\ {a_2}_2\\ b_2\end{array}\\ -{\Omega }_{2}M\begin{array}{c}{a_1}_1\\ {a_2}_1\\ b_1\end{array}-\begin{array}{c}\frac{-g_1{a}_1^2}{K}\\ 0\\ 0\end{array}\begin{array}{c}{a_1}_1\\ {a_2}_1\\ b_1\end{array}\begin{array}{c}{a_1}_2\\ {a_2}_2\\ b_2\end{array}-\begin{array}{c}\frac{-g_1{a}_1{a}_2}{K}\\ 0\\ 0\end{array}.\end{array}$$

An example of first order $ϵ$ is considered. Since the initial linear operator of the system is $L\mathfrak{c}$, then ${{a_1}_1,{a_2}_1,b_1}^T$ is the linear combination of the eigenvectors corresponding to the zero eigenvalue.$$\begin{array}{}\text{(45)}& \begin{array}{c}{a}_1\\ a_2\\ b\end{array}={\displaystyle \sum }_{i=1}^3\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}{e}^{\mathrm{i}{\mathrm{k}}_{i}r}+\text{Complex\hspace{0.17em}conjugate},\end{array}$$and have deduced that $x_i=H{y}_i$. Assuming $y_i=1$ and $z_i=0$, let $x_i=H$.$$\begin{array}{}\text{(46)}& \begin{array}{c}{a_1}_1\\ {a_2}_1\\ b_1\end{array}={\displaystyle \sum }_{i=1}^3\begin{array}{c}H\\ 1\\ 0\end{array}{W}_i{e}^{\mathrm{i}{\mathrm{k}}_{i}r}+\text{Complex\hspace{0.17em}conjugate}.\end{array}$$

Let us now consider the ${ϵ}^2$ case. By Fredholm’s solubility condition, the vector function of the preceding equation must be orthogonal to the operator $L_c^{+}$’s zero eigenvectors.

$L_c^{+}$’s zero eigenvectors are as follows:$$\begin{array}{}\text{(47)}& \begin{array}{c}1\\ H\\ 0\end{array}{e}^{-\mathrm{i}{\mathrm{k}}_{i}r}+\text{Complex\hspace{0.17em}conjugate}.\end{array}$$

The orthogonality criteria implies that$$\begin{array}{}\text{(48)}& \begin{array}{ccc}1& H& 0\end{array}{e}^{-\mathrm{i}{\mathrm{k}}_{i}r}{L}_c\begin{array}{c}{a_1}_2\\ {a_2}_2\\ b_2\end{array}=0.\end{array}$$

By only investigating at $e^{-\mathrm{i}{\mathrm{k}}_{1}r}$ in the following, we get another case by modifying the subscripts as follows:$$\begin{array}{}\text{(49)}& B+H\frac{\partial W_1}{\partial T_1}=h_{11}H+h_{12}+h_{21}H+h_{22}B{W}_1-Bh+Hh\frac{-g_1{a}_1^2}{K}\overline{W_2}\overline{W_3}.\end{array}$$

The same procedures can be used to achieve the same results as follows:$$\begin{array}{c}\text{(50)}& \begin{array}{c}{a_1}_2\\ {a_2}_2\\ b_2\end{array}=\begin{array}{c}{b}_0\\ h_0\\ g_0\end{array}+\begin{array}{c}{b}_i\\ h_i\\ g_i\end{array}{e}^{\mathrm{i}{\mathrm{k}}_{i}r}+\begin{array}{c}{b}_{\text{ii}}\\ h_{\text{ii}}\\ g_{\text{ii}}\end{array}{e}^{i2k_{i}r}+\begin{array}{c}{b}_{12}\\ h_{12}\\ g_{12}\end{array}{e}^{i{k}_1-k_{2}r}+\begin{array}{c}{b}_{23}\\ h_{23}\\ g_{23}\end{array}{e}^{i{k}_2-k_{3}r}\\ +\begin{array}{c}{b}_{31}\\ h_{31}\\ g_{31}\end{array}{e}^{i{k}_3-k_{1}r}+\text{Complex\hspace{0.17em}conjugate}.\end{array}$$

For ${ϵ}^3$ case,$$\begin{array}{c}\text{(51)}& B+H\frac{\partial z}{\partial T_1}+\frac{\partial z}{\partial T_2}={\Omega }_1{z}_1+{\Omega }_2{w}_1{h}_{11}H+h_{12}+h_{21}H+h_{22}B-\text{Bh}+\text{Hh}\frac{-g_1{a}_1^2}{K}\overline{Z_2{W}_3}+\overline{Z_3{W}_2}\\ -h\text{BH}+\mathrm{h}{\mathrm{H}}^2+\text{hB}+\text{hH}\frac{g_2{a}_2^2}{K}{W}_1^2+W_2^2+W_3^2{W}_1.\end{array}$$

After simplification, we derive the coefficient expressions as ${\Omega }_0=B+H$, ${\Omega }^{\ast }={\Omega }_{\mathfrak{c}}-\Omega /{\Omega }_{\mathfrak{c}}$, $f=-\text{hB}+\text{hH}-g_1{a}_1^2/K$, and $g=h\text{BH}+\mathrm{h}{\mathrm{H}}^2+\text{hB}+\text{hH}{g}_2{a}_2^2/K$. As a result, the amplitude equation is as follows:$$\begin{array}{c}\text{(52)}& {\Omega }_0\frac{\partial W_1}{\partial t}={\Omega }^{\ast }{W}_1+f\overline{W_2{W}_3}-g_1{W_1}^2+g_2{W_2}^2+{W_3}^2{W}_1,{\Omega }_0\frac{\partial W_2}{\partial t}\hspace{1em}={\Omega }^{\ast }{W}_2+f\overline{W_1{W}_3}-g_1{W_2}^2+g_2{W_1}^2+{W_3}^2{W}_2,\\ {\Omega }_0\frac{\partial W_3}{\partial t}\hspace{1em}={\Omega }^{\ast }{W}_3+f\overline{W_2{W}_1}-g{W_3}^2+g_2{W_2}^2+{W_1}^2{W}_3.\end{array}$$

The equation for the amplitude can be written as follows:$$\begin{array}{}\text{(53)}& W_j={\varsigma }_j\exp i{\iota }_j,\end{array}$$where the phase is $\varpi =2\pi /3$ and the mode ${\varsigma }_j$ can be defined as ${\varsigma }_j=W_j$.

From (53),$$\begin{array}{c}\text{(54)}& \frac{\partial \varpi }{\partial t}={\displaystyle \sum }_{j=1}^3\frac{\partial {\varpi }_j}{\partial t}=\frac{1}{i}{\displaystyle \sum }_{j=1}^3\frac{1}{A_j}\frac{\partial A_j}{\partial t}=-\frac{f}{{\Omega }_0}{\varsigma }_2^2{\varsigma }_3^2+{\varsigma }_3^2{\varsigma }_1^2+{\varsigma }_1^2{\varsigma }_2^2{{\varsigma }_1{\varsigma }_2{\varsigma }_3}^{-1}i\cos \varpi +\sin \varpi \\ +i{\Omega }_{0}3{\Omega }^{\ast }-g_1+2g_2{{W_1}_2{W}_2}^2+{W_3}^2.\end{array}$$

Considering the real part, we get the following:$$\begin{array}{}\text{(55)}& {\Omega }_0\frac{\partial \varsigma }{\partial t}=-f{\varsigma }_2^2{\varsigma }_3^2+{\varsigma }_3^2{\varsigma }_1^2+{\varsigma }_1^2{\varsigma }_2^2{{\varsigma }_1{\varsigma }_2{\varsigma }_3}^{-1}\sin \varpi .\end{array}$$

We derive the following mathematical relations by substituting the value of (53) in (52) and taking the real part of the equations into account as follows:$$\begin{array}{c}\text{(56)}& {\Omega }_0\frac{\partial {\varsigma }_1}{\partial t}={\Omega }^{\ast }{\varpi }_1+f{\varsigma }_2{\varsigma }_3\cos \varpi -g_1{\varsigma }_1^3-g_2{\varsigma }_2^2+{\varsigma }_3^2{\varsigma }_1,{\Omega }_0\frac{\partial {\varsigma }_2}{\partial t}\hspace{1em}={\Omega }^{\ast }{\varpi }_2+f{\varsigma }_3{\varsigma }_1\cos \varpi -g_1{\varsigma }_2^3-g_2{\varsigma }_2^3+{\varsigma }_1^2{\varsigma }_2,\\ {\Omega }_0\frac{\partial {\varsigma }_3}{\partial t}\hspace{1em}={\Omega }^{\ast }{\varpi }_3+f{\varsigma }_1{\varsigma }_2\cos \varpi -g_1{\varsigma }_3^3-g_2{\varsigma }_1^2+{\varsigma }_2^2{\varsigma }_3.\end{array}$$

When $\varpi =0$ and $\varpi =\pi$, the value of $f$ is positive and negative, respectively. If we consider $f>0$, then the system’s pattern will be stable and defined as $F_0$, if we consider $f<0$, then the system’s pattern will be stable and defined as $F_{\pi }$. As a result of (56), the following system of equations is obtained as follows:$$\begin{array}{c}\text{(57)}& {\Omega }_0\frac{\partial {\varsigma }_1}{\partial t}={\Omega }^{\ast }{\varpi }_1+f{\varsigma }_2{\varsigma }_3-g_1{\varsigma }_1^3-g_2{\varsigma }_2^2+{\varsigma }_3^2{\varsigma }_1,{\Omega }_0\frac{\partial {\varsigma }_2}{\partial t}\hspace{1em}={\Omega }^{\ast }{\varpi }_2+f{\varsigma }_3{\varsigma }_1-g_1{\varsigma }_2^3-g_2{\varsigma }_2^3+{\varsigma }_1^2{\varsigma }_2,\\ {\Omega }_0\frac{\partial {\varsigma }_3}{\partial t}\hspace{1em}={\Omega }^{\ast }{\varpi }_3+f{\varsigma }_1{\varsigma }_2-g_1{\varsigma }_3^3-g_2{\varsigma }_1^2+{\varsigma }_2^2{\varsigma }_3.\end{array}$$

We examine the stability properties of the amplitude equations using (57).

Theorem 6.

The prerequisite for stationary state stability is ${\Omega }^{\ast }<{\Omega }_0^{\ast }$ and the condition for the unstable is ${\Omega }^{\ast }>{\Omega }_0^{\ast }$, where ${\Omega }_0^{\ast }=0$.

Proof.

Consider the stationary state to be (O), and from (57), the prerequisite for stability is given by ${\varsigma }_1={\varsigma }_2={\varsigma }_3=0$ for ${\Omega }^{\ast }<{\Omega }_0^{\ast }=0$ and for ${\Omega }^{\ast }>{\Omega }_0^{\ast }$ it is unstable.

Theorem 7.

The prerequisite for stability for stripe patterns is ${\Omega }^{\ast }>{\Omega }_m^{\ast }$ and unstable if ${\Omega }^{\ast }<{\Omega }_m^{\ast }$ where ${\Omega }_s=f^2{g}_1/{g_1-g_2}^2$.

Proof.

Let us denote the stripe pattern as (M) and it is given by ${\varsigma }_1=\sqrt{{\Omega }^{\ast }/g_1}$ and ${\varsigma }_2={\varsigma }_3=0$.

Consider ${\varsigma }_j={\varsigma }_m+\theta {\varsigma }_j\text{\hspace{0.17em}and\hspace{0.17em}}j=\mathrm{1,2,3}$.

Substitute ${\varsigma }_1={\varsigma }_m+\theta {\varsigma }_1,{\varsigma }_2=\theta {\varsigma }_2,{\mathrm{a}\mathrm{n}\mathrm{d}\varsigma }_3=\theta {\varsigma }_3$ in (57) and replace the steady state condition. The following are the mathematical relationships:$$\begin{array}{c}\text{(58)}& {\Omega }_0\frac{\partial }{\partial t}\theta {\varsigma }_1={\Omega }^{\ast }-3g_1{\varsigma }_m^2\theta {\varsigma }_1,{\Omega }_0\frac{\partial }{\partial t}\theta {\varsigma }_2={\Omega }^{\ast }-2g_2{\varsigma }_m^2\theta {\varsigma }_2+f{\varsigma }_m\theta {\varsigma }_3,\\ {\Omega }_0\frac{\partial }{\partial t}\theta {\varsigma }_3={\Omega }^{\ast }-2g_2{\varsigma }_m^2\theta {\varsigma }_3+f{\varsigma }_m\theta {\varsigma }_2.\end{array}$$

The characteristic equation for stripe pattern stability is deduced from (58) as follows:$$\begin{array}{}\text{(59)}& {\Omega }^{\ast }-3g_1{\varsigma }_m^2-\lambda {{\Omega }^{\ast }-2g_2{\varsigma }_m^2-\lambda }^2-f^2{\varsigma }_m^2=0,\end{array}$$where $\lambda$ is the eigenvalue.

Substituting ${\varsigma }_m=\sqrt{{\Omega }^{\ast }/g_1}$ in (64), we get the following:$$\begin{array}{}\text{(60)}& -2{\Omega }^{\ast }-\lambda {{\Omega }^{\ast }-\frac{g_2}{g_1}{\Omega }^{\ast }-\lambda }^2-f^2\frac{{\Omega }^{\ast }}{g_1}=0.\end{array}$$

We have ${\lambda }_1=-2{\Omega }^{\ast }$ and ${\lambda }_{\mathrm{2,3}}={\Omega }^{\ast }1-g_1/g_2+-f\sqrt{{\Omega }^{\ast }/g_1}$.

If the eigenvalues are negative, the stripe patterns are stable as follows:$$\begin{array}{}\text{(61)}& {\Omega }^{\ast }>\frac{f^2{g}_1}{{g_1-g_2}^2}={\Omega }_m^{\ast }.\end{array}$$

Thus, the prerequisite for stability for stripe patterns is ${\Omega }^{\ast }>{\Omega }_s^{\ast }$ and unstable if ${\Omega }^{\ast }<{\Omega }_m^{\ast }$.

Theorem 8.

The hexagonal pattern is stable if ${\Omega }^{\ast }<{\Omega }_{f_2}^{\ast }$ and unstable if ${\Omega }^{\ast }>{\Omega }_{f_2}^{\ast }$, where ${\Omega }_{f_2}^{\ast }=f^{2}2g_1+g_2/{g_1-g_2}^2$.

Proof.

Substituting ${\varsigma }_1={\varsigma }_2={\varsigma }_3={\varsigma }_h$ in (57), we get the following:$$\begin{array}{}\text{(62)}& {\varsigma }_f^{\pm }=\frac{f\pm \sqrt{f^2+4g_1+2g_2{\Omega }^{\ast }}}{2g_1+2g_2}.\end{array}$$

On the basis of the real value, the value of ${\varsigma }_f^{\pm }$ discriminant must be positive and it becomes the following:$$\begin{array}{}\text{(63)}& {\Omega }^{\ast }>\frac{-f^2}{4g_1+2g_2}{\Omega }_{f_1}^{\ast }.\end{array}$$

Substitute ${\varsigma }_1={\varsigma }_m+\theta {\varsigma }_1\mathrm{a}\mathrm{n}\mathrm{d}{\varsigma }_2=\theta {\varsigma }_2,{\varsigma }_3=\theta {\varsigma }_3$ in (57) and replace the steady state condition. The following are the mathematical relationships:$$\begin{array}{c}\text{(64)}& {\Omega }_0\frac{\partial }{\partial t}\theta {\varsigma }_1=w_1\theta {\varsigma }_1+w_2\theta {\varsigma }_2+w_2\theta {\varsigma }_3,{\Omega }_0\frac{\partial }{\partial t}\theta {\varsigma }_2=w_2\theta {\varsigma }_1+w_1\theta {\varsigma }_2+w_2\theta {\varsigma }_3,\\ {\Omega }_0\frac{\partial }{\partial t}\theta {\varsigma }_3=w_2\theta {\varsigma }_1+w_2\theta {\varsigma }_2+w_1\theta {\varsigma }_3,\end{array}$$where $w_1={\Omega }^{\ast }-3g_1+2g_2{\varsigma }_f^2$ and $w_2=f{\varsigma }_f-2g_2{\varsigma }_f^2$.

From (64), the characteristic equation is as follows:$$\begin{array}{}\text{(65)}& {\lambda }^3-3w_1{\lambda }^2+3g_1^2-g_2^2\lambda +3g_1{g}_2^2-g_1^3-2g_2^3=0.\end{array}$$

After solving the above cubic equation, the following eigenvalues are obtained:$$\begin{array}{}\text{(66)}& {\lambda }_1=w_1+2w_2{\mathrm{a}\mathrm{n}\mathrm{d}\lambda }_{\mathrm{2,3}}=w_1-w_2.\end{array}$$

Let us examine eigenvalues’ stability conditions, say ${\varsigma }_f^{+}$ and ${\varsigma }_f^{-}$.

Case 1.

${\varsigma }_f={\varsigma }_f^{+}$:

Substituting ${\varsigma }_f^{+}$ from (62) in ${\lambda }_1$ defined as ${\lambda }_1=w_1+2w_2={\Omega }^{\ast }-3g_1+2g_2{\varsigma }_f^{+2}+2f{\varsigma }_f^{+}$, we get the following:$$\begin{array}{}\text{(67)}& {\lambda }_1=-\frac{4{\Omega }^{\ast }{g}_1+2g_2+f^{2}f\sqrt{f^2+4{\Omega }^{\ast }{g}_1+2g_2}}{g_1+2g_2}<0.\end{array}$$

Assume ${\lambda }_2={\lambda }_3<0$.$$\begin{array}{c}\text{(68)}& {\lambda }_2={\lambda }_3={\Omega }^{\ast }-3g_1{\varsigma }_f^{+2}-f{\varsigma }_f^{+}=\frac{4{\Omega }^{\ast }{g_1+2g_2}^2-6g_1{f}^2-12g_1{\Omega }^{\ast }{g}_1+2g_2-2h^2{g}_1+2g_2-4f2g_1+g_2\sqrt{g_1+2g_2}}{4{g_1+2g_2}^2}<0.\end{array}$$

On simplification, we get the following equation:$$\begin{array}{}\text{(69)}& {\Omega }^{\ast }=\frac{f^{2}2g_1+g_2}{{g_1-g_2}^2}.\end{array}$$

Thus, the hexagonal pattern is stable if ${\Omega }^{\ast }<{\Omega }_{f_2}^{\ast }$.

Case 2.

${\varsigma }_f={\varsigma }_f^{-}$:

Substitute ${\varsigma }_f^{-}$ from (62) in ${\lambda }_1$ defined as ${\lambda }_1=w_1+2w_2={\Omega }^{\ast }-3g_1+2g_2{\varsigma }_f^{-2}+2f{\varsigma }_f^{-}$, we get the following:$$\begin{array}{}\text{(70)}& {\lambda }_1=-\frac{4{\Omega }^{\ast }{g}_1+2g_2-f^2+f\sqrt{f^2+4{\Omega }^{\ast }{g}_1+2g_2}}{g_1+2g_2}>0.\end{array}$$

On comparing (67) and (70), ${\lambda }_1$ is positive. Now, calculating ${\lambda }_2$ and ${\lambda }_3$, we get the following:$$\begin{array}{c}\text{(71)}& {\lambda }_2={\lambda }_3={\Omega }^{\ast }-3g_1{\varsigma }_f^{-2}-f{\varsigma }_f^{-}=\frac{4{\Omega }^{\ast }{g_1+2g_2}^2-6g_1{f}^2-12g_1{\Omega }^{\ast }{g}_1+2g_2-2h^2{g}_1+2g_2+4f2g_1+g_2\sqrt{g_1+2g_2}}{4{g_1+2g_2}^2}<0.\end{array}$$

Again comparing (68) and (71), we get (71) as positive.

As a result, when ${\varsigma }_f={\varsigma }_f^{-}$, all eigenvalues turn positive, which assures the instability of hexagonal patterns.

Now, the system can be simulated numerically in a 350 $\times$ 350 2D square lattice.

6. Numerical Computation

Analytical studies are never completed unless the results are numerically validated. System (2) is simulated by the ODE solver, and the results are displayed with parameter values. Various numerical results are shown to validate the analytical stability analysis presented in the preceding sections. The analyses have been carried out with positive parameter values. As shown in the images below, we use a variety of parameter values to acquire a better understanding of the dynamics of system (2). The numbers specified for the parameters are not based on field data, and they are purely hypothetical parameters meant to depict the system’s dynamic behaviour.

The density plot graph depicts the interplay of prey and predator populations at various geographical regions. Overlapping density plots of both species allow for the identification of locations with more strong predator-prey interactions. Places with high predator density and low prey density may indicate significant predation pressure, whereas places having high prey density and low predator density may indicate prey refuge or favourable breeding grounds.

Let us examine model (2) with the following parameter values $g_1=1.3,{\mathbb{K}}_1=5,\sigma =0.027,\nu =0.08,g_2=0.856,{\mathbb{K}}_2=5,\phi =0.002,\eta =0.09,{\mu }_1=0.1,{\mu }_2=12,\text{and\hspace{0.17em}}\mathrm{\varrho }=1.2$. From Figures 1 and 2, it is clear that both populations survive for the long term. In addition, Figure 1 demonstrates that the periodic solution emerges for the first prey, second prey, and predator. Consider the same set of parameters as above and slightly change the values of the intrinsic growth rate such that $g_1=1.5$ and $g_2=0.5$. From Figure 2, it is noticeable that all three species oscillate at a certain time interval then it becomes stable. Figures 1 and 2 are qualitatively equivalent. The outcomes of the remaining instances with instances II and III parameters are qualitatively equivalent. The sole distinction is the duration of the fluctuation period. It appears that being close to steady states is critical, even if the steady state is unstable. Only then can the species cohabit in an oscillating fashion and prevent a large fluctuation (or blooming). It is concluded that prey 1, prey 2, and predator populations coexist simultaneously.

[figure(s) omitted; refer to PDF]

Now, we assume the following parameters $c_1=1.5,{\mathbb{K}}_1=15,\alpha =0.027,\gamma =0.08,c_2=1.9,{\mathbb{K}}_2=15,\beta =0.002,\delta =0.09,w_1=0.1,w_2=1.2,and\rho =1.2$. Figure 3, shows that the population persists and are stable.

[figure(s) omitted; refer to PDF]

For Figures 4 and 5, let the initial conditions be $x_0=1.2,y_0=1.8,{andz}_0=2.4$, and all other constants are 1. From Figure 4, it is noticed that even when ${\mathbb{K}}_1$ is set to a peak value, prey 1 persists at zero level. That is, the prey population vanishes for a large value of ${\mathbb{K}}_1$, while the other two populations oscillate at time $t$. Similarly, Figure 5 shows that for the carrying capacity ${\mathbb{K}}_1=20$, prey 1 remains at zero level and vanishes. A periodic solution occurs between prey 2 and the predator. The dynamics of the predator species hunting the prey species are also visible in Figure 5(c). The solution consists of multiple layers within the support and generates periodic patterns.

[figure(s) omitted; refer to PDF]

Now, if we decrease the value of the conversion coefficient ${\mu }_2$, it leads to the extinction of the predator, while the prey population remains constant, which is shown in Figure 6. Similarly, from Figure 7, when the growth rate is increased, the prey population disappears while the predator population endures and remains stable. Only circular patterns exist within the solution’s support in Figure 7(c). It is also possible to observe that the predator population is chasing the prey population which leads to the extinction of prey.

[figure(s) omitted; refer to PDF]

7. Pattern Formation

This section examines simulations of system (2) computationally, proving that diffusion creates spatial patterns. Turing patterns can appear in ecological systems as spatially different zones of high and low population densities, resulting in elaborate patterns that resemble stripes, dots, or other complex structures. The appearance of Turing patterns in prey-predator models emphasises the relevance of spatial dynamics and reveals a more complete knowledge of ecological interactions. The emergence of patterns for specific parameter values and diffusion coefficients of (19) is studied, demonstrating how diffusion can cause a stable steady state to become unstable. The simulation proceeds by randomly or predefined distributing predator and prey populations over the geographical region. The system of differential equations is then numerically solved, which advances the populations in time. The diffusion parameters contribute to the expansion of populations at each time step, allowing predators and prey to travel and distribute throughout the domain.

Patterns develop in the population distributions as the experiment proceeds. Depending on the model’s unique characteristics, these patterns might include spatial clusters, waves, or spiral structures. Patterns occur as a result of the intricate interplay between predator-prey interactions, diffusion processes, and the system’s underlying nonlinear dynamics. Contour plots, showing the number of predators and prey at different places, can be used to visualise the simulated patterns. This visualisation gives insights on the population’s geographical organisation and structure throughout time.

We provide here numerical demonstrations to further clarify the prior theoretical results using various sets of parameters. We analyse model (2) by considering the following values, $\begin{array}{c}{g}_1=0.002,{\mathbb{K}}_1=1,\sigma =0.005,\nu =0.001,g_2=0.0001,\\ {\mathbb{K}}_2=1,\phi =0.01,\eta =0.02,{\mu }_1=0.001,{\mu }_2=0.0001,\text{and\hspace{0.17em}}\\ \mathrm{\varrho }=0.00011\end{array}$. When the diffusion parameters $d_1$, $d_2$, and $d_3$ are increased from 0.01 to 50, four basic dynamics such as chaos, intermittent chaos, smooth oscillatory state, and stationary behaviour emerge which are shown in Figures 8– Figure 10.

[figure(s) omitted; refer to PDF]

Throughout the prey-predator concept, a pattern is formed when both prey and predator diffuse at different times in three-dimensional space. Figure 11 depicts the coexistence of spot and stripe patterns at intervals of 0, 0.1, 0.2, and 0.8.

[figure(s) omitted; refer to PDF]

Figure 12 illustrates that predator death rates influence the pattern of spatial distribution. We get a pattern for $t=500$ and choose diffusion coefficients $d_1$ = 5 and $d_2$ = 0.2. There are blue and brown spot patterns. In Figure 12(a), blue spots indicate high population density and brown spots indicate low population density, whereas in Figure 12(b), both populations are equal.

[figure(s) omitted; refer to PDF]

We used Turing patterns in the spatial domain for interpreting system (2). As time passes, the resulting pattern structure will tend to stabilize. The wave number modifies the structure of Turing patterns in homotopy series solutions, which is sensitive to the starting solution. In short, as time goes by in the range of Turing pattern parameters, the system will produce striped patterns, spot patterns, and both striped and spot patterns at the same time.

8. Conclusion

In the present work, we have constructed and explored a predator-prey model characterized by a system of differential equations incorporating reaction and diffusion. Our investigation encompassed an in-depth examination of Turing patterns and spatial pattern formation within the model. We conducted an extensive numerical analysis of the Turing system, investigating how pattern generation is affected by varying parameters. The outcomes included the creation of intricate spatiotemporal patterns, showcasing the spatial complexity attainable in reaction-diffusion systems, including chaotic patterns. We further introduced an innovative sinc function interpolation method applicable to three-species predator-prey systems (PPS) with intricate behaviours. In addition, we have used Taylor’s expansion method to obtain the amplitude equation for the reaction-diffusion system, which accurately described Turing patterns near the bifurcation point. Through simulation results, we successfully demonstrated the effectiveness of this novel approach.

The investigation of Turing patterns in diffusive prey-predator models provides important insights into the spatial organisation and stability of ecological systems. Researchers can get a better grasp of the variables impacting species coexistence, population dynamics, and ecological stability by studying the underlying mechanisms and circumstances that give birth to these patterns. Furthermore, these discoveries have practical applications in sectors such as conservation biology, pest control, and biodiversity protection, where knowing spatial dynamics is critical for making ethical choices. Thus, studying Turing patterns in diffusive prey-predator models is an important area of study in ecological systems. These models improve our knowledge of complicated ecological processes and contribute to the development of sustainable management methods by accounting for spatial dynamics. Further research into the mechanics and consequences of Turing patterns in predator-prey dynamics will surely contribute to a better understanding of the complex interplay between species interactions and geographical variability.

Acknowledgments

The authors would like to thank the Deanship of Scientific Research, Qassim University for funding publication of this project.

References

[1] N. B. Sharmila, C. Gunasundari, "Mathematical analysis of prey predator models with Holling type I functional responses and time delay," Commun. Math. Biol. Neurosci., vol. 2023, 2023.

[2] B. Dubey, "A prey-predator model with a reserved area," Nonlinear Analysis Modelling and Control, vol. 12 no. 4, pp. 479-494, DOI: 10.15388/NA.2007.12.4.14679, 2007.

[3] C. M. Roberts, N. V. Polunin, "Marine Reserves: Simple Solutions to Managing Complex Fisheries," Ambio, vol. 22 no. 6, pp. 363-368, 1993.

[4] M. S. Arif, K. Abodayeh, A. Ejaz, "On the stability of the diffusive and non-diffusive predator-prey system with consuming resources and disease in prey species," Mathematical Biosciences and Engineering, vol. 20 no. 3, pp. 5066-5093, DOI: 10.3934/mbe.2023235, 2023.

[5] N. B. Sharmila, C. Gunasundari, "Travelling wave solutions for a diffusive prey-predator model with one predator and two preys," International Journal of Apllied Mathematics, vol. 35 no. 5, pp. 46-67, DOI: 10.12732/ijam.v35i5.3, 2022.

[6] N. B. Sharmila, C. Gunasundari, "Stability analysis of a fractional order prey-predator model with disease in prey," Mathematica Applicanda, vol. 50 no. 2, pp. 287-302, DOI: 10.14708/ma.v50i2.7142, 2023.

[7] M. Aakash, C. Gunasundari, "Effect of partially and fully vaccinated individuals in some regions of India: a mathematical study on COVID-19 outbreak," Communications in Mathematical Biology and Neuroscience, vol. 2023, 2023.

[8] G. Santhosh Kumar, C. Gunasundari, "Turing instability of a diffusive predator-prey model along with an Allee effect on a predator," Communications in Mathematical Biology and Neuroscience, vol. 2022, 2022.

[9] M. Aakash, C. Gunasundari, Q. M. Al-Mdallal, "Mathematical modeling and simulation of SEIR model for COVID-19 outbreak: a case study of Trivandrum," Frontiers in Applied Mathematics and Statistics, vol. 9,DOI: 10.3389/fams.2023.1124897, 2023.

[10] M. Shoaib Arif, K. Abodayeh, A. Ejaz, "Computational modeling of reaction-diffusion COVID-19 model having isolated compartment," Computer Modeling in Engineering and Sciences, vol. 135 no. 2, pp. 1719-1743, DOI: 10.32604/cmes.2022.022235, 2023.

[11] E. Phan Van Long, "Sufficient condition for synchronization in complete networks of reaction-diffusion equations of hindmarsh-rose type with linear coupling," IAENG International Journal of Applied Mathematics, vol. 52 no. 2, pp. 315-319, 2022.

[12] P. Zheng, "On a two-species competitive predator-prey system with density-dependent diffusion," Mathematical Biosciences and Engineering, vol. 19 no. 12, pp. 13421-13457, DOI: 10.3934/mbe.2022628, 2022.

[13] A. Ejaz, Y. Nawaz, M. Shoaib Arif, D. S Mashat, K. Abodayeh, "Stability analysis of predator-prey system with consuming resource and disease in predator species," Computer Modeling in Engineering and Sciences, vol. 132 no. 2, pp. 489-506, DOI: 10.32604/cmes.2022.019440, 2022.

[14] M. S. Arif, K. Abodayeh, A. Ejaz, "Stability analysis of fractional-order predator-prey system with consuming food resource," Axioms, vol. 12 no. 1,DOI: 10.3390/axioms12010064, 2023.

[15] L. N. Guin, P. K. Mandal, "Effect of prey refuge on spatiotemporal dynamics of the reaction–diffusion system," Computers & Mathematics with Applications, vol. 68 no. 10, pp. 1325-1340, DOI: 10.1016/j.camwa.2014.08.025, 2014.

[16] L. N. Guin, B. Mondal, S. Chakravarty, "Existence of spatiotemporal patterns in the reaction-diffusion predator-prey model incorporating prey refuge," International Journal of Biomathematics, vol. 09 no. 06,DOI: 10.1142/s1793524516500856, 2016.

[17] L. N. Guin, S. Acharya, "Dynamic behaviour of a reaction-diffusion predator-prey model with both refuge and harvesting," Nonlinear Dynamics, vol. 88 no. 2, pp. 1501-1533, DOI: 10.1007/s11071-016-3326-8, 2017.

[18] R. Han, L. N. Guin, B. Dai, "Consequences of refuge and diffusion in a spatiotemporal predator-prey model," Nonlinear Analysis: Real World Applications, vol. 60,DOI: 10.1016/j.nonrwa.2021.103311, 2021.

[19] L. N. Guin, S. Djilali, S. Chakravarty, "Cross-diffusion-driven instability in an interacting species model with prey refuge," Chaos, Solitons & Fractals, vol. 153 no. 1,DOI: 10.1016/j.chaos.2021.111501, 2021.

[20] R. Han, L. N. Guin, S. Acharya, "Complex dynamics in a reaction-cross-diffusion model with refuge depending on predator-prey encounters," The European Physical Journal Plus, vol. 137 no. 1,DOI: 10.1140/epjp/s13360-022-02358-7, 2022.

[21] R. Han, G. Mandal, L. N. Guin, S. Chakravarty, "Dynamical response of a reaction-diffusion predator-prey system with cooperative hunting and prey refuge," Journal of Statistical Mechanics: Theory and Experiment, vol. 2022 no. 10,DOI: 10.1088/1742-5468/ac946d, 2022.

[22] L. N. Guin, R. Murmu, H. Baek, K.-H. Kim, "Dynamical analysis of a Beddington-DeAngelis interacting species system with harvesting," Mathematical Problems in Engineering, vol. 2020,DOI: 10.1155/2020/7596394, 2020.

[23] L. N. Guin, P. J. Pal, J. Alzahrani, N. Ali, K. Sarkar, S. Djilali, A. Zeb, I. Khan, S. M. Eldin, "Influence of Allee effect on the spatiotemporal behavior of a diffusive predator-prey model with Crowley-Martin type response function," Scientific Reports, vol. 13, pp. 4710-4719, DOI: 10.1038/s41598-023-28419-0, 2023.

[24] A. M. Turing, "The chemical basis of morphogenesis," Philosophical Transactions of the Royal Society of London,Series A B, vol. 237 no. 641, pp. 37-72, 1952.

[25] L. A. Segel, J. L. Jackson, "Dissipative structure: an explanation and an ecological example," Journal of Theoretical Biology, vol. 37 no. 3, pp. 545-559, DOI: 10.1016/0022-5193(72)90090-2, 1972.

[26] A. Gierer, H. Meinhardt, "A theory of biological pattern formation," Kybernetik, vol. 12 no. 1, pp. 30-39, DOI: 10.1007/bf00289234, 1972.

[27] S. A. Levin, L. A. Segel, "Hypothesis for origin of planktonic patchiness," Nature, vol. 259 no. 5545,DOI: 10.1038/259659a0, 1976.

[28] K. A. J. White, C. A. Gilligan, "Spatial heterogeneity in three-species, plant-parasite-hyperparasite, systems," Philosophical Transactions of the Royal Society of London- Series B: Biological Sciences, vol. 353 no. 1368, pp. 543-557, DOI: 10.1098/rstb.1998.0226, 1998.

[29] W. Chen, R. Peng, "Stationary patterns created by cross-diffusion for the competitor-competitor-mutualist model," Journal of Mathematical Analysis and Applications, vol. 291 no. 2, pp. 550-564, DOI: 10.1016/j.jmaa.2003.11.015, 2004.

[30] N. A. Mbroh, J. B. Munyakazi, "A fitted operator finite difference method of lines for singularly perturbed parabolic convection-diffusion problems," Mathematics and Computers in Simulation, vol. 165, pp. 156-171, DOI: 10.1016/j.matcom.2019.03.007, 2019.

[31] N. Dhiman, M. Tamsir, "A collocation technique based on modified form of trigonometric cubic B-spline basis functions for Fisher’s reaction–diffusion equation," Multidiscipline Modeling in Materials and Structures, vol. 14 no. 5, pp. 923-939, DOI: 10.1108/mmms-12-2017-0150, 2018.

[32] C. M. Elliott, T. Ranner, "Evolving surface finite element method for the Cahn–Hilliard equation," Numerische Mathematik, vol. 129 no. 3, pp. 483-534, DOI: 10.1007/s00211-014-0644-y, 2015.

[33] H. Wu, X. F. Han, "Discontinuous Galerkin spectral element method for a class of nonlinear reaction–diffusion equations," Journal of Shanxi Normal University (Philosophy and Social Sciences edition), vol. 20, 2014.

[34] J. H. He, "Homotopy perturbation method: a new nonlinear analytical technique," Applied Mathematics and Computation, vol. 135 no. 1, pp. 73-79, DOI: 10.1016/s0096-3003(01)00312-5, 2003.

[35] M. Ivan Azis, "Numerical solution for unsteady anisotropic diffusion-convection-reaction equation of exponentially varying coefficients and compressible flow," IAENG International Journal of Applied Mathematics, vol. 52 no. 2, pp. 327-335, 2022.

[36] J. Lund, K. L. Bowers, Sinc Methods for Quadrature and Differential Equations, 1992.

[37] F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, 1993.

[38] D. Dai, X. Lv, Y. Wang, "Numerical simulation for a class of predator–prey system with homogeneous Neumann boundary condition based on a sinc function interpolation method," Boundary Value Problems, vol. 2020 no. 1,DOI: 10.1186/s13661-020-01402-8, 2020.

[39] X. Zhang, Efficient Solution of Differential Equation Based on MATLAB: Spectral Method Principle and Implementation, 2015.

[40] Q. Zheng, J. Shen, "Pattern formation in the FitzHugh–nagumo model," Computers & Mathematics with Applications, vol. 70 no. 5, pp. 1082-1097, DOI: 10.1016/j.camwa.2015.06.031, 2015.