Dynamics analysis of a diffusional immunosuppressive infection model with Beddington-DeAngelis functional response

1.   Introduction

It is well known that optimal control problems play a very important role in the fields of science and engineering. In the operation of physical and economic processes, optimal control problems have a variety of applications. Therefore, highly effective numerical methods are key to the successful application of the optimal control problem in practice. The finite element method is an important method for solving optimal control problems and has been extensively studied in the literature. Many researchers have made various contributions on this topic. A systematic introduction to the finite element method for partial differential equations (PDEs) and optimal control problems can be found in ^[1,2]. For example, a priori error estimates of finite element approximation were established for the optimal control problems governed by linear elliptic and parabolic state equations, see ^[3,4]. Using adaptive finite element method to obtain posterior error estimation; see ^[5,6]. Furthermore, some superconvergence results have been established by applying recovery techniques, see ^[7,8].

The two-grid method based on two finite element spaces on one coarse and one fine grid was first proposed by Xu ^[9,10,11]. It is combined with other numerical methods to solve many partial differential equations, e.g., nonlinear elliptic problems ^[12], nonlinear parabolic equations ^[13], eigenvalue problems ^[14,15,16] and fractional differential equations ^[17].

Many real applications, such as heat conduction control of storage materials, population dynamics control and wave control problems governed by integro-differential equations, need to consider optimal control problems governed by elliptic integral equations and parabolic integro-differential equations. More and more experts and scholars began to pay attention to the numerical simulation of these optimal control problems. In ^[18], the authors analyzed the finite element method for optimal control problems governed by integral equations and integro-differential equations. In ^[19], the authors considered the error estimates of expanded mixed methods for optimal control problems governed by hyperbolic integro-differential equations. As far as we know, there is no research on a two-grid finite element method for parabolic integro-differential control problems in the existing literature.

In this paper, we design a two-grid scheme of fully discrete finite element approximation for optimal control problems governed by parabolic integro-differential equations. It is shown that when the coarse and fine mesh sizes satisfy $h = H^2$, the two-grid method achieves the same convergence property as the finite element method. We are interested in the following optimal control problems:

where $\Omega$ is a bounded domain in $\mathbb{R}^{2}$ and $J = (0, T]$. Let $K$ be a closed convex set in $U = L^2(J; L^2(\Omega))$, $f\in L^{2}(J; L^{2}(\Omega))$, $y_{d}\in H^{1}(J; L^{2}(\Omega))$ and $y_{0}\in H^{1}(\Omega)$. $K$ is a set defined by

$A = A(x) = (a_{ij}(x))$ is a symmetric matrix function with $a_{ij}(x)\in W^{1, \infty}(\Omega)$, which satisfies the ellipticity condition

Moreover, $B(t, s) = B(x, t, s)$ is also a $2\times2$ matrix; assume that there exists a positive constant $M$ such that

In this paper, we adopt the standard notation $W^{m, p}(\Omega)$ for Sobolev spaces on $\Omega$ with a norm $\|\cdot\|_{m, p}$ given by $\| v \|_{m, p}^{p} = \sum\limits_{|\alpha|\leq m}\| D^\alpha v\|_{L^{p}(\Omega)}^{p}$, as well as a semi-norm $|\cdot|_{m, p}$ given by $| v|_{m, p}^{p} = \sum\limits_{|\alpha| = m}\| D^\alpha v\|_{L^{p}(\Omega)}^{p}$. We set $W_0^{m, p}(\Omega) = \{v\in W^{m, p}(\Omega): v|_{\partial \Omega} = 0\}$. For $p = 2$, we denote $H^m(\Omega) = W^{m, 2}(\Omega), \ H_0^m(\Omega) = W_0^{m, 2}(\Omega)$, and $\|\cdot\|_{m} = \|\cdot\|_{m, 2}, \ \|\cdot\| = \|\cdot\|_{0, 2}$.

We denote by $L^s(J; W^{m, p}(\Omega))$ the Banach space of all $L^s$ integrable functions from $J$ into $W^{m, p}(\Omega)$ with the norm $\|v\|_{L^s(J; W^{m, p}(\Omega))} = \Big(\int_0^T||v||_{W^{m, p}(\Omega)}^sdt\Big)^{\frac{1}{s}}\ \rm{for}\ s\in [1, \infty)$ and the standard modification for $s = \infty$. For simplicity of presentation, we denote $\| v\|_{L^s(J; W^{m, p}(\Omega))}$ by $\| v\|_{L^s(W^{m, p})}$. Similarly, one can define the spaces $H^1(J; W^{m, p}(\Omega))$ and $C^k(J; W^{m, p}(\Omega))$. In addition $C$ denotes a general positive constant independent of $h$ and $\Delta t$, where $h$ is the spatial mesh size and $\Delta t$ is a time step.

The outline of this paper is as follows. In Section 2, we first construct a fully discrete finite element approximation scheme for the optimal control problems (1.1)–(1.4) and give its equivalent optimality conditions. In Section 3, we derive a priori error estimates for all variables, and then analyze the global superconvergence by using the recovery techniques. In Section 4, we present a two-grid scheme and discuss its convergence. In Section 5, we present a numerical example to verify the validity of the two-grid method.

2.   Fully discrete finite element scheme

In this section, we shall construct a fully discrete finite element approximation scheme for the control problems (1.1)–(1.4). For sake of simplicity, we take the state space $Q = L^{2}(J; V)$ and $V = H_0^{1}(\Omega)$.

We recast (1.1)–(1.4) in the following weak form: find $(y, u)\in Q\times K$ such that

where $(\cdot, \cdot)$ is the inner product of $L^{2}(\Omega)$.

Since the objective functional is convex, it follows from ^[2] that the optimal control problems (2.1)–(2.3) have a unique solution $(y, u)$, and that $(y, u)$ is the solution of (2.1)–(2.3) if and only if there is a co-state $p\in Q$ such that $(y, p, u)$ satisfies the following optimality conditions:

As in ^[20], the inequality (Eq 2.8) can be expressed as

where $\bar{p} = \frac{\int_\Omega pdx}{\int_\Omega dx}$ denotes the integral average on $\Omega$ of the function $p$.

Let $T_{h}$ denote a regular triangulation of the polygonal domain $\Omega$, $h_{\tau}$ denote the diameter of $\tau$ and $h = \max\limits_{\tau \in T_{h}}{h_{\tau}}$. Let $V_h\subset V$ be defined by the following finite element space:

And the approximated space of control is given by

Set $K_h = U_h\cap K.$

Before the fully discrete finite element scheme is given, we introduce some projection operators. First, we define the Ritz-Volterra projection ^[21] $R_h:\ V\rightarrow V_h$, which satisfies the following: for any $y, p\in V$

Next, we define the standard $L^2$-orthogonal projection ^[22] $Q_h:\ L^{2}(\Omega) \rightarrow U_h$, which satisfies the following: for any $\phi\in L^{2}(\Omega)$

At last, we define the element average operator ^[7] $\pi_h: L^{2}(\Omega)\rightarrow U_h$ by

We have the approximation property

We now consider the fully discrete finite element approximation for the control problem. Let $\Delta t > 0$, $N = T/\Delta t\in \mathbb{Z}$ and $t_n = n \Delta t$, $n\in \mathbb{Z}$. Also, let

Like in ^[23], we define for $1\leq s\leq \infty$ and $s = \infty$, the discrete time dependent norms

where $l = 0$ for the control variable $u$ and the state variable $y$, and $l = 1$ for the co-state variable $p$.

Then the fully discrete approximation scheme is to find $(y_h^n, u_h^n)\in V_h\times K_h$, $n = 1, 2, \cdots, N$, such that

Again, we can see that the above optimal control problem has a unique solution $(y_h^n, u_h^n)$, and that $(y_h^n, u_h^n)\in V_h\times K_{h}$ is the solution of (2.20)–(2.22) if and only if there is a co-state $p_h^{n-1}\in V_h$ such that $(y_h^n, p_h^{n-1}, u_h^n)$ satisfies the following optimality conditions:

Similarly, employing the projection (2.9), the optimal condition (2.27) can be rewritten as follows:

where $\overline{p_h^{n-1}} = \frac{\int_{\Omega}p_h^{n-1}}{\int_{\Omega}1}$.

In the rest of the paper, we shall use some intermediate variables. For any control function $\tilde{u}\in K$ satisfies the following:

3.   A priori error estimates and superconvergence

In this section, we will discuss a priori error estimates and superconvergence of the fully discrete case for the state variable, the co-state variable and the control variable. In order to do it, we need the following lemmas.

Lemma 3.1. Let $(y_h^n(u), p_h^{n-1}(u))$ be the solution of (2.29)–(2.32) with $\tilde{u} = u$ and $(y, p)$ be the solution of (2.4)–(2.8). Assume that the exact solution $(y, p)$ has enough regularities for our purpose. Then, for $\Delta t$ small enough and $1\leq n\leq N$, we have

Proof. For convenience, let

Taking $t = t_n$ in (2.4), subtracting (2.29) from (2.4) and then using (2.12), we have

Choosing $v_h = dt\xi_{y}^n$ in (3.3), we get

Notice that

Multiplying $\Delta t$ and summing over $n$ from $1$ to $l$ $(1\leq l\leq N)$ on both sides of (3.4), and by using (3.5) and $\xi_y^0 = 0$, we find that

Now, we estimate the right-hand terms of (3.6). For $A_1$, from the results given in ^[24], we have

For $A_2$, using (2.13), the Hölder inequality and the Cauchy inequality, we have

At last, for $A_3$, it follows from the Cauchy inequality, Cauchy mean value theorem and assumptions on $A$ and $B$ that

where

and

where $t_{i}^*$ is located between $t_{i-1}$ and $t_{i}$, and we also used

From (3.7)–(3.9), we have

Adding $\sum_{n = 1}^l\|\nabla \xi_{y}^n\|^2\Delta t$ to both sides of (3.10), by use of the assumption on $A$ and discrete Gronwall's inequality, we have

Using (2.13), the Poincare inequality and the triangle inequality, we get

Taking $t = t_{n-1}$ in (2.6), subtracting (2.31) from (2.6) and then using (2.14), we have

Choosing $q_h = -dt\xi_{p}^{n}$ in (3.13), multiplying by $\Delta t$ and summing over $n$ from $l+1$ to $N\ (0\leq l\leq N-1)$ on both sides of (3.13), since $\xi_p^N = 0$, we find that

Notice that

Now, we estimate the right-hand terms of (3.14). Similar to (3.7), we have

For $B_2$, using (2.15) and the Cauchy inequality, we have

For $B_3$, applying the same estimates as $A_3$, we conclude that

where

For $B_4$, using the Cauchy inequality and the smoothness of $y$ and $y_d$, we have

Combining (3.16)–(3.19), we have

By adding $\sum_{n = l+1}^N\|\nabla \xi_{p}^n\|^2\Delta t$ to both sides of (3.20) and applying the assumption on $A$, discrete Gronwall's inequality and (3.12), we conclude that

Using (2.15) and the triangle inequality, we get

we have completed the proof of the Lemma 3.1.

Lemma 3.2. Choose $\tilde{u}^n = Q_h u^n$ and $\tilde{u}^n = u^n$ in (2.29)–(2.32) respectively. Then, for $\Delta t$ small enough and $1\leq n\leq N$, we have

Proof. For convenience, let

Taking $\tilde{u}^n = u^n$ and $\tilde{u}^n = Q_hu^n$ in (2.29), we easily get

By choosing $v_h = dt\lambda_y^n$ in (3.24), multiplying by $\Delta t$ and summing over $n$ from $1$ to $l$ $(1\leq l\leq N)$ on both sides of (3.24), we find that

where we use (2.17) and the assumption on $B$; additionally, $\theta^n$ is located between $t_{n}$ and $t_{n+1}$.

Add $\sum_{n = 1}^l\|\nabla \lambda_{y}^n\|^2\Delta t$ to both sides of (3.25); then for sufficiently small $\Delta t$, combining (3.25) and the discrete Gronwall inequality, we have

Similar to (3.24), we have

By choosing $q_h = -dt\lambda_p^n$ in (3.27), multiplying by $\Delta t$ and summing over $n$ from $l+1$ to $N$ $(0\leq l\leq N-1)$ on both sides of (3.27), combining (3.26) and Poincare inequality gives

Add $\sum_{n = l+1}^N\|\nabla \lambda_{p}^{n-1}\|^2\Delta t$ to both sides of (3.28); then for sufficiently small $\Delta t$, applying the discrete Gronwall inequality and the assumptions on $A$ and $B$, we have

Using the stability analysis as in Lemma 3.2 yields Lemma 3.3.

Lemma 3.3. Let $(y_{h}^n, p_{h}^n)$ and $(y_{h}^n(Q_hu), p_{h}^n(Q_hu))$ be the discrete solutions of (2.29)$-$(2.32) with $\tilde{u}^{n} = u_h^{n}$ and $\tilde{u}^{n} = Q_hu^{n}$, respectively. Then, for $\Delta t$ small enough and $1\leq n\leq N$, we have

Next, we derive the following inequality.

Lemma 3.4. Choose $\tilde{u}^n = Q_h u^n$ and $\tilde{u}^n = u_h^n$ in (2.29)$-$(2.32) respectively. Then, we have

Proof. For $n = 0, 1, \ldots, N$, let

From (2.29)–(2.32), we have

Notice that

and

By choosing $v_h = -r_{p}^{n-1}$ in (3.32), $q_h = r_{y}^{n}$ in (3.33), and then multiplying the two resulting equations by $\Delta t$ and summing it over $n$ from $1$ to $N$, we have

which completes the proof of the lemma.

Lemma 3.5. Let $u$ be the solution of (2.4)–(2.8) and $u_{h}^{n}$ be the solution of (2.23)–(2.27). Assume that all of the conditions in Lemmas 3.1–3.4 are valid. Then, for $\Delta t$ small enough and $1\leq n\leq N$, we have

Proof. Take $\tilde{u} = u_h^{n}$ in (2.8) and $\tilde{u}_h = Q_h u^{n}$ in (2.27) to get the following two inequalities:

and

Note that $u_h^{n}-u^{n} = u_h^{n}-Q_hu^{n}+Q_hu^{n}-u^{n}$. Adding the two inequalities (3.36) and (3.37), we have

Thus, by (3.38), (2.16), (2.8) and Lemma 3.4, we find that

It follows from the Cauchy inequality, Lemma 3.1, Lemma 3.2 and Poincare's inequality that

Substituting the estimates for $F_{1}$–$F_{3}$ into (3.39), we derive (3.35).

Using (3.11), (3.21), Lemmas 3.2–3.5 and the triangle inequality, we derive the following superconvergence for the state variable.

Lemma 3.6. Let $u$ be the solution of (2.4)–(2.8) and $u_{h}^{n}$ be the solution of (2.23)–(2.27). Assume that all of the conditions in Lemmas 3.1–3.5 are valid. Then, for $\Delta t$ small enough and $1\leq n\leq N$, we have

Now, the main result of this section is given in the following theorem.

Theorem 3.1. Let $(y, p, u)$ and $(y_{h}^{n}, p_{h}^{n-1}, u_{h}^n)$ be the solutions of (2.4)–(2.8) and (2.23)–(2.27), respectively. Assume that $y$, $p$ and $u$ have enough regularities for our purpose; then, for $\Delta t$ small enough and $1\leq n\leq N$, we have

Proof. The proof of the theorem can be completed by using Lemmas 3.1–3.5, (2.17) and the triangle inequality.

To provide the global superconvergence for the control and state, we use the recovery techniques on uniform meshes. Let us construct the recovery operators $P_h$ and $G_h$. Let $P_h v$ be a continuous piecewise linear function (without the zero boundary constraint). The value of $P_h v$ on the nodes are defined by a least squares argument on element patches surrounding the nodes; the details can be found in ^[25,26].

We construct the gradient recovery operator $G_h v = (P_h v_x, P_h v_y)$ for the gradients of $y$ and $p$. In the piecewise linear case, it is noted to be the same as the Z-Z gradient recovery (see ^[25,26]). We construct the discrete co-state with the admissible set

Now, we can derive the global superconvergence result for the control variable and state variable.

Theorem 3.2. Let $u$ and $u_{h}^{n}$ be the solutions of (2.4)–(2.8) and (2.29)–(2.32), respectively. Assume that all of the conditions in Lemmas 3.1–3.5 are valid. Then we have

Proof. Using (2.9), (3.47) and Theorem 3.1, we have

Theorem 3.3. Let $(y, p)$ and $(y_{h}^{n}, p_{h}^{n-1})$ be the solutions of (2.4)–(2.8) and (2.29)–(2.32), respectively. Assume that all of the conditions in Lemmas 3.1–3.5 are valid. Then we have

Proof. Notice that

It follows from Lemma 3.6 that

It can be proved by the standard interpolation error estimate technique (see ^[1]) that

Therefore, it follows from (3.52) and (3.53) that

Similarly, it can be proved that

Therefore, we complete the proof.

4.   Two-grid scheme

In this section, we will present a two-grid scheme and analyze a priori error estimates. Now, we present our two-grid algorithm which has the following two steps:

$\bf{Step}$ $\bf{1}$. On the coarse grid $T_H$, find $(y_H^n, p_H^{n-1}, u_H^n)\in V_H^{2}\times K_{H}$ that satisfies the following optimality conditions:

$\bf{Step}$ $\bf{2}$. On the fine grid $T_h$, find $(\overline{\widetilde{y}}_h^n, \overline{\widetilde{p}}_h^{n-1}, \overline{\widetilde{u}}_h^n)\in V_h^{2}\times K_{h}$ such that

Combining Theorem 3.1 and the stability estimates, we easily get the following results.

Theorem 4.1. Let $(y, p, u)$ and $(\overline{\widetilde{y}}_h^n, \overline{\widetilde{p}}_h^n, \overline{\widetilde{u}}_h^n)$ be the solutions of (2.4)–(2.8) and (4.1)–(4.10), respectively. Assume that $y$, $y_{d}$, $p$, $p_d$ and $u$ have enough regularities for our purpose; then, for $\Delta t$ small enough and $1\leq n\leq N$, we have

Proof. For convenience, let

Taking $t = t_{n}$ in (2.4), subtracting (4.6) from (2.4) and then using (2.12), we have

Selecting $v_h = dte_y^n$ in (4.13), multiplying by $\Delta t$ and summing over $n$ from $1$ to $l$ $(1\leq l\leq N)$ on both sides of (4.13), we find that

Similar to Lemma 3.1, it is easy to show that

Similar to $A_3$, we find that

For $I_4$, using Theorem 3.2, we have

Combining (4.15)–(4.17), the discrete Gronwall inequality, the triangle inequality and (2.13), we get

By taking $t = t_{n-1}$ in (2.6), subtracting (4.8) from (2.6) and using (2.12), we have

By selecting $q_h = -dte_p^n$ in (4.19), multiplying by $\Delta t$ and summing over $n$ from $l+1$ to $N$ $(0\leq l\leq N-1)$ on both sides of (4.19), we find that using (2.15), (4.18) and the triangle inequality, similar to (3.14), gives

Note that

Using (2.19), (4.20) and the mean value theorem, we have

which completes the proof.

5.   Numerical experiments

In this section, we present the following numerical experiment to verify the theoretical results. We consider the following two-dimensional parabolic integro-differential optimal control problems

subject to

where $\Omega = (0, 1)^2$.

We applied a piecewise linear finite element method for the state variable $y$ and co-state variable $p$. The stopping criterion of the finite element method was chosen to be the abstract error of control variable $u$ between two adjacent iterates less than a prescribed tolerance, i.e.,

where $\epsilon = 10^{-5}$ was used in our numerical tests. For the linear system of equations, we used the algebraic multigrid method with tolerance $10^{-9}$.

The numerical experiments were conducted on a desktop computer with a 2.6 GHz 4-core Intel i7-6700HQ CPU and 8 GB 2133 MHz DDR4 memory. The MATLAB finite element package iFEM was used for the implementation ^[27].

Example: We chose the following source function $f$ and the desired state $y_d$ as

such that the exact solutions for $y$, $p$, $u$ are respectively,

In order to see the convergence order with respect to time step size $\triangle t$ and mesh size $h$, we choose $\triangle t = h$ or $\triangle t = h_2$ with $h = \frac{1}{4}, \frac{1}{16}, \frac{1}{64}$. To see the convergence order of the two-grid method, we choose the coarse and fine mesh size pairs $(\frac{1}{2}, \frac{1}{4}), (\frac{1}{4}, \frac{1}{16}), (\frac{1}{8}, \frac{1}{64})$. Let us use $y^h, p^h$ and $u^h$ as two-grid solutions in the following tables. In Tables 1 and 2, we let $\triangle t = h_2$ and present the errors of the finite element method and two-grid method for $y$ and $p$ in the $L^2$-norm. Next, in Tables 3 and 4, we set $\triangle t = h$ and show the errors of the two methods for $y$ and $p$ in the $H^1$-norm and $u$ in the $L^2$-norm. We can see that the two-grid method maintains the same convergence order as the finite element method. Moreover, we also display the computing times of the finite element method and the two-grid method in these tables. By comparison, we find that the two-grid method is more effective for solving the optimal control problems (1.1)–(1.4).

6.   Conclusions

In this paper, we presented a two-grid finite element scheme for linear parabolic integro-differential control problems (1.1)–(1.4). A priori error estimates for the two-grid method and finite element method have been derived. We have used recovery operators to prove the superconvergence results. These results seem to be new in the literature. In our future work, we will investigate a posteriori error estimates. Furthermore, we shall consider a priori error estimates and a posteriori error estimates for optimal control problems governed by hyperbolic integro-differential equations.

Use of AI tools declaration

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

Conflict of interest

The authors declare there is no conflict of interest.