Stability and Hopf bifurcation of the coexistence equilibrium for a differential-algebraic biological economic system with predator harvesting

1.   Introduction

In this manuscript, we detail the development of a new stabilizer-free weak Galerkin (WG) finite element method of any polynomial order in 2D and 3D, on triangular and tetrahedral meshes respectively, for obtaining the solutions of the stationary Stokes equations: Find unknown functions ${\bf{u}}$ (velocity) and $p$ (pressure) such that

where the viscosity $\mu > 0$, and the domain $\Omega$ is a polygon or a polyhedron in $\mathbb{R}^d\; (d = 2, 3)$.

The new pressure-robust, stabilizer-free WG method has a new variational formulation. Find ${\bf{u}}_h\in V_h$ and $p_h\in W_h$ such that

where $\nabla_w$ denotes weak gradients to be defined in (2.3) and (2.4), and $V_h$ and $W_h$ are $P_k$ and $P_{k-1}$ WG finite element spaces, to be defined in (2.1) and (2.2), respectively. Some studies on the WG methods are as follows: the parabolic equation ^[1,2,3], with two-order superconvergence on triangular meshes ^[4], the convection-diffusion equation ^[5,6], 4th order problems ^[7,8], with the conforming discontinuous Galerkin formulation ^[9,10,11], the Navier-Stokes equations ^[12,13,14], with the discrete maximum principle ^[15], the second order elliptic equations ^[16,17], with the energy conservation ^[18], the Darcy flow ^[19,20], the Oseen equations ^[21], the Stokes equations with pressure-robustness ^[22,23,24], the Maxwell equations ^[25], the Cahn-Hilliard equation ^[26], the div-curl equations ^[27], with adaptive refinements ^[28,29], the biharmonic equation ^[30,31], the biharmonic equation with continuous finite elements^[32,33], the Stokes equations with $H({\operatorname{{div}}})$ elements ^[34,35], with two-order superconvergence under the CDG formulation ^[36,37], and with two-order superconvergence for the Stokes equations ^[38].

We note that the moment equation (1.1) is not tested by applying $H({\operatorname{{div}}}, \Omega)$ functions in (1.4), different from most other pressure-robust methods such as those detailed in ^[22,23,24], but is tested by using discontinuous polynomials. The pointwise divergence-free property (pressure-robustness) is achieved by introducing pressure face variables, i.e., $p_h = \{p_0, p_b\}$ where $p_b|_e\in P_k(e)$ on every face edge/triangle $e$ of a mesh $\mathcal T_h$ (see (2.2) below). The method was previously applied before in hybridized discontinuous Galerkin methods in ^[39].

It is shown that the WG finite element pair, $V_h-W_h$, is inf-sup stable. The method is shown pressure robust, i.e., both errors of the velocity and the pressure are independent of the pressure $p$ and the viscosity $\mu$ in (1.1). This is necessary when $\mu$ is small. However, very few methods can achieve this. Under the inf-sup stability, we shall prove quasi-optimal approximation for the velocity in an $H^1$-norm and in the $L^2$ norm. The quasi-optimal convergence for pressure in the $L^2$ norm is also proved. The theory is numerically verified by applying varying degrees of WG finite elements to both triangular and tetrahedral meshes.

2.   The WG finite element method

Let ${ \mathcal T}_h$ be a mesh on the domain $\Omega$, consisting of conforming shape-regular triangles or shape-regular tetrahedrons. Here we let $h_T$ be the element diameter of $T\in {{\mathcal T}}_h$, and we let mesh size $h = \max_{T\in{{\mathcal T}}_h} h_T$.

For $k\ge 1$ and given ${{\mathcal T}}_h$, the finite element space for velocity is defined by

and the finite element space for pressure is defined by

where $P_k(e)$ and $P_k(T)$ denote the space of polynomials of degree $k$ or less on the edge/triangle $e$ and triangle/tetrahedron $T$ respectively, $(\cdot, \cdot)_{\mathcal T_h} = \sum_{T\in \mathcal T_h} (\cdot, \cdot)_T$ and $\langle\cdot, \cdot \rangle_{\partial \mathcal T_h} = \sum_{T\in \mathcal T_h} \langle \cdot, \cdot\rangle_{\partial T}$.

For a function ${\bf{v}}\in V_h$, the $(k+1)$-degree weak gradient $\nabla_w {\bf{v}}$ is a piecewise polynomial on the mesh $\mathcal T_h$, $\nabla_w {\bf{v}}|_T \in [P_{k+1}(T)]^{d\times d}$, such that

For a function $q = \{q_0, q_b\}\in W_h$, its weak gradient $\nabla_w q$ is defined as a piecewise vector-valued polynomial such that for each $T\in{{\mathcal T}}_h$, $\nabla_w q \in [P_k(T)]^d$ satisfies

We denote by $\Pi_k$ the local/element-wise $L^2$-orthogonal projection onto $[P_k(T)]^j$ where $j = 1, \; d, \; d \; \times\; d$ and $T\in{{\mathcal T}}_h$. Let $\Pi_k^b$ be a generic edge/face-wise defined $L^2$ projection onto $[P_k(e)]^j$ for $e\in{{\partial T}}$. Define ${\bf{Q}}_h {\bf{u}} = \{\Pi_k {\bf{u}}, \Pi_{k+1}^b {\bf{u}}\}\in V_h$ and $Q_h p = \{\Pi_{k-1}p, \Pi_k^bp\}\in W_h$.

Lemma 2.1. Let $\boldsymbol\phi\in [H_0^1(\Omega)]^d$ and $\psi\in H^1(\Omega)$, then, for $T\in{{\mathcal T}}_h$

Proof. Using (2.3) and integration by parts, we obtain the following for any $\tau\in [P_{k+1}(T)]^{d\times d}$:

This implies that (2.5) holds. The scalar version, (2.6), is proved in the same manner.

For any function $\varphi\in H^1(T)$, the following trace inequality holds true:

We define two semi-norms $||| {\bf{v}}|||$ and $\| {\bf{v}}\|_{1, h}$ for any ${\bf{v}}\in V_h$:

We also define two semi-norms $||| q|||$ and $\|q\|_{1, h}$ for any $q\in W_h$:

In fact, $\| {\bf{v}}\|_{1, h}$ is a norm in $V_h$ and $||| \cdot |||$ is also a norm in $V_h$, as they have been proved in ^[4] by applying the norm equivalence as follows:

and

Lemma 2.2. The following inf-sup conditions hold, for all $q = \{q_0, q_b\}\in W_h$ and ${\bf{v}} = \{ {\bf{v}}_0, {\bf{v}}_b\}\in V_h$:

and

where $\beta > 0$ is independent of $h$ and $\mathcal T_h$.

Proof. For any given $q = \{q_0, q_b\}\in W_h$, it is known that there exists a function $\tilde {\bf{v}}\in H_0({\operatorname{{div}}}, \Omega)$ and $\tilde {\bf{v}}|_T\in [P_k(T)]^d$ (see, e.g., [40, (7.4)–(7.6)]) such that

where

Let

where $T_1$ and $T_2$ are the two elements on the two sides of edge/triangle $e$. For such specially defined ${\bf{v}}$, we have

It follows from (2.13) and (2.17) that

By (2.4), we have

Combining (2.16), (2.19) and (2.18) implies that

which implies (2.14).

Next we shall derive (2.15). For any ${\bf{v}} = \{ {\bf{v}}_0, {{\mathbf 0}}\}\in V_h$ and $\tau\in [P_{k+1}(T)]^{d\times d}$, we have the following by (2.3), (2.7) and the inverse inequalities:

which implies that

It follows from (2.20) that for any ${\bf{v}} = \{ {\bf{v}}_0, {{\mathbf 0}}\}\in V_h$

Then we have

which implies (2.15). This completes the proof of the lemma.

Lemma 2.3. There is a unique solution for the WG finite element equations (1.4) and (1.5).

Proof. We only need to show that zero is the unique solution of (1.4) and (1.5) if ${{\bf{f}}} = {{\mathbf 0}}$. We let ${{\bf{f}}} = {{\mathbf 0}}$ and ${\bf{v}} = {\bf{u}}_h$ in (1.4) and $q = p_h$ in (1.5). By summing the two equations, we get

It implies that $\nabla_w {\bf{u}}_h = 0$ on $T$. By (2.12), we also obtain that $\| {\bf{u}}_h\|_{1, h} = 0$. Thus, ${\bf{u}}_h = {{\mathbf 0}}$.

Since ${\bf{u}}_h = {{\mathbf 0}}$ and ${{\bf{f}}} = {{\mathbf 0}}$, (1.4) is reduced to $({\bf{v}}_0, \nabla_wp_h) = 0$ for any ${\bf{v}} = \{ {\bf{v}}_0, {\bf{v}}_b\}\in V_h$. Then the inf-sup conditions (2.15) and (2.14) imply that $\nabla_wp_h = 0$ and $p_0 = 0$. By (2.13), we have that $\|p_h\|_{1, h} = 0$ and $p_0 = p_b = 0$ on ${{\partial T}}$. We have proved the lemma.

3.   Error equations

To derive the equations that the errors satisfy, we introduce ${\bf{e}}_h = {\bf{Q}}_h {\bf{u}}- {\bf{u}}_h$ and $\epsilon_h = Q_hp-p_h$.

Lemma 3.1. [41, Theorem 3.1] For $\tau\in [H^{k+2}(\Omega)]^d$, a quasi-projection $\Pi_h$ can be defined such that $\Pi_h\tau\in [H({\operatorname{{div}}}, \Omega)]^d$, $\Pi_h\tau|_T\in [P_{k+1}(T)]^{d\times d}$ and for ${\bf{v}}_0\in [P_k(T)]^d$,

where $|\cdot|_{k+2}$ is the semi-$H^{k+2}$ Sobolev norm on the space.

Lemma 3.2. [41, Theorem 3.1] Let $\tau\in H^{k+1}(\Omega)$. A quasi-projection $\pi_h$ can be defined such that $\pi_h\tau\in H({\operatorname{{div}}}, \Omega)$, $\pi_h\tau|_T\in [P_{k}(T)]^d$ and for $q_0\in P_{k-1}(T)$,

where $|\cdot|_{k+1}$ is the semi-$H^{k+1}$ Sobolev norm on the space.

Lemma 3.3. For any ${\bf{v}}\in V_h$ and $q\in W_h$, the following error equations hold true:

where

Proof. First, we test (1.1) by applying ${\bf{v}}_0$ with ${\bf{v}} = \{ {\bf{v}}_0, {\bf{v}}_b\}\in V_h$ to obtain

It follows from (3.2) and (2.5) that

It follows from (2.6) that

We substitute (3.12) and (3.13) into (3.11) to obtain

We subtract (3.14) from (1.4) to get

Multiplying (1.2) by $q = \{q_0, q_b\}\in W_h$, by applying (3.2), it follows that

which implies that

The difference between (3.17) and (1.5) implies (3.8). We have proved the lemma.

4.   Quasi-optimal finite element solutions

We shall first prove the optimal order error estimates of the $|||\cdot|||$ norm for the velocity ${\bf{u}}_h$, and of the $L^2$ norm for the pressure $p_h$.

Lemma 4.1. Let ${\bf{u}}\in [H^{k+1}(\Omega)]^d$, ${\bf{v}}\in V_h$ and $q\in W_h$. The following estimates hold:

where $e_1(\cdot, \cdot)$ and $e_2(\cdot, \cdot)$ have been defined in (3.9) and (3.10), respectively.

Proof. By the Cauchy-Schwarz inequality, the definitions of $\Pi_{k+1}$ and $\Pi_h$, we compute

Similarly, we have

We have proved the lemma.

Theorem 4.1. Let $({\bf{u}}, p)\in \big([H^{k+1}(\Omega)]^d \cap [H^1_0(\Omega)]^d \big) \times \big(H^k(\Omega)\cap L^2_0(\Omega)\big)$ be the solutions of (1.1)–(1.3). Let $({\bf{u}}_h, p_h)\in V_h \times W_h$ be the solutions of (1.4) and (1.5). Then, the following error estimates hold true:

Proof. It follows from (3.7) that for any ${\bf{v}} = \{ {\bf{v}}_0, {\bf{v}}_b\}\in V_h$, by (4.1), we have

Then applying the estimate (4.5) and (2.15) yields

By letting ${\bf{v}} = {\bf{e}}_h$ in (3.7) and $q = \epsilon_h$ in (3.8), and by using (3.8), we have

It then follows from (4.1), (4.2) and (4.6) that

which implies (4.3). The estimate (4.4) follows from (4.5), (4.3) and the inf-sup conditions (2.14) and (2.15). We have proved the theorem.

We shall derive next the optimal-order convergence for velocity in the $L^2$ norm by using the duality argument. Recall that ${\bf{e}}_h = \{ {\bf{e}}_0, {\bf{e}}_b\} = {\bf{Q}}_h {\bf{u}}- {\bf{u}}_h$ and $\epsilon_h = Q_hp-p_h$. Consider the problem of seeking $({{\boldsymbol{\psi}}}, \xi)$ such that

Assume that the duality problem given by (4.7)–(4.9) has the ${{\bf H}}^{2}(\Omega)\times H^1(\Omega)$-regularity property that the solution $({{\boldsymbol{\psi}}}, \xi)\in {{\bf H}}^{2}(\Omega)\times H^1(\Omega)$ and the following a priori estimate holds true:

We need the following lemma first.

Lemma 4.2. For any ${\bf{v}}\in V_h$ and $q\in W_h$, the following equations hold true:

where ${{\boldsymbol{\psi}}}$ and $\xi$ are defined in (4.7), and

Proof. Testing (4.7) by applying ${\bf{v}}_0$ with ${\bf{v}} = \{ {\bf{v}}_0, {\bf{v}}_b\}\in V_h$ gives

By performing integration by parts, and setting $\langle\nabla {{\boldsymbol{\psi}}}\cdot{{\bf{n}}}, {\bf{v}}_b\rangle_{\partial{{\mathcal T}}_h} = 0$, we derive

By integration by parts, and given(2.3) and (2.5), we have

Combining (4.14) and (4.15) gives

By applying the definition of $\Pi_k$, (2.6), and (3.13), we obtain

Combining (4.16) and (4.17) with (4.13) yields (4.11).

Testing (4.8) by applying $q_0$ with $q = \{q_0, q_b\}\in W_h$ gives

By applying integration by parts, we obtain

Combining (4.18) and (4.19) yields

We have proved the lemma.

By the same argument as that for (4.16), (3.7) has another form, i.e.,

Letting ${\bf{v}} = {\bf{Q}}_h{{\boldsymbol{\psi}}}$ and $q = Q_h\xi$ in (4.21) and (4.22), we obtain

Letting ${\bf{v}} = {\bf{e}}_h$ and $q = \epsilon_h$ in (4.11) and (4.12), we have

By applying (4.26), (4.23) becomes

Theorem 4.2. Let $({\bf{u}}, p)\in \big([H^{k+1}(\Omega)]^d \cap [H^1_0(\Omega)]^d \big) \times \big(H^k(\Omega)\cap L^2_0(\Omega)\big)$ be the solutions of (1.1)–(1.3). Let $({\bf{u}}_h, p_h)\in V_h \times W_h$ denote the unique solutions of (1.4) and (1.5). With the condition (4.10), the following error bound holds:

Proof. Letting ${\bf{v}} = {\bf{e}}_h$ in (4.11) yields

By applying (4.27) and (4.24), (4.29) becomes

Next we shall estimate all of the terms on the right hand side of (4.30). Using the Cauchy-Schwarz inequality, the trace inequality (2.7), and the definition of $\Pi_{k+1}$, we obtain

Similarly, we have

Using the Cauchy-Schwarz inequality and the trace inequalities, applying (2.12) and (4.3), we obtain

By (4.4), we have

Substituting all the four bounds above into (4.30), we get

By applying this inequality and the regularity condition (4.10), (4.28) is proved.

5.   Numerical test

In the first example in 2D, we have chosen the domain $\Omega = (0, 1)\; \times\; (0, 1)$ for the Stokes equations (1.1)–(1.3). We have chosen an ${{\mathbf f}}$ (depending on $\mu$) in (1.1) such that the exact solution of (1.1)–(1.3) is as follows (independent of $\mu$):

We computed the solution (5.1) on triangular grids shown, as in Figure 1 for the $P_k$-WG/$P_{k-1}$-WG mixed finite elements for $k = 1, 2, 3, 4$, and $5$. The results are listed in Tables 1–5. As shown, the optimal order of convergence has been achieved for all solutions in all norms. From the data, we can see the method is pressure-robust and the error is independent of the viscosity $\mu$.

We note that for some high level grids the computer round-off error was found to exceed the truncation error when $\mu = 10^{-6}$, as described in Tables 3–5.

We computed the 2D solution (5.1) again for slightly perturbed triangular meshes, as illustrated in Figure 2 by employing the $P_k$-WG/$P_{k-1}$-WG mixed finite elements for $k = 1, 2, 3, 4$ and $5$. The computational results are listed in Tables 6–10. The quasi-optimal convergence has been achieved for all solutions in all norms.

In the third test, we performed 3D numerical computation on domain $\Omega = (0, 1)\times(0, 1)\times(0, 1)$. We chose an ${{\mathbf f}}$ in (1.1) such that we would have the following exact solution

The 3D meshes are illustrated in Figure 3. The computational results are listed in Tables 11–13.

Use of AI tools declaration

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

Acknowledgments

Yan Yang is supported in part by the Program of Sichuan National Applied Mathematics Center, No. 2023-KFJJ-01-001.

Conflict of interest

The authors declare there are no conflicts of interest.