
Citation: Geraldo Ceni Coelho, Glaci Benvenuti-Ferreira, Jorge Schirmer, Osório A. Lucchese. Survival, growth and seed mass in a mixed tree species planting for Atlantic Forest restoration[J]. AIMS Environmental Science, 2016, 3(3): 382-394. doi: 10.3934/environsci.2016.3.382
[1] | Vincenzo Capasso, Sebastian AniȚa . The interplay between models and public health policies: Regional control for a class of spatially structured epidemics (think globally, act locally). Mathematical Biosciences and Engineering, 2018, 15(1): 1-20. doi: 10.3934/mbe.2018001 |
[2] | Dongmei Li, Tana Guo, Yajing Xu . The effects of impulsive toxicant input on a single-species population in a small polluted environment. Mathematical Biosciences and Engineering, 2019, 16(6): 8179-8194. doi: 10.3934/mbe.2019413 |
[3] | Rajalakshmi Manoharan, Reenu Rani, Ali Moussaoui . Predator-prey dynamics with refuge, alternate food, and harvesting strategies in a patchy habitat. Mathematical Biosciences and Engineering, 2025, 22(4): 810-845. doi: 10.3934/mbe.2025029 |
[4] | Zeyan Yue, Sheng Wang . Dynamics of a stochastic hybrid delay food chain model with jumps in an impulsive polluted environment. Mathematical Biosciences and Engineering, 2024, 21(1): 186-213. doi: 10.3934/mbe.2024009 |
[5] | Sheng Wang, Lijuan Dong, Zeyan Yue . Optimal harvesting strategy for stochastic hybrid delay Lotka-Volterra systems with Lévy noise in a polluted environment. Mathematical Biosciences and Engineering, 2023, 20(4): 6084-6109. doi: 10.3934/mbe.2023263 |
[6] | Martin Bohner, Sabrina Streipert . Optimal harvesting policy for the Beverton--Holt model. Mathematical Biosciences and Engineering, 2016, 13(4): 673-695. doi: 10.3934/mbe.2016014 |
[7] | Yuanpei Xia, Weisong Zhou, Zhichun Yang . Global analysis and optimal harvesting for a hybrid stochastic phytoplankton-zooplankton-fish model with distributed delays. Mathematical Biosciences and Engineering, 2020, 17(5): 6149-6180. doi: 10.3934/mbe.2020326 |
[8] | Fangfang Zhu, Xinzhu Meng, Tonghua Zhang . Optimal harvesting of a competitive n-species stochastic model with delayed diffusions. Mathematical Biosciences and Engineering, 2019, 16(3): 1554-1574. doi: 10.3934/mbe.2019074 |
[9] | Ali Moussaoui, Pierre Auger, Christophe Lett . Optimal number of sites in multi-site fisheries with fish stock dependent migrations. Mathematical Biosciences and Engineering, 2011, 8(3): 769-783. doi: 10.3934/mbe.2011.8.769 |
[10] | Maoxiang Wang, Fenglan Hu, Meng Xu, Zhipeng Qiu . Keep, break and breakout in food chains with two and three species. Mathematical Biosciences and Engineering, 2021, 18(1): 817-836. doi: 10.3934/mbe.2021043 |
In today's world of industrial pollution, toxicants are pervading the air, ecological problems have become increasingly prominent, and environmental pollution has become a major problem. SARS, Ebola virus, AIV, H1N1 influenza, and COVID-19 are threatening the ecological balance as well as the survival of human beings and other creatures. It is necessary to study the effects of toxicants on the ecosystem. Hallam et al. proposed using a dynamic methodology to examine ecotoxicology in [1,2,3]. They established a model of the interaction between toxicants and population, and provided sufficient conditions for the persistence and extinction of a population stressed by a toxicant. Researchers have been studying ecotoxicology since the 1980s, and a large amount of literature has been devoted to problems in the area. Luo et al. [4,5,6] studied a new age-dependent model of toxicant population in an environment with a small capacity for toxicants. The threshold between persistence in the mean and extinction was obtained for each species in a polluted environment [7,8]. The results of thresholding in [9] were then extended to a stochastic Lotka–Volterra cooperative model for n-species.
The effects of environmental pollution on biological population, the dynamical behavioral analysis of ecosystem models, and the control problem have attracted the attention of many scholars. A large number of ecological studies have shown that differences in individual size structures have a more important effect on population development than those in the age structure. This kind of model has achieved remarkable results through theory, numerical calculations, and experimental methods. Among them, the maximum principle and bang-bang structure of the optimal control were established in [10] for a size-structured forestry model with the benefit of carbon sequestration. In [11,12,13], several authors investigated size-structured population models with separable mortality rates. He and Liu [14] studied an optimal birth control problem for a size-structured population model that takes fertility as the control variable. Later, they studied a nonlinear egg-juvenile-adult model in which eggs and juveniles were structured by age while adults were structured by size [15]. In [16], Liu et al. proposed a unified size-structured PDE model for the growth of metastatic tumors and provided several numerical examples. In [17], Mansal et al. studied the dynamics of the predator and the prey in a fishery model by using the fractional-order derivative. In [18], Thiao et al. discussed the fractional optimal economic control of a continuous game theory model described by the fractional-order derivative. Optimal harvesting problems for a size-structured population model were analyzed in [19,20]. However, most of these studies have focused on a single species, and few have examined interactions among species, especially predator-prey models of two species. In addition, due to the influence of seasonal changes and other factors, the living environment of populations often undergoes periodic changes. Research on optimal harvesting problems dependent on the model of individual size in a periodic environment has been reported in [21,22]. In this study, we establish mathematical population models to consider the size factor to render them reasonable.
Few studies to date have examined optimal control problems of size-dependent population models and periodic effects in a polluted environment. In order to bridge this gap, we discuss optimal harvesting for a periodic, n-dimensional food chain model dependent on size structure in a polluted environment, and detail a simulation and an example based on the population density.
The remainder of this paper is organized as follows: Section 2 describes a population model with size structure in a polluted environment. Its well-posedness is proved in Section 3 and its optimality conditions are established in Section 4. In Section 5, we discuss the existence and uniqueness of the optimal control pair. Section 6 is devoted to approximating the T-periodic solution by a numerical algorithm. A discussion of the results and the conclusions of this study are provided in Section 7.
In [1,2,3], Hallam et al. proposed the following dynamic population model with toxicant effects:
{dxdt=x[r0−r1C0−fx],dC0dt=kCE−gC0−mC0,dCEdt=−k1CEx+g1C0x−hCE+u, | (2.1) |
where x=x(t) is the population biomass at time t; C0=C0(t) is the concentration of toxicants in the organism at time t; CE=CE(t) is the concentration of toxicants in the environment of the population at time t. The exogenous rate of input of toxicants into the environment was represented by u. They investigated the persistence and extinction of a population in a polluted environment.
Luo et al. [23] studied optimal birth control for the following age-dependent n-dimensional food chain model:
{∂p1∂t+∂p1∂a=−μ1(a,t)p1−λ1(a,t)P2(t)p1,∂pi∂t+∂pi∂a=−μi(a,t)pi+λ2i−2(a,t)Pi−1(t)pi−λ2i−1(a,t)Pi+1(t)pi,∂pn∂t+∂pn∂a=−μn(a,t)pn+λ2n−2(a,t)Pn−1(t)pn,pi(0,t)=βi(t)∫a2a1mi(a,t)pi(a,t)da,pi(a,0)=pi0(a),Pi(t)=∫a+0pi(a,t)da,(a,t)∈Q, | (2.2) |
where Q=(0,a+)×(0,+∞), [a1,a2] is the fertility interval. pi(a,t) represents the density of the ith population of age a at time t, and a+ is the life expectancy of individuals. The control variable βi(t) is the average fertility of the ith population. The vitality rate μi(a,t) denotes the average mortality of the ith population. λk(a,t) is the interaction coefficient, and mi(a,t) is the ratio of females in the ith population. Maximum principles for the control problems with free terminal states, infinite horizons and target sets have been derived.
By combining (2.1) and (2.2), we propose the following periodic, n-dimensional food chain model with size structure in a polluted environment:
{∂p1∂t+∂(V1(x,t)p1)∂x=f1(x,t)−μ1(x,c10(t))p1−λ1(x,t)P2(t)p1−u1(x,t)p1,∂pi∂t+∂(Vi(x,t)pi)∂x=fi(x,t)−μi(x,ci0(t))pi+λ2i−2(x,t)Pi−1(t)pi−λ2i−1(x,t)Pi+1(t)pi−ui(x,t)pi,i=2,3,…n−1,∂pn∂t+∂(Vn(x,t)pn)∂x=fn(x,t)−μn(x,cn0(t))pn+λ2n−2(x,t)Pn−1(t)pn−un(x,t)pn,dci0dt=k1ce(t)−g1ci0(t)−mci0(t),dcedt=−k2ce(t)n∑i=1Pi(t)+g2n∑i=1ci0(t)Pi(t)−h1ce(t)+v(t),Vi(0,t)pi(0,t)=∫l0βi(x,ci0(t))pi(x,t)dx,pi(x,t)=pi(x,t+T),Pi(t)=∫l0pi(x,t)dx,(x,t)∈Q, | (2.3) |
where Q=(0,l)×R+,l∈R+ is the maximal size of an individual in the population, T∈R+ is the period of evolution of the population. The meanings of the other parameters are as follows: pi(x,t): the density of the ith population of size x at time t; ci0(t): the concentration of toxicants in the ith population; ce(t): the concentration of toxicants in the environment; Vi(x,t): the growth rate of the ith population; μi(x,ci0(t)),βi(x,ci0(t)): the mortality and fertility rates of the ith population, respectively; λk(x,t): the interaction coefficient (k=1,2,…,2n−2); v(t): the input rate of exogenous toxicants; Pi(t): total number of individuals in the ith population; fi(x,t): the immigration rate of the ith population; ui(x,t): function of the harvesting efforts of the ith population. k1,g1,m,h1,k2, and g2 are nonnegative constants. For the sake of convenience, if there is no special description, i=1,2,…,n. The model represents the population dynamics, and couples toxicants with the population. Moreover, it is more realistic for the fertility and mortality rates to depend on the size and the concentration of toxicants in an organism.
This paper considers the following objective functional:
max{J(u,v):u=(u1(x,t),u2(x,t),…,un(x,t),v=v(t),(u,v)∈Ω)}, | (2.4) |
where
J(u,v)=n∑i=1∫T0∫l0wi(x,t)ui(x,t)pi(x,t)dxdt−12n∑i=1∫T0∫l0ciu2i(x,t)dxdt−12∫T0cn+1v2(t)dt, |
represents the revenue obtained from harvesting less the costs of harvesting and curbing environmental pollution. The weight function wi(x,t) is the selling price factor of an individual belonging to the ith population. The positive constants ci and cn+1 are the cost factors of the ith harvested population and the curbing of environmental pollution, respectively. Therefore, the objective function represents the total profit from the harvested populations during period T. To utilize the resources of the species over a long time, we must develop them rationally and manage them scientifically. We should not only consider the maximization of current economic benefits, but should also consider the ecological balance to ensure the maximization of long-term economic benefits. Our optimal control pair (u∗,v∗) in Ω satisfies J(u∗,v∗)=max(u,v)∈ΩJ(u,v). The admissible control set Ω is as follows:
Ω={(u,v)∈[L∞T(Q)]n×L∞T(R+):0≤ui(x,t)≤Ni,a.e.(x,t)∈Q,0≤v0≤v(t)≤v1,a.e.t∈R+}, |
where
L∞T(Q)={η∈L∞(Q):η(x,t)=η(x,t+T),a.e.(x,t)∈Q}, |
L∞T(R+)={η∈L∞(R+):η(t)=η(t+T),a.e.t∈∈R+}. |
This paper makes the following assumptions:
(A1) Vi:[0,l)×R+→R+ are bounded continuous functions, Vi(x,t)>0 and Vi(x,t)=Vi(x,t+T) for (x,t)∈Q, limx↑lVi(x,t)=0, and Vi(0,t)=1 for t∈R+. There are Lipschitz constants LVi such that
|Vi(x1,t)−Vi(x2,t)|≤LVi|x1−x2|forx1,x2∈[0,l],t∈R+. |
(A2) 0≤βi(x,ci0(t))=βi(x,ci0(t+T))≤¯βi,¯βi are constants.
(A3) {μi(x,ci0(t))=μ0(x)+¯μi(x,ci0(t))a.e.(x,t)∈Q,whereμ0∈L1loc([0,l)),μ0(x)≥0a.e.x∈[0,l),∫l0μ0(x)dx=+∞,¯μi∈L∞(Q),¯μi(x,ci0(t))≥0and¯μi(x,ci0(t))=¯μi(x,ci0(t+T))a.e.(x,t)∈Q.
(A4) fi∈L∞(Q),0≤fi(x,t)=fi(x,t+T),0≤wi(x,t)≤wi(x,t+T)≤¯wi,0≤λi(x,t)≤¯λi,Pi(t) ≤M0,¯wi, ¯λi and M0 are constants.
(A5) ∣βi(x,c1i0(t))−βi(x,c2i0(t))∣≤Lβ∣c1i0(t)−c2i0(t)∣,∣μi(x,c1i0(t))−μi(x,c2i0(t))∣ ≤Lμ∣c1i0(t)−c2i0(t)∣.
(A6) g1≤k1≤g1+m,v1≤h1. (see [24])
This section discusses the existence and uniqueness of the solution to the state system (2.3). For convenience of research, the following definitions are introduced:
Definition 3.1 For i=1,2,…,n, the unique solution x=φi(t;t0,xi0) of the initial value problem x′(t)=Vi(x,t),x(t0)=xi0 is said to be a characteristic curve of (2.3) through (t0,xi0). In particular, let zi(t):=φi(t;0,0) be the characteristic curve through (0, 0) in the x−t plane.
Definition 3.2 Suppose g(x) is an essentially bounded measurable function on E; let
‖g‖∞=infE0⊂E(supE−E0|g(x)|), |
where the infimum is taken for all zero sets E0 in E that make f(x) a bounded function on E−E0. This is also denoted by
Esssupx∈E|g(x)|. |
For any point (x,t)∈[0,l)×[0,T] such that x≤zi(t), define the initial time τ:=τ(t,x); then, φi(t;τ,0)=x⇔φi(τ;t,x)=0. The solution of (2.3) is
pi(x,t)=pi(0,t−z−1i(x))Πi(x;x,t)+∫x0fi(r,φ−1i(r;t,x))Vi(r,φ−1i(r;t,x))Πi(x;x,t)Πi(r;x,t)dr, | (3.1) |
where
Π1(s;x,t)=exp{−∫s0μ1(r,c10(φ−11(r;t,x)))+λ1(r,φ−11(r;t,x))P2(φ−11(r;t,x))V1(r,φ−11(r;t,x))=+u1(r,φ−11(r;t,x))+V1x(r,φ−11(r;t,x))V1(r,φ−11(r;t,x))dr}, |
Πi(s;x,t)=exp{−∫s0μi(r,ci0(φ−1i(r;t,x)))−λ2i−2(r,φ−1i(r;t,x))Pi−1(φ−1i(r;t,x))Vi(r,φ−1i(r;t,x))=+λ2i−1(r,φ−1i(r;t,x))Pi+1(φ−1i(r;t,x))+ui(r,φ−1i(r;t,x))+Vix(r,φ−1i(r;t,x))Vi(r,φ−1i(r;t,x))dr}i=2,3,…,n−1, |
Πn(s;x,t)=exp{−∫s0μn(r,cn0(φ−1n(r;t,x)))−λ2n−2(r,φ−1n(r;t,x))Pn−1(φ−1n(r;t,x))Vn(r,φ−1n(r;t,x))=+un(r,φ−1n(r;t,x))+Vnx(r,φ−1n(r;t,x))Vn(r,φ−1n(r;t,x))dr}. |
ci0(t)=ci0(0)exp{−(g1+m)t}+k1∫t0ce(σ)exp{(σ−t)(g1+m)}dσ. | (3.2) |
ce(t)=ce(0)exp{−∫t0[k2n∑i=1Pi(τ)+h1]dτ}+∫t0[g2n∑i=1ci0(σ)Pi(σ)+v(σ)]⋅exp{∫σt[k2n∑i=1Pi(τ)+h1]dτ}dσ. | (3.3) |
Theorem 3.1 Assume that (A1)−(A6) hold; then, the hybrid system (2.3) has a nonnegative and unique solution (p1(x,t),p2(x,t),…,pn(x,t),c10(t),c20(t),…,cn0(t),ce(t)), such that
(i)(pi(x,t),ci0(t),ce(t))∈L∞(Q)×L∞(0,T)×L∞(0,T).
(ii)0≤ci0(t)≤1,0≤ce(t)≤1,∀t∈(0,T),0≤pi(x,t),∫l0pi(x,t)dx≤M,∀(x,t)∈Q.
Proof. Without loss of generality, we assume that ui(x,t)≡0. p(x,t)=(p1(x,t),p2(x,t),…,pn(x,t)), c0(t)=(c10(t),c20(t),…,cn0(t)), X=[L∞T(R+,L1(0,l))]n×[L∞(R+)]n×L∞(R+). Then, the state space is defined as follows:
Y={(p,c0,ce)∈X∣pi(x,t)≥0a.e.inQ,∫l0pi(x,t)dx≤M,0≤ci0(t)≤1,0≤ce(t)≤1.}, |
where
M:=max{¯β1T‖f1(⋅,⋅)‖L1(Q)exp{¯β1T}+‖f1(⋅,⋅)‖L1(Q),¯βiTexp{¯λ2i−2M0T}‖fi(⋅,⋅)‖L1(Q)=⋅exp{¯βiTexp{¯λ2i−2M0T}}+exp{¯λ2i−2M0T}‖fi(⋅,⋅)‖L1(Q),i=2,3,…,n}. |
Define a mapping
G:Y→X,G(p,c0,ce)=(G1(p,c0,ce),G2(p,c0,ce),…,G2n+1(p,c0,ce)). |
where
Gi(p,c0,ce)(x,t)=pi(0,t−z−1i(x))Πi(x;x,t)+∫x0fi(r,φ−1i(r;t,x))Vi(r,φ−1i(r;t,x))Πi(x;x,t)Πi(r;x,t)dr,i=1,2,…,n. |
Gj(p,c0,ce)(t)=cj0(0)exp{−(g1+m)t}+k1∫t0ce(σ)exp{(σ−t)(g1+m)}dσ,j=n+1,n+2,…,2n. |
G2n+1(p,c0,ce)(t)=ce(0)exp{−∫t0[k2n∑i=1Pi(τ)+h1]dτ}+∫t0[g2n∑i=1ci0(σ)Pi(σ)+v(σ)]⋅exp{∫σt[k2n∑i=1Pi(τ)+h1]dτ}dσ. |
By assumption (A1), we have Vi(0,t)=1. Let bi(t)=pi(0,t), then, noting φ−1i=t−z−1i(x), we have
b1(t)=V1(0,t)p1(0,t)=∫l0β1(x,c10(t))p1(x,t)dx=∫l0β1(x,c10(t))b1(φ−11(0;t,x))dx+∫l0β1(x,c10(t))∫x0f1(r,φ−11(r;t,x))V1(r,φ−11(r;t,x))drdx=∫l0β1(x,c10(t))b1(t−z−11(x))dx+∫l0β1(x,c10(t))∫tφ−11(0;t,x)f1(φ1(σ;t,x),σ)dσdx≤¯β1∫l0b1(t−z−11(x))dx+¯β1∫l0∫t0f1(φ1(σ;t,x),σ)dσdx≤¯β1∫t0b1(σ)dσ+¯β1‖f1(⋅,⋅)‖L1(Q). |
It follows from Bellman's lemma that
b1(t)≤¯β1‖f1(⋅,⋅)‖L1(Q)exp{¯β1T}. |
Then we can see that
∫l0∣G1(p,c0,ce)∣(x,t)dx=∫l0p1(0,φ−11(0;t,x))Π1(x;x,t)dx+∫l0∫x0f1(r,φ−11(r;t,x))V1(r,φ−11(r;t,x))Π1(x;x,t)Π1(r;x,t)drdx≤∫l0p1(0,φ−11(0;t,x))dx+∫l0∫x0f1(r,φ−11(r;t,x))V1(r,φ−11(r;t,x))drdx=∫l0b1(t−z−11)dx+∫l0∫tφ−11(0;t,x)f1(φ1(σ;t,x),σ)dσdx≤∫t0b1(σ)dσ+‖f1(⋅,⋅)‖L1(Q)≤ˉβ1T‖f1(⋅,⋅)‖L1(Q)exp{ˉβ1T}+‖f1(⋅,⋅)‖L1(Q). |
Similarly, we have
∫l0∣Gi(p,c0,ce)∣(x,t)dx≤¯βiTexp{¯λ2i−2M0T}‖fi(⋅,⋅)‖L1(Q)exp{¯βiTexp{¯λ2i−2M0T}}+exp{¯λ2i−2M0T}‖fi(⋅,⋅)‖L1(Q),i=2,3,…,n. |
It follows that G is a mapping from Y to Y. We now discuss the compressibility of the mapping G.
G1(p,c0,ce)(x,t)=b1(φ−11(0;t,x))E1(φ−11(0;t,x);x,t,P2(t))+∫lφ−11(0;t,x)f1(φ1(σ;t,x),σ)E1(σ;x,t,P2(t))dσ, |
where
E1(r;x,t,P2(t))=exp{−∫trμ1(φ1(σ;t,x),c10(σ))+λ1(φ1(σ;t,x),σ)P2(t)+V1x(φ1(σ;t,x),σ)dσ}. |
Then, we have
∫l0∣G1(p1,c10,c1e)−G1(p2,c20,c2e)∣dx=∫l0∣b11(φ−11(0;t,x))E11(φ−11(0;t,x);x,t,P12(t))−b21(φ−11(0;t,x))E21(φ−11(0;t,x);x,t,P22(t))∣dx+∫l0∫tφ−11(0;t,x)∣f1(φ1(σ;t,x),σ)E11(σ;x,t,P12(t))−f1(φ1(σ;t,x),σ)E21(σ;x,t,P22(t))∣dσdx≤∫l0∣b11(φ−11(0;t,x))−b21(φ−11(0;t,x))∣dx+∫l0b11(φ−11(0;t,x))∫tφ−11(0;t,x)∣μ1(φ1(σ;t,x),c110(σ))−μ1(φ1(σ;t,x),c210(σ))∣dσdx+∫l0b11(φ−11(0;t,x))∫tφ−11(0;t,x)¯λ1∣P12(t)−P22(t)∣dσdx+∫l0∫tφ−11(0;t,x)f1(φ1(σ;t,x),σ)∫tσ∣μ1(φ1(σ;t,x),c110(σ))−μ1(φ1(σ;t,x),c210(σ))∣drdσdx+∫l0∫tφ−11(0;t,x)f1(φ1(σ;t,x),σ)∫tσ∣P12(t)−P22(t)∣drdσdx≤∫t0∣b11(σ)−b21(σ)∣dσ+MLμ∫t0∣c110(σ)−c210(σ)∣dσ+¯λ1MT∫l0∣p12(x,σ)−p22(x,σ)∣dx+lLμ‖f1(⋅,⋅)‖L1(Q)∫t0∣c110(σ)−c210(σ)∣dσ+lT‖f1(⋅,⋅)‖L1(Q)∫l0∣p12(x,σ)−p22(x,σ)∣dx≤¯β1∫t0∫l0∣p11(x,σ)−p21(x,σ)∣dxdσ+[MLμ+lLμ‖f1(⋅,⋅)‖L1(Q)+MLβ]∫t0∣c110(σ)−c210(σ)∣dσ+[¯λ1MT+lT‖f1(⋅,⋅)‖L1(Q)]∫l0∣p12(x,σ)−p22(x,σ)∣dx. |
Similarly, we have
∫l0∣Gi(p1,c10,c1e)−Gi(p2,c20,c2e)∣dx≤exp{¯λ2i−2M0T}{¯βi∫t0∫l0∣p1i(x,σ)−p2i(x,σ)∣dxdσ=+[MLμ+lLμ‖f(⋅,⋅)‖L1(Q)+MLβ]∫t0∣c1i0(σ)−c2i0(σ)∣dσ=+[¯λ2i−1MT+lT‖f(⋅,⋅)‖L1(Q)]∫l0∣p1i+1(x,σ)−p2i+1(x,σ)∣dx},i=2,3,…,n−1. |
∫l0∣Gn(p1,c10,c1e)−Gn(p2,c20,c2e)∣dx≤exp{¯λ2n−2M0T}{¯βn∫t0∫l0∣p1n(x,σ)−p2n(x,σ)∣dxdσ=+[MLμ+lLμ‖f(⋅,⋅)‖L1(Q)+MLβ]∫t0∣c1n0(σ)−c2n0(σ)∣dσ}. |
Thus, we have
∫l0∣Gi(p1,c10,c1e)−Gi(p2,c20,c2e)∣dx,≤M1(n∑i=1∫t0∫l0∣p1i(x,σ)−p2i(x,σ)∣dxdσ+∫t0∣c1i0(σ)−c2i0(σ)∣dσ), | (3.4) |
where
M1=max{¯β1,[1+¯λ1MT+lT‖f(⋅,⋅)‖L1(Q)]¯β2exp{¯λ2M0T},…,=[1+¯λ2n−3MT+lT‖f(⋅,⋅)‖L1(Q)]¯βnexp{¯λ2n−2M0T},[MLμ+lLμ‖f(⋅,⋅)‖L1(Q)+MLβ]}. |
By (3.2) and (3.3), we obtain
∣Gj(p1,c10,c1e)−Gj(p2,c20,c2e)∣(t)(j=n+1,n+2,…,2n)≤M2∫t0∣c1e(σ)−c2e(σ)∣dσ, | (3.5) |
where M2=k1.
∣G2n+1(p1,c10,c1e)−G2n+1(p2,c20,c2e)∣(t)≤M3(n∑i=1∫t0∫l0∣p1i(x,σ)−p2i(x,σ)∣dxdσ+∫t0∣c1i0(σ)−c2i0(σ)∣dσ), | (3.6) |
where M3=max{k2+g2+k2h1T+k2g2MT,g2M}.
We now show that the mapping G has a unique fixed point. Due to the periodicity of elements in the set Y, we consider the case [0,T] only. We define an equivalent norm in L∞(0,T) as follows:
‖(p,c0,ce)‖∗=Esssupt∈[0,T]e−λt{n∑i=1∫l0∣pi(x,t)∣dx+n∑i=1∣ci0(t)∣+∣ce(t)∣}, |
where λ>0 is large enough. Then, we have
‖G(p1,c10,c1e)−G(p2,c20,c2e)‖∗≤M4Esssupt∈[0,T]e−λt∫t0{n∑i=1∫l0∣p1i(x,σ)−p2i(x,σ)∣dx=+n∑i=1∣c1i0(σ)−c2i0(σ)∣+∣c1e(σ)−c2e(σ)∣}dσ≤M4Esssupt∈[0,T]e−λt∫t0eλσ{e−λσ[n∑i=1∫l0∣p1i(x,σ)−p2i(x,σ)∣dx=+n∑i=1∣c1i0(σ)−c2i0(σ)∣+∣c1e(σ)−c2e(σ)∣]}dσ≤M4‖(p1−p2,c10−c20,c1e−c2e)‖∗Esssupt∈[0,T]e−λt∫t0eλσdσ≤M4λ‖(p1−p2,c10−c20,c1e−c2e)‖∗, |
where M4=max{M1,M2,M3}. Thus, choosing λ>M4 yields G as a strict contraction on the space of (Y,‖⋅‖∗). By the Banach fixed point theory, G has a unique fixed point, which is the solution of the system (2.3).
Theorem 3.2 If T is small enough, then there are constants Kj(T) with limT→0Kj(T)>0,j=1,2, such that
n∑i=1‖p1i−p2i‖L∞(0,T;L1(0,l))+n∑i=1‖c1i0−c2i0‖L∞(0,T)+‖c1e−c2e‖L∞(0,T)≤K1(T)T(n∑i=1‖u1i−u2i‖L∞(0,T;L1(0,l))+‖v1−v2‖L∞(0,T)). | (3.7) |
n∑i=1‖p1i−p2i‖L1(Q)+n∑i=1‖c1i0−c2i0‖L1(0,T)+‖c1e−c2e‖L1(0,T)≤K2(T)T(n∑i=1‖u1i−u2i‖L1(Q)+‖v1−v2‖L1(0,T)). | (3.8) |
Proof. Let (pj,cj0,cje) be the solution of (2.3) corresponding to (uj,vj),j=1,2, it follows from (3.1) that
∫l0∣p11(x,t)−p21(x,t)∣dx=∫l0∣b11(φ−11(0;t,x))Π11(x;x,t)−b21(φ−11(0;t,x))Π21(x;x,t)∣dx+∫l0∫x0|f1(r,φ−11(r;t,x))V1(r,φ−11(r;t,x))Π11(x;x,t)Π11(r;x,t)−f1(r,φ−11(r;t,x))V1(r,φ−11(r;t,x))Π21(x;x,t)Π21(r;x,t)|drdx≤∫l0∣b11(φ−11(0;t,x))−b11(φ−11(0;t,x))∣dx+∫l0b11(φ−11(0;t,x))[∫tφ−11(0;t,x)∣μ1(φ1(r;t,x),c110(r))−μ1(φ1(r;t,x),c210(r))∣dr]dx+∫l0b11(φ−11(0;t,x))[∫tφ−11(0;t,x)∣λ1(φ1(r;t,x),r)(P12(t)−P22(t))∣dr]dx+∫l0b11(φ−11(0;t,x))[∫tφ−11(0;t,x)∣u11(φ1(r;t,x),r)−u21(φ1(r;t,x),r)∣dr]dx+∫l0∫tφ−11(0;t,x)f1(φ1(σ;t,x),σ)[∫tσ∣μ1(φ1(σ;t,x),c110(σ))−μ1(φ1(σ;t,x),c210(σ))∣dr]dσdx+∫l0∫tφ−11(0;t,x)f1(φ1(σ;t,x),σ)[∫tσ∣λ1(φ1(σ;t,x),σ)(P12(t)−P22(t))∣dr]dσdx+∫l0∫tφ−11(0;t,x)f1(φ1(σ;t,x),σ)[∫tσ∣u11(φ1(σ;t,x),σ)−u21(φ1(σ;t,x),σ)∣dr]dσdx≤∫t0∣b11(σ)−b21(σ)∣dx+¯β1MLμ∫t0∣c110(σ)−c210(σ)∣dσ+¯β1¯λ1MT∫l0∣p12(x,σ)−p12(x,σ)∣dx+¯β1M∫t0∫l0∣u11(φ1(σ;t,x),σ)−u21(φ1(σ;t,x),σ)∣dxdσ+lLμ‖f1(⋅,⋅)‖L1(Q)∫t0∣c110(σ)−c210(σ)∣dσ+¯λ1lT‖f1(⋅,⋅)‖L1(Q)∫l0∣p12(x,σ)−p22(x,σ)∣dx+‖f1(⋅,⋅)‖L1(Q)∫t0∫l0∣u11(φ1(σ;t,x),σ)−u21(φ1(σ;t,x),σ)∣dxdσ≤¯β1∫t0∫l0∣p11(x,σ)−p21(x,σ)∣dxdσ+¯λ1[MT+lT‖f1(⋅,⋅)‖L1(Q)]∫l0∣p12(x,σ)−p22(x,σ)∣dx+[¯β1MLμ+lLμ‖f1(⋅,⋅)‖L1(Q)+MLβ]∫t0∣c110(σ)−c210(σ)∣dσ+[¯β1M+‖f1(⋅,⋅)‖L1(Q)]∫t0∫l0∣u11(φ1(σ;t,x),σ)−u21(φ1(σ;t,x),σ)∣dxdσ. | (3.9) |
Similarly, we can see that
\begin{array}{l} \int_{0}^{l}\mid p_{i}^{1}(x,t)-p_{i}^{2}(x,t)\mid{\mathrm{d}}x\\ \leq\overline{\beta}_{i}\exp\left\{\overline{\lambda}_{2i-2}M_{0}T\right\}\int_{0}^{t}\int_{0}^{l}\mid p_{i}^{1}(x,\sigma)-p_{i}^{2}(x,\sigma)\mid{\mathrm{d}}x{\mathrm{d}}\sigma\\ +\exp\left\{\overline{\lambda}_{2i-2}M_{0}T\right\}\left[\overline{\beta}_{i}ML_{\mu}+lL_{\mu}\|f_{i}(\cdot,\cdot)\|_{L^{1}(Q)} +ML_{\beta}\right]\\\int_{0}^{t}\mid c_{i0}^{1}(\sigma)-c_{i0}^{2}(\sigma)\mid{\mathrm{d}}\sigma\\ +\overline{\lambda}_{i}\exp\left\{\overline{\lambda}_{2i-2}M_{0}T\right\}[MT+lT\|f_{i}(\cdot,\cdot)\|_{L^{1}(Q)}]\int_{0}^{l}\mid p_{i+1}^{1}(x,\sigma)-p_{i+1}^{2}(x,\sigma)\mid{\mathrm{d}}x\\ +\exp\left\{\overline{\lambda}_{2i-2}M_{0}T\right\}\left[\overline{\beta}_{i}M+\|f_{i}(\cdot,\cdot)\|_{L^{1}(Q)}\right]\int_{0}^{t}\int_{0}^{l}\mid u_{i}^{1}(\varphi_{i}(\sigma;t,x),\sigma)\\-u_{i}^{2}(\varphi_{i}(\sigma;t,x),\sigma)\mid{\mathrm{d}}x{\mathrm{d}}\sigma,\\ \; i = 2,3,\ldots,n-1. \end{array} | (3.10) |
\begin{array}{l} \int_{0}^{l}\mid p_{n}^{1}(x,t)-p_{n}^{2}(x,t)\mid{\mathrm{d}}x\\ \leq\; \overline{\beta}_{n}\exp\left\{\overline{\lambda}_{2n-2}M_{0}T\right\}\int_{0}^{t}\int_{0}^{l}\mid p_{n}^{1}(x,\sigma)-p_{n}^{2}(x,\sigma)\mid{\mathrm{d}}x{\mathrm{d}}\sigma\\ +\exp\left\{\overline{\lambda}_{2n-2}M_{0}T\right\}\left[\overline{\beta}_{n}ML_{\mu}+lL_{\mu}\|f_{n}(\cdot,\cdot)\|_{L^{1}(Q)} +ML_{\beta}\right]\int_{0}^{t}\mid c_{n0}^{1}(\sigma)-c_{n0}^{2}(\sigma)\mid{\mathrm{d}}\sigma\\ +\exp\left\{\overline{\lambda}_{2n-2}M_{0}T\right\}\left[\overline{\beta}_{n}M+\|f_{n}(\cdot,\cdot)\|_{L^{1}(Q)}\right]\\\int_{0}^{t}\int_{0}^{l}\mid u_{n}^{1}(\varphi_{n}(\sigma;t,x),\sigma)-u_{n}^{2}(\varphi_{n}(\sigma;t,x),\sigma)\mid{\mathrm{d}}x{\mathrm{d}}\sigma. \end{array} | (3.11) |
From (3.2) and (3.3), we obtain
\begin{align} \mid c_{i0}^{1}(t)-c_{i0}^{2}(t)\mid\leq k_{1}\int_{0}^{t}\mid c_{e}^{1}(s)-c_{e}^{2}(s)\mid{\mathrm{d}}s. \end{align} | (3.12) |
\begin{align} \mid c_{e}^{1}(t)-c_{e}^{2}(t)\mid\leq& (k_{2}+g_{2}+k_{2}h_{1}T+k_{2}g_{2}MT)\sum\limits_{i = 1}^{n}\int_{0}^{l}\int_{0}^{t}\mid p_{i}^{1}(x,s)-p_{i}^{2}(x,s)\mid{\mathrm{d}}s{\mathrm{d}}x\\ \phantom{ = \;\;}&+g_{2}M\sum\limits_{i = 1}^{n}\int_{0}^{t}\mid c_{i0}^{1}(s)-c_{i0}^{2}(s)\mid{\mathrm{d}}s +\int_{0}^{t}\mid v^{1}(s)-v^{2}(s)\mid{\mathrm{d}}s. \end{align} | (3.13) |
Finally, from (3.9)–(3.13), we can deduce (3.7), where K_{1}(T) is a constant depending on the bounds of parameters in (2.3). In addition, combined with (3.9)–(3.13), there exists K_{2}(T) , which enables us to obtain (3.8) of the standard norm in L^{1} space.
In this section, we employ tangent-normal cone techniques in nonlinear functional analysis to deduce the necessary conditions for the optimal control pair.
Theorem 4.1 If (u^{\ast}, v^{\ast}) is an optimal control pair and (p^{\ast}, c_{0}^{\ast}, c_{e}^{\ast}) is the corresponding optimal state, then
\begin{align} u_{i}^{\ast}(x,t) = \mathcal{F}_{i}\left(\frac{[w_{i}(x,t)-\xi_{i}(x,t)]p_{i}^{\ast}(x,t)}{c_{i}}\right), \end{align} | (4.1) |
\begin{align} v^{\ast}(t) = \mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}(t)}{c_{n+1}}\right), \end{align} | (4.2) |
in which \mathcal{F}_{j} are given by
\begin{align} \mathcal{F}_{j}(\eta) = \left\{ \begin{array}{ll} 0,\; \; \; \; \; \; \; \; \eta < 0,\\ \eta,\; \; \; \; \; \; \; \; 0\leq \eta\leq N_{j},\; \; j = 1,2,\ldots,n+1,\\ N_{j},\; \; \; \; \; \; \eta > N_{j},\\ \end{array} \right. \end{align} | (4.3) |
and (\xi_{1}, \xi_{2}, \ldots, \xi_{2n+1}) is the solution of the following adjoint system corresponding to (u^{\ast}, v^{\ast}) :
\begin{align} \left\{ \begin{array}{ll} \frac{\partial \xi_{1}}{\partial t}+V_{1}\frac{\partial \xi_{1}}{\partial x} = [\mu_{1}(x,c_{10}^{\ast}(t))+\lambda_{1}P_{2}^{\ast}(t)+u_{1}^{\ast}]\xi_{1} \\+[k_{2}c_{e}^{\ast}(t)-g_{2}c_{10}^{\ast}(t)]\xi_{2n+1}-\xi_{1}(0,t)\beta_{1}(x,c_{10}^{\ast}(t))\\-\int_{0}^{l}(\lambda_{2}p_{2}^{\ast}\xi_{2})(x,t){\mathrm{d}}x+w_{1}u_{1}^{\ast},\\ \frac{\partial \xi_{i}}{\partial t}+V_{i}\frac{\partial \xi_{i}}{\partial x} = [\mu_{i}(x,c_{i0}^{\ast}(t))-\lambda_{2i-2}P_{i-1}^{\ast}(t)\\+\lambda_{2i-1}P_{i+1}^{\ast}(t) +u_{i}^{\ast}]\xi_{i}+[k_{2}c_{e}^{\ast}(t)-g_{2}c_{i0}^{\ast}(t)]\xi_{2n+1}\\ -\xi_{i}(0,t)\beta_{i}(x,c_{i0}^{\ast}(t)) +\int_{0}^{l}(\lambda_{2i-3}p_{i-1}^{\ast}\xi_{i-1}-\lambda_{2i}p_{i+1}^{\ast}\xi_{i+1})(x,t){\mathrm{d}}x+w_{i}u_{i}^{\ast},\\ i = 2,3,\ldots n-1,\\ \frac{\partial \xi_{n}}{\partial t}+V_{n}\frac{\partial \xi_{n}}{\partial x} = [\mu_{n}(x,c_{n0}(t))-\lambda_{2n-2}P_{n-1}^{\ast}(t)+u_{n}^{\ast}]\xi_{n} +[k_{2}c_{e}^{\ast}(t)-g_{2}c_{n0}^{\ast}(t)]\xi_{2n+1}\\ -\xi_{n}(0,t)\beta_{n}(x,c_{n0}^{\ast}(t)) +\int_{0}^{l}(\lambda_{2n-3}p_{n-1}^{\ast}\xi_{n-1})(x,t){\mathrm{d}}x+w_{n}u_{n}^{\ast},\\ \frac{{\mathrm{d}}\xi_{n+i}}{{\mathrm{d}}t} = \int_{0}^{l}\frac{\partial\mu_{i}(x,c_{i0}^{\ast}(t))}{\partial c_{i0}}p_{i}^{\ast}\xi_{i}{\mathrm{d}}t+(g_{1}+m)\xi_{n+i}\\-g_{2}P_{i}^{\ast}(t)\xi_{2n+1} -\xi_{i}(0,t)\int_{0}^{l}\frac{\partial\beta_{i}(x,c_{i0}^{\ast}(t))}{\partial c_{i0}}p_{i}^{\ast}{\mathrm{d}}x,\\ \frac{{\mathrm{d}}\xi_{2n+1}}{{\mathrm{d}}t} = -k_{1}\sum\limits_{i = 1}^{n}\xi_{i+n}+\left[k_{2}\sum\limits_{i = 1}^{n}P_{i}^{\ast}(t)+h_{1}\right]\xi_{2n+1},\\ \xi_{i}(l,t) = 0,\xi_{i}(x,t) = \xi_{i}(x,t+T),\\ \xi_{j}(T) = 0,\; \; \; \; j = n+1,n+2,\ldots,2n+1.\\ \end{array} \right. \end{align} | (4.4) |
Proof. The existence of a unique, bounded solution to the system (4.4) can be treated in the same manner as that for (2.3). For any given (\nu_{1}, \nu_{2})\in\mathcal{T}_{\Omega}(u^{\ast}, v^{\ast}) (the tangent cone of \Omega at (u^{\ast}, v^{\ast}) ), \nu_{1} = (\nu_{11}, \nu_{21}, \ldots, \nu_{n1}) , such that (u^{\ast}+\varepsilon\nu_{1}, v^{\ast}+\varepsilon\nu_{2})\in\Omega , where \varepsilon > 0 is small enough. Then, from J(u^{\ast}+\varepsilon\nu_{1}, v^{\ast}+\varepsilon\nu_{2})\leq J(u^{\ast}, v^{\ast}) , we derive
\begin{align*} &\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}w_{i}(u_{i}^{\ast}+\varepsilon \nu_{i1})p_{i}^{\varepsilon}{\mathrm{d}}x{\mathrm{d}}t -\frac{1}{2}\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}c_{i}(u_{i}^{\ast}+\varepsilon\nu_{i1})^{2}{\mathrm{d}}x{\mathrm{d}}t -\frac{1}{2}\int_{0}^{T}c_{n+1}(v^{\ast}+\varepsilon\nu_{2})^{2}{\mathrm{d}}t\\ \leq&\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}w_{i}u_{i}^{\ast}p_{i}^{\ast}{\mathrm{d}}x{\mathrm{d}}t -\frac{1}{2}\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}c_{i}u_{i}^{\ast^{2}}{\mathrm{d}}x{\mathrm{d}}t -\frac{1}{2}\int_{0}^{T}c_{n+1}v^{\ast^{2}}{\mathrm{d}}t, \end{align*} |
that is
\begin{align*} &\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}w_{i}u_{i}^{\ast}\frac{p_{i}^{\varepsilon}-p_{i}^{\ast}}{\varepsilon}{\mathrm{d}}x{\mathrm{d}}t +\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}w_{i}\nu_{i1}p_{i}^{\varepsilon}{\mathrm{d}}x{\mathrm{d}}t\\ \leq&\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}c_{i}u_{i}^{\ast}\nu_{i1}{\mathrm{d}}x{\mathrm{d}}t +\int_{0}^{T}c_{n+1}v^{\ast}\nu_{2}{\mathrm{d}}t+\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}c_{i}\varepsilon\nu_{i1}^{2}{\mathrm{d}}x{\mathrm{d}}t +\int_{0}^{T}c_{n+1}\varepsilon\nu_{2}^{2}{\mathrm{d}}t, \end{align*} |
as \varepsilon\rightarrow0^{+} ,
\begin{align} \sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}w_{i}(u_{i}^{\ast}z_{i}+\nu_{i1}p_{i}^{\ast}) {\mathrm{d}}x{\mathrm{d}}t-\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}c_{i}u_{i}^{\ast}\nu_{i1}{\mathrm{d}}x{\mathrm{d}}t -\int_{0}^{T}c_{n+1}v^{\ast}\nu_{2}{\mathrm{d}}t\leq0, \end{align} | (4.5) |
where
\begin{align*} \frac{1}{\varepsilon}(p_{i}^{\varepsilon}-p_{i}^{\ast})\rightarrow z_{i}, \frac{1}{\varepsilon}(c_{i0}^{\varepsilon}-c_{i0}^{\ast})\rightarrow z_{i+n}, \frac{1}{\varepsilon}(c_{e}^{\varepsilon}-c_{e}^{\ast})\rightarrow z_{2n+1},\; \; \; \; \mbox{as}\; \; \varepsilon\rightarrow0. \end{align*} |
By Theorem 3.2 we see that z_{1}, z_{2}, \ldots, z_{2n+1} makes sense ([25], p.18). (p^{\varepsilon}, c_{0}^{\varepsilon}, c_{e}^{\varepsilon}) is the state corresponding to (u^{\ast}+\varepsilon\nu_{1}, v^{\ast}+\varepsilon\nu_{2}) . It follows from (2.3) that ( z_{1} , z_{2}, \ldots , z_{2n+1} ) satisfies
\begin{align} \left\{ \begin{array}{ll} \frac{\partial z_{1}}{\partial t}+V_{1}\frac{\partial z_{1}}{\partial x} = -[\mu_{1}(x,c_{10}^{\ast}(t))+\lambda_{1}P_{2}^{\ast}(t)+V_{1x}+u_{1}^{\ast}]z_{1}\\-\lambda_{1}Z_{2}(t)p_{1}^{\ast} -\frac{\partial\mu_{1}(x,c_{10}^{\ast}(t))}{\partial c_{10}}p_{1}^{\ast}z_{n+1}-v_{11}p_{1}^{\ast},\\ \frac{\partial z_{i}}{\partial t}+V_{i}\frac{\partial z_{i}}{\partial x} = -[\mu_{i}(x,c_{i0}(t))\\-\lambda_{2i-2}P_{i-1}^{\ast}(t)+\lambda_{2i-1}P_{i+1}^{\ast}(t) +V_{ix}+u_{i}^{\ast}]z_{i}+\lambda_{2i-2}Z_{i-1}(t)p_{i}^{\ast}\\ -\lambda_{2i-1}(x,t)Z_{i+1}(t)p_{i}^{\ast}\\-\frac{\partial\mu_{i}(x,c_{i0}^{\ast}(t))}{\partial c_{i0}}p_{i}^{\ast}z_{n+i}-v_{i1}p_{i}^{\ast},\; \; \; \; i = 2,3,\ldots n-1,\\ \frac{\partial z_{n}}{\partial t}+V_{n}\frac{\partial z_{n}}{\partial x} = -[\mu_{n}(x,c_{n0}(t))-\lambda_{2n-2}P_{n-1}^{\ast}(t)+V_{nx}+u_{n}^{\ast}]z_{n}+\lambda_{2n-2}Z_{n-1}(t)p_{n}^{\ast}\\ -\frac{\partial\mu_{n}(x,c_{n0}^{\ast}(t))}{\partial c_{n0}}p_{n}^{\ast}z_{2n}-v_{n}p_{n}^{\ast},\\ \frac{{\mathrm{d}}z_{n+i}}{{\mathrm{d}}t} = k_{1}z_{2n+1}-g_{1}z_{n+i}-mz_{n+i},\\ \frac{{\mathrm{d}}z_{2n+1}}{{\mathrm{d}}t} = -k_{2}c_{e}^{\ast}(t)\sum\limits_{i = 1}^{n}Z_{i}(t) +g_{2}\sum\limits_{i = 1}^{n}(c_{i0}^{\ast}Z_{i}(t)+z_{n+i}P_{i}^{\ast}(t)) \\-\left[k_{2}\sum\limits_{i = 1}^{n}P_{i}^{\ast}(t)+h_{1}\right]z_{2n+1}+v_{2},\\ V_{i}(0,t)z_{i}(0,t) = \int_{0}^{l}\beta_{i}(x,c_{i0}^{\ast}(t))z_{i}{\mathrm{d}}x +\int_{0}^{l}\frac{\partial\beta_{i}(a,c_{i0}^{\ast}(t))}{\partial c_{i0}}p_{i}^{\ast}z_{n+i}{\mathrm{d}}x,\\ z_{i}(x,t) = z_{i}(x,t+T),\; z_{n+i}(0) = z_{2n+1}(0) = 0,\\ P_{i}^{\ast}(t) = \int_{0}^{l}p_{i}^{\ast}(x,t){\mathrm{d}}x,\; Z_{i}(t) = \int_{0}^{l}z_{i}(x,t){\mathrm{d}}x.\\ \end{array} \right. \end{align} | (4.6) |
We multiply the first five equations in (4.6) by \xi_{i}(i = 1, 2, \ldots, 2n+1) respectively, and integrate on Q and (0, T) . By using (4.4), we have
\begin{align} \sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}w_{i}u_{i}^{\ast}z_{i}{\mathrm{d}}x{\mathrm{d}}t = -\sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}\nu_{i1}\xi_{i}p_{i}^{\ast}{\mathrm{d}}x{\mathrm{d}}t+\int_{0}^{T}\nu_{2}\xi_{2n+1}{\mathrm{d}}t. \end{align} | (4.7) |
Substituting (4.7) and (4.5) gives
\begin{align*} \sum\limits_{i = 1}^{n}\int_{0}^{T}\int_{0}^{l}[(w_{i}-\xi_{i})p_{i}^{\ast}-c_{i}u_{i}^{\ast}] \nu_{i1}{\mathrm{d}}x{\mathrm{d}}t+\int_{0}^{T}(-c_{n+1}v^{\ast}+\xi_{2n+1})\nu_{2}{\mathrm{d}}t\leq0, \end{align*} |
for any (\nu_{1}, \nu_{2})\in\mathcal{T}_{\Omega}(u^{\ast}, v^{\ast}) . Consequently, the structure of normal cone tells us that [(w_{i}-\xi_{i})p_{i}^{\ast}-c_{i}u_{i}^{\ast}, -c_{n+1}v^{\ast}+\xi_{2n+1}]\in \mathcal{N}_{\Omega}(u^{\ast}, v^{\ast}) , (the normal cone of \Omega at (u^{\ast}, v^{\ast}) ), which gives the desired result.
Theorem 4.2 If T is small enough, then there is a constant K_{3} such that
\begin{align} &\sum\limits_{i = 1}^{n}\|\xi_{i}^{1}-\xi_{i}^{2}\|_{L^{\infty}(Q)}+\sum\limits_{q = n+1}^{2n}\|\xi_{q}^{1}-\xi_{q}^{2}\|_{L^{\infty}(0,T)} +\|\xi_{2n+1}^{1}-\xi_{2n+1}^{2}\|_{L^{\infty}(0,T)}\\ \leq& K_{3}(T)T\left(\sum\limits_{i = 1}^{n}\|u_{i}^{1}-u_{i}^{2}\|_{L^{\infty}(Q)}+\|v^{1}-v^{2}\|_{L^{\infty}(0,T)}\right). \end{align} | (4.8) |
The proof process of Theorem 4.2 is similar to that of Theorem 3.2, and is omitted here.
In order to show that there exists a unique optimal control pair by means of the Ekeland variational principle, we embed the functional \widetilde{J}(u, v) into [L^{1}(Q)]^{n}\times L^{1}(0, T) . We define
\begin{align*} \widetilde{J}(u,v) = \left\{ \begin{array}{ll} J(u,v),\; \; \; \; (u,v)\in\Omega,\\ -\infty,\; \; \; \; \; \; \; \; \mbox{otherwise}.\\ \end{array} \right. \end{align*} |
Lemma 5.1([6], Lemma 4.1) \; \widetilde{J}(u, v) is upper semi-continuous in [L^{1}(Q)]^{n}\times L^{1}(0, T) .
Theorem 5.1 If T is sufficiently small, there exists one and only one optimal control pair (u^{\ast}, v^{\ast}) , which is in the feedback, and is determined by (4.1)–(4.4) and (2.3), where C_{1} and C_{2} are the supremum of |p_{i}| and |\xi_{j}|, i = 1, 2, \ldots, n, j = 1, 2, \ldots, 2n+1 , respectively.
Proof. Define the mapping \mathcal{L}:\Omega\rightarrow \Omega as follows:
\begin{align*} \mathcal{L}(u,v)& = \mathcal{F}\left(\frac{(w_{1}-\xi_{1})p_{1}}{c_{1}}, \frac{(w_{2}-\xi_{2})p_{2}}{c_{2}},\ldots, \frac{(w_{n}-\xi_{n})p_{n}}{c_{n}},\frac{\xi_{2n+1}}{c_{n+1}}\right)\\ & = \left(\mathcal{F}_{1}\left(\frac{(w_{1}-\xi_{1})p_{1}}{c_{1}}\right), \mathcal{F}_{2}\left(\frac{(w_{2}-\xi_{1})p_{2}}{c_{2}}\right),\ldots,\right. \left.\mathcal{F}_{n}\left(\frac{(w_{n}-\xi_{n})p_{n}}{c_{n}}\right), \mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}}{c_{n+1}}\right)\right), \end{align*} |
where (p, c_{0}, c_{e}) and (\xi_{1}, \xi_{2}, \ldots, \xi_{2n+1}) are the state and adjoint state, respectively, corresponding to the control (u, v) . We show that \mathcal{L} admits a unique fixed point, which maximizes the functional \widetilde{J} .
From Lemma 5.1 and the Ekeland variational principle ([26], p.180), for any given \varepsilon > 0 , there exists (u^{\varepsilon}, v^{\varepsilon})\in\Omega\subset[L^{1}(Q)]^{n}\times L^{1}(0, T) such that
\begin{align} \widetilde{J}(u^{\varepsilon},v^{\varepsilon})\geq\sup\limits_{(u,v)\in\Omega}\widetilde{J}(u,v)-\varepsilon, \end{align} | (5.1) |
\begin{align} \widetilde{J}(u^{\varepsilon},v^{\varepsilon})\geq\sup\limits_{(u,v)\in\Omega}\left\{\widetilde{J}(u,v)- \sqrt{\varepsilon}\left(\sum\limits_{i = 1}^{n}\|u_{i}^{\varepsilon}-u_{i}\|_{L^{1}(Q)}+ \|v^{\varepsilon}-v\|_{L^{1}(0,T)}\right)\right\}, \end{align} | (5.2) |
Thus, the perturbed functional
\begin{align*} \widetilde{J}_{\varepsilon}(u,v) = \widetilde{J}(u,v)-\sqrt{\varepsilon}\left(\sum\limits_{i = 1}^{n}\|u_{i}^{\varepsilon}-u_{i}\|_{L^{1}(Q)}+\|v^{\varepsilon}-v\|_{L^{1}(0,T)}\right), \end{align*} |
attains its supremum at (u^{\varepsilon}, v^{\varepsilon}) . Then, we argue as in Theorem 4.1:
\begin{align} (u^{\varepsilon},v^{\varepsilon}) = \mathcal{L}(u^{\varepsilon},v^{\varepsilon}) = &\left(\mathcal{F}_{1}\left(\frac{(w_{1}-\xi_{1}^{\varepsilon}) p_{1}^{\varepsilon}+\sqrt{\varepsilon}\theta_{1}^{\varepsilon}}{c_{1}}\right),\mathcal{F}_{2}\left(\frac{(w_{2}-\xi_{2}^{\varepsilon}) p_{2}^{\varepsilon}+\sqrt{\varepsilon}\theta_{2}^{\varepsilon}}{c_{2}}\right),\ldots,\right.\\ \phantom{ = \;\;}&\left.\mathcal{F}_{n}\left(\frac{(w_{n}-\xi_{n}^{\varepsilon})p_{n}^{\varepsilon}+\sqrt{\varepsilon}\theta_{n}^{\varepsilon}}{c_{n}}\right), \mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}^{\varepsilon}+\sqrt{\varepsilon}\theta_{n+1}^{\varepsilon}}{c_{n+1}}\right)\right), \end{align} | (5.3) |
where (p^{\varepsilon}, c_{0}^{\varepsilon}, c_{e}^{\varepsilon}) and (\xi_{1}^{\varepsilon}, \xi_{2}^{\varepsilon}, \ldots, \xi_{2n+1}^{\varepsilon}) are the state and adjoint state, respcetively, corresponding to the control (u^{\varepsilon}, v^{\varepsilon}) , \theta_{1}^{\varepsilon}, \theta_{2}^{\varepsilon}, \ldots, \theta_{n}^{\varepsilon}\in L^{\infty}(Q) , and \theta_{n+1}^{\varepsilon}\in L^{\infty}(0, T) with \mid\theta_{i}^{\varepsilon}\mid\leq1, i = 1, 2, \ldots, n+1 .
First, we show that \mathcal{L} has only one fixed point. Let (p^{j}, c_{0}^{j}, c_{e}^{j}) and (\xi_{1}^{j}, \xi_{2}^{j}, \ldots, \xi_{2n+1}^{j}) be the state and adjoint state corresponding to the control (u^{j}, v^{j}), j = 1, 2 . By (3.7) and (4.8), we have
\begin{align*} &\|\mathcal{L}(u^{1},v^{1})-\mathcal{L}(u^{2},v^{2})\|_{\infty}\\ = &\sum\limits_{i = 1}^{n}\left\|\mathcal{F}_{i}\left(\frac{(w_{i}-\xi_{i}^{1})p_{i}^{1}}{c_{i}}\right) -\mathcal{F}_{i}\left(\frac{(w_{i}-\xi_{i}^{2})p_{i}^{2}}{c_{i}}\right)\right\|_{L^{\infty}(Q)} +\left\|\mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}^{1}}{c_{n+1}}\right)-\mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}^{1}}{c_{n+1}}\right)\right\|_{L^{\infty}(0,T)}\\ \leq&\sum\limits_{i = 1}^{n}\left\|\frac{(w_{i}-\xi_{i}^{1})p_{i}^{1}}{c_{i}}- \frac{(w_{i}-\xi_{i}^{1})p_{i}^{1}}{c_{i}}\right\|_{L^{\infty}(Q)}+\left\|\frac{\xi_{2n+1}^{1}}{c_{n+1}}-\frac{\xi_{2n+1}^{2}}{c_{n+1}}\right\|_{L^{\infty}(0,T)}\\ \leq&\sum\limits_{i = 1}^{n}\left\|\frac{w_{i}(p_{i}^{1}-p_{i}^{2})}{c_{i}}+\frac{\mid\xi_{i}^{1}\mid(p_{i}^{1}-p_{i}^{2})}{c_{i}} +\frac{\mid p_{i}^{2}\mid(\xi_{i}^{1}-\xi_{i}^{2})}{c_{i}}\right\|_{L^{\infty}(Q)} +\left\|\frac{\xi_{2n+1}^{1}-\xi_{2n+1}^{2}}{c_{n+1}}\right\|_{L^{\infty}(0,T)}\\ \leq& T\left(\sum\limits_{i = 1}^{n}\frac{1}{c_{i}}(\overline{w}_{i}K_{1}+C_{2}K_{1}+C_{1}K_{3})+K_{3}\right) \cdot\left(\sum\limits_{i = 1}^{n}\|u_{i}^{1}-u_{i}^{2}\|_{L^{\infty}(Q)}+\|v^{1}-v^{2}\|_{L^{\infty}(0,T)}\right). \end{align*} |
Clearly, \mathcal{L} is a contraction if T is sufficiently small. Hence, \mathcal{L} has a unique fixed point (u^{\ast}, v^{\ast}) .
Next, we prove (u^{\varepsilon}, v^{\varepsilon})\rightarrow (u^{\ast}, v^{\ast}) as \varepsilon\rightarrow0^{+} . The relations (4.1), (4.2), and (5.3) lead to
\begin{align*} &\|\mathcal{L}(u^{\varepsilon},v^{\varepsilon})-(u^{\varepsilon},v^{\varepsilon}\|_{\infty}\\ = &\left\|\left(\mathcal{F}_{1}\left(\frac{(w_{1}-\xi_{1}^{\varepsilon})p_{1}^{\varepsilon}}{c_{1}}\right), \mathcal{F}_{2}\left(\frac{(w_{2}-\xi_{2}^{\varepsilon})p_{2}^{\varepsilon}}{c_{2}}\right),\ldots, \mathcal{F}_{n}\left(\frac{(w_{n}-\xi_{n}^{\varepsilon})p_{n}^{\varepsilon}}{c_{n}}\right)\right., \mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}^{\varepsilon}}{c_{n+1}}\right)\right)\\ &-\left(\mathcal{F}_{1}\left(\frac{(w_{1}-\xi_{1}^{\varepsilon})p_{1}^{\varepsilon}}{c_{1}} +\frac{\sqrt{\varepsilon}\theta_{1}^{\varepsilon}}{c_{1}}\right)\right., \mathcal{F}_{2}\left(\frac{(w_{2}-\xi_{2}^{\varepsilon})p_{2}^{\varepsilon}}{c_{2}}+\frac{\sqrt{\varepsilon}\theta_{2}^{\varepsilon}}{c_{2}}\right),\ldots,\\ &\; \mathcal{F}_{n}\left(\frac{(w_{n}-\xi_{n}^{\varepsilon})p_{n}^{\varepsilon}}{c_{n}}+\frac{\sqrt{\varepsilon}\theta_{n}^{\varepsilon}}{c_{n}}\right), \left.\left.\mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}^{\varepsilon}}{c_{n+1}}+\frac{\sqrt{\varepsilon}\theta_{n+1}^{\varepsilon}}{c_{n+1}}\right)\right)\right\|_{\infty}\\ \leq&\sum\limits_{i = 1}^{n}\left\|\left(\mathcal{F}_{i}\left(\frac{(w_{i}-\xi_{i}^{\varepsilon}) p_{i}^{\varepsilon}}{c_{i}}\right)-\mathcal{F}_{i}\left(\frac{(w_{i}-\xi_{i}^{\varepsilon})p_{i}^{\varepsilon}}{c_{i}}+ \frac{\sqrt{\varepsilon}\theta_{i}^{\varepsilon}}{c_{i}}\right)\right)\right\|_{L^{\infty}(Q)}\\ &+\left\| \mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}^{\varepsilon}}{c_{n+1}}\right)-\mathcal{F}_{n+1}\left(\frac{\xi_{2n+1}^{\varepsilon}}{c_{n+1}} +\frac{\sqrt{\varepsilon}\theta_{n+1}^{\varepsilon}}{c_{n+1}}\right)\right\|_{L^{\infty}(0,T)}\\ = &\sum\limits_{i = 1}^{n}\left\|\frac{(w_{i}-\xi_{i}^{\varepsilon})p_{i}^{\varepsilon}}{c_{i}}- \frac{(w_{i}-\xi_{i}^{\varepsilon})p_{i}^{\varepsilon}}{c_{i}}-\frac{\sqrt{\varepsilon}\theta_{i}^{\varepsilon}}{c_{i}}\right\|_{L^{\infty}(Q)} +\left\|\frac{\xi_{2n+1}^{\varepsilon}}{c_{n+1}}-\frac{\xi_{2n+1}^{\varepsilon}}{c_{n+1}}-\frac{\sqrt{\varepsilon}\theta_{n+1}^{\varepsilon}}{c_{n+1}}\right\|_{L^{\infty}(0,T)}\\ \leq&\sqrt{\varepsilon}\sum\limits_{i = 1}^{n}\frac{\|\theta_{i}^{\varepsilon}(x,t)\|_{L^{\infty}(Q)}}{c_{i}} +\sqrt{\varepsilon}\frac{\|\theta_{n+1}^{\varepsilon}(t)\|_{L^{\infty}(0,T)}}{c_{n+1}}\\ \leq&\sqrt{\varepsilon}\sum\limits_{i = 1}^{n+1}\frac{1}{c_{i}}, \end{align*} |
it is easy to derive that
\|(u^{\ast},v^{\ast})-(v^{\varepsilon},v^{\varepsilon})\|_{\infty}\\ \leq\|\mathcal{L}((u^{\ast},v^{\ast}))-\mathcal{L}(u^{\varepsilon},v^{\varepsilon})\|_{\infty}+\|\mathcal{L}(u^{\varepsilon},v^{\varepsilon})-(u^{\varepsilon},v^{\varepsilon})\|_{\infty}\\ \leq T\left(\sum\limits_{i = 1}^{n}\frac{1}{c_{i}}(\overline{w}_{i}K_{1}+C_{2}K_{1}+C_{1}K_{3})+K_{3}\right) \cdot\left(\sum\limits_{i = 1}^{n}\|u_{i}^{\ast}-u_{i}^{\varepsilon}\|_{L^{\infty}(Q)} +\|v^{\ast}-v^{\varepsilon}\|_{L^{\infty}(0,T)}\right)\\+\sqrt{\varepsilon}\sum\limits_{i = 1}^{n+1}\frac{1}{c_{i}}. |
So, if T is small enough, the following result holds:
\begin{align*} &\sum\limits_{i = 1}^{n}\|u_{i}^{\ast}-u_{i}^{\varepsilon}\|_{L^{\infty}(Q)}+\|v^{\ast}-v^{\varepsilon}\|_{L^{\infty}(0,T)} \leq\frac{\sqrt{\varepsilon}\sum\limits_{i = 1}^{n+1}\frac{1}{c_{i}}}{1-T\left(\sum\limits_{i = 1}^{n}\frac{1}{c_{i}}(\overline{w}_{i}K_{1}+C_{2}K_{1}+C_{1}K_{3})+K_{3}\right)}. \end{align*} |
which gives the desired result.
Finally, passing to the limit \varepsilon\rightarrow0^{+} in the inequality of (5.2) and using Lemma 5.1 yield \widetilde{J}(u^{\ast}, v^{\ast})\geq\limsup_{(u, v)\in\Omega}\widetilde{J}(u, v), which finishes the proof.
Our goal is to obtain a numerical approximation for the nonnegative T-periodic solution of the system (2.3). We discuss the problem of evolution of a single species in a polluted environment. Suppose the computational domain \tilde{Q} = [0, l]\times[0, \tilde{T}] is divided into an J\times N mesh with the spatial step size h = \frac{l}{J} = 0.01 in the x direction and time step size \tau = \frac{\tilde{T}}{N} = 0.03 . The grid points (x_{j}, t_{n}) are defined by
\begin{equation} x_{j} = jh,\; \; j = 0,1,2,\ldots,J;\nonumber \end{equation} |
\begin{equation} t_{n} = n\tau,\; \; n = 0,1,2,\ldots,N,\nonumber \end{equation} |
where J and N are two integers. The notation p_{j}^{n} denotes the solution p(jh, n\tau) of the difference equation.
Based on the above analysis, the finite difference scheme of system (2.3) can be written as follows:
\begin{equation} \frac{p_{j}^{n}-p_{j}^{n-1}}{\tau}+V\frac{p_{j}^{n}-p_{j-1}^{n}}{h}+V_{x}p_{j}^{n}+\mu p_{j}^{n}-f_{j}^{n} = 0, \end{equation} | (6.1) |
where j = 1, 2, \ldots, J; n = 1, 2, \ldots, N. It follows from (6.1) that
\begin{equation} -\lambda V p_{j-1}^{n}+[1+\lambda V+\tau(V_{x}+\mu)]p_{j}^{n} = p_{j}^{n-1}+\tau f_{j}^{n}, \end{equation} | (6.2) |
where \lambda = \frac{\tau}{h} .
The boundary and initial conditions can be discretized as
\begin{align} \left\{ \begin{array}{ll} p_{j}^{0} = p_{0j},\\ p_{0}^{n} = \sum\limits_{j = 1}^{J}\beta_{j}p_{j}h. \end{array} \right. \end{align} | (6.3) |
From (6.2) and (6.3), we have
\begin{equation} {\bf{AP}}^{n} = {\bf{P}}^{n-1}+\tau{\bf{F}}, \end{equation} | (6.4) |
where
\begin{align*} {\bf{A}} = \left[ \begin{array}{ccccc} 1+\lambda V+\tau(V_{x}+\mu)-\lambda V\beta h & -\lambda V\beta h & \ldots & -\lambda V\beta h &-\lambda V\beta h\\ -\lambda V & 1+\lambda V+\tau(V_{x}+\mu) & \\ & \ddots & \ddots & \ddots \\ & & & -\lambda V & 1+\lambda V+\tau(V_{x}+\mu) \\ & & & & -\lambda V \\ \end{array} \right], \end{align*} |
\begin{equation} {\bf{P}}^{n} = (p_{1}^{n},p_{2}^{n},\ldots,p_{J}^{n})^{{\mathrm{T}}},\; \; {\bf{F}} = (f_{1}^{n},f_{2}^{n},\ldots,f_{J}^{n})^{{\mathrm{T}}}.\nonumber \end{equation} |
Note that {\bf{A}} is an upper triangular matrix, so the nonlinear algebraic equation (6.4) have solutions.
We first prove the stability of the implicit difference scheme (6.1) by using the Fourier method. We suppose that f\equiv0 . The domain of the function defined on the grid point is extended according to the usual method, that is, when x\in(x_{j}-\frac{h}{2}, x_{j}+\frac{h}{2}) , p^{n}(x) = p_{j}^{n} ; then, we have
\begin{equation} -\lambda V p^{n}(x-h)+[1+\lambda V+\tau(V_{x}+\mu)]p^{n}(x) = p^{n-1}(x),x\in R. \end{equation} | (6.5) |
Taking the Fourier transform of both sides of (6.5), we obtain
\begin{equation} -\lambda V \hat{U}^{n}(k)e^{-{\mathrm{i}}kh}+[1+\lambda V+\tau(V_{x}+\mu)]\hat{U}^{n}(k) = \hat{U}^{n-1}(k). \end{equation} | (6.6) |
By (6.6), we get
\begin{equation} \hat{U}^{n}(k) = \frac{\hat{U}^{n-1}(k)}{-\lambda Ve^{-{\mathrm{i}}kh}+1+\lambda V+\tau(V_{x}+\mu)}.\nonumber \end{equation} |
Thus, the growth factor is as follows:
\begin{align*} G(\tau,k)& = \frac{1}{-\lambda Ve^{-{\mathrm{i}}kh}+1+\lambda V+\tau(V_{x}+\mu)}\\ & = \frac{1}{1+\lambda V(1-\cos kh)+\tau(V_{x}+\mu)+\lambda V{\mathrm{i}}\sin kh}. \end{align*} |
Then, we have
\begin{align*} |G(\tau,k)|^{2}& = \frac{1}{[1+\lambda V(1-\cos kh)+\tau(V_{x}+\mu)]^{2}+\lambda^{2}V^{2}\sin^{2}kh}\\ & = \frac{1}{[1+2\lambda V\sin^{2}\frac{kh}{2}+\tau(V_{x}+\mu)]^{2}+4\lambda^{2}V^{2}\sin^{2}\frac{kh}{2}(1-\sin^{2}\frac{kh}{2})}. \end{align*} |
If V_{x}+\mu > 0 , there is |G(\tau, k)|\leq1 , that is, the von Neumann condition is satisfied, the difference scheme (6.1) is stable under condition V_{x}+\mu > 0 .
Next, we analyze the convergence of the difference scheme (6.1). Note that the difference scheme is compatible, we thus use the Lax equivalence theorem. The difference scheme must be convergent and have first-order accuracy.
We chose the following parameters:
\left\{\begin{array}{ll} \beta(x,c_{0}(t)) = 100x^{2}(1-x)(1+\sin(\pi x))\left|\sin\frac{2\pi c_{0}(t)}{T}\right|,\\ \mu(x,c_{0}(t) = e^{-4x}(1-x)^{-1.4}(2+\cos\frac{2\pi c_{0}(t)}{T}),\\ V(t,x) = 1-x,\; f(x,t) = 2+(1+x)\sin(\frac{2\pi t}{T}),\\ p_{0}(x) = e^{x},\; u(x,t) = 0,\; x = 1,\; T = \frac{1}{3},\; \tilde{T} = 3. \end{array}\right. |
In this paper, we used the backward difference scheme and chasing method, and (6.4) was solved through programming. The fertility rate, mortality rate, and immigration rate were T-periodic and were all greater than zero, which is consistent with the assumptions. We considered T = \frac{1}{3} . Their graphs are given in Figures 1–3, respectively. The fertility rate was the highest when the size was half and the mortality rate was the highest when the size was the maximum, which conformed to the empirical situation. Therefore, the selection of parameters \beta , \mu , and f was reasonable.
The graph of the numerical solution p is given in Figure 4. Over time, solution p showed T-periodic changes. We take the numerical solution of system (2.3), corresponding to an arbitrary positive initial datum p_{0} , on some interval [kT, (k+1)T] , where k is large enough. We can then get the periodic solution of system (2.3) by extending the numerical solution p . Over time, the density of the population increased first and then decreased. With respect to size, the density of the population gradually decreased, and finally tended to be stable. The maximum value was reached at x = 0 . On such an interval with a sufficiently large k , the numerical solution was already stable. During computation we found that any positive initial datum p_{0} was appropriate for use.
In this paper, we considered the problem of optimal harvesting for a periodic n -dimensional food chain model dependent on the size structure, that combined toxicants with the population. In the foregoing, we have discussed the existence and uniqueness of a nonnegative solution of the state system. The necessary optimality conditions were provided, and the existence of the optimal policy was investigated. Some numerical results were finally given. The results implied that the solution of (2.3) always maintains the pattern of increasing periodically, and any positive initial datum p_{0} is appropriate. From Theorem 4.1, the optimal strategy had a bang-bang structure and provided threshold conditions for the optimal control problems (2.3) and (2.4). The bang-bang structure of solutions is much more common in optimal population management.
We now comment on the differences between the methods and results of this paper as well as closely related work. The existence of the optimal strategy was proved by compactness and the extreme sequence in [13]. In this paper, the corresponding problem was treated by using Ekeland variational principle. The authors [10] established the maximum principle and bang-bang structure for optimal control but paid no attention to existence and uniqueness results. The existence of the optimal harvesting rate in the harvesting problem was only given in [11]. Furthermore, if V_{i}(x, t) = 1 for Q = (0, l) \times R_{+}, i = 1, 2\ldots, n , the state system degenerates into an age-structured model, and our results cover the corresponding results [4,5,6].
Note that the individual price factor w_{i}(x, t) plays an important role in the structure of the optimal controller (4.1). However, as we do not have a clear biological meaning for the solutions \xi_{i} (i = 1, 2, \ldots, 2n+1) of the adjoint system (4.4), it is difficult to give a precise explanation of the threshold conditions (4.1) and (4.2). In specific applications, the optimal population density and optimal policy are calculated by combining the state system and the adjoint system. In this way, we enable the objective functional (2.4) to maximize the total profit during period T . This is a challenging problem, and future work in the area should address it. In addition, inspired by [17,18], we intend to consider the problem of fractional optimal control problem of the population model in future research.
The authors thank the referees for their valuable comments and suggestions on the original manuscript that helped improve its quality. The work was supported by the National Natural Science Foundation of China under grant 11561041.
None of the authors has a conflict of interest in the publication of this paper.
[1] |
Joly CA, Metzger JP, Tabarelli M (2014) Experiences from the Brazilian Atlantic Forest: ecological findings and conservation initiatives. New Phytol 204: 459-473. doi: 10.1111/nph.12989
![]() |
[2] | Ribeiro MC, Metzger JP, Martensen AC, et al. (2009) The Brazilian Atlantic Forest: how much is left, and how is the remaining forest distributed? Implications for conservation. Biol Conserv 142: 1141-1153. |
[3] |
Rodrigues RR, Lima RAF, Gandolfi S, et al. (2009) On the restoration of high diversity forests: 30 years of experience in the Brazilian Atlantic Forest. Biol Conserv 142: 1242-1251. doi: 10.1016/j.biocon.2008.12.008
![]() |
[4] | Coelho GC (2010) Restauração florestal em pequenas propriedades: desafios e oportunidades. In: Hüller A (Ed.) Gestão Ambiental nos Municípios: Instrumentos e Experiências na Administração Pública. Santo Ângelo: FURI, 195-215. |
[5] |
Rodrigues RR, Gandolfi S, Nave AG, et al. (2011) Large-scale ecological restoration of high-diversity tropical forests in SE Brazil. For Ecol Manage 261: 1605-1613. doi: 10.1016/j.foreco.2010.07.005
![]() |
[6] |
Sampaio AB, Holl KD, Scariot A (2007) Does restoration enhance regeneration of seasonal deciduous forests in pastures in central Brazil? Rest Ecol 15: 462-471. doi: 10.1111/j.1526-100X.2007.00242.x
![]() |
[7] |
Engel VL, Parrotta JA (2001) An evaluation of direct seeding for reforestation of degraded lands in central São Paulo State, Brazil. For Ecol Manage 152: 169-181. doi: 10.1016/S0378-1127(00)00600-9
![]() |
[8] | Kageyama PY, Castro CFA, Carpanezzi AA 1989. Implantação de matas ciliares: estratégia para auxiliar a sucessão secundária. In: Barbosa L M (ed.) Simpósio sobre mata ciliar. Fundação Cargill, Campinas, Brazil, 130-143. |
[9] |
Camargo JLC, Ferraz IDK, Imakawa AM (2002). Rehabilitation of degraded areas of Central Amazonia using direct sowing of forest tree seeds. Rest Ecol 10: 636-644. doi: 10.1046/j.1526-100X.2002.01044.x
![]() |
[10] | Uhl C, Buschbacher R, Serrão EAS (1988). Abandoned pastures in Eastern Amazonia. I. Patterns of plant succession. J Ecol 76: 663- 681. |
[11] | Pompéia S (2005) Recuperação da vegetação da Serra do Mar em áreas afetadas pela poluição atmosférica de Cubatão: uma análise histórica. In: Galvão APM, Porfírio-Da-Silva V (Eds.). Restauração Florestal—Fundamentos e Estudos de Caso. Colombo: Embrapa Florestas, 119-143. |
[12] |
Campoe OC, Stape JL, Mendes JCT (2010) Can intensive management accelerate the restoration of Brazil’s Atlantic forests? For Ecol Manage 259: 1808-1814. doi: 10.1016/j.foreco.2009.06.026
![]() |
[13] | Barbosa LM, Barbosa JM, Barbosa KC, et al. (2003) Recuperação florestal com espécies nativas no Estado de São Paulo: pesquisas apontam mudanças necessárias. Florestar estat 6: 28-34. |
[14] |
Martínez-Garza C, Howe HF (2003) Restoring tropical diversity: beating the time tax on species loss. J Appl Ecol 40: 423-429. doi: 10.1046/j.1365-2664.2003.00819.x
![]() |
[15] |
Newmaster SG, Bell FW, Roosenboom CR, et al. (2006) Restoration of floral diversity through plantations on abandoned agricultural land. Can J For Res 36: 1218-1235. doi: 10.1139/x06-021
![]() |
[16] |
Trindade DFV, Coelho GC (2012) Woody species recruitment under monospecific plantations of pioneer trees - facilitation or inhibition? iForest 5: 1-5. doi: 10.3832/ifor0601-009
![]() |
[17] |
Souza FM, Batista JLF (2004) Restoration of seasonal semideciduous forests in Brazil: influence of age and restoration design on forest structure. For Ecol Manage 191: 185-200. doi: 10.1016/j.foreco.2003.12.006
![]() |
[18] |
Denslow JS (1987) Tropical rainforest gaps and tree species diversity. Ann Rev Ecol Syst 18: 431-451. doi: 10.1146/annurev.es.18.110187.002243
![]() |
[19] |
Swaine MD, Whitmore TC (1988) On the definition of ecological species groups in tropical rain forests. Vegetatio 75: 81-86. doi: 10.1007/BF00044629
![]() |
[20] | Budowski G (1965) Distribution of American rain forest species in the light of successional process. Turrialba 15: 40-42. |
[21] |
Denslow JS (1980) Gap partioning among tropical rain forest trees. Biotropica 12: 47-55. doi: 10.2307/2388156
![]() |
[22] | Liebsch D, Marques MC, Goldenberg R (2008) How long does the Atlantic Rain Forest take to recover after a disturbance? Changes in species composition and ecological features during secondary succession. Biol Conserv 141: 1717-1725. |
[23] | Ferretti AR, Kageyama PY, Árbocz GF, et al. (1995) Classificação das espécies arbóreas em grupos ecofisiológicos para revegetação com nativas no estado de São Paulo. Florestar estat 3: 73-77. |
[24] | Coelho GC, Rigo MS, Libardoni JB, et al. (2011) Understory structure in two successional stage of Semi-deciduous Seasonal Forest remnant of Southern Brazil. Biota Neotrop 11: 63-74. |
[25] |
Parrotta JA, Turnbull JW, Jones N (1997) Catalyzing native forest regeneration on degraded tropical lands. For Ecol Manage 99: 1-7. doi: 10.1016/S0378-1127(97)00190-4
![]() |
[26] |
Holl KD, Aide TM (2011) When and where to actively restore ecosystems? For Ecol Manage 261: 1558-1563. doi: 10.1016/j.foreco.2010.07.004
![]() |
[27] | Kageyama PY, Gandara FB (2000) Recuperação de áreas ciliares. In: Rodrigues R R, Leitão Filho H F (Eds). Mata ciliares: uma abordagem multidisciplinar. São Paulo: EDUSP/FAPESP. 249-269. |
[28] | Hooper E, Condit R, Legendre P (2002) Responses of 20 native tree species to reforestation strategies for abandoned farmland in Panama. Ecol Appl 12: 1626-1641. |
[29] |
St-Denis A, Messier C, Kneeshaw D (2013) Seed size, the only factor positively affecting direct seeding success in an abandoned field in Quebec, Canada. Forests 4: 500-516. doi: 10.3390/f4020500
![]() |
[30] |
Westoby M, Falster DS, Moles AT, et al. (2002) Plant ecological strategies: some leading dimensions of variation between species. Ann Rev Ecol Syst 33: 125-159. doi: 10.1146/annurev.ecolsys.33.010802.150452
![]() |
[31] |
Howe HF, Richter WM (1982) Effects of seed size on seedling size in Virola surinamensis: a within and between tree analysis. Oecologia 53: 347-351. doi: 10.1007/BF00389011
![]() |
[32] |
Green PT, Juniper PA (2004) Seed mass, seedling herbivory and the reserve effect in tropical rainforest seedlings. Funct Ecol 18: 539-547. doi: 10.1111/j.0269-8463.2004.00881.x
![]() |
[33] |
Osunkoya OO, Ash JE, Hopkins MS, et al. (1992) Factors affecting survival of tree seedlings in North Queensland rainforests. Oecologia 91: 569-578. doi: 10.1007/BF00650333
![]() |
[34] |
Osunkoya OO, Ash JE, Hopkins MS, et al. (1994) Influence of seed size and seedling ecological attributes on shade-tolerance of rain-forest tree species in northern Queensland. J Ecol 82: 149-163. doi: 10.2307/2261394
![]() |
[35] |
Baraloto C, Forget PM, Goldberg DE (2005) Seed mass, seedling size and neotropical tree seedling Establishment. J. Ecol 93: 1156-1166. doi: 10.1111/j.1365-2745.2005.01041.x
![]() |
[36] | Poorter L, Wright SJ, Paz H, et al. (2008) Are functional traits good predictors of demographic rates? Evidence from five neotropical forests. Ecology 89: 1908-1920. |
[37] |
Cole RJ, Holl KD, Keene CL, et al. (2011) Direct seeding of late-successional trees to restore tropical montane Forest. For Ecol Manage 261: 1590-1597. doi: 10.1016/j.foreco.2010.06.038
![]() |
[38] | Benvenuti-Ferreira G, Coelho GC, Schirmer J, et al. ((2009) Dendrometry and litterfall of neotropical pioneer and early secondary tree species. Biota Neotrop 9: 65-71. |
[39] | Carvalho PER (2003-2014). Espécies arbóreas brasileiras (Vols. 1-5). Brasília: Embrapa Informação Tecnológica. |
[40] | Alvares CA, Stape JL, Sentelhas PC, et al. (2013) Köppen’s climate classification map for Brazil. Met Zeit 22: 711-728. |
[41] |
Ramírez-Marcial N, González-Espinosa M, Camacho-Cruz A, et al. (2010) Forest restoration in Lagunas de Montebello National Park, Chiapas, Mexico. Ecol Rest 28: 354-360. doi: 10.3368/er.28.3.354
![]() |
[42] |
Martínez-Garza VP, Ricker M, Campos A, et al. (2005) Restoring tropical biodiversity: leaf traits predict growth and survival of late-successional trees in early-successional environments. For Ecol Manage 217: 365-379. doi: 10.1016/j.foreco.2005.07.001
![]() |
[43] |
Siddique I, Engel VL, Parrotta JA, et al. (2008) Dominance of legume trees alters nutrient relations in mixed species forest restoration plantings within seven years. Biogeochemistry 88: 89-101. doi: 10.1007/s10533-008-9196-5
![]() |
[44] |
Martínez-Garza C, Bongers F, Poorter L (2013) Are functional traits good predictors of species performance in restoration plantings in tropical abandoned pastures? For Ecol Manage 303: 35-45. doi: 10.1016/j.foreco.2013.03.046
![]() |
[45] |
Moles AT, Westoby M (2006) Seed size and plant strategy across the whole life cycle. Oikos 113: 91-105. doi: 10.1111/j.0030-1299.2006.14194.x
![]() |
[46] |
Tunjai P, Elliot S (2012) Effects of seed traits on the success of direct seeding for restoring southern Thailand’s lowland evergreen forest ecosystem. New for 43: 319-333. doi: 10.1007/s11056-011-9283-7
![]() |
[47] |
Flores O, Hérault B, Delcamp M, et al. (2014) Functional traits help predict post-disturbance demography of tropical trees. PloS one 9: e105022. doi: 10.1371/journal.pone.0105022
![]() |
[48] | Melo FPL, Lemire D, Tabarelli M (2007) Extirpation of large-seeded seedlings from the edge of a large Brazilian Atlantic forest fragment. Écoscience 14: 124-129. |
[49] |
Easdale TA, Healey JR (2009) Resource-use-related traits correlate with population turnover rates, but not stem diameter growth rates, in 29 subtropical montane tree species. Persp Plant Ecol Evol Syst 11: 203-218. doi: 10.1016/j.ppees.2009.03.001
![]() |
[50] | Stehmann JR, Forzza RC, Salino A, et al. (2009) Plantas da Floresta Atlântica. Rio de Janeiro: Jardim Botânico do Rio de Janeiro. |
[51] | Oliveira-Filho AT, Budke AT, Jarenkow JC, et al. (2015) Delving into the variations in tree species composition and richness across South American subtropical Atlantic and Pampean forests. J Plant Ecol 8: 242-260. |
1. | Juan Liu, Jie Hu, Peter Yuen, Fuzhong Li, A Seasonally Competitive M-Prey and N-Predator Impulsive System Modeled by General Functional Response for Integrated Pest Management, 2022, 10, 2227-7390, 2687, 10.3390/math10152687 | |
2. | Kareem T. Elgindy, Numerical solution of nonlinear periodic optimal control problems using a Fourier integral pseudospectral method, 2024, 144, 09591524, 103326, 10.1016/j.jprocont.2024.103326 | |
3. | Fenghui Qin, Zhanping Wang, Optimal Control of Forest Population System with Size Structure in Polluted Environment, 2025, 206, 0022-3239, 10.1007/s10957-025-02687-4 | |
4. | Fenghui Qin, Zhanping Wang, Optimal boundary control for a competitive population system with size structure and time delay in a polluted environment, 2025, 18, 1793-5245, 10.1142/S1793524524500141 |