Stability and bifurcation of a delayed prey-predator eco-epidemiological model with the impact of media

1.   Introduction

With the development of computer technology, the finite element method follows closely in the early 1960s. The research on the reliability and validity of finite element analysis promoted the development of the finite element method ^{[3,4,5,25,27]}. Applying finite element methods, the emergence of errors has captured the attention of scholars. One of the main sources of errors is the error caused by the discretisation of the model, while researchers dissect various aspects of finite element analysis. In addition, the generation of finite element mesh is also a concern for scholars. In the early conventional FEA, scholars usually used experience, intuition or even guesses to generate meshes to simple judge whether the approximation results are reasonable or not. If it is not reasonable, the grid needs to be redesigned whenever necessary for the efficiency of the analysis and the reliability of the results. Therefore, the emergence of the adaptive finite element method is because the computer, after checking the current conditions, decides whether the solution is accurate enough to meet the determination needs according to the error information obtained during the adjustment process.

To the best of our knowledge, the adaptive finite element method is a numerical method that can automatically adjust the algorithm to improve the process of solving ^[2]. An appropriate mesh can greatly reduce the errors arising from the discretisation of the finite element approximation process during the replication. According to the current situation, solutions of optimal control problems for nonlinear systems are usually not available. Besides the complexity and diversity of nonlinear equations, it is also very practical to use the adaptive finite element method for solving nonlinear equations.

Adaptive finite element methods have been widely and successfully applied in various linear optimal control problems, for example, Eriksson and Johnson proposed that the adaptive finite element algorithm produces a series of successively refined meshes in which the final mesh satisfies a given error tolerance ^[11]. Gaevaskaya and Hope et al. proposed an adaptive finite element method for a class of distributed optimal control problems with control constraints is analysed by applying the reliability and discrete local efficiency of the posterior estimator and the quasi-orthogonality property as basic tools ^[13]. Braess and Carstensen et al. Investigate the residual jump contributions of a standard explicit residual-based a posteriori error estimator for an adaptive finite element method ^[1]. Hu and Xu et al. performed research on the convergence and optimality of the adaptive nonconforming linear element method for the stokes problem ^[17].

However, a large number of literature show that the researches on adaptive finite element method for nonlinear optimal control problems have not reached its peak yet.

Recently, for instance, Wriggers and Scherf propose an adaptive finite element technique for nonlinear contact problems ^[29]. Eriksson and Johnson consider adaptive finite element methods for parabolic problems to a class of nonlinear scalar problems, the authors obtain posteriori error estimates and design corresponding adaptive algorithms ^[12]. Nowadays, the problem of nonlinear optimal control problems, similar to big data, is the focus of scholars worldwide. Hence it is worthy of investigating the adaptive finite element method for such nonlinear problems.

Many scholars have also studied the prior error estimates of the bilinear optimal control problem. For example, Yang, Demlow, and Dobrowolski et al. investigated the prior error estimates and superconvergence of optimal control problems for bilinear models and give the optimal $L^2$-norm error estimates and the almost optimal $L^\infty$-norm estimates about the state and co-state variables ^[7,8,31]. Lu investigated a second-order parabolic bilinear optimal control problems and provided a priori error estimates for the finite element solutions of the state equations describing the system ^[24]. Shen and Yang et al. investigated a quadratic optimal control problem governed by a linear hyperbolic integro differential equation and its finite element approximation, a priori estimates have been carried out using the standard functional analysis techniques, and the existence and regularity of the solution are provided by using these estimates. At the same time, some scholars have analysed the posteriori error estimates of the finite element method for bilinear optimal control problems ^[15,26]. Lu, Chen, and Leng et al. discussed the discretisation of Raviart-Thomas mixed finite element for general bilinear optimal control problems, a posteriori error estimates are derived for both the coupled state and the control solutions ^[18,23]. Although bilinear optimal control problems are frequently met in applications, they are much more difficult to handle than linear or nonlinear cases. There is little work on nonlinear optimal control problems.

In this paper, we focus on nonlinear optimal control problem with integral control constraints where we deal with the control via adopting piecewise constant discetization while applying continuous piecewise linear discretization for the state and the co-state, respectively. Then a posteriori error estimates is gained. For the convergence and the quasi-optimality, we prove them relying on quasi-orthogonality and discrete local upper bound. Based on the mild assumption to the initial grids, we obtain the proof of convergence and quasi-optimality by means of the solution operator of nonlinear elliptic equations. Finally, some numerical simulations are provided support for our theoretical analysis.

Here are some notations will be used in this paper. Let $\Omega$ be a bounded Lipschitz domain in $\mathbb{R}^2$ and $\partial\Omega$ denote the boundary of $\Omega$. We use the standard notation $W^{m, q}(\omega)$ with norm $\|\cdot\|_{m, q, \omega}$ and seminorm $|\cdot|_{m, q, \omega}$ to express the standard Sobolev space for $\omega\subset\Omega$. Moreover, we will omit the subscription if $\omega = \Omega$. For $q = 2$, we denote $W^{m, 2}(\Omega)$ by $H^m(\Omega)$ and $\|\cdot\|_m = \|\cdot\|_{m, 2}$. Also for $m = 0 \ \mathrm{and}\ q = 2$, we denote $W^{0, 2}(\omega) = L^2(\omega)$ and $\|\cdot\|_{0, 2, \omega} = \|\cdot\|_{0, \omega}$. Additionally, we observe that $H^1_0(\Omega) = \{v\in H^1(\Omega):v = 0\ \mathrm{on}\ \partial\Omega\}$. Let $\mathcal{T}_{h_0}$ be the initial partition of $\bar{\Omega}$ into disjoint triangles. By newest-vertex bisections for $\mathcal{T}_{h_0}$, we can obtain a class $\mathbb{T}$ of conforming partitions. For $\mathcal{T}_h, \tilde{\mathcal{T}}_h \in \mathbb{T}$, we use $\mathcal{T}_h\subset\mathcal{\tilde{T}}_h$ to indicate that $\tilde{\mathcal{T}}_h$ is a refinement of $\mathcal{T}_h$ and $h_{T} = |T|^{1/2}$, $T$ is the partition diameter. In addition $(\cdot, \cdot)$ denotes the $L^2$ inner product. Beyond that, let $C$ be a constant which independent of grids size, then we use $A\approx B$ to represent $cA\leq B\leq CA$.

The rest of the paper is organized as follows. In Section 2, we introduce the optimal control problem of our interest and obtain a posteriori error estimates. The relevant algorithms are introduced in Section 3. In Section 4, we use quasi-orthogonality and discrete local upper bounds to prove the convergence of the adaptive finite element method, as well as quasi-optimality in Section 5. We provide an adaptive finite element algorithm and some numerical simulations to verify our theoretical analysis in Section 6. Finally, we summarize the results of this paper and develop a plan for future work.

2.   A posteriori error estimates for finite element method

In this paper we mainly enter into meaningful discussions with the following nonlinear optimal control problem governed by nonlinear elliptic equations:

where $y$ is the state variable, $u$ is the control variable, $f$ is a function of the control variable, $\alpha$ is a constant greater than zero, $y_{d}\in L^2(\Omega), \ U_{ad} = \{v:v\in L^2(\Omega), \ \int_{\Omega}v \ dx\geq 0\}$ is a closed convex subset of $U = L^2(\Omega)$ and $\phi(\cdot)\in W^{2, \infty}(-R, R)$ for any $R > 0, \ \phi'(y)\in L^2(\Omega)$ for any $y\in H^1(\Omega), \ \phi'\geq0$. Let $V = H_{0}^1(\Omega)$, we give the weak formulation to deal with state equation, namely, find $y\in V$ such that

where

Then the nonlinear optimal control problem can be restated as follows

It is well known ^[14,20] that the nonlinear optimal control problem has at least one solution $(y, u)$, and that if a pair $(y, u)$ is the solution of the nonlinear optimal control problem, then there is a co-state $p\in V$ such that the triplet $(y, p, u)$ satisfies the following optimality conditions:

Due to coercivity of $a(\cdot, \cdot)$, we define a linear operator $S:L^2(\Omega)\rightarrow H_{0}^1(\Omega)$ such that $S(f+u) = y$ and let $S^*$ be the adjoint of $S$ such that $S^*(y-y_{d}) = p$ where $V_h$ is the continuous piecewise linear finite element space with respect to the partition $\mathcal{T}_h\in \mathbb{T}$. For $\mathcal{T}_h\in \mathbb{T}$, we define $U^h$ as the piecewise constant finite element space with respect to $\mathcal{T}_h$. Let $U_{ad}^h = \{v_h\in U^h:\int_{\Omega}v_hdx\geq 0\}$. Then we derive the standard finite element discretization for the nonlinear optimal control problem as follows:

Similarly the nonlinear optimal control problem (2.6)–(2.7) has a solution $(y_h, u_h)$, and that if a pair $(y_h, u_h)\in V^h\times U^h$ is the solution of (2.6)–(2.7), then there is a co-state $p_h \in V_{h}$ such that the triplet $(y_{h}, p_{h}, u_{h})$ satisfies the following optimality conditions:

Here we define some error indicators we largely and frequently use in this paper where $\eta(\cdot)$ are error indicators and $osc(\cdot)$ represent the data oscillations. For $\mathcal{T}_h\in \mathbb{T}, \ T\in \mathcal{T}_h$, we define

where $u_{h}\in U_{ad}^{h}$, $y_{h}, p_{h}\in V_{h}$, and where $f_{T}$ is $L^2$-projection of $f$ onto piecewise constant space on $T$ and $f_T = \frac{\int_Tf}{|T|}$. For $\omega\subset \mathcal{T}_h$, we have

by which $\eta_{2, \mathcal{T}_h}^2(u_{h}, y_{h}, \omega), \ \eta_{3, \mathcal{T}_h}^2(u_{h}, y_{h}, \omega)$ and $osc^2_{\mathcal{T}_h}(y_{h}-y_{d}, \omega)$ can be denoted similarly.

2.1.   A posteriori error estimates

Also, a reliable and effective a posteriori error estimates for the nonlinear optimal control problem (2.1)–(2.2) which is presented in next Theorem.

Theorem 2.1. For $\mathcal{T}_h\in \mathbb{T}$, let $(y, p, u)$ be the exactsolution of problem (2.3), (2.4) and (2.5) and $(y_{h}, p_{h}, u_{h})$ be the solution of problem (2.8), (2.9) and (2.10) with respect to $\mathcal{T}_h$. Then thereexist constants $c$ and $C$ such that

and

Proof. Step 1. We seek functions $y^h, p^h\in V$ satisfying the following auxiliary problems

It follows that

With the help of the proof of Lemma 3.4 in ^[14] and the lemma 3.4 in ^[21], we have

where $\eta^2_{1, \mathcal{T}_h}(p_{h}, \mathcal{T}_h) = \sum\limits_{T\in\mathcal{T}_h}\eta^2_{1, \mathcal{T}_h}(p_{h}, T)$, and

where $\pi_{h}$ is the $L^2$-projection operator onto piecewise constant space on $\mathcal{T}_h$, $C(\sigma)$ is a universal constant, which depends on $\sigma$, $\sigma$ is an arbitrary positive number, and $C$ is a general universal constant, which can include $C(\sigma)$. Obviously, from (2.15) and (2.16) there holds

Let $e^y = y^{h}-y_{h}$, and $e_I^y = \hat{\pi}_{h}e^y$, where $\hat{\pi}_{h}$ is the average interpolation operator defined in Lemma 3.2 of ^[14], then we can obtain

where $\sigma$ is an arbitrary positive number. The definition of $\eta_2$ will be given later on. Then let $\sigma = \frac{c}{2C}$, we have

Similarly, let $e^p = p^{h}-p_h$ and $e_I^p$ be the average interpolation of $e^p$, then we can get

in which we apply the embedding theorem $\|v\|_{0, 4}\leq C\|v\|_1$ (see ^[6]) and the property: $\|p_{h}\|_1\leq C$, $C(\sigma)\|y^{h}-y_{h}\|_0^2+C(\sigma)\|p_{h}\|_1^2 \|y^{h}-y_{h}\|_0^2\leq C(\sigma)\|y^{h}-y_{h}\|_0^2+C(\sigma)C \|y^{h}-y_{h}\|_0^2\leq C\|y^{h}-y_{h}\|_0^2$. The definition of $\eta_3$ will be given later on. Absolutely, we can obtain

Hence we gain

By the triangular inequality we obtain that

In connection with what we discussed above, we have

Step 2. Now we are in the position to get the lower bound. Consulting to the proof of Lemma 3.6 in ^[14] and Lemma 3.4 in ^[21], we conclude that

Proof. It is easily seen that

Since $u+p = \max(0, \bar{p}) = const$, hence

such that

This completes the proof.

According to ^[28], we empoly the standard bubble function technique to estimate error indicators $\eta_{2, \mathcal{T}_h}(u_h, y_h, \mathcal{T}_h)$ and $\eta_{3, \mathcal{T}_h}(y_h, p_h, \mathcal{T}_h)$. Similar to Chapter 7.2 in ^[22], there exists polynomials $w_T\in H_0^1(T)$ and $w_{\partial T}\in H_0^1({\partial T\backslash\partial\Omega})$ such that

and apparently

By using the Schwarz inequality, it follows from (2.18), (2.20), and (2.21) that

where $\sigma$ is an arbitrary positive number and $\phi(\cdot)\in W^{2, \infty}(\Omega)$ has been used. Then let $\sigma = \frac{1}{2C}$ and we have

Next, it follows from (2.19), (2.22), and (2.23) that

where $\phi(\cdot)\in W^{2, \infty}(\Omega)$ has been used. Absolutely, we can deduce that

Then we have

In connection with (2.24) and (2.25) we are easy to gain

It can also be deduced that

Above-mentioned results tell the proof of Theorem 2.1 is accomplished.

3.   The algorithms

In this section, we introduce two related algorithms as follows:

Algorithm 3.1. Adaptive finite element algorithm for nonlinear optimal controlproblems:

(0) Given an initial grids $\mathcal{T}_{h_0}$ and construct finiteelement space $U_{ad}^{h_0}$ and $V_{h_0}$. Select marking parameter $0 < \theta\leq 1$ and set $k: = 0$.

(1) Solve the discrete nonlinear optimal control problem (2.8)–(2.10), then obtain approximate solution $(u_{h_k}, y_{h_k}, p_{h_k})$ with respect to $\mathcal{T}_{h_k}$.

(2) Compute the local error estimator $\eta_{\mathcal{T}_{h_k}}(T)$ for all $T\in\mathcal{T}_{h_k}$.

(3) Select a minimal subset $\mathcal{M}_{h_k}$ of $\mathcal{T}_{h_k}$ such that

where $\eta_{\mathcal{T}_{h_k}}^2(\omega) = \eta_{1, \mathcal{T}_h}^2(p_h, \omega)+\eta_{2, \mathcal{T}_h}^2(u_{h}, y_{h}, \omega)+\eta_{3, \mathcal{T}_h}^2(y_{h}, p_{h}, \omega)$ for all $\omega\subset\mathcal{T}_{h_k}$.

(4) Refine $\mathcal{M}_{h_k}$ by bisecting $b\geq 1$ times inpassing from $\mathcal{T}_{h_k}$ to $\mathcal{T}_{h_{k+1}}$ andgenerally additional elements are refined in the process in order toensure that $\mathcal{T}_{h_{k+1}}$ is conforming.

(5) Solve the discrete nonlinear optimal control problem (2.8)–(2.10), then obtain approximate solution $(u_{h_{k+1}}, y_{h_{k+1}}, p_{h_{k+1}})$ with respect to $\mathcal{T}_{h_{k+1}}$.

(6) Set $k = k+1$ and go to step (2).

Algorithm 3.2. Given an initial control $u_h^0\in U_{ad}^{h}$, then seek $(y_h^k, p_h^k, u_h^k)$ such that

for $k = 1, 2, \cdots$, and apparently

where $\mathcal{P}_h$ is the $L^2$-projection from $L^2(\Omega)$ to $U^h$ and $\bar{p}_h^k = \frac{\int_\Omega p_h^k}{|\Omega|}$.

4.   Convergence analysis for adaptive finite element method

In this section we first consider the nonlinear elliptic equations as follows:

where $f\in L^2(\Omega)$. We introduce the quantity $\mathcal{J}(h)$ in view of the idea in ^[30] as follows:

where $S$ is the solution operator for nonlinear elliptic equations. Obviously, $\mathcal{J}(h)\ll 1$ for $h\in (0, h_{0})$ if $h_{0}\ll 1$. Hereinafter there holds the Lemma about the quantity referring to ^[10,30].

Lemma 4.1. For each $f\in L^2(\Omega)$, there exists a constant $C$ such that

and

where $S_{h}$ is the discretesolution operator for nonlinear elliptic equations.

Here, Lemma 4.1 is a preparation for Lemma 4.4.

4.1.   Local perturbation property

According to ^[19], the local perturbation property plays an important role for the proof of the convergence. It impels us to combine the sum of the error estimates.

Lemma 4.2. For $\mathcal{T}_h\in \mathbb{T}, \ T\in \mathcal{T}_h$, let $u_{h_1}, u_{h_2}\in U_{ad}^h, \ y_{h_1}, y_{h_2}, p_{h_1}, p_{h_2}\in V_{h}$, we have

Proof. Step 1. According to ^[16,21] we have

By adopting the inverse estimates and (4.8) we obtain

where $\omega_{T}$ denotes the patch of elements that share an edge with $T$. By the definition of $\eta_{1, \mathcal{T}_h}(p_h, T)$ and (4.10) we can deduce that

Then we have

where the proof of (4.4) is finished.

Step 2. By the definition of $\eta_{2, \mathcal{T}_h}(u_h, y_h, T)$ and (4.9), we can deduce that

where $\phi(\cdot)\in W^{2, \infty}(\Omega)$ has been used and the proof of (4.5) is finished.

Step 3. By the definition of $\eta_{3, \mathcal{T}_h}(y_h, p_h, T)$ and (4.10) we can derive that

where $\phi(\cdot)\in W^{2, \infty}(\Omega)$ has been used and the proof of (4.6) is finished. \\Step 4. Similarly, we have

In brief, Lemma 4.2 is proved.

4.2.   Error reduction

The authors in ^[25] demonstrate an error reduction provided the current errors are larger than the desired errors, that is to say, the errors may not be reduced in the process of coarse grids refinement before introducing a node of the refined grids inside each marked element while D$\ddot{\mathrm{o}}$rfler proves a similar result assumption ^[9].

Lemma 4.3. Let $\mathcal{T}_h\subset\tilde{\mathcal{T}_h}$ for $\mathcal{T}_h, \tilde{\mathcal{T}_h}\in\mathbb{T}.\mathcal{M}_h\subset\mathcal{T}_h$ denotes the set of elements whichare marked from $\mathcal{T}_h$ to $\tilde{\mathcal{T}}_h$. Then for $u_h\in U_{ad}^h, \ \tilde{u}_h \in U_{ad}^{\tilde{h}}, \ y_h, p_h\in V_{h}, \ \tilde{y}_h, \tilde{p}_h\in V_{\tilde{h}}$ and any $\delta, \delta_{1}\in(0, 1]$, we have

and

and

and

where $\lambda = 1-2^{-\frac{b}{2}}, \ h_{0} = \max\limits_{{T\in\mathcal{T}_{h_{0}}}}h_{T}$ and $\mathcal{R}_h$ denotes the set ofelements which are refined from $\mathcal{T}_h$ to $\tilde{\mathcal{T}_h}$.

Proof. Step 1. Applying the Young's inequality with parameter $\delta_{1}$ and (4.4) we get

Note that $T$ will be bisected at least one time for the element $T\in \mathcal{R}_h\subset\mathcal{T}_h$, then we have

For $T\in \mathcal{T}_h\backslash\mathcal{R}_h$, we gain

In connection with the above estimates we demonstrate that

which illustrates (4.11) has been proved.

Step 2. Employing the Young's inequality with parameter $\delta$ and (4.5) we obtain

where it is similar to (4.15). Then let $\tilde{\mathcal{T}}_{h_{T'}} = \{T\in \tilde{\mathcal{T}_h}:T\subset T'\}$ where $T'\in\mathcal{M}_h$ is a marked element. For arbitrary $p_h\in V_{\mathcal{T}_h}\subset V_{\tilde{\mathcal{T}_h}}$, we find the jump $[\nabla p_h] = 0$ on the interior sides of $\cup\tilde{\mathcal{T}}_{h_{T'}}$. Suppose $b$ is the number of bisections, we can deduce that

due to refinement by bisection, then we obtain

It is easy to find that

for any $T\in\mathcal{T}_h\backslash\mathcal{M}_h$. In connection with the above estimates we expound that

As has been said, we can achieve (4.13) connecting with the above estimates to which the proof of (4.12) is similar.

Step 3. For arbitrary $T\in\mathcal{T}_h\cap\tilde{\mathcal{T}_h}$ via using (4.7) and the Young's inequality we obtain that

Obviously, (4.14) can be got by summing the above inequality over $T\in\mathcal{T}_h\cap\tilde{\mathcal{T}_h}$. To sum up, Lemma 4.3 is proved.

4.3.   Quasi-orthogonality

As to the proof of the convergence, one of the main obstacle is that there do not have the orthogonality while it is vital to prove the convergence. Thus getting back to the second place we transfer proof of the quasi-orthogonality. The latter is popularly adopted in the adaptive mixed and the nonconforming adaptive finite element methods ^[19]. Apparently it is true for the following basic relationships with $\mathcal{T}_{h_k}, \mathcal{T}_{h_{k+1}}\in\mathbb{T}$ and $\mathcal{T}_{h_k}\subset\mathcal{T}_{h_{k+1}}$,

where $(u, y, p)$ are the solution of (2.3)–(2.5), $(u_{h_k}, y_{h_k}, p_{h_k})$ and $(u_{h_{k+1}}, y_{h_{k+1}}, p_{h_{k+1}})$ are the solution of (2.8)–(2.10) with respect to $\mathcal{T}_{h_k}$ and $\mathcal{T}_{h_{k+1}}$, respectively. Accordingly we have the quasi-orthogonality below.

Lemma 4.4. For $\mathcal{T}_{h_k}, \mathcal{T}_{h_{k+1}}\in\mathbb{T}$ and $\mathcal{T}_{h_k}\subset\mathcal{T}_{h_{k+1}}$, we have

and

and

Proof. Step 1. It follows from Lemma 4.3 in ^[19] that for $U_{ad}^{h_k}\subset U_{ad}^{h_{k+1}}$ we have

where $\mathcal{R}_h$ is the set of elements which are refined from $\mathcal{T}_{h_k}$ to $\mathcal{T}_{h_{k+1}}$.

For the right-hand first item of (4.22) we let $\zeta_h\in H_{0}^1(\Omega)$ be the solution of the following problem based on Lemma 4.1

Hence we can get

From the proof of Lemma 3.3 in ^[19], we infer that

and

and

Similar to Lemma 4.6 in ^[4], we derive that

In connection with the above estimates we conclude that

For the third term at the right end of the above inequality, we let $\varphi_h\in H_{0}^1(\Omega)$ be the solution of the following problem

According to the standard duality theory, we can deduce that

where $\varphi_{h_k}$ is the standard finite element estimate of $\varphi_h$ with respect to $V_{\mathcal{T}_{h_k}}$. So we have

As mentioned above, we can get

In combination with (4.22) and above inequality, we deduct that

It is easy to derive the desired result (4.19) with the help of (4.16) and (4.23).

Step 2. Our task now is to prove (4.20) and so is (4.21). Obviously we have

By using the Cauchy inequality, we obtain

It is easy to derive the desired result (4.20) with the assistance of (4.17) and (4.23)–(4.25).

4.4.   Convergence

For $\mathcal{T}_{h_k}\in \mathbb{T}$, we will denote $U_{ad}^{h_{k}}, \ V_{h_{k}}$ and the solution $(u_{h_{k}}, y_{h_{k}}, p_{h_{k}})$ of (2.8)–(2.10) with respect to $\mathcal{T}_{h_{k}}$ by $U_{ad}^{h_k}, \ V_{h_k}$, and $(u_{h_k}, y_{h_k}, p_{h_k})$ and we define some notations before we prove the convergence of Algorithm 3.1 as follows:

for $\omega\subset\Omega$.

Theorem 4.1. Let $(\mathcal{T}_{h_k}, U_{ad}^{h_k}, V_{h_k}, u_{h_k}, y_{h_k}, p_{h_k})$ be the sequence of grids, finite element spaces and discretesolutions produced by the Algorithm 3.1. Then there existconstants $\gamma_{1} > 0, \ \gamma_{2} > 0$ and $\alpha\in (0, 1)$, onlydepending on the shape regularity of initial grids $\mathcal{T}_{h_0}, \ b$, and the marking parameter $\theta\in(0, 1]$, such that

provided $h_{0}\ll 1$.

Proof. We get the following results from Theorem 2.1, Lemma 4.3 and Lemma 4.4,

where $\mathcal{R}_{h_k}$ is the set of elements which are refined from $\mathcal{T}_{h_k}$ to $\mathcal{T}_{h_{k+1}}$. Applying the upper bound in Theorem 2.1, (4.29) can be simplified into

Multiplying (4.28) and (4.31) with $\tilde{\gamma_{2}}$ and $\tilde{\gamma_{1}}$, respectively, and adding the results to (4.30) yields

If $\tilde{\gamma_{1}}$ is such that

and one chooses $\tilde{\gamma_{2}}$ such that

then we have

where

By using the marking strategy in Algorithm 3.1 and the upper bound in Theorem 2.1 to arrive at

where $\beta\in(0, 1)$. Then we deduce that

where

As long as $\delta < 1$ and $\beta$ is small enough, $\alpha_{1}\in(0, 1)$ can be guaranteed. To facilitate judgment, we transfer the following adjustments to the above formula:

It is absolutely clear that

if $h_{0}\ll 1$ and $\delta_{1}$ is sufficiently small. Then consider (4.32) to deduce that

which can say

if $h_{0}\ll 1$ and $\delta$ is sufficiently small. Therefore if choose $\alpha = \max\{\alpha_{1}, \alpha_{2}, \alpha_{3}\}$, we can derive the expected results.

Theorem 4.2. Let $(\mathcal{T}_{h_k}, U_{ad}^{h_k}, V_{h_k}, u_{h_k}, y_{h_k}, p_{h_k})$ be the sequence of grids, finite element spaces and discretesolutions produced by the Algorithm 3.1 and the conditionsof Theorem 4.1 keep. Then we have

Proof. It is obviously true combining Theorem 2.1 and Theorem 4.1.

5.   Quasi-optimality for adaptive finite element algorithm

In this section we consider the quasi-optimality for the adaptive finite element method. Firstly we give the notations interpretation. For $\mathcal{T}_h, \mathcal{T}_{h_1}, \mathcal{T}_{h_2}\in \mathbb{T}$, let $\#\mathcal{T}_h$ be the number of elements in $\mathcal{T}_h$, and $\mathcal{T}_{h_1}\oplus\mathcal{T}_{h_2}$ be the smallest common conforming refinement of $\mathcal{T}_{h_1}$ and $\mathcal{T}_{h_2}$ and satisfies ^[4,27]

According to ^[19], we need defining a function approximation class

where

and

5.1.   Localized upper bound

To illustrate the quasi-optimality of the adaptive finite element method, we need a local upper bound on the distance between nested solutions ^[4], since the error of this method can only be estimated by using the indicators of refined elements without a buffer layer.

Lemma 5.1. For $\mathcal{T}_h, \tilde{\mathcal{T}}_h\in \mathbb{T}$ and $\mathcal{T}_h\subset\tilde{\mathcal{T}}_h$, let $\mathcal{R}_h$ bethe set of refined elements from $\mathcal{T}_h$ to $\tilde{\mathcal{T}}_h$. Let $(u_{h}, y_{h}, p_{h})$ and $(\tilde{u}, \tilde{y}, \tilde{p})$ be the solutions of (2.8)–(2.10) with respect to $\mathcal{T}_h$ and $\tilde{\mathcal{T}}_h$, respectively. Then there exists a constant $C$, depending on the shape regularity of initial grids $\mathcal{T}_{h_{0}}$ and $b$ such that

where

Proof. From (4.23) of Lemma 4.4, we have

By Lemma 4.6 in ^[4], we deduce that

The corresponding result is attained by the similar method above

In connection with (5.3), (5.4) and (5.5) we can derive (5.2).

5.2.   Dörfler property

Dörfler introduced a crucial marking and proved strict energy error reduction for the Laplacian provided the initial grids $\mathcal{T}_{h_0}$ satisfying a mild assumption ^[9]. If the sum of errors satisfy suitable error reductions, the error indicators on the coarse grids must satisfy a Dörfler property on the refined one ^[19].

Lemma 5.2. Assume that the marking parameter $\theta\in(0, \theta^*)$, where

For $\mathcal{T}_h, \tilde{\mathcal{T}}_h\in\mathbb{T}$ and $\mathcal{T}_h\subset\tilde{\mathcal{T}}_h$, let $(u_{h}, y_{h}, p_{h})$ and $(\tilde{u}, \tilde{y}, \tilde{p})$ be thesolutions of (2.8)–(2.10) with respect to $\mathcal{T}_h$ and $\tilde{\mathcal{T}}_h$, respectively. If

is satisfied for $\mu: = \frac{1}{2}\Big(1-\frac{\theta}{\theta^*}\Big)$. Then, the set $\mathcal{R}_h$ of elements which are refined from $\mathcal{T}_h$ to $\tilde{\mathcal{T}}_h$ satisfies the Dörfler property

where

for $\omega\subset\mathcal{T}_h$ and $e_{\tilde{\mathcal{T}}_h}^2$, $osc_{\tilde{\mathcal{T}}_h}^2(\tilde{\mathcal{T}}_h)$ similarly defined.

Proof. By the lower bound in Theorem 2.1 and (5.6) to obtain that

It is well-known that there exists the fundamental relationships:

Hence we can get the following result form Lemma 5.1

For $T\in \mathcal{T}_h\cap\tilde{\mathcal{T}}_h$, we can get the following result from (4.14) of Lemma 4.3

According to Remark 2.1 in ^[4], we get the following result as the indicator $\eta_{\mathcal{T}_h}(T)$ dominates oscillation $osc_{\mathcal{T}_h}(T)$

for all $T\in\mathcal{R}_h$. Then in connection with (5.7), (5.8) and (5.9) one obtains

where

Lemma 5.3. Let $(u, y, p)$ and $(\mathcal{T}_{h_k}, U_{ad}^{h_k}, V_{h_k}, u_{h_k}, y_{h_k}, p_{h_k})$ be the solution of (2.3)–(2.5) and the sequence ofgrids, finite element spaces and discrete solutions produced by Algorithm 3.1, respectively. Assume that the markingparameter $\theta$ satisfies the condition in Lemma 5.2, then the following estimate is valid

if $(u, y, p, y_{d}, f)\in \mathcal{A}^s$.

Proof. Let $\epsilon^2: = \mu N^{-1}(e_{\mathcal{T}_{k}}^2+osc_{\mathcal{T}_{k}}^2(\mathcal{T}_{k}))$, where $N$ shall be produced in the proof of (5.13) and $\mu$ is defined in Lemma 5.2. Because of $(u, y, p, y_{d}, f)\in \mathcal{A}^s$, there exists a $\mathcal{T}_{h_\epsilon}\in\mathbb{T}$ and a $(u_{h_\epsilon}, y_{h_\epsilon}, p_{h_\epsilon})\in U_{ad}^{h_{\epsilon}}\times V_{h_{\epsilon}}\times V_{h_{\epsilon}}$ such that

and

Let $(u_{h_*}, y_{h_*}, p_{h_*})$ be the solution of (2.8)–(2.10) with respect to $\mathcal{T}_{h_*}$, where $\mathcal{T}_{h_*} = \mathcal{T}_{h_\epsilon}\oplus\mathcal{T}_{h_k}$ is the smallest common refinement of $\mathcal{T}_{h_\epsilon}$ and $\mathcal{T}_{h_k}$. In the following we will prove the following inequality firstly

where

According to the principle of adding one item and subtracting one item, then we have

With the help of the Young's inequality one obtains

and so are $a(y-y_{h_*}, y_{h_*}-y_{h_\epsilon})$ and $a(p-p_{h_*}, p_{h_*}-p_{h_\epsilon})$. Hence in connection with what we get above we deduce that

From Remark 2.1 in ^[4] and (4.14) in Lemma 4.3 with $\mathcal{T}_h = \tilde{\mathcal{T}}_h = \mathcal{T}_{h_*}, \ y = y_{h_*}$ and $\tilde{y} = y_{h_\epsilon}$, we obtain that

For any $T'\in\mathcal{T}_{h_\epsilon}$, let $\mathcal{T}_{h_{T'}}: = \{T\in \mathcal{T}_{h_*}:T\in T'\}$. From the proof of Lemma 4.3 in ^[19], we derive that

And then we can get

Combining (5.14)–(5.16) to obtain (5.13) and using (5.12) and the definition of $\epsilon^2$, we have

It is true for the following result from Lemma 5.2

Combining (5.11), (5.17) and the definition of $\epsilon^2$ to derive the desired result (5.10).

5.3.   Quasi-optimality

The following consequence is the result of previous estimates where the truth is that the number of elements is a dwindle for the errors. Namely, if a given adaptive method is used to approximate the exact solution at a certain convergence rate, the iteratively constructed grids sequence will achieve this rate until a constant factor.

Theorem 5.1. Let $(u, y, p)$ and $(\mathcal{T}_{h_k}, U_{ad}^{h_k}, V_{h_k}, u_{h_k}, y_{h_k}, p_{h_k})$ be the solution of (2.3)–(2.5) and the sequence ofgrids, finite element spaces and discrete solutions produced by Algorithm 3.1, respectively. Assume that $\mathcal{T}_{h_0}$ satisfies the condition (b) of Section 4 in ^[27]. Let $(u, y, p, y_{d}, f)\in \mathcal{A}_{s}$, then we have

provided $h_{0}\ll 1$.

Proof. It follows from Theorem 2.1 and Lemma 5.3 that

From Lemma 2.3 in ^[4], we have

With the assistance of Lemma 5.3 and (5.20) to gain

where

Then it follows from (5.19), (5.21) and Theorem 4.1 that

which tells the proof of Theorem 5.1.

6.   Numerical experiments

In this section, we firstly present an adaptive finite element and then give the adaptive iteration method where the purpose is to provide empirical analysis for our theory.

Example 1. We consider the nonlinear optimal control problem governed by nonlinear elliptic equations subject to the state equation

where we choose $\alpha = 1$ and $\Omega = [0, 1]\times[0, 1]$ and apparently exact solution

By simple calculation we have $\int_\Omega pdx = -\frac{4}{\pi}$ which satisfies $u\in U_{ad}$.

We choose 15 adaptive loops for Example 1, then we plot the profiles of the exact state and the numerical state on adaptively refined grids with $\theta = 0.5$ in Figure 1. It is easy to see that the solution is smooth, but we can find larger gradients in some regions, hence comparing with uniform refinement, adaptive finite element method can provide smaller error. In Figure 2, we provide the triangle refined grids after 6 and 12 adaptive iterations of Algorithm 3.1.

In Figure 3, we plot the convergence history for the errors where the left is adaptive refinement $(\theta = 0.5)$ and the right is uniform refinement $(\theta = 1)$. We can find it very intuitively when we provide the optimal convergence rate slope $-1$ via adopting the linear finite elements, the error reduction can be observed.

Example 2. We consider the same nonlinear optimal control problem as Example 1 with $\alpha = 0.1$, $\Omega = (-1, 1)\times(0, 1)\cup(-1, 0)\times(-1, 0]$ and apparently exact solution

where $u\in U_{ad}$ can be guaranteed.

Comparing with Example 1, we provide some plots concerning with 21 adaptive loops for Example 2 with $\theta = 0.5$. In Figure 4, it is more easy to say that the solution is smooth while the large gradients can be found in some regions illustrating that adaptive refinement can obtain smaller errors than the uniformly refinement where we offer the error estimate graphs to explain. In Figure 5, the left plot tells us the error estimates on adaptive refinement $(\theta = 0.5)$ and the right shows the error estimates on uniform refinement $(\theta = 1)$. With the slope $-1$ being the optimal convergence rate expected, we see the error reduction from Figure 5. Meanwhile we can also find that the convergence order of the total-error and the error estimate indicators are approaching to straight line slope $-1$ in which they are roughly parallel where it is showed that the posteriori error estimates we obtained in Section 2 are reliable.

In Figure 6, we show the adaptive grids after 9 and 19 adaptive iterations for Example 2 of 21 adaptive loops with $\theta = 0.5$. We can find that the grids are concentrated on the regions where the solutions have larger gradients. Only can we note that reduced orders are observed for the uniform refinement because of the singularity of the solutions.

7.   Conclusions and future expectation

In this paper, we first study the adaptive finite element method for nonlinear optimal control problems and give the corresponding adaptive algorithm. To evaluate the adaptive finite element method, we obtain the a posteriori error estimates for the nonlinear elliptic equations with upper and lower bound convergence and optimality, which are also important indicators for evaluating the algorithms. Therefore, we prove that the sum of the posterior errors of the control, state and covariance variables are convergent, as shown in Theorem 4.2. Based on the local upper bound, we prove the quasi-optimality of the proposed adaptive algorithm, see Section 5. To verify our theoretical analysis, we finally provide some numerical simulations. In previous research papers, the finite element methods for linear optimal control problems were studied. Our innovation is to extend the method of linear optimal control problems to a series of nonlinear optimal control problems.

There are a lot of problems which can not be tackled, such as the $L^2-L^2$ posteriori error estimates for nonlinear elliptic equations as well as the convergence and quasi-optimality for nonlinear parabolic equations. Furthermore, we note that the analysis in this paper can be generalized to common nonlinear parabolic problems and boundary problems, and we will work on these problems.

Acknowledgments

This work is supported by National Science Foundation of China (11201510), National Social Science Fund of China (19BGL190), Chongqing Research Program of Basic Research and Frontier Technology (cstc2019jcyj-msxmX0280), General project of Chongqing Natural Science Foundation (cstc2021jcyj-msxmX0949), Scientific and Technological Research Program of Chongqing Municipal Education Commission (KJZD-K20200120), Chongqing Key Laboratory of Water Environment Evolution and Pollution Control in Three Gorges Reservoir Area (WEPKL2018YB-04), Science and technology research project of Chongqing Education Commission (KJQN202101212), Research Center for Sustainable Development of Three Gorges Reservoir Area (2019sxxyjd07), Guangdong Basic and Applied Basic Research Foundation of Joint Fund Project (2021A1515111048), and Guangdong Province Characteristic Innovation Project (2021WTSCX120).

Conflict of interest

The authors declare that they have no competing interests.