Mamadou Birba
Laboratoire LAMI, Université Joseph KI-ZERBO, Ouaga, Burkina Faso.
DOI: 10.4236/am.2026.179033   PDF    HTML   XML   2 Downloads   33 Views   Citations

Abstract

This work is devoted to the study of a class of nonlinear degenerate parabolic systems arising from the modelling of two age-structured populations interacting within a one-dimensional heterogeneous spatial environment. The system is characterized by nonlinear degenerate boundary conditions and nonlocal renewal laws. The presence of cross-diffusion operators $\text{div}\left(v{\partial }_{x}u\right)$ and $\text{div}\left(u{\partial }_{x}v\right)$ leads to coupling between the equations and, simultaneously, strong degeneracy. Consequently, we propose a method based on regularizing the system using a family of highly regular functions. We then establish the existence of global weak solutions using the Faedo-Galerkin method, with a finite element basis. Uniform estimates and compactness arguments subsequently allow us to perform a double limiting procedure. We finally show that the solution of the regularized system converges to that of the original degenerate system, while preserving in particular the nonlocal boundary conditions and the initial data.

Keywords

Cross-Diffusion System, Degenerate System, Weak Solutions, Faedo-Galerkin, Degenerate Boundary Conditions

Share and Cite:

1. Introduction

The mathematical modeling of structured populations is a fundamental area of mathematical biology and theoretical ecology. Since the pioneering work of Kermack and McKendrick [1] and Von Foerster [2], age-structured equations have demonstrated their power to describe population dynamics where the age of individuals influences their demographic characteristics (mortality rates, fertility, migration). When these populations evolve in a spatially heterogeneous environment, the introduction of diffusion operators makes it possible to capture the phenomena of dispersal and territorial colonization.

In many real biological systems, several populations interact in complex ways: competition for resources, predator-prey models [3] [4]. For mutualism, or epidemiological interactions, see for example [5] [6]. These interactions often give rise to nonlinear coupling in terms of diffusion, reflecting phenomena such as the tendency to avoid local competition or to follow density gradients of other species. In this paper, we are interested in the existence of solutions of the following nonlinear system

$\{\begin{array}{ll}\frac{\partial y}{\partial t}+\frac{\partial y}{\partial a}-\frac{\partial }{\partial x}\left(z\frac{\partial y}{\partial x}\right)+{\mu }_1\left(t,a\right)y=0,\hfill & \left(t,a,x\right)\in Q,\hfill \\ \frac{\partial z}{\partial t}+\frac{\partial z}{\partial a}-\frac{\partial }{\partial x}\left(y\frac{\partial z}{\partial x}\right)+{\mu }_2\left(t,a\right)z=0,\hfill & \left(t,a,x\right)\in Q,\hfill \\ y\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_1\left(a,t\right)y\left(t,a,x\right)\text{d}a},\hfill & \left(t,x\right)\in Q_T,\hfill \\ z\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_2\left(a,t\right)z\left(t,a,x\right)\text{d}a},\hfill & \left(t,x\right)\in Q_T,\hfill \\ y\left(0,a,x\right)=y_0\left(a,x\right),z\left(0,a,x\right)=z_0\left(a,x\right),\hfill & \left(a,x\right)\in Q_A,\hfill \\ z\frac{\partial y}{\partial x}\left(t,a,0\right)=z\frac{\partial y}{\partial x}\left(t,a,L\right)=0,\hfill & \left(t,a\right)\in \Sigma ,\hfill \\ y\frac{\partial z}{\partial x}\left(t,a,0\right)=y\frac{\partial z}{\partial x}\left(t,a,L\right)=0,\hfill & \left(t,a\right)\in \Sigma ,\hfill \end{array}$ (1)

where $y\left(t,a,x\right)$ and $z\left(t,a,x\right)$ are respectively the densities of two populations of age $a$ , at time $t$ and position $x$ , and are assumed to be nonnegative, $T>0$ , $A$ is the maximum age, $Q=\left(0,T\right)\times \left(0,A\right)\times \left(0,L\right)$ , $Q_T=\left(0,T\right)\times \left(0,L\right)$ , $Q_A=\left(0,A\right)\times \left(0,L\right)$ and $\Sigma =\left(0,T\right)\times \left(0,A\right)$ . The main novelty and difficulty of our system lies in the terms $\frac{\partial }{\partial x}\left(z\frac{\partial y}{\partial x}\right)$ and $\frac{\partial }{\partial x}\left(y\frac{\partial z}{\partial x}\right)$ . This is a cross-diffusion type coupling, which makes the system strongly nonlinear and degenerate, presenting substantial analytical challenges while corresponding to relevant biological situations, such as cell populations interacting in a biological tissue or competing species in an ecosystem.

Specifically, system (1) exhibits several characteristics that make it a mathematically rich and demanding subject of study. First, there is the dual spatiotemporal structuring. The unknowns $y\left(t,a,x\right)$ and $z\left(t,a,x\right)$ depend simultaneously on time $t$ , age $a$ , and spatial position $x$ . This triple dependence necessitates the use of anisotropic function spaces and complicates obtaining a priori estimates. Second, the principal operators exhibit degenerate nonlinearity. The terms $\text{div}\left(z{\partial }_{x}y\right)$ and $\text{div}\left(y{\partial }_{x}z\right)$ are nonlinear, nonmonotonic, and potentially degenerate when $y$ or $z$ vanishes. This degeneracy prevents the direct application of classical theories to parabolic equations. Third, we have atypical boundary conditions. The conditions $z\frac{\partial y}{\partial x}=y\frac{\partial z}{\partial x}=0$ on the boundary are non-standard and degenerate: they provide no information about the flux when one of the densities vanishes at the boundary. This situation requires a carefully chosen weak formulation. We also have non-local birth conditions. The conditions

$y\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_1\left(a,t\right)y\left(t,a,x\right)\text{d}a},z\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_2\left(a,t\right)z\left(t,a,x\right)\text{d}a},$

introduce a non-local coupling between all ages at each time step, typical of renewal models in population dynamics. Finally, there is a strong coupling between the equations. The diffusion of each population is controlled by the density of the other population, creating a bilateral interdependence that complicates the analysis of existence and uniqueness.

Our work lies at the intersection of several active research fields: with regard to structured age equations: we generalize classical Gurtin-MacCamy [7] models by introducing nonlinear and coupled spatial diffusion. While the work of Busenberg and Iannelli [8] and Anita [9] primarily deals with linear diffusion or nonlinear reactions, our model incorporates diffusion operators that depend on the unknowns in a cross-referencing manner.

With regard to reaction-diffusion systems: our system differs from Lotka-Volterra type systems with diffusion studied, for example, by Conway, Hoff, and Smoller [10], by the presence of age structuring and by the specific form of the diffusion terms.

Regarding degenerate equations: the potentially degenerate nature of our system brings it closer to porous medium equations [11] or cross-diffusion systems [12]. However, the combination with age structuring is novel and requires adapted techniques. Regarding non-standard boundary conditions: the conditions $z\frac{\partial y}{\partial x}=y\frac{\partial z}{\partial x}=0$ exhibit a formal similarity to “zero flux” type conditions but with additional degeneracy when the densities vanish on the boundary. Moreover, for general degenerate cross-diffusion systems without age-structure, the question of well-posedness has been addressed by Choquet et al. [13]. Their work establishes existence and uniqueness of entropy solutions for a broad class of such systems using techniques such as the doubling of variables method. However, the presence of age-structuring in our model (1) introduces an additional layer of complexity, requiring a specific approximation (Faedo-Galerkin) and a priori estimates that are not covered by their framework. Our result complements theirs by providing an existence proof for the coupled age-space degenerate problem.

Introduced by Shigesada et al. [14] with a system known as SKT, cross-diffusion models aim to describe the interaction between species that tend to avoid each other. Cross-diffusion describes how the population flow of a given species is affected by the presence of other species. This type of model has been researched by numerous authors, including: Levin and Segel [15]; Okubo and Levin [16]; Mimura and Murray [17]; Mimura and Kawasaki [18]; Mimura and Yamaguti [19]; Bendahmane et al. [20]; to name just a few. Regarding age-structured models, the pioneering work of Gurtin and MacCamy [21] can be cited.

The remainder of this paper is organized as follows. In Section 0, we introduce a non-degenerate approximation of the original system by regularizing the crossdiffusion terms with a family of smooth functions $f_{\epsilon }$ . We then prove the existence of weak solutions to this approximate system using the Faedo-Galerkin method, a fixed point argument, and uniform a priori estimates. In Section, we establish uniform bounds with respect to the regularization parameter $\epsilon$ and prove the positivity of the solutions. Finally, in Section, we pass to the limit $\epsilon \to 0$ to recover a weak solution of the original degenerate cross-diffusion system (1). The main existence result is stated in Theorem 6.

2. The Well Posedness of Auxiliary System

To overcome the difficulty of higher-order derivative coupling, we introduce a specific approximation. The system (1) is a degenerate parabolic system and therefore difficult to analyze in its original form. One of the difficulties in analyzing of system lies the managing of the coupling of higher-order derivatives. To overcome this difficulty, we introduce a specific approximation based on the following function.

Let $\epsilon >0$ and $f_{\epsilon }\in C^{\infty }\left(ℝ\right)$ be a function satisfying

$\{\begin{array}{ll}\left|f_{\epsilon }\left(s\right)-f_{\epsilon }\left(s^{\prime }\right)\right|\le \left|s-s^{\prime }\right|,\hfill & \forall s,s^{\prime }\in ℝ,\hfill \\ \epsilon \le f_{\epsilon }\left(s\right)\le C_{\epsilon }\left(1+\epsilon \right),\hfill & \forall s\in ℝ,\hfill \\ \underset{\epsilon \to 0^{+}}{\lim }{f}_{\epsilon }\left(s\right)=s,\hfill & \forall s\in {ℝ}^{,}\hfill \end{array}$ (2)

where $C_{\epsilon }>0$ is a constant dependent of $\epsilon$ .

For example, we can take

$f_{\epsilon }\left(s\right)=\frac{s}{1+\epsilon s}+\epsilon ,\forall s\in ℝ,$

which satisfies all the conditions 2.

2.1. The Non-Degenerate Approximate System

In this subsection, we introduce a specific approximation of the system (1) and then prove the existence of weak solutions of the approximate problem by using the Faedo-Galerkin method.

Consider the approximate system

$\{\begin{array}{ll}\frac{\partial y^{\epsilon }}{\partial t}+\frac{\partial y^{\epsilon }}{\partial a}-\frac{\partial }{\partial x}\left(f_{\epsilon }\left(z^{\epsilon }\right)\frac{\partial y^{\epsilon }}{\partial x}\right)+{\mu }_1\left(t,a\right)y^{\epsilon }=0,\hfill & \left(t,a,x\right)\in Q,\hfill \\ \frac{\partial z^{\epsilon }}{\partial t}+\frac{\partial z^{\epsilon }}{\partial a}-\frac{\partial }{\partial x}\left(f_{\epsilon }\left(y^{\epsilon }\right)\frac{\partial z^{\epsilon }}{\partial x}\right)+{\mu }_2\left(t,a\right)z^{\epsilon }=0,\hfill & \left(t,a,x\right)\in Q,\hfill \\ y^{\epsilon }\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_1\left(a,t\right)y^{\epsilon }\left(t,a,x\right)\text{d}a},\hfill & \left(t,x\right)\in Q_T,\hfill \\ z^{\epsilon }\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_2\left(a,t\right)z^{\epsilon }\left(t,a,x\right)\text{d}a},\hfill & \left(t,x\right)\in Q_T,\hfill \\ y^{\epsilon }\left(0,a,x\right)=y_0\left(a,x\right),z^{\epsilon }\left(0,a,x\right)=z_0\left(a,x\right),\hfill & \left(a,x\right)\in Q_A,\hfill \\ f_{\epsilon }\left(z^{\epsilon }\right)\frac{\partial y^{\epsilon }}{\partial x}\left(t,a,0\right)=f_{\epsilon }\left(z^{\epsilon }\right)\frac{\partial y^{\epsilon }}{\partial x}\left(t,a,L\right)=0,\hfill & \left(t,a\right)\in \Sigma ,\hfill \\ f_{\epsilon }\left(y^{\epsilon }\right)\frac{\partial z^{\epsilon }}{\partial x}\left(t,a,0\right)=f_{\epsilon }\left(y^{\epsilon }\right)\frac{\partial z^{\epsilon }}{\partial x}\left(t,a,L\right)=0,\hfill & \left(t,a\right)\in \Sigma .\hfill \end{array}$ (3)

Remark. The function $f_{\epsilon }$ guarantees uniform parabolic regularity of the regularized system, preserves the physical behaviour for positive densities, and provides $C^{\infty }\left({ℝ}^{*}+\right)$ regularity for the analysis.

Definition 1. A pair $\left(y^{\epsilon },z^{\epsilon }\right)$ is called a weak solution of the regularized system (3) if:

$y^{\epsilon },z^{\epsilon }\in L^{\infty }\left(0,T;L^2\left(0,A;L^2\left(\Omega \right)\right)\right)\cap L^2\left(0,T;L^2\left(0,A;H^1\left(\Omega \right)\right)\right)$

with

$\left({\partial }_t+{\partial }_a\right)y^{\epsilon },\left({\partial }_t+{\partial }_a\right)z^{\epsilon }\in L^2\left(0,T;L^2\left(0,A;{\left(H^1\left(\Omega \right)\right)}^{*}\right)\right),$

and for all ${\varphi }_1,{\varphi }_2\in C_c^{\infty }\left(\left[0,T\right)\times \left[0,A\right)\times \Omega \right)$ ,

$\{\begin{array}{l}{\displaystyle {\int }_Q\left[\left({\partial }_t+{\partial }_a\right)y^{\epsilon }{\varphi }_1+f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }{\partial }_x{\varphi }_1+{\mu }_1{y}^{\epsilon }{\varphi }_1\right]\text{d}t\text{d}a\text{d}x}=0,\\ {\displaystyle {\int }_Q\left[\left({\partial }_t+{\partial }_a\right)z^{\epsilon }{\varphi }_2+f_{\epsilon }\left(y^{\epsilon }\right){\partial }_x{z}^{\epsilon }{\partial }_x{\varphi }_2+{\mu }_2{z}^{\epsilon }{\varphi }_2\right]\text{d}t\text{d}a\text{d}x}=0,\end{array}$ (4)

with initial conditions

$y^{\epsilon }\left(0,a,x\right)=y_0\left(a,x\right),z^{\epsilon }\left(0,a,x\right)=z_0\left(a,x\right)$

in the sense of traces in $L^2\left(\text{Ω}\times \left(0,A\right)\right)$ , and renewal conditions

$y^{\epsilon }\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)y^{\epsilon }\left(t,a,x\right)\text{d}a},z^{\epsilon }\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_2\left(t,a\right)z^{\epsilon }\left(t,a,x\right)\text{d}a}$

in the sense of traces in $L^2\left(\left(0,T\right)\times \Omega \right)$ .

We will include a similar definition for the original degenerate system (with $f_{\epsilon }$ replaced by the identity).

2.2. Existence of Solution to the Approximate System

We will now show the existence of solution of the approximate problem (3) using the Faedo-Galerkin method.

Hypotheses 1. The following hypotheses hold

• ${\beta }_i\in L^{\infty }\left(\left(0,T\right)\times \left(0,A\right)\right)$ , ${\beta }_i\left(t,a\right)\ge 0$ a.e., ${\mu }_i\in L_{\text{loc}}^{\infty }\left(\left[0,T\right]\times \left[0,A\right)\right)$ a.e.;

• $y_0,z_0\in L^{\infty }\left(\left(0,A\right)\times \Omega \right)$ are nonnegative.

Theorem 2. Under the assumptions on the data 1 and under the conditions 2, for each fixed $\epsilon >0$ , the system (3) admits a weak solution $\left(y^{\epsilon },z^{\epsilon }\right)$ .

Proof of Theorem 2. The proof proceeds in several steps.

Step 1: Finite-dimensional approximation.

We use a ${\mathbb{P}}^1$ finite element basis ${\left\{{\varphi }_i\right\}}_{i=1}^n$ on a uniform mesh of $\left(0,L\right)$ , orthonormalized in $L^2\left(0,L\right)$ . Let $M_{ij}={\displaystyle {\int }_0^L{\varphi }_i{\varphi }_j\text{d}x}$ and factor $M=L{L}^{\text{T}}$ (Cholesky). Define

$e_k\left(x\right)={\displaystyle \sum }_{j=1}^n{\left(L^{-\text{T}}\right)}_{jk}{\varphi }_j\left(x\right).$ (5)

The family ${\left\{e_k\right\}}_{k=1}^n$ is orthonormal in $L^2\left(0,L\right)$ and satisfies

${\Vert e_k\Vert }_{L^{\infty }\left(0,L\right)}\le C,{\Vert {\partial }_x{e}_k\Vert }_{L^{\infty }\left(0,L\right)}\le C,$ (6)

where $C$ depends on $n$ (via the inverse mesh size).

We seek approximations

$y_n\left(t,a,x\right)={\displaystyle \sum }_{k=1}^n{c}_{n,k}^y\left(t,a\right)e_k\left(x\right),z_n\left(t,a,x\right)={\displaystyle \sum }_{k=1}^n{c}_{n,k}^z\left(t,a\right)e_k\left(x\right),$ (7)

with initial conditions

$c_{n,k}^y\left(0,a\right)={\langle y_0\left(a,\cdot \right),e_k\rangle }_{L^2\left(\Omega \right)},c_{n,k}^z\left(0,a\right)={\langle z_0\left(a,\cdot \right),e_k\rangle }_{L^2\left(\Omega \right)}.$ (8)

The birth conditions become

$c_{n,k}^y\left(t,0\right)={\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)c_{n,k}^y\left(t,a\right)\text{d}a},c_{n,k}^z\left(t,0\right)={\displaystyle {\int }_0^A{\beta }_2\left(t,a\right)c_{n,k}^z\left(t,a\right)\text{d}a}.$ (9)

The weak formulation of (4) with test function $e_m$ yields

$\left(\frac{\partial }{\partial t}+\frac{\partial }{\partial a}\right)c_m^y\left(t,a\right)=g_m^1\left(t,a\right)\text{and}\left(\frac{\partial }{\partial t}+\frac{\partial }{\partial a}\right)c_m^z\left(t,a\right)=g_m^2\left(t,a\right),$ (10)

where

$\{\begin{array}{l}{g}_m^1=-{\displaystyle {\int }_0^L\left(f_{\epsilon }\left(z_n\right){\partial }_x{y}_n{\partial }_x{e}_m+{\mu }_1{y}_n{e}_m\right)\text{d}x},\\ g_m^2=-{\displaystyle {\int }_0^L\left(f_{\epsilon }\left(y_n\right){\partial }_x{z}_n{\partial }_x{e}_m+{\mu }_2{z}_n{e}_m\right)\text{d}x}.\end{array}$ (11)

Note that by integration along the characteristic lines, see [9], the unique solution of (10) associated to the conditions (8) and (9) is given by

$c_{n,m}^y\left(t,a\right)=\{\begin{array}{ll}{c}_{n,m}^y\left(t-a,0\right)+{\displaystyle {\int }_0^a{g}_m^1\left(t-a+s,s\right)\text{d}s},\hfill & \text{if}a<t,\hfill \\ c_m^1\left(a-t\right)+{\displaystyle {\int }_0^t{g}_m^1\left(s,a-t+s\right)\text{d}s},\hfill & \text{if}a\ge t,\hfill \end{array}$ (12)

and the other solution is

$c_{n,m}^z\left(t,a\right)=\{\begin{array}{ll}{c}_{n,m}^z\left(t-a,0\right)+{\displaystyle {\int }_0^a{g}_m^2\left(t-a+s,s\right)\text{d}s},\hfill & \text{if}a<t,\hfill \\ c_m^2\left(a-t\right)+{\displaystyle {\int }_0^t{g}_m^2\left(s,a-t+s\right)\text{d}s},\hfill & \text{if}a\ge t.\hfill \end{array}$ (13)

Now, define the ball

$B_X\left(0,R\right)=\left\{\left(c^y,c^z\right)\in X:{\Vert \left(c^y,c^z\right)\Vert }_{X,\lambda }\le R\right\}.$ (14)

We need the following result which allows the control of $L^{\infty }$ as a function of age.

Lemma 3 ($L^{\infty }$ control in age). Let $\left(c^y,c^z\right)\in B_X\left(0,R\right)$ with $y_0,z_0\in L^{\infty }\left(\left(0,A\right)\times \Omega \right)$ nonnegative. Then there exists a constant $C>0$ , depending on $R$ , $T$ , $A$ , ${\Vert {\beta }_1\Vert }_{L^{\infty }}$ , ${\Vert {\beta }_2\Vert }_{L^{\infty }}$ , ${\Vert y_0\Vert }_{L^{\infty }}$ and ${\Vert z_0\Vert }_{L^{\infty }}$ , such that for all $t\in \left(0,T\right)$ and all $a\in \left(0,A\right)$ ,

${\left|c^y\left(t,a\right)\right|}_n+{\left|c^z\left(t,a\right)\right|}_n\le C.$ (15)

Proof. Let $t\in \left(0,T\right)$ be fixed. By the characteristic representation (see [4]), we have

$c_m^y\left(t,a\right)=\{\begin{array}{ll}{\langle y_0\left(a-t,\cdot \right),e_m\rangle }_{L^2\left(\Omega \right)},\hfill & a>t,\hfill \\ {\displaystyle {\int }_0^A{\beta }_1\left(t-a+s,s\right)c_m^y\left(t-a+s,s\right)\text{d}s},\hfill & a\le t.\hfill \end{array}$ (16)

For $a>t$ , since $y_0\in L^{\infty }\left(\left(0,A\right)\times \Omega \right)$ , we have

${\left|c^y\left(t,a\right)\right|}_n^2={\displaystyle \sum }_{m=1}^n{\left|{\langle y_0\left(a-t,\cdot \right),e_m\rangle }_{L^2\left(\Omega \right)}\right|}^2\le {\Vert y_0\left(a-t,\cdot \right)\Vert }_{L^2\left(\Omega \right)}^2\le \left|\Omega \right|{\Vert y_0\Vert }_{L^{\infty }\left(\left(0,A\right)\times \Omega \right)}^2.$ (17)

For $a\le t$ , using Cauchy-Schwarz inequality:

$\begin{array}{c}{\left|c^y\left(t,a\right)\right|}_n^2={\displaystyle \sum }_{m=1}^n{\left|{\displaystyle {\int }_0^A{\beta }_1\left(t-a+s,s\right)c_m^y\left(t-a+s,s\right)\text{d}s}\right|}^2\\ \le {\Vert {\beta }_1\Vert }_{L^{\infty }\left(\left(0,T\right)\times \left(0,A\right)\right)}^{2}A{\displaystyle {\int }_0^A{\left|c^y\left(t-a+s,s\right)\right|}_n^2\text{d}s}.\end{array}$ (18)

By integrating over $a\in \left(0,t\right)$ and using Gronwall’s inequality (see [23]), we obtain

$\underset{a\in \left(0,A\right)}{\sup }{\left|c^y\left(t,a\right)\right|}_n^2\le C\left(R\right)\left({\Vert y_0\Vert }_{L^{\infty }\left(\left(0,A\right)\times \Omega \right)}^2+1\right),$ (19)

where $C\left(R\right)$ depends on $R$ , $T$ , $A$ and ${\Vert {\beta }_1\Vert }_{L^{\infty }}$ .

The same argument applies to $c^z$ , yielding

$\underset{a\in \left(0,A\right)}{\sup }{\left|c^z\left(t,a\right)\right|}_n^2\le C\left(R\right)\left({\Vert z_0\Vert }_{L^{\infty }\left(\left(0,A\right)\times \Omega \right)}^2+1\right).$ (20)

Combining (19) and (20) gives the desired result (15). □

Step 2: Fixed point formulation.

Define the space

$X=L^2{\left(0,T;L^2\left(0,A\right)\right)}^n\times L^2{\left(0,T;L^2\left(0,A\right)\right)}^n,X_1=L^2{\left(0,T\right)}^n\times L^2{\left(0,T\right)}^n,$

with weighted norm for $\left(c^y,c^z\right)\in X_1$

${\Vert \left(c^y,c^z\right)\Vert }_{X,\lambda }^2={\displaystyle \sum }_{m=1}^n{\displaystyle {\int }_0^T{\text{e}}^{-\lambda t}\left({\Vert c_m^y\left(t\right)\Vert }_{L^2\left(0,A\right)}^2+{\Vert c_m^z\left(t\right)\Vert }_{L^2\left(0,A\right)}^2\right)\text{d}t},$

and for $\left(b^1,b^2\right)\in X_1$

${\Vert \left(b^1,b^2\right)\Vert }_{X_1,\lambda }^2={\displaystyle \sum }_{m=1}^n{\displaystyle {\int }_0^T{\text{e}}^{-\lambda t}\left({\left|b_m^1\left(t\right)\right|}^2+{\left|b_m^2\left(t\right)\right|}^2\right)\text{d}t}.$

$ℱ:X\to X,\left(c^y,c^z\right)\mapsto ℱ\left(c^y,c^z\right)=\mathcal{S}\left(b\left[c^y,c^z\right],g\left[c^y,c^z\right]\right),$

where $\mathcal{S}$ is the solution of the continuous linear operator of the transport Equation (10) with

$\{\begin{array}{l}b\left[c^y,c^z\right]=\left({\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)c^y\left(t,a\right)\text{d}a},{\displaystyle {\int }_0^A{\beta }_2\left(t,a\right)c^z\left(t,a\right)\text{d}a}\right),\\ g\left[c^y,c^z\right]=\left(g_m^1\left[c^y,c^z\right],g_m^2\left[c^y,c^z\right]\right).\end{array}$

Showing the existence of a solution of the system 10 associated with conditions (8) and (9) amounts to finding $\left(c^y,c^z\right)\in X$ such that $\left(c^y,c^z\right)=ℱ\left(c^y,c^z\right)$ and this pair is a fixed point of $ℱ$ .

Preservation of the ball. Let $\left(c^y,c^z\right)\in B_X\left(0,R\right)$ . Using the characteristic representation, we obtain

${\Vert ℱ\left(c^y,c^z\right)\Vert }_{X,\lambda }^2\le C\left({\Vert c^y\Vert }_{X,\lambda }^2+{\Vert c^z\Vert }_{X,\lambda }^2+{\Vert y_0\Vert }_{L^{\infty }\left(Q_A\right)}^2+{\Vert z_0\Vert }_{L^{\infty }\left(Q_A\right)}^2\right).$ (21)

Thus, choosing $R$ sufficiently large such that $R^2\ge 2C\left(R^2+{\Vert y_0\Vert }_{L^{\infty }\left(Q_A\right)}^2+{\Vert z_0\Vert }_{L^{\infty }\left(Q_A\right)}^2\right)$ , we have $ℱ\left(B_X\left(0,R\right)\right)\subset B_X\left(0,R\right)$ .

For simplicity, let us denote by $\delta c^y=c^y-{\tilde{c}}^y$ and $\delta c^z=c^z-{\tilde{c}}^z$ . By linearity of $\mathcal{S}$ , we have

$ℱ\left(c^y,c^z\right)-ℱ\left({\tilde{c}}^y,{\tilde{c}}^z\right)=\mathcal{S}\left(\delta b,\delta g\right),$

with $\delta b=b\left[c^y,c^z\right]-b\left[{\tilde{c}}^y,{\tilde{c}}^z\right]$ , and $\delta g=g\left[c^y,c^z\right]-g\left[{\tilde{c}}^y,{\tilde{c}}^z\right]$ .

By the continuity of $\mathcal{S}$ , there exists $C_0>0$ such that

${\Vert \mathcal{S}\left(\delta b,\delta g\right)\Vert }_{X,\lambda }\le \frac{C_0}{\sqrt{\lambda }}\left({\Vert \delta b\Vert }_{X_1,\lambda }+{\Vert \delta g\Vert }_{X_2,\lambda }\right).$ (22)

Estimation of ${\Vert \delta b\Vert }_{X_1,\lambda }$ . By Cauchy-Schwarz inequality, from

$\delta b_m^1\left(t\right)={\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)\delta c_m^y\left(t,a\right)\text{d}a}.$

Hence,

${\left|\delta b_m^1\left(t\right)\right|}^2\le {\Vert {\beta }_1\Vert }_{L^{\infty }}^{2}A{\Vert \delta c_m^y\left(t\right)\Vert }_{L^2\left(0,A\right)}^2.$ (23)

Multiplying (23) by ${\text{e}}^{-\lambda t}$ and by integrating over $\left[0,T\right]$ we find

${\Vert \delta b_m^1\Vert }_{\lambda }^2\le C_{{\beta }_1}{\displaystyle {\int }_0^T{\text{e}}^{-\lambda t}{\Vert \delta c_m^y\left(t\right)\Vert }_{L^2\left(0,A\right))}^2\text{d}t},\text{with}{C}_{{\beta }_1}={\Vert {\beta }_1\Vert }_{L^{\infty }}^{2}A.$ (24)

Likewise for $\delta b_m^2$ , we find

${\Vert \delta b_m^2\left(t\right)\Vert }_{\lambda }^2\le C_{{\beta }_2}{\displaystyle {\int }_0^T{\text{e}}^{-\lambda t}{\Vert \delta c_m^z\left(t\right)\Vert }_{L^2\left(0,A\right))}^2\text{d}t},\text{with}{C}_{{\beta }_2}={\Vert {\beta }_2\Vert }_{L^{\infty }}^{2}A.$ (25)

Summing (24) and (25) over $m$ , we get

${\Vert \delta b\Vert }_{X_1,\lambda }^2\le C_{\beta }{\Vert \left(\delta c^y,\delta c^z\right)\Vert }_{X,\lambda }^2,\text{with}{C}_{\beta }=\text{max}\left(C_{{\beta }_1},C_{{\beta }_2}\right).$ (26)

Estimation of ${\Vert \delta g\Vert }_{X,\lambda }$ .

Upper bound of $\delta g_1^m$ . Let us write

$\delta g_1^m=-{\displaystyle {\int }_0^L\left[f_{\epsilon }\left(z_n\right){\partial }_x{y}_n-f_{\epsilon }\left({\tilde{z}}_n\right){\partial }_x{\tilde{y}}_n\right]{\partial }_x{e}_m\text{d}x}-{\displaystyle {\int }_0^L{\mu }_1\delta y_n{e}_m\text{d}x},$

and

$f_{\epsilon }\left(z_n\right){\partial }_x{y}_n-f_{\epsilon }\left({\tilde{z}}_n\right){\partial }_x{\tilde{y}}_n=\left[f_{\epsilon }\left(z_n\right)-f_{\epsilon }\left({\tilde{z}}_n\right)\right]{\partial }_x{y}_n+f_{\epsilon }\left({\tilde{z}}_n\right){\partial }_x\left(\delta y_n\right).$

Thus,

$\begin{array}{c}\left|\delta g_1^m\left(t,a\right)\right|\le {\displaystyle {\int }_0^L\left|f_{\epsilon }\left(z_n\right)-f_{\epsilon }\left({\tilde{z}}_n\right)\right|\left|\left|{\partial }_x{y}_n\right|\right|{\partial }_x{e}_m\text{d}x}\\ +{\displaystyle {\int }_0^L\left|f_{\epsilon }\left({\tilde{z}}_n\right)\right|\left|{\partial }_x\left(\delta y_n\right)\right|\left|{\partial }_x{e}_m\right|\text{d}x}+{\Vert {\mu }_1\Vert }_{L^{\infty }}{\displaystyle {\int }_0^L\left|\delta y_n\right|\left|e_m\right|\text{d}x}.\end{array}$(27)

From (2), $f_{\epsilon }$ is globally Lipschitz and so

$\left|f_{\epsilon }\left(z_n\right)-f_{\epsilon }\left({\tilde{z}}_n\right)\right|\le \left|z_n-{\tilde{z}}_n\right|=\left|\delta z_n\right|.$ (28)

Using the inequality (6), we have

$\left|{\partial }_x{y}_n\right|\le {\displaystyle \sum }_{m=1}^n\left|c_k^y\left(t,a\right)\right|\left|{\partial }_x{e}_m\left(x\right)\right|\le C{\displaystyle \sum }_{m=1}^n\left|c_m^y\left(t,a\right)\right|.$ (29)

From (28) and (29) the first term of (27) can be bounded above as follows

${\displaystyle {\int }_0^L\left|f_{\epsilon }\left(z_n\right)-f_{\epsilon }\left({\tilde{z}}_n\right)\right|\left|{\partial }_x{y}_n\right|\left|{\partial }_x{e}_m\right|\text{d}x}\le C^{2}nL{\left({\displaystyle \sum }_{m=1}^n{\left|\delta c_m^z\left(t,a\right)\right|}^2\right)}^{1/2}{\left({\displaystyle \sum }_{m=1}^n{\left|c_m^y\left(t,a\right)\right|}^2\right)}^{1/2}.$(30)

We also have

$\left|f_{\epsilon }\left({\tilde{z}}_n\right)\right|\le 1+\epsilon \text{and}\left|{\partial }_x\left(\delta y_n\right)\right|\le CL{\displaystyle \sum }_{m=1}^n\left|\delta c_m^y\left(t,a\right)\right|$(31)

From (31) and the second term of (27) can be bounded above as follows

${\displaystyle {\int }_0^L\left|f_{\epsilon }\left({\tilde{z}}_n\right)\right|\left|{\partial }_x\left(\delta y_n\right)\right|\left|{\partial }_x{e}_m\right|\text{d}x}\le C^{2}L\sqrt{n}{\left({\displaystyle \sum }_{m=1}^n{\left|\delta c_m^y\left(t,a\right)\right|}^2\right)}^{1/2}.$(32)

For the last term of (27), with

$\left|\delta y_n\right|\le C{\displaystyle \sum }_{m=1}^n\left|\delta c_m^y\left(t,a\right)\right|,$

we obtain

${\Vert {\mu }_1\Vert }_{L^{\infty }}{\displaystyle {\int }_0^L\left|\delta y_n\right|\left|e_m\right|\text{d}x}\le {\Vert {\mu }_1\Vert }_{L^{\infty }}{C}^2\sqrt{n}\left|\Omega \right|{\left({\displaystyle \sum }_{m=1}^n{\left|\delta c_m^y\left(t,a\right)\right|}^2\right)}^{1/2}.$ (33)

Notice that

${\Vert c^y\Vert }_X^2={\displaystyle \sum }_{m=1}^n{\left|c_m^y\right|}_{L^2\left(0,T,L^2\left(0,A\right)\right)}^2\text{and}{\Vert c^z\Vert }_X^2={\displaystyle \sum }_{m=1}^n{\left|c_m^z\right|}_{L^2\left(0,T,L^2\left(0,A\right)\right)}^2,$

${\left|c^y\left(t,a\right)\right|}_n^2={\displaystyle \sum }_{m=1}^n{\left|c_m^y\left(t,a\right)\right|}^2\text{and}{\left|c^z\left(t,a\right)\right|}_n^2={\displaystyle \sum }_{m=1}^n{\left|c_m^z\left(t,a\right)\right|}^2.$

By summing inequalities (30), (32) and (33), inequality (27) gives

$\left|\delta g_1^m\left(t,a\right)\right|\le C_1{\left|c^y\left(t,a\right)\right|}_n{\left|\delta c^z\left(t,a\right)\right|}_n+C_2{\left|\delta c^y\left(t,a\right)\right|}_n,$ (34)

with $C_1,C_2>0$ depending only on $M,n,\left|\Omega \right|,{\Vert {\mu }_1\Vert }_{L^{\infty }}$ and $\epsilon$ .

Integrating (34) over $a\in \left[0,A\right]$ and using ${\left(y+z\right)}^2\le 2\left(y^2+z^2\right)$ , we find

${\Vert \delta g_1^m\left(t\right)\Vert }_{L^2\left(0,A\right)}^2\le C_3{\displaystyle {\int }_0^A\left({\Vert c^y\left(t,a\right)\Vert }_n^2{\Vert \delta c^z\left(t,a\right)\Vert }_n^2+{\Vert \delta c^y\left(t,a\right)\Vert }_n^2\right)\text{d}a},$ (35)

with $C_3=2\max \left(C_1^2,C_2^2\right)$ .

By Hölder’s inequality and Lemma 3, there exist $R_1>0$ such that ${\Vert c^y\left(t\right)\Vert }_{L^{\infty }\left(0,A\right)}\le R_1$ for all $t\in \left[0,T\right]$ , (35) becomes

${\Vert \delta g_1^m\left(t\right)\Vert }_{L^2\left(0,A\right)}^2\le C_1\left(R\right)\left({\Vert \delta c^z\left(t\right)\Vert }_{L^2\left(0,A\right)}^2+{\Vert \delta c^y\left(t\right)\Vert }_{L^2\left(0,A\right)}^2\right),$ (36)

with $C_1\left(R\right)=C_3\left(2n{R}_1^2+1\right)$ .

Similarly, there exists a positive constant $C_2\left(R\right)$ such that

${\Vert \delta g_2^m\left(t\right)\Vert }_{L^2\left(0,A\right)}^2\le C_2\left(R\right)\left({\Vert \delta c^z\left(t\right)\Vert }_{L^2\left(0,A\right)}^2+{\Vert \delta c^y\left(t\right)\Vert }_{L^2\left(0,A\right)}^2\right).$ (37)

With the weighted norm and summing over $m$ in (36) and (37), we obtain

${\Vert \delta g\Vert }_{X,\lambda }^2\le \text{max}\left(C_1\left(R\right),C_1\left(R\right)\right){\Vert \left(\delta c^y,\delta c^z\right)\Vert }_{X,\lambda }^2.$ (38)

Finally from (26) and (38), (22) leads to

${\Vert ℱ\left(c^y,c^z\right)-ℱ\left({\tilde{c}}^y,{\tilde{c}}^z\right)\Vert }_{X,\lambda }^2\le \frac{2C_0^2\max \left(C_{\beta },C\left(R\right)\right)}{\lambda }{\Vert \left(c^y,c^z\right)-\left({\tilde{c}}^y,{\tilde{c}}^z\right)\Vert }_{X,\lambda }^2,$ (39)

with, $C\left(R\right)=\text{max}\left(C_1\left(R\right),C_1\left(R\right)\right)$ .

Now, choosing $\lambda >2C_0^2\max \left(C_{\beta },C\left(R\right)\right)$ , we show that $ℱ$ is a contraction on $B_X\left(0,R\right)$ for all $T>0$ . By Banach’s fixed point theorem, there exists a unique fixed point $\left(c^y,c^z\right)\in B_X\left(0,R\right)$ , which gives the existence of the Galerkin solution.

Step 3: A priori estimates independent of $n$ .

The global existence of the Faedo-Galerkin approximation solution is obtained by deriving independent a priori estimates of $n$ for the sequence of functions ${\left(y_n\right)}_n$ et ${\left(z_n\right)}_n$ in various Banach spaces. Thus, we have the following lemma.

Lemma 4. Assume that assumptions (1) hold. If $y_0,z_0\in L^{\infty }\left(\left(0,A\right)\times \Omega \right)$ , then there exist constants $c_1$ , $c_2$ , and $c_3$ independent of $n$ such that

${\Vert y_n^{\epsilon }\Vert }_{\mathscr{H}}+{\Vert z_n^{\epsilon }\Vert }_{\mathscr{H}}\le c_1,$ (40)

${\Vert {\partial }_x{y}_n^{\epsilon }\Vert }_{{\mathscr{H}}_1}+{\Vert {\partial }_x{z}_n^{\epsilon }\Vert }_{{\mathscr{H}}_1}\le c_2,$ (41)

${\Vert \left({\partial }_t+{\partial }_a\right)y_n^{\epsilon }\Vert }_{{\mathscr{H}}_2}+{\Vert \left({\partial }_t+{\partial }_a\right)z_n^{\epsilon }\Vert }_{{\mathscr{H}}_2}\le c_3,$ (42)

with $\mathscr{H}=L^{\infty }\left(\left(0,T\right)\times \left(0,A\right);L^2\left(0,L\right)\right)$ , ${\mathscr{H}}_1=L^2\left(\left(0,T\right)\times \left(0,A\right);L^2\left(0,L\right)\right)$ , ${\mathscr{H}}_2=L^{\infty }\left(\left(0,T\right)\times \left(0,A\right);{\left(H^1\left(0,L\right)\right)}^{\text{*}}\right)$ .

Proof. Proof of Lemma 4 Consider the following weak system

$\{\begin{array}{l}{\displaystyle {\int }_0^L\left[\left({\partial }_t+{\partial }_a\right)y_n^{\epsilon }\cdot e_m+f_{\epsilon }\left(z_n^{\epsilon }\right){\partial }_x{y}_n^{\epsilon }{\partial }_x{e}_m+{\mu }_1{y}_n^{\epsilon }{e}_m\right]\text{d}x}=0\\ {\displaystyle {\int }_0^L\left[\left({\partial }_t+{\partial }_a\right)z_n^{\epsilon }\cdot e_m+f_{\epsilon }\left(y_n^{\epsilon }\right){\partial }_x{z}_n^{\epsilon }{\partial }_x{e}_m+{\mu }_2{z}_n^{\epsilon }{e}_m\right]\text{d}x}=0.\end{array}$ (43)

Multiplying the first and second equations of the system 43 by $y_n^{\epsilon }$ and $z_n^{\epsilon }$ respectively, then the Faedo-Galerkin solutions satisfy for each fixed $t\ge 0$ and $a\ge 0$ , the following weak formulations:

$\{\begin{array}{l}{\displaystyle {\int }_0^L\left(\frac{\text{d}}{2\text{d}t}{\left(y_n^{\epsilon }\right)}^2+\frac{\text{d}}{2\text{d}a}{\left(y_n^{\epsilon }\right)}^2+f_{\epsilon }\left(z_n^{\epsilon }\right){\left|{\partial }_x{y}_n^{\epsilon }\right|}^2\right)\text{d}x}=-{\displaystyle {\int }_0^L{\mu }_1{\left(y_n^{\epsilon }\right)}^2\text{d}x},\\ {\displaystyle {\int }_0^L\left(\frac{\text{d}}{2\text{d}t}{\left(z_n^{\epsilon }\right)}^2+\frac{\text{d}}{2\text{d}a}{\left(z_n^{\epsilon }\right)}^2+f_{\epsilon }\left(y_n^{\epsilon }\right){\left|{\partial }_x{z}_n^{\epsilon }\right|}^2\right)\text{d}x}=-{\displaystyle {\int }_0^L{\mu }_2{\left(z_n^{\epsilon }\right)}^2\text{d}x}.\end{array}$ (44)

Integrating over $\left[0,T\right]\times \left[0,A\right]$ and summing the two equations of (44), we obtain

$\begin{array}{l}A\left(t\right)+2{\displaystyle {\int }_Q\left[f_{\epsilon }\left(z_n^{\epsilon }\right){\left|{\partial }_x{y}_n^{\epsilon }\right|}^2+f_{\epsilon }\left(y_n^{\epsilon }\right){\left|{\partial }_x{z}_n^{\epsilon }\right|}^2+{\mu }_1{\left(y_n^{\epsilon }\right)}^2+{\mu }_2{\left(z_n^{\epsilon }\right)}^2\right]\text{d}x\text{d}\tau \text{d}a}\\ =A\left(0\right)+{\displaystyle {\int }_{Q_T}\left[{\left(y_n^{\epsilon }\right)}^2+{\left(z_n^{\epsilon }\right)}^2\right]\left(\tau ,0,x\right)\text{d}x\text{d}\tau },\end{array}$ (45)

with $A\left(t\right)={\displaystyle {\int }_{Q_A}\left[{\left(y_n^{\epsilon }\right)}^2+{\left(z_n^{\epsilon }\right)}^2\right]\left(t,s,x\right)\text{d}x\text{d}s}$ .

By the boundary condition and the Cauchy-Schwarz inequality, we have

$\begin{array}{l}{\displaystyle {\int }_0^L{\left|y_n^{\epsilon }\left(t,0,x\right)\right|}^2\text{d}x}+{\displaystyle {\int }_0^L{\left|z_n^{\epsilon }\left(t,0,x\right)\right|}^2\text{d}x}\\ \le C{\displaystyle {\int }_{Q_A}{\left|y_n^{\epsilon }\left(t,a,x\right)\right|}^2\text{d}a\text{d}x}+C{\displaystyle {\int }_{Q_A}{\left|z_n^{\epsilon }\left(t,a,x\right)\right|}^2\text{d}a\text{d}x}\end{array}$ (46)

and

${\displaystyle {\int }_Q\left({\mu }_1{\left(y_n^{\epsilon }\right)}^2+{\mu }_2{\left(z_n^{\epsilon }\right)}^2\right)}\le C_1{\displaystyle {\int }_Q\left({\left|y_n^{\epsilon }\left(\tau ,a,x\right)\right|}^2+{\left|z_n^{\epsilon }\left(\tau ,a,x\right)\right|}^2\right)},$ (47)

with $C=A{\max }_{i\in \left\{1,2\right\}}\left({\Vert {\beta }_i\Vert }_{\infty }^2\right)$ and $C_1={\text{max}}_{i\in \left\{1,2\right\}}\left({\Vert {\mu }_i\Vert }_{\infty }\right)$ .

Moreover, we have

$\begin{array}{c}A\left(0\right)={\displaystyle {\int }_{Q_A}\left({\left(y_n^{\epsilon }\right)}^2+{\left(z_n^{\epsilon }\right)}^2\right)\left(0,s,x\right)\text{d}x\text{d}s}\\ \le AL\left({\Vert y_0\Vert }_{L^{\infty }\left(\left(0,A\right)\times \Omega \right)}^2+{\Vert z_0\Vert }_{L^2\left(\left(0,A\right)\times \Omega \right)}^2\right)=C_2.\end{array}$

Using (46), inequality (45) becomes

$A\left(t\right)\le A\left(0\right)+C{\displaystyle {\int }_0^{t}A\left(\tau \right)\text{d}\tau }.$ (48)

From (48), We apply Gronwall’s lemma to obtain

$A\left(t\right)\le A\left(0\right){\text{e}}^{Ct}\Rightarrow {\displaystyle {\int }_0^{T}A\left(t\right)\text{d}t}\le \frac{A\left(0\right)}{C}\left({\text{e}}^{CT}-1\right).$ (49)

Hence the result (40) with $c_1=\frac{C_2}{C}\left({\text{e}}^{CT}-1\right)$ .

From (45), (2) and (49), we have

$2\epsilon {\displaystyle {\int }_Q\left({\left|{\partial }_x{y}_n^{\epsilon }\right|}^2+{\left|{\partial }_x{z}_n^{\epsilon }\right|}^2\right)}\le {c^{\prime }}_2.$ (50)

Hence the result (41) with $c_2=\frac{{c^{\prime }}_2}{2\epsilon }$ .

Let $\varphi \in H^1\left(0,L\right)$ . From (44), we find

$\left|{\displaystyle {\int }_0^L\left({\partial }_t+{\partial }_a\right)y_n^{\epsilon }\cdot \varphi \text{d}x}\right|\le {\displaystyle {\int }_0^L\left(\left|f_{\epsilon }\left(z_n^{\epsilon }\right)\right|\left|{\partial }_x{y}_n^{\epsilon }\right|\left|{\partial }_x\varphi \right|+{\Vert {\mu }_1\Vert }_{\infty }\left|y_n^{\epsilon }\right|\Vert \varphi \Vert \right)\text{d}x}$

From (2), we obtain

$\begin{array}{c}\left|{\displaystyle {\int }_0^L\left({\partial }_t+{\partial }_a\right)y_n^{\epsilon }\cdot {\varphi }_1\text{d}x}\right|\le \left(1+\epsilon \right){\Vert {\partial }_x{y}_n\Vert }_{L^2\left(0,L\right)}{\Vert {\partial }_x\varphi \Vert }_{L^2\left(0,L\right)}\\ +{\Vert {\mu }_1\Vert }_{L^{\infty }}{\Vert y_n\Vert }_{L^2\left(0,L\right)}{\Vert \varphi \Vert }_{L^2\left(0,L\right)}\end{array}$

Using ${\Vert \varphi \Vert }_{L^2\left(0,L\right)}\le {\Vert \varphi \Vert }_{H^1\left(0,L\right)}$ and ${\Vert {\partial }_x\varphi \Vert }_{L^2\left(0,L\right)}\le {\Vert \varphi \Vert }_{H^1\left(0,L\right)}$ , we have

${\Vert \left({\partial }_t+{\partial }_a\right)y_n\Vert }_{{\left(H^1\left(0,L\right)\right)}^{\text{*}}}\le \left(1+\epsilon \right){\Vert {\partial }_x{y}_n\Vert }_{L^2\left(0,L\right)}+{\Vert {\mu }_1\Vert }_{L^{\infty }}{\Vert y_n\Vert }_{L^2\left(0,L\right)}.$

Integrating at $t,a$ with the $L^{\infty }$ norm, and thanks to the inequalities (40) and (41), we obtain (42). This completes the proof of Lemma 4. □

Step 4: Passage to the limit $n\to \infty$ .

By Lemma 4, the sequences ${\left(y_n^{\epsilon }\right)}_{n\ge 1}$ and ${\left(z_n^{\epsilon }\right)}_{n\ge 1}$ are uniformly bounded in $L^{\infty }\left(\left(0,T\right)\times \left(0,A\right);L^2\left(\Omega \right)\right)\cap L^2\left(\left(0,T\right)\times \left(0,A\right);H^1\left(\Omega \right)\right)$ and the sequence ${\left(\left({\partial }_t+{\partial }_a\right)y_n^{\epsilon }\right)}_{n\ge 1}$ is uniformly bounded in $L^2\left(\left(0,T\right)\times \left(0,A\right);{\left(H^1\left(\Omega \right)\right)}^{*}\right)$ independently of $n$ . By the Aubin-Lions lemma, we extract subsequences such that

$\left(y_n^{\epsilon },z_n^{\epsilon }\right)\to \left(y^{\epsilon },z^{\epsilon }\right)\text{strongly}\text{in}{\left[L^2\left(Q\right)\right]}^2.$

Moreover, from (2), $f_{\epsilon }$ is continuous and adding (40), we have following convergences

$\{\begin{array}{l}{f}_{\epsilon }\left(z_n^{\epsilon }\right)\to f_{\epsilon }\left(z^{\epsilon }\right)\text{strongly}\text{in}{L}^2\left(\left(0,T\right)\times \left(0,A\right);L^2\left(\Omega \right)\right),\\ {\partial }_x{y}_n^{\epsilon }\rightharpoonup {\partial }_x{y}^{\epsilon }\text{weakly}\text{in}{L}^2\left(\left(0,T\right)\times \left(0,A\right);L^2\left(\Omega \right)\right),\\ \left({\partial }_t+{\partial }_a\right)y_n^{\epsilon }\rightharpoonup \left({\partial }_t+{\partial }_a\right)y^{\epsilon }\text{weakly}\text{in}{L}^2\left(\left(0,T\right)\times \left(0,A\right);{\left(H^1\left(\Omega \right)\right)}^{*}\right).\end{array}$

Thus, $f_{\epsilon }\left(z_n^{\epsilon }\right){\partial }_x{y}_n^{\epsilon }$ converges weakly to $f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }$ in $L^2\left(\left(0,T\right)\times \left(0,A\right);L^2\left(\Omega \right)\right)$ .

Now, using test functions ${\varphi }_1\in {\mathcal{C}}_c^{\infty }\left(\left[0,T\right)\times \left[0,A\right)\times \Omega \right)$ and by passing to the limit in the following weak formulation

${\displaystyle {\int }_Q\left(\left({\partial }_t+{\partial }_a\right)y_n^{\epsilon }{\varphi }_1+f_{\epsilon }\left(z_n^{\epsilon }\right){\partial }_x{y}_n^{\epsilon }{\partial }_x{\varphi }_1+{\mu }_1{y}_n^{\epsilon }{\varphi }_1\right)\text{d}t\text{d}a\text{d}x}=0,$ (51)

we find the following convergences

$y_n^{\epsilon }{\varphi }_1\to y^{\epsilon }{\varphi }_1,\text{strongly}\text{in}{L}^2\left(\left(0,T\right)\times \left(0,A\right);L^2\left(\Omega \right)\right),$

$f_{\epsilon }\left(z_n^{\epsilon }\right){\partial }_x{y}_n^{\epsilon }{\partial }_x{\varphi }_1\rightharpoonup f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }{\partial }_x{\varphi }_1,\text{weakly}\text{in}{L}^2\left(\left(0,T\right)\times \left(0,A\right);L^2\left(\Omega \right)\right),$

$\left({\partial }_t+{\partial }_a\right)y_n^{\epsilon }{\varphi }_1\rightharpoonup \left({\partial }_t+{\partial }_a\right)y^{\epsilon }{\varphi }_1,\text{weakly}\text{in}{L}^2\left(\left(0,T\right)\times \left(0,A\right);{\left(H^1\left(\Omega \right)\right)}^{\text{*}}\right).$

The same procedure is used for the equation concerning $z_n$ .

This shows that for all ${\varphi }_1,{\varphi }_2\in {\mathcal{C}}_c^{\infty }\left(\left[0,T\right)\times \left[0,A\right)\times \Omega \right)$ , $\left(y,z\right)$ is a weak solution to the following system,

$\{\begin{array}{l}{\displaystyle {\int }_Q\left(\left[\left({\partial }_t+{\partial }_a\right)y^{\epsilon }\right]{\varphi }_1+f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }{\partial }_x{\varphi }_1+{\mu }_1{y}^{\epsilon }{\varphi }_1\right)\text{d}t\text{d}a\text{d}x}=0,\\ {\displaystyle {\int }_Q\left(\left[\left({\partial }_t+{\partial }_a\right)z^{\epsilon }\right]{\varphi }_2+f_{\epsilon }\left(y^{\epsilon }\right){\partial }_x{z}^{\epsilon }{\partial }_x{\varphi }_2+{\mu }_2{z}^{\epsilon }{\varphi }_2\right)\text{d}t\text{d}a\text{d}x}=0.\end{array}$ (52)

Now we show that the limit functions $\left(y,z\right)$ satisfy the initial conditions and birth conditions. We will adapt a standard argument given in [23].

Let $N\ge 1$ be a fixed integer. We consider as test functions, ${\varphi }_1,{\varphi }_2$ of the form $\forall \left(t,a,x\right)\in Q$ , ${\varphi }_i\left(t,a,x\right)={\displaystyle {\sum }_{k=1}^N{c}_k^i\left(t,a\right)e_k\left(x\right)}$ , where $c_k^i\in {\mathcal{C}}^1\left(\left[0,T\right]\times \left[0,A\right]\right)$ such that ${\varphi }_i\left(T,\cdot ,\cdot \right)={\varphi }_i\left(\cdot ,A,\cdot \right)=0$ .

By integrating (51) and the first equation of the system (52) by part, we obtain the following equations respectively

$\begin{array}{l}-{\displaystyle {\int }_Q\left(\left[\left({\partial }_t+{\partial }_a\right){\varphi }_1\right]y_n^{\epsilon }+f_{\epsilon }\left(z_n^{\epsilon }\right){\partial }_x{y}_n^{\epsilon }{\partial }_x{\varphi }_1+{\mu }_1{y}_n^{\epsilon }{\varphi }_1\right)\text{d}t\text{d}a\text{d}x}\\ ={\displaystyle {\int }_{Q_A}{y}_n^{\epsilon }\left(0,a,x\right){\varphi }_1\left(0,a,x\right)\text{d}a\text{d}x}+{\displaystyle {\int }_{Q_T}{y}_n^{\epsilon }\left(t,0,x\right){\varphi }_1\left(t,0,x\right)\text{d}t\text{d}x},\end{array}$ (53)

and

$\begin{array}{l}-{\displaystyle {\int }_Q\left(\left[\left({\partial }_t+{\partial }_a\right){\varphi }_1\right]y^{\epsilon }+f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }{\partial }_x{\varphi }_1+{\mu }_1{y}^{\epsilon }{\varphi }_1\right)\text{d}t\text{d}a\text{d}x}\\ ={\displaystyle {\int }_{Q_A}{y}_0\left(a,x\right){\varphi }_1\left(0,a,x\right)\text{d}a\text{d}x}+{\displaystyle {\int }_{Q_T}{y}^{\epsilon }\left(t,0,x\right){\varphi }_1\left(t,0,x\right)\text{d}t\text{d}x}.\end{array}$ (54)

By passing to the limit when $n$ tends to infinity ($n\to +\infty$ ) in (53), we obtain for all ${\varphi }_1\in C_c^{\infty }\left(\left[0,T\right)\times \left[0,A\right)\times \Omega \right)$

$\begin{array}{l}-{\displaystyle {\int }_Q\left(\left[\left({\partial }_t+{\partial }_a\right){\varphi }_1\right]y^{\epsilon }+f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }{\partial }_x{\varphi }_1+{\mu }_1{y}^{\epsilon }{\varphi }_1\right)\text{d}t\text{d}a\text{d}x}\\ ={\displaystyle {\int }_{Q_A}{y}_0\left(a,x\right){\varphi }_1\left(0,a,x\right)\text{d}a\text{d}x}+{\displaystyle {\int }_{Q_T}\left({\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)y^{\epsilon }\left(t,a,x\right)\text{d}a}\right){\varphi }_1\left(t,0,x\right)\text{d}t\text{d}x},\end{array}$ (55)

Comparing (54) and (55), we have this leads to

$y^{\epsilon }\left(t,0,x\right)=y_0\left(a,x\right)\text{and}{y}^{\epsilon }\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)y^{\epsilon }\left(t,a,x\right)\text{d}a}.$

The same procedure is used for the equation concerning $z_n$ , to have

$z^{\epsilon }\left(t,0,x\right)=z_0\left(a,x\right)\text{and}{z}^{\epsilon }\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)z^{\epsilon }\left(t,a,x\right)\text{d}a}.$

Moreover, since $f_{\epsilon }\left(z_n^{\epsilon }\right)\ge \epsilon >0$ and $f_{\epsilon }\left(y^{\epsilon }\right)\ge \epsilon >0$ , these conditions are well-defined and imply ${\partial }_x{y}_n^{\epsilon }={\partial }_x{z}_n^{\epsilon }=0$ at the boundary.

As $n\to +\infty$ , using the convergences $f_{\epsilon }\left(z_n^{\epsilon }\right)\to f_{\epsilon }\left(z^{\epsilon }\right)$ , $f_{\epsilon }\left(y_n^{\epsilon }\right)\to f_{\epsilon }\left(y^{\epsilon }\right)$ , ${\partial }_x{y}_n^{\epsilon }\rightharpoonup {\partial }_x{y}^{\epsilon }$ , and ${\partial }_x{z}_n^{\epsilon }\rightharpoonup {\partial }_x{z}^{\epsilon }$ , we recover the original degenerate conditions:

$f_{\epsilon }\left(z^{\epsilon }\right){\partial }_{x}y=0,f_{\epsilon }\left(y^{\epsilon }\right){\partial }_{x}z=0\text{on}\Sigma ,$ (56)

in the weak sense Which completes the proof of Theorem 2. □

3. Uniform Estimates with Respect to $\epsilon$ and Positivity

Lemma 5. Assume that assumptions (1) hold. If $y_0,z_0\in L^{\infty }\left(\left(0,A\right)\times \Omega \right)$ are positives, then the solution $\left(y^{\epsilon },z^{\epsilon }\right)$ to the system (3) is positive. Moreover, there exist constants $c_4,c_5,c_6$ independent of $\epsilon$ such that

${\Vert y^{\epsilon }\Vert }_{\mathscr{H}}+{\Vert z^{\epsilon }\Vert }_{\mathscr{H}}\le c_4,$ (57)

${\displaystyle {\int }_Q{f}_{\epsilon }\left(z^{\epsilon }\right){\left({\partial }_x{y}^{\epsilon }\right)}^2\text{d}x\text{d}a\text{d}t}+{\displaystyle {\int }_Q{f}_{\epsilon }\left(y^{\epsilon }\right){\left({\partial }_x{z}^{\epsilon }\right)}^2\text{d}x\text{d}a\text{d}t}\le c_5,$ (58)

${\Vert \left({\partial }_t+{\partial }_a\right)y^{\epsilon }\Vert }_{{\mathscr{H}}_3}+{\Vert \left({\partial }_t+{\partial }_a\right)z^{\epsilon }\Vert }_{{\mathscr{H}}_3}\le c_6,$ (59)

with ${\mathscr{H}}_3=L^2\left(\left(0,T\right)\times \left(0,A\right),{\left(H^1\left(0,L\right)\right)}^{\text{*}}\right)$ .

Proof of Lemma 5. Let’s show the positivity of the weak solution.

Let $\left(y^{\epsilon },z^{\epsilon }\right)$ be a weak solution of (4) obtained by Galerkin’s method. We consider the negative part $y^{-}=\text{min}\left(y^{\epsilon },0\right)$ and $z^{-}=\text{min}\left(z^{\epsilon },0\right)$ . We test the first equation of (4) with $\varphi =-y^{-}$ (which belongs to the feasible region). We obtain

$\frac{\text{d}}{2\text{d}t}{\displaystyle {\int }_0^L{\left|y^{-}\right|}^2\text{d}x}+\frac{\text{d}}{2\text{d}a}{\displaystyle {\int }_0^L{\left|y^{-}\right|}^2\text{d}x}=-{\displaystyle {\int }_0^L{f}_{\epsilon }\left(z^{\epsilon }\right){\left|{\partial }_x{y}^{-}\right|}^2\text{d}x}-{\displaystyle {\int }_0^L{\mu }_1{\left|y^{-}\right|}^2\text{d}x}.$ (60)

By a simple calculation we find the following

$y^{\epsilon }{y}^{-}=-{\left|y^{-}\right|}^2,{\partial }_x{y}^{\epsilon }\cdot {\partial }_x{y}^{-}=-{\left|{\partial }_x{y}^{-}\right|}^2\text{and}{f}_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }\cdot {\partial }_x{y}^{-}\le 0.$ (61)

Since the functions ${\mu }_1$ and ${\mu }_2$ are non-negative, by using (61), the inequality (60) becomes

$\frac{\text{d}}{2\text{d}t}{\displaystyle {\int }_0^L{\left|y^{-}\right|}^2\text{d}x}+\frac{\text{d}}{2\text{d}a}{\displaystyle {\int }_0^L{\left|y^{-}\right|}^2\text{d}x}\le 0$ (62)

By integrating (62) over $\left(0,T\right)\times \left(0,A\right)$ , we get

$\begin{array}{l}{\displaystyle {\int }_Q{\left|y^{-}\left(t,a,x\right)\right|}^2\text{d}x\text{d}a}+{\displaystyle {\int }_Q{\left|y^{-}\left(t,a,x\right)\right|}^2\text{d}x\text{d}t}\\ \le {\displaystyle {\int }_Q{\left|y^{-}\left(0,s,x\right)\right|}^2\text{d}x\text{d}a}+{\displaystyle {\int }_Q{\left|y^{-}\left(t,0,x\right)\right|}^2\text{d}x\text{d}t},\end{array}$

Taking into account the positivity of the initial data and the birth data, we have

${\displaystyle {\int }_Q{\left|y^{-}\left(t,s,x\right)\right|}^2\text{d}x\text{d}a}+{\displaystyle {\int }_Q\left|y^{-}\left(t,a,x\right)\right|{|}^2\text{d}x\text{d}t}\le 0$

Hence $y^{-}\equiv 0$ a.e. in $Q$ . The same holds for $z^{\epsilon }$ .

Let’s now show inequalities (57), (59), and (58).

For (57). We multiply the first equation of (4) by $y^{\epsilon }$ and integrate over $x$

$\frac{1}{2}{\partial }_t{\Vert y^{\epsilon }\Vert }_{L_x^2}^2+\frac{1}{2}{\partial }_a{\Vert y^{\epsilon }\Vert }_{L_x^2}^2+{\displaystyle {\int }_0^L{f}_{\epsilon }\left(z^{\epsilon }\right){\left({\partial }_x{y}^{\epsilon }\right)}^2\text{d}x}+{\displaystyle {\int }_0^L{\mu }_1{\left(y^{\epsilon }\right)}^2\text{d}x}=0.$ (63)

Then we integrate over $a\in \left(0,A\right)$ and $t\in \left(0,T\right)$ . Since $f_{\epsilon }\left(z^{\epsilon }\right)\ge \epsilon$ and ${\mu }_1\ge 0$ , we have

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_0^A{\Vert y^{\epsilon }\Vert }_{L^2\left(0,L\right)}^2\text{d}a}+{\left[{\Vert y^{\epsilon }\Vert }_{L^2\left(0,L\right)}^2\right]}_{a=0}^{a=A}+2\epsilon {\displaystyle {\int }_0^A{\Vert {\partial }_x{y}^{\epsilon }\Vert }_{L^2\left(0,L\right)}^2\text{d}a}\le 0.$ (64)

The boundary term at $a=A$ is the maximum age of the two populations, therefore $y^{\epsilon }\left(t,A,x\right)=0$ . The term at $a=0$ is

${\Vert y^{\epsilon }\left(t,0,x\right)\Vert }_{L^2\left(0,L\right)}^2={\displaystyle {\int }_0^L{\left({\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)y^{\epsilon }\left(t,a,x\right)\text{d}a}\right)}^2\text{d}x}.$ (65)

By Cauchy-Schwarz and the assumption ${\beta }_1\in L^{\infty }\left(\left(0,T\right)\times \left(0,A\right)\right)$ , (65) leads to

${\Vert y^{\epsilon }\left(t,0,x\right)\Vert }_{L^2\left(0,L\right)}^2\le A{\Vert {\beta }_1\left(t,\cdot \right)\Vert }_{L^{\infty }\left(0,A\right)}^2{\displaystyle {\int }_0^A{\Vert y^{\epsilon }\left(t,a,\cdot \right)\Vert }_{L^2\left(0,L\right)}^2\text{d}a}.$ (66)

We set $E\left(t\right)={\displaystyle {\int }_0^A{\Vert y^{\epsilon }\left(t,a,\cdot \right)\Vert }_{L^2\left(0,L\right)}^2\text{d}a}$ . Then by the inequality (67), (64) becomes

$\frac{1}{2}\frac{\text{d}E}{\text{d}t}\le CE\left(t\right),\text{with}C=A{\Vert {\beta }_1\Vert }_{L^{\infty }\left(\left(0,T\right)\times \left(0,A\right)\right)}^2.$ (67)

We integrate in time and use Gronwall’s lemma to obtain

$E\left(t\right)\le E\left(0\right){\text{e}}^{\left(2C+1\right)t},$ (68)

with $E\left(0\right)={\displaystyle {\int }_0^A{\Vert y_0\Vert }_{L^2\left(0,L\right)}^2\text{d}a}\le AL{\Vert y_0\Vert }_{L^{\infty }\left(Q_A\right)}^2$ .

Thus, we obtain a uniform bound at $\epsilon$ for ${\Vert y^{\epsilon }\Vert }_{\mathscr{H}}$ . The same holds for $z^{\epsilon }$ . Hence (57).

Estimate (58) comes directly from the energy equality

$\begin{array}{l}{\Vert y^{\epsilon }\left(T\right)\Vert }_{L^2\left(Q_A\right)}^2+{\Vert z^{\epsilon }\left(T\right)\Vert }_{L^2\left(Q_A\right)}^2+2{\displaystyle {\int }_Q\left(f_{\epsilon }\left(z^{\epsilon }\right){\left({\partial }_x{y}^{\epsilon }\right)}^2+f_{\epsilon }\left(y^{\epsilon }\right){\left({\partial }_x{z}^{\epsilon }\right)}^2\right)}\\ ={\Vert y^{\epsilon }\left(\cdot ,0,\cdot \right)\Vert }_{L^2\left(Q_T\right)}^2+{\Vert z^{\epsilon }\left(\cdot ,0,\cdot \right)\Vert }_{L^2\left(Q_T\right)}^2+AL\left({\Vert y_0\Vert }_{L^{\infty }\left(Q_A\right)}^2+{\Vert z_0\Vert }_{L^{\infty }\left(Q_A\right)}^2\right).\end{array}$ (69)

From (68), the inequality (69) leads to

$\begin{array}{l}\frac{1}{2}{\Vert y^{\epsilon }\left(T\right)\Vert }_{L^2\left(Q_A\right)}^2+\frac{1}{2}{\Vert z^{\epsilon }\left(T\right)\Vert }_{L^2\left(Q_A\right)}^2+{\displaystyle {\int }_Q{f}_{\epsilon }\left(z^{\epsilon }\right){\left({\partial }_x{y}^{\epsilon }\right)}^2+f_{\epsilon }\left(y^{\epsilon }\right){\left({\partial }_x{z}^{\epsilon }\right)}^2}\\ \le \frac{AL}{2}\left({\text{e}}^{\left(2C+1\right)T}+1\right)\left({\Vert y_0\Vert }_{L^{\infty }\left(Q_A\right)}^2+{\Vert z_0\Vert }_{L^{\infty }\left(Q_A\right)}^2\right).\end{array}$ (70)

This gives us (58), with $c_5=\frac{AL}{2}\left({\text{e}}^{\left(2C+1\right)T}+1\right)\left({\Vert y_0\Vert }_{L^{\infty }\left(Q_A\right)}^2+{\Vert z_0\Vert }_{L^{\infty }\left(Q_A\right)}^2\right)$ .

To prove the last inequality (59), we use duality. Let $\psi \in H^1\left(0,L\right)$ with ${\Vert \psi \Vert }_{H^1}\le 1$ .

We replace $\varphi$ with $\psi$ in the first equation of (4). Thus, we obtain

$\langle \left({\partial }_t+{\partial }_a\right)y^{\epsilon },\psi \rangle =-{\displaystyle {\int }_{\Omega }{f}_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }{\partial }_x\psi \text{d}x}-{\displaystyle {\int }_{\Omega }{\mu }_1{y}^{\epsilon }\psi \text{d}x}.$ (71)

According to (58) and the Cauchy-Schwarz inequality, (71) gives

$\left|\langle \left({\partial }_t+{\partial }_a\right)y^{\epsilon },\psi \rangle \right|\le \sqrt{2}{\left({\displaystyle {\int }_{\Omega }}{f}_{\epsilon }\left(z^{\epsilon }\right){({\partial }_x{y}^{\epsilon })}^{2}dx\right)}^{1/2}+\parallel {\mu }_1{\parallel }_{L^{\infty }}\parallel y^{\epsilon }{\parallel }_{L^2\left(\text{Ω}\right)}.$ (72)

Taking the supremum over $\psi$ , we obtain

${\Vert \left({\partial }_t+{\partial }_a\right)y^{\epsilon }\left(t,a,\cdot \right)\Vert }_{{\left(H^1\right)}^{\text{*}}}\le \sqrt{2}{\left({\displaystyle {\int }_{\Omega }{f}_{\epsilon }\left(z^{\epsilon }\right){\left({\partial }_x{y}^{\epsilon }\right)}^2\text{d}x}\right)}^{1/2}+{\Vert {\mu }_1\Vert }_{L^{\infty }}{\Vert y^{\epsilon }\Vert }_{L^2\left(\Omega \right)}.$ (73)

From ${\left(a+b\right)}^2\le 2a^2+2b^2$ and integrating over $\left(0,T\right)\times \left(0,A\right)$ , inequality (73) leads to

$\begin{array}{c}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^A{\Vert \left({\partial }_t+{\partial }_a\right)y^{\epsilon }\Vert }_{{\left(H^1\right)}^{*}}^2\text{d}a\text{d}t}}\le 4{\displaystyle {\int }_Q{f}_{\epsilon }\left(z^{\epsilon }\right){\left({\partial }_x{y}^{\epsilon }\right)}^2\text{d}x\text{d}a\text{d}t}\\ +2{\Vert {\mu }_1\Vert }_{L^{\infty }}^2{\displaystyle {\int }_0^T{\displaystyle {\int }_0^A{\Vert y^{\epsilon }\Vert }_{L^2\left(\Omega \right)}^2\text{d}a\text{d}t}}.\end{array}$ (74)

According to (57) and (58), we have

${\displaystyle {\int }_0^T{\displaystyle {\int }_0^A{\Vert \left({\partial }_t+{\partial }_a\right)y^{\epsilon }\Vert }_{{\left(H^1\right)}^{\text{*}}}^2\text{d}a\text{d}t}}\le c_6^2,$ (75)

with $c_6^2=4C+2{\Vert {\mu }_1\Vert }_{L^{\infty }\left(\left(0,T\right)\times \left(0,A\right)\right)}^{2}TA{c}_4^2$ .

This completes the proof of Lemma 5. □

4. Existence of Weak Solutions to the Original Degenerate System

We now pass to the limit $\epsilon \to 0$ to recover a solution of the original degenerate system (1).

Theorem 6. Under Assumption 1, the original system (1) admits at least one weak solution $\left(y,z\right)$ .

Proof of Theorem 6. From Lemma 5, the families ${\left\{y^{\epsilon }\right\}}_{\epsilon >0}$ and ${\left\{z^{\epsilon }\right\}}_{\epsilon >0}$ are uniformly bounded in $\mathscr{H}$ and their material derivatives ar bounded in ${\mathscr{H}}_3$ . By the Aubin-Lions lemma (see [22] [24]), there exist subsequences (still denoted $y^{\epsilon }$ and $z^{\epsilon }$ ) such that

$\left(y^{\epsilon },z^{\epsilon }\right)\to \left(y,z\right),\text{strongly}\text{in}{\left[L^2\left(Q\right)\right]}^2.$ (76)

Moreover, ${\partial }_x{y}^{\epsilon }\rightharpoonup {\partial }_{x}y$ weakly in $L^2\left(Q\right)$ and $\left({\partial }_t+{\partial }_a\right)y^{\epsilon }\rightharpoonup \left({\partial }_t+{\partial }_a\right)y$ weakly in $L^2\left(\left(0,T\right)\times \left(0,A\right);{\left(H^1\left(\Omega \right)\right)}^{*}\right)$ .

Since $f_{\epsilon }\left(z^{\epsilon }\right)\to z$ a.e. and $f_{\epsilon }\left(z^{\epsilon }\right)$ is uniformly bounded, we have $f_{\epsilon }\left(z^{\epsilon }\right)\to z$ strongly in $L^2\left(Q\right)$ . Consequently,

$f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }\rightharpoonup z {\partial }_{x}y\text{weakly}\text{in}{L}^2\left(Q\right).$

By passing to the limit when $\epsilon$ tends to zero ($\epsilon \to 0$ ) in the first equation of the weak formulation (4), we obtain for all ${\varphi }_1\in C_c^{\infty }\left(\left[0,T\right)\times \left[0,A\right)\times \Omega \right)$

${\displaystyle {\int }_Q\left(\left[\left({\partial }_t+{\partial }_a\right)y\right]{\varphi }_1+z{\partial }_{x}y{\partial }_x{\varphi }_1+{\mu }_{1}y{\varphi }_1\right)\text{d}t\text{d}a\text{d}x}=0.$ (77)

Proceeding in the same way with the sequence in $z^{\epsilon }$ , we find the following result

${\displaystyle {\int }_Q\left(\left[\left({\partial }_t+{\partial }_a\right)z\right]{\varphi }_2+y{\partial }_{x}z{\partial }_x{\varphi }_2+{\mu }_{2}z{\varphi }_2\right)\text{d}t\text{d}a\text{d}x}=0,$ (78)

with ${\varphi }_2\in C_c^{\infty }\left(\left[0,T\right)\times \left[0,A\right)\times \Omega \right)$ .

To recover the initial and birth conditions, we use the integrated formulation. For any $\varphi \in C^1\left(\left[0,T\right]\times \left[0,A\right];H^1\left(\Omega \right)\right)$ with $\varphi \left(T,a,x\right)=\varphi \left(t,A,x\right)=0$ , integration by parts in the weak formulation (4) for $y^{\epsilon }$ yields

$\begin{array}{l}-{\displaystyle {\int }_Q\left(y^{\epsilon }\left({\partial }_t+{\partial }_a\right){\varphi }_1+f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }{\partial }_x{\varphi }_1+{\mu }_1{y}^{\epsilon }{\varphi }_1\right)\text{d}t\text{d}a\text{d}x}\\ ={\displaystyle {\int }_{Q_A}{y}^{\epsilon }\left(0,a,x\right){\varphi }_1\left(0,a,x\right)\text{d}a\text{d}x}+{\displaystyle {\int }_{Q_T}{y}^{\epsilon }\left(t,0,x\right){\varphi }_1\left(t,0,x\right)\text{d}t\text{d}x}.\end{array}$ (79)

Passing to the limit $\epsilon \to 0$ in (79) and using the strong convergence of $y^{\epsilon }$ and $f_{\epsilon }\left(z^{\epsilon }\right){\partial }_x{y}^{\epsilon }\rightharpoonup z{\partial }_{x}y$ , we get

$\begin{array}{l}-{\displaystyle {\int }_{Q}y\left({\partial }_t+{\partial }_a\right){\varphi }_1\text{d}t\text{d}a\text{d}x}-{\displaystyle {\int }_{Q}z{\partial }_{x}y{\partial }_x{\varphi }_1\text{d}t\text{d}a\text{d}x}-{\displaystyle {\int }_Q{\mu }_{1}y{\varphi }_1\text{d}t\text{d}a\text{d}x}\\ ={\displaystyle {\int }_{Q_A}y\left(0,a,x\right){\varphi }_1\left(0,a,x\right)\text{d}a\text{d}x}+{\displaystyle {\int }_{Q_T}y\left(t,0,x\right){\varphi }_1\left(t,0,x\right)\text{d}t\text{d}x}.\end{array}$ (80)

Integrating by parts in the limit Equation (77) gives the same expression, hence by uniqueness of the limit we identify

$y\left(t,0,x\right)={\displaystyle {\int }_0^A{\beta }_1\left(t,a\right)y\left(t,a,x\right)\text{d}a},y\left(0,a,x\right)=y_0\left(a,x\right).$ (81)

The same holds for $z$ .

Moreover since $f_{\epsilon }\left(z^{\epsilon }\right)\ge \epsilon >0$ and $f_{\epsilon }\left(y^{\epsilon }\right)\ge \epsilon >0$ , these conditions are well-defined and imply ${\partial }_x{y}^{\epsilon }={\partial }_x{z}^{\epsilon }=0$ at the boundary.

As $\epsilon \to 0$ , using the convergences $f_{\epsilon }\left(z^{\epsilon }\right)\to z$ , $f_{\epsilon }\left(y^{\epsilon }\right)\to y$ , ${\partial }_x{y}^{\epsilon }\rightharpoonup {\partial }_{x}y$ , and ${\partial }_x{z}^{\epsilon }\rightharpoonup {\partial }_{x}z$ , we recover the original degenerate conditions

$z{\partial }_{x}y=0,y{\partial }_{x}z=0\text{on}\Sigma ,$

in the weak sense. This completes the proof of Theorem 6. □

5. Conclusion and Perspectives

In this work, we have studied a class of nonlinear degenerate parabolic systems modeling the dynamics of two age-structured populations interacting within a one-dimensional heterogeneous spatial environment. The main difficulty of the model lies in the presence of cross-diffusion operators of the form $\text{div}\left(v{\partial }_{x}u\right)$ and $\text{div}\left(u{\partial }_{x}v\right)$ , which induce a strong coupling between the equations and, simultaneously, a double degeneracy whenever one of the population densities vanishes.

To overcome these difficulties, we introduced a regularization strategy based on the family of functions $f_{ϵ}$ satisfying the conditions (2.1), which ensures uniform parabolic regularity while preserving the physical behavior of the solutions for positive densities. Using the Faedo-Galerkin method with a finite element basis, we established the existence of global weak solutions to the regularized system (2.2). The choice of a finite element basis, rather than the classical eigenfunctions of the Laplacian, was crucial to handle the cross-diffusion structure and to derive the necessary compactness estimates.

Uniform a priori estimates, independent of both the Galerkin dimension $n$ and the regularization parameter $ϵ$ , were rigorously established. These estimates, combined with the Aubin-Lions compactness lemma, allowed us to perform a double limiting procedure: first as $n\to \infty$ to obtain solutions of the regularized system, and then as $ϵ\to 0$ to recover weak solutions of the original degenerate system (1.1). We also proved the positivity of the solutions, which is biologically essential since the unknowns represent population densities.

Our results extend classical Gurtin-MacCamy type models by incorporating nonlinear and coupled spatial diffusion, and complement existing works on reaction-diffusion systems by including age structuring and nonlocal renewal laws. The methods developed in this article provide a robust framework for studying similar degenerate cross-diffusion systems arising in mathematical biology and ecology.

Several directions for future research emerge naturally from this work:

1) Long-time behavior: The asymptotic behavior of the solutions as $t\to \infty$ deserves a thorough investigation. In particular, one could study the existence and stability of steady states, the extinction or persistence of populations, and the possible emergence of spatial patterns. Such results would have significant biological implications for understanding the long-term dynamics of interacting populations.

2) Multi-dimensional extensions: The extension of our results to higher spatial dimensions ($d\ge 2$ ) is a natural and important direction. However, this extension presents additional technical difficulties, particularly in the compactness arguments and in the treatment of the boundary conditions. The use of more sophisticated tools, such as the compensated compactness method or the theory of Young measures, may be necessary.

3) More general cross-diffusion structures: Our method could be adapted to study more general cross-diffusion systems of the form

${\partial }_{t}u+{\partial }_{a}u-\Delta \left(\Phi \left(u,v\right)\right)+{\mu }_{1}u=0,{\partial }_{t}v+{\partial }_{a}v-\Delta \left(\Psi \left(u,v\right)\right)+{\mu }_{2}v=0,$

where Φ and Ψ are suitable nonlinear functions. This includes, for example, the SKT model [14] and other chemotaxis-type systems.

4) Stochastic perturbations: Environmental fluctuations are ubiquitous in biological systems. Incorporating stochastic perturbations into the model, such as multiplicative noise, and studying the resulting stochastic partial differential equations would provide a more realistic description of population dynamics.

5) Numerical simulations: The development of efficient and accurate numerical schemes for the degenerate system is an important direction for applications. The finite element framework developed in this work could serve as a starting point for designing numerical methods, provided that the stability and convergence of the schemes are rigorously established.

In summary, this work provides a solid mathematical foundation for the analysis of degenerate cross-diffusion systems in age-structured population dynamics. We hope that the techniques developed herein will stimulate further research in this fascinating and interdisciplinary field, bridging the gap between mathematical analysis, numerical methods, and biological applications.

Acknowledgements

Sincere thanks to the members of JAMP for their professional performance, and special thanks to managing editor Hellen XU for a rare attitude of high quality.