Research on pattern dynamics of a class of predator-prey model with interval biological coefficients for capture

1.   Introduction

With the development of human civilization, maintaining ecological balance has become a major challenge facing mankind. The predator-prey model is a mathematical model for studying the interaction between predators and prey in ecosystems ^[1,2]. It can be used to analyze the impact of environmental changes or the introduction of new species. At present, about $28\%$ of the assessed species in the world are at risk of extinction ^[3]. In response to this situation, the implementation of a reasonable policy for the development of biological resources can protect populations from possible extinction ^[4].

Yan ^[5] studied the global stability of a delayed diffusion predator-prey model with Michaelis-Menten-type prey capture. Ou ^[6] obtained the parameter conditions for the stability and bifurcation of two delayed predator-prey systems. Yao ^[7] studied the dynamics of a Leslie-Gower predator-prey system with proportionally dependent Holling IV functional response and a constant prey capture rate. Zhang ^[8] applied a delayed predator-prey model with non-constant mortality and a constant prey capture rate to investigate the effect of time delay on equilibrium stability. Chen ^[9] studied the pattern dynamics of a harvested predator-prey model with no-flux boundary conditions. Cui ^[10] proposed a new Lotka-Volterra commensal symbiosis system accompanying delay. Therefore, motivated by the above discussion, we consider the following predator-prey model with a harvesting term:

where $u\left(x, t \right)$ and $v\left(x, t \right)$ represent the densities of the prey and predators at location $x$ and time $t$, respectively. $\varDelta$ is the Laplacian operator. $d_1$ and $d_2$ represent the diffusion rates of prey $u$ and predator $v$, respectively. $a$ and $\beta$ are two positive constants. $bu$ implies harvesting of the prey population, and $\alpha v$ denotes harvesting of the predator population.

It can be seen from the existing literature that the parameters in most biological models are considered to be accurate. However, due to factors such as climate change, natural disasters, and deforestation, most measurement processes and initial data may have errors ^[11]. Therefore, models with imprecise parameters are more realistic. A model with imprecise biological parameters can be reflected by random, fuzzy, and interval methods. In fuzzy methods, imprecise parameters are represented in the form of fuzzy sets or fuzzy numbers. Nowadays there are many studies on fuzzy parameters of dynamical systems ^[12,13]. In the stochastic method, imprecise parameters are represented by random variables with appropriate probability distribution ^[14,15]. Different from random and fuzzy methods, Pal ^[16] first proposed the concept of an interval to describe the imprecise parameters of ecological models. This method is easier and more effective than the first two methods. Since then, many scholars have used an interval method to study ecological models with inaccurate parameters. Pal uses the exponential form of interval parameters, and Ramezanadeh uses the linear form of interval parameters ^[17,18].

Fractional derivatives have memory and nonlocality. In biological systems, memory refers to the ability of the system to retain the information of past events and use it to influence future behavior ^{[19,20,21,22]}. Due to the fact that most populations have long-term memory ^[23,24], integer order population systems ignore the influence of memory. Therefore, we consider a fractional order predator-prey model. There are many ways to solve fractional differential equations, for example, the fractional reduced differential transform method ^[25], the Adomian decomposition sumudu transform method ^[26], the Fourier spectral method ^{[27,28,29,30]}, and the piecewise reproducing kernel method ^[31]. In this paper, the Euler discrete method is used. It can describe the ecological process of the reaction-diffusion equation on a long time scale and has the advantage of fast calculation speed.

The main contributions of this paper are as follows:

$1)$ Due to other factors, most measurement processes and initial data may have errors. In this paper, a more realistic interval parameter model is used.

$2)$ The Turing pattern is theoretically classified and verified by numerical methods. Some different results are obtained.

$3)$ We numerically prove that the dispersal rate of the prey population will suppress the spatiotemporal chaos of the model.

The rest of this article is organized as follows: In Section 2, we establish a model and explore the positivity and uniqueness of solutions for models without diffusion terms. In Section 3, we discuss the stability of the model and Hopf bifurcation. In Section 4, the Turing instability of the model is discussed. In Section 5, weak nonlinear analysis is used to derive the amplitude equation. In Section 6, we conduct numerical simulations. Conclusions are given in Section 7.

2.   Model and preliminaries

2.1.   Model formation

To establish our model, we introduce the following two definitions.

Definition 2.1. (^[16,18,32]) An interval number $A$ is defined by $A = \left[ \acute{a}, \grave{a} \right] = \left\{ y|\acute{a}\le y\le \grave{a}, y\in \mathbb{R} \right\}$. Moreover, each real number $a\in \mathbb{R}$ can be represented by $[a, a]$. Let $A = \left[ \acute{a}, \grave{a} \right]$ and $B = \left[ \acute{a}, \grave{a} \right]$. Define the following algorithms:

$1) \; A+B = \left[ \acute{a}, \grave{a} \right] +[ \acute{b}, \grave{b}] = [\acute{a}+\acute{b}, \grave{a}+\grave{b} ]$ for $\acute{a}+\acute{b} > 0$;

$2) \; A-B = \left[ \acute{a}, \grave{a} \right] -[ \acute{b}, \grave{b}] = [\acute{a}-\acute{b}, \grave{a}-\grave{b} ]$ for $\acute{a}-\acute{b} > 0$;

$3) \; \rho A = \rho \left[ \acute{a}, \grave{a} \right] = \left[ \rho \acute{a}, \rho \grave{a} \right]$ for $\rho \ge 0$;

$4) \; \rho A = \rho \left[ \acute{a}, \grave{a} \right] = \left[ \rho \grave{a}, \rho \acute{a} \right]$ for $\rho < 0$.

Definition 2.2. (^[16,18,32]) Let $a > 0$ and $b > 0$ have interval $[a, b]$. The interval-valued function is presented as $f\left(r \right) = a^{1-r}b^r$ for $r\in \left[ 0, 1 \right]$.

Some authors have found that the dynamic behavior of a fractional model ^{[33,34,35,36]} is much more complex than that of the corresponding integer model. According to Definition 2.1, we consider the following fractional predator-prey model with interval biological coefficient

where $\acute{a} > 0, \acute{b} > 0, \acute{\alpha} > 0$, and $\acute{\beta} > 0$. Now, there are many types of fractional derivatives ^{[34,35,36,37,38]}. Here, $D_{t}^{\alpha}$ represents Caputo fractional differentiation. It has the advantage of relatively simple calculation, which is defined as follows:

Using Definition 2.2, we get

where $0\le r\le 1$.

2.2.   Positivity and uniqueness

In this section, we prove the positivity and uniqueness of the solution of the fractional order model without a diffusion term.

The non-diffusion version of model (2.3) is as follows:

Theorem 2.3. All solutions of the model system (2.4) are nonnegative.

Proof. For a similar proof of Theorem (2.3), readers are referred to ^[19]. □

Theorem 2.4. The fractional system (2.4) has a unique solution under any nonnegative initial conditions.

Proof. According to the method proposed in ^[39,40,41], we define the following operator:

Let

with

Using the Banach fixed point theorem, we can obtain the uniform norm:

The Picardis operator is as follows:

This is defined as:

where $X\left(t \right) = \left[ u\left(t \right), v\left(t \right) \right] ^T, X_0\left(t \right) = \left[ u_0\left(t \right), v_0\left(t \right) \right] ^T$, and $F\left(t, X\left(t \right) \right) = \left[ f_1\left(t, u \right), f_2\left(t, v \right) \right] ^T$.

We assume that the solution of the model is bounded in a time period:

We can get:

with $N = \max \left\{ N_1, N_2 \right\}, b = \max \left\{ b_1, b_2 \right\}$, and $a < \frac{b}{N}$.

Since $F$ is a contraction and $\beta < 1$, we obtain $a\beta < 1$, that is, the defined operator $O$ is also a contraction. Therefore, the uniqueness proof of the system solution is complete. □

3.   Stability and Hopf bifurcation analysis

Here, we determine the equilibrium point of the system (2.4). By analyzing the stability of the equilibrium point and the Hopf bifurcation, the conditions under which different states of the system appear are given.

3.1.   Stability analysis

Before determining the stability of the equilibrium point, we first give the stability criterion of the fractional differential system.

Theorem 3.1. (^[42,43]) Consider a fractional differential system

Let $x_*$ be an equilibrium point and the $\lambda _i, \left(i = 1, 2, \cdots, n \right)$ are eigenvalues of Jacobian matrix $J = \frac{\partial f}{\partial x_*}$.

$1)$ The equilibrium point $x_*$ is asymptotically stable if and only if

$2)$ The equilibrium point $x_*$ is stable if and only if

$3)$ The equilibrium point $x_*$ is unstable if and only if

Definition 3.2. (^[44]) The roots of the equation $f\left(t, x\left(t \right) \right) = 0$ are called the equilibria of fractional differential system

where $x\left(t \right) = \left(x_1\left(t \right), x_2\left(t \right), \cdots, \; x_n\left(t \right) \right) ^T\in \mathbb{R}^n, \; f\left(t, x\left(t \right) \right) \in R^n$, and $D_{t}^{\alpha}x\left(t \right) = \left(D_{t}^{\alpha_1}x_1\left(t \right), \; D_{t}^{\alpha_2}x_2\left(t \right), \cdots, \; D_{t}^{\alpha_n}x_n\left(t \right) \right) ^T, \; \alpha _i\in \mathbb{R}^+, i = 1, 2, \cdots, n$.

From a biological perspective, we are only interested in the positive equilibrium point. Obtain the equilibrium point by solving the following system of equations.

Denote by $f\left(u, v \right) = 0$ and $g\left(u, v \right) = 0$. The positive equilibrium point $E_* = \left(u_*, v_* \right)$ is obtained, where

and

We can obtain the Jacobi matrix for system (3.6) at the equilibrium point $E_*$ as follows:

As such, the characteristic equation at equilibrium point $E_*$ is read as follows:

where

Roots of the characteristic equations are

Through Theorem (3.1), we draw the following conclusions:

Theorem 3.3. (^[45]) The stability of equilibrium point $E_*$ is determined by $tr_0, \; det _0$, and $\alpha$.

If $tr_{0}^{2}-4\det _0\ge 0$, then,

$1)$ The equilibrium point $E_*$ is asymptotically stable if and only if $tr_0\leq0$ and $det_0 > 0$;

$2)$ The equilibrium point $E_*$ is unstable if and only if $tr_0 > 0$ or $det_0 < 0$.

If $tr_{0}^{2}-4\det _0 < 0$, then,

$1)$ The equilibrium point $E_*$ is stable if and only if $\alpha \frac{\pi}{2} < \Big{|}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) \Big{|}$;

$2)$ The equilibrium point $E_*$ is unstable if and only if $\alpha \frac{\pi}{2} > \Big{|}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) \Big{|}$.

Proof. The eigenvalues are real when $tr_{0}^{2}-4\det _0\ge 0$.

For $tr_0 = 0$, $\lambda _{1, 2} = \pm i\sqrt{\det _0}$ is obtained, therefore $\left| arg\left(\lambda _{1, 2} \right) \right| = \frac{\pi}{2} > \alpha \frac{\pi}{2}$ implies $E_*$ is asymptotically stable; the eigenvalues are negative real when $tr_0 < 0$ and $det_0 > 0$, so $\left| arg\left(\lambda _{1, 2} \right) \right| = \pi > \alpha \frac{\pi}{2}$ implies $E_*$ is asymptotically stable.

For $tr_0 > 0$ and $det_0 > 0$, both the eigenvalues are positive real, hence $\left| arg\left(\lambda _{1, 2} \right) \right| = 0 < \alpha \frac{\pi}{2}$ implies $E_*$ is unstable; when $det_0 < 0$, the two eigenvalues are real numbers with opposite signs, so there exists $\left| arg\left(\lambda _1 \right) \right| = 0 < \alpha \frac{\pi}{2}$ which implies that $E_*$ is unstable.

The two eigenvalues are now complex conjugate when $tr_{0}^{2}-4\det _0 < 0$.

Therefore, $E_*$ is stable if $\alpha \frac{\pi}{2} < \Big{|}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right)\Big{|}$ and is unstable for $\alpha \frac{\pi}{2} > \Big{|}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) \Big{|}$. □

3.2.   Hopf bifurcation analysis

When $tr_0 = 0$ and $det_0 > 0$, the system (2.4) with $\alpha = 1$ loses stability through Hopf bifurcation. Since the stability of system (2.4) is affected by the fractional derivative, the fractional derivative can be regarded as a parameter of Hopf bifurcation. In the following, we establish the conditions for the Hopf bifurcation of system (3.1) around $E_*$ at parameter $\alpha = \alpha_h$ ^[19,46]:

$1)$ The Jacobian matrix at the equilibrium point $E_*$ has a pair of complex conjugate eigenvalues $\lambda _{1, 2} = a_i+ib_i$ which become purely imaginary at $\alpha = \alpha_h$;

$2) \; m\left(\alpha _h \right) = 0$ where $m\left(\alpha \right) = \alpha \frac{\pi}{2}-\underset{1\le i\le 2}{\min}|arg\left(\lambda _i \right) |$;

$3) \; \left. \frac{\partial m\left(\alpha \right)}{\partial \alpha} \right|_{\alpha = \alpha _h}\ne 0$.

Now, we prove that $E_*$ has Hopf bifurcation when $\alpha$ goes through $\alpha_h$.

Theorem 3.4. Suppose that the equilibrium point $E_*$ is unstable when $tr_{0}^{2}-4\det _0 < 0$ and $tr_0 > 0$. The fractional parameter $\alpha$ passes through the critical value $\alpha_h$, and the system (2.4) undergoes Hopf bifurcation near $E_*$, where

Proof. For $tr_{0}^{2}-4\det _0 < 0$ and $tr_0 > 0$, the eigenvalues are complex conjugates with a positive real part. Hence,

and $\alpha \frac{\pi}{2} > \Big{|}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) \Big{|}$ for some $\alpha$. Let $\alpha_h \frac{\pi}{2} = \Big{|}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right) \Big{|}$, get $\alpha _h = \frac{2}{\pi}\tan ^{-1}\left(\frac{\sqrt{4\det _0-tr_{0}^{2}}}{tr_0} \right)$. Moreover, $\left. \frac{\partial m\left(\alpha \right)}{\partial \alpha} \right|_{\alpha = \alpha _h} = \frac{\pi}{2}\ne 0$. Therefore, all Hopf conditions satisfy. □

4.   Turing instability

In this section, we present the Turing instability condition of the system (2.3) at $\alpha = 1$.

Perturbate the equilibrium point with $u = u_*+\tilde{u}, v = v_*+\tilde{v}$, substitute it into system (2.3), expand it through the Taylor series, and remove higher-order terms to obtain the linear perturbation equation

where

and $J$ is a Jacobian matrix about $E_*$. For convenience, we still denote $\tilde{u}$ and $\tilde{v}$ as $u$ and $v$.

Expanding the perturbation variables in Fourier space and substituting $U = \left(\begin{array}{c} c_{k}^{1}\\ c_{k}^{2}\\ \end{array} \right) e^{\lambda t+ikr}$ into the perturbation Eq (4.1) yields the characteristic equation

where $\lambda$ is the growth rate, $k$ is the wave number, $r$ is the spatial vector, and $c_{k}^{1}, c_{k}^{2}$ are constants.

Solve characteristic Eq (4.3), and obtain the following dispersion relationship:

where

The solution of characteristic Eq (4.4) is in the following form:

In order to explore the existence conditions of Turing instability at $k\neq 0$, we should ensure that $tr_k < 0$ and $det _k < 0$. When $\acute{a}^r\grave{a}^{1-r} > \acute{\beta}^{1-r}\grave{\beta}^r$, $tr_k < 0$ is easy to satisfy. In order to ensure the occurrence of $det _k < 0$, the condition of marginal stability $\min \left(\det _{k_{c}^{2}} \right) = 0$ should be satisfied. Here $k_{c}^{2} = \frac{a_{11}d_2+a_{22}d_1}{2d_1d_2}$ is the minimum value of $det_k$ with respect to $k_{c}^{2}$.

From $\min \left(\det _{k_{c}^{2}} \right) = 0$, we can obtain:

Let $H\left(d_1 \right) = a_{22}^{2}d_{1}^{2}+2d_2\left(a_{11}a_{22}-2det_0 \right) d_1+a_{11}^{2}d_{2}^{2}$. If there are $H\left(d_1 \right) = 0$ and $H\left(d_1 \right) = -\min \left(\det _{k_{c}^{2}} \right)$, then there must be two roots, $d_{1}^{+}$ and $d_{1}^{-}$, where

Theorem 4.1. Suppose $0\le r\le 1, d_1 > 0, d_2 > 0$, and $\acute{a}^r\grave{a}^{1-r} > \acute{\beta}^{1-r}\grave{\beta}^r$ are valid.

$1)$ The equilibrium point $E_*$ is asymptotically stable if and only if $d_{+}^{-} < d_1 < d_{1}^{+}$;

$2)$ The equilibrium point $E_*$ is unstable if and only if $d_1 > d_{+}^{-}$ or $d_1 < d_{+}^{-}$;

$3)$ Turing bifurcation occurs at $d_1 = d_{+}^{-}$ or $d_1 = d_{+}^{-}$, and the critical wave number is $k_{c}^{2} = \sqrt{\frac{\det _0}{d_{1}^{+}d_2}}$ or $k_{c}^{2} = \sqrt{\frac{\det _0}{d_{1}^{-}d_2}}$.

Proof. The eigenvalues are negative real when $d_{+}^{-} < d_1 < d_{1}^{+}$, so $\left| arg\left(\lambda _k \right) \right| = \pi > \alpha \frac{\pi}{2}$ implies $E_*$ is asymptotically stable. When $d_1 > d_{+}^{-}$ or $d_1 < d_{+}^{-}$, the two eigenvalues are real numbers with opposite signs, so there exists $\left| arg\left(\lambda _k \right) \right| = 0 < \alpha \frac{\pi}{2}$ which implies that $E_*$ is unstable. From $\min \left(\det \left(k_{c}^{2} \right) \right) = 0$, we have $k_{c}^{2} = \sqrt{\frac{\det _0}{d_{1}^{+}d_2}}$ or $k_{c}^{2} = \sqrt{\frac{\det _0}{d_{1}^{-}d_2}}$. □

Remark 4.2. Take $\left[ \acute{a}, \grave{a} \right] = \left[ 1.31, 3.3 \right], \; [ \acute{b}, \grave{b} ] = \left[ 0.01, 0.04 \right], \; \left[ \acute{\alpha}, \grave{\alpha} \right] = \left[ 0.5, 1.5 \right], \; [ \acute{\beta}, \grave{\beta}] = \left[ 1.05, 1.5 \right], \; r = 1,$ and $d_2 = 1.24$. We have drawn the stable region of equilibrium point $E_*$ on the plane when $d_1 > 0, d_2 > 0$. According to Theorem 4.1, the stable region and the unstable region are distinguished. The critical value $d_1 = 0.1458$ was obtained through fixed parameters in Figure 1(a). Furthermore, we draw a graph $k$ about the wave number and the real part $\lambda$ of the eigenvalue as shown in Figure 1(b).

5.   Weakly nonlinear analysis

In this section, we will use weak nonlinear analysis to calculate the amplitude equation near the Turing instability threshold $d_1 = d_{1}^{c}$. Write system (2.3) in the following form:

where $L$ is a linear operator and $N$ is a nonlinear operator

and

with

We only consider the behavior of the control parameter near the bifurcation point, so the control parameter $d_1$ can be expanded as follows:

where $\varepsilon$ is a small parameter. At the same time, the variable $U$ and the nonlinear term $N$ are expanded according to this small parameter:

with

and

The linear operator $L$ can be decomposed into

where

We set $T_0 = t, T_1 = \varepsilon t, T_2 = \varepsilon ^2t$, and $T_3 = \varepsilon ^3t$, then the partial derivative of time can be written as follows:

Substitute formulas (5.8)–(5.15) into Eq (5.1), and the left side of the equation becomes:

and the right side of the equation becomes:

Comparing the order of $\varepsilon$ on both sides of the equation, the following three cases are obtained:

They are discussed separately below. For $\mathcal{O}\left(\varepsilon \right)$:

That is, $\left(\begin{array}{c} u_1\\ v_1\\ \end{array} \right)$ is a linear combination of eigenvectors corresponding to eigenvalues of $0$. Therefore,

and the general solution of Eq (5.19) can be written as:

where $\phi = -\frac{a_{22}+d_2k_{c}^{2}}{a_{21}}$, $|k_j| = k_c$, $k_{c}^{2} = \sqrt{\frac{\det _0}{d_1d_2}}$, and $k_j$ is the amplitude about the mode of $e^{-ik_jr}$. For $\mathcal{O}\left(\varepsilon ^2 \right)$:

According to the de Fredholm solvability condition, the vector function at the right end of Eq (5.25) must be orthogonal to the zero eigenvalue of $L_{c}^{+}$ for this equation to have a nontrivial solution.

The zero eigenvector as:

with $\varphi = -\frac{a_{12}}{a_{22}+d_2k_{c}^{2}}$. According to the orthogonal condition of Eq (5.20), we have

where $P_{u}^{j}$ and $P_{v}^{j}$ are the coefficients corresponding to $e^{ik_jr}$ in $P_u$ and $P_v$. The system of equations related to amplitude $A_j$ obtained from Eq (5.28) is:

where $h_1 = \frac{f_{uu}}{2}\phi ^2+f_{uv}\phi +\frac{f_{vv}}{2}$, and $h_2 = \frac{g_{uu}}{2}\phi ^2+g_{uv}\phi +\frac{g_{vv}}{2}$. Introducing a second-order disturbance term as:

Substitute formulas (5.24) and (5.30) into Eq (5.20), and we have

with

For $\mathcal{O}\left(\varepsilon ^3 \right)$:

According to the orthogonal condition of Eq (5.21), we have

Direct calculation produces the amplitude equation:

where

Suppose that the perturbation of amplitude $G$ under $\varepsilon$ is as follows:

Then, from formulas (5.15), (5.29), (5.36), and (5.38), we can derive

with

Since each amplitude $A_j = \rho _je^{i\psi _j}\left(j = 1, 2, 3 \right)$ in Eq (5.39) can be decomposed into mode $\rho _j = |A_j|$ and phase angle $\psi_j$, substituting $A_j$ into Eq (5.39) to separate the real and imaginary parts yields the following equation:

where $\psi = \psi _1+\psi _2+\psi _3$. We can infer from Eq (5.41) that the solution to the equation is stable when $h > 0, \psi = 0$ and $h < 0, \psi = \pi$. Eq (5.41) has the following solutions:

$1)$ Stationary state:

stable when $\mu < \mu _2 = 0$, and unstable when $\mu > \mu _2 = 0$.

$2)$ Strip pattern:

stable when $\mu > \mu _3 = \frac{h^2g_1}{\left(g_2-g_1 \right) ^2}$, and unstable when $\mu < \mu _3 = \frac{h^2g_1}{\left(g_2-g_1 \right) ^2}$.

$3)$ Hexagon pattern:

When $\mu > \mu _1 = \frac{-h^2}{4\left(g_1+2g_2 \right)}$ is satisfied, there exists

When $\mu < \mu _4 = \frac{\left(2g_1+g_2 \right) h^2}{\left(g_2-g_1 \right) ^2}$, $\rho ^+ = \frac{|h|+\sqrt{h^2+4\left(g_1+2g_2 \right) \mu}}{2\left(g_1+2g_2 \right)}$ is stable, and $\rho ^- = \frac{|h|-\sqrt{h^2+4\left(g_1+2g_2 \right) \mu}}{2\left(g_1+2g_2 \right)}$ is always unstable.

$4)$ Mixed state:

When $\mu > \mu _3 = \frac{h^2g_1}{\left(g_2-g_1 \right) ^2}$ is satisfied, there exists

It is always unstable with $g_1 < g _2$.

6.   Numerical simulation

In this section, we use the Euler discrete method for numerical simulation in two-dimensional space $\varOmega = \left[ 0, L_x \right] \times \left[ 0, L_y \right]$. Choose $L_x = 200, L_y = 200, t = 1000$, time step $\varDelta t = 0.9$, and space step $\varDelta h = 2$. We define $u_{pq}^{n} = u\left(x_p, y_q, n\varDelta t \right)$ and $v_{pq}^{n} = v\left(x_p, y_q, n\varDelta t \right)$ where $p, q = 1, 2, \cdots, \frac{L_x}{\varDelta h}$. System (2.3) is discretized by the Euler method as follows:

where

The parameters in system (2.3) are selected as follows:

and then, we obtain

The initial data are as follows:

The numerical simulation results show that there is a mixed structure solution under this set of parameters, and there are spots and stripe patterns in the graphics, as shown in Figure 2.

Under the same parameters, changing $r = 0.6$ yields $E_* = \left(0.21788, 0.14732 \right), \; \mu = -0.31413, \; \mu _3 = -4.52393, \; g_1 = -15680.19532$, and $g_2 = -20512.90109$. The numerical results show that there are also mixed structure solutions under this set of parameters, as seen in Figure 3.

Although we only change the interval variable $r$ in the graph, it can still be seen that the interval variable $r$ will affect the positive equilibrium point $E_*$ and the critical value $d_{1}^{c}$ of the Turing instability. Therefore, these two spatial patterns are slightly different and also prove the correctness of the theory.

Select the following symmetric initial conditions:

Applying symmetric initial conditions, numerical simulations are performed with varying the parameters $r, t$, and $d_1$ while other parameters remain constant, and the results show the existence of a symmetric hybrid structure solution, as seen in Figure 4. The figure shows the detailed pattern evolution of symmetric mixed patterns under different parameters. At first, circular patterns begin to appear under the initial conditions. As time $t$ progresses, the bounded domain is gradually destroyed and the circular pattern decomposes into striped and spotted patterns. It is found that as the diffusion rate $d_1$ of the prey population increases, the spatiotemporal chaos of the model is gradually suppressed.

7.   Conclusions

In this paper, the predator-prey model with an interval biological coefficient was analyzed theoretically and simulated numerically. We proved the existence and uniqueness of the model solution, discussed the stability of the positive equilibrium point, studied the Hopf bifurcation around the equilibrium point related to the fractional parameter $\alpha$, and discussed the Turing instability of the model at the starting point $d_1 = d_{1}^{-}$. It was found that the fractional $\alpha$ and the diffusion term $d_1$ played an important role in controlling the existence of the Hopf bifurcation and Turing instability. Then, the amplitude equation near the threshold of the Turing instability was given by using the weak nonlinear analysis method, and different mode selections were classified by using the amplitude equation. Finally, through numerical simulation, we showed the symmetric and asymmetric patterns of the model, and found that the diffusion rate of the prey population inhibited the spatiotemporal chaos of the model.

Author contributions

Conceptualization: Xiao-Long Gao and Hao-Lu Zhang; Methodology and software: Xiao-Yu Li; Data curation, formal analysis, and funding acquisition: Xiao-Long Gao and Xiao-Yu Li; Writing-original draft and writing-review and editing: Xiao-Long Gao, Hao-Lu Zhang, and Xiao-Yu Li. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

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

Acknowledgements

This paper is supported by the doctoral research start-up Fund of Inner Mongolia University of Technology (DC2300001252).

Conflicts of interest

The authors declare that there are no conflicts of interest regarding the publication of this article.