# Spectral properties of a semi-discretized McKendrick equation

Laurent Attias, Vincent Siess, Stéphane Labbé

<figure><img src="image_1.png" /></figure>

## To cite this version:

Laurent Attias, Vincent Siess, Stéphane Labbé. Spectral properties of a semi-discretized McKendrick equation. 2025. hal-04698907v2

HAL Id: hal-04698907

https://hal.science/hal-04698907v2

Preprint submitted on 17 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Type of article

Spectral properties of a semi-discretized McKendrick equation

Laurent Attias¹*, Stéphane Labbé¹ and Vincent Siess²

¹ Sorbonne Université, CNRS, Université Paris Cité, Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France

² Dassault Systèmes, 10 rue Marcel Dassault, 78140 Vélizy-Villacoublay, France

* Correspondence: laurent.attias@sorbonne-universite.fr.

**Abstract:** In this paper, we perform a numerical study on the eigenvalues of the operator of a semi-discretized age-structured predator-prey model. The age-structuration is based on an agile modeling framework for population dynamics. We also prove a theoretical result on the distribution of the eigenvalues of a semi-discretized McKendrick equation. We show that, under regularity assumptions, and after a rescaling step, the eigenvalues converge to the unit circle as the number of age groups goes to infinity. A sufficient condition to get this convergence consists in having piecewise-constant birth and death rates; by using random values of birth and death rates, we also show numerically that having not enough regular values invalidates the result. We complete this study with a fine analysis of the properties of the original McKendrick equation, particularly highlighting that the eigenvectors form a complete basis of the function space L²(]0, A[), where A is the maximal age considered in our population.

**Keywords:** Population dynamics, McKendrick equation, predator-prey model, spectral problems, linear stability

1. Introduction

The McKendrick equation and the predator-prey model are well known population dynamics models. The McKendrick-von Foerster equation first appeared in articles dealing with epidemiological and disease modeling, but dealt with demographic modeling [1, 2]. This equation is classically written

$$ \partial_t\phi(t, a) + \partial_a\phi(t, a) + \mu(a)\phi(t, a) = 0, $$
where $\phi$ is the population density, $a$ the age and $\mu$ the mortality rate, and completed with a nonlocal boundary condition

$$\phi(t, 0) = \int_{0}^{+\infty} b(a)\phi(t, a) \, da,$$

where $b$ is the natality rate, and an initial condition

$$\phi(0, a) = \phi_0(a).$$

The mortality and natality rates are nonnegative functions. Moreover $b$ is usually zero out of a fertility interval $[a_0, a_1]$, with $0 < a_0 < a_1$.

The McKendrick equation constitutes an important tool in demography, and is often rediscovered by mathematicians studying nonlinear population dynamics problems [3, 4].

Several approaches exist to study structured population dynamics. Among them, let us cite integral formulation [5], difference equations [6, 7], partial differential equations [1, 8] and projection matrix [9, 10]. On the fundamental PDE approach, some preliminary observations have been carried on the spectrum of the operator involved in the integral formulation. As a compact operator, it has a countable spectrum, and one can express the solution as a superposition of eigenfunctions. In [11], authors compare the McKendrick PDE to the abovementioned approaches and derive solutions under suitable hypotheses.

More recently, the introduction of nonlinear models breathed new life into age-structured populations study. In [12], using nonlinear integral Volterra equations, authors showed the existence, uniqueness and convergence to an equilibrium of solutions of a nonlinear Sharpe-Lotka-McKendrick model. This work paved the way to breakthroughs on nonlinear models, as much on theoretical developments as on biological applications (see [13] for an exhaustive approach). The growing mathematical complexity of nonlinearities in age-structured models led to developments of new tools and methods, such as linear and nonlinear operator semigroups in Banach spaces, with a functional analysis approach [14].

When a population's description goes beyond the assumptions of [11], it is often necessary to introduce new methods for modeling and analyzing such complex structured population models. For a complete view on structured populations, we refer to [15]. We also mention the works of [16, 17] for an approach of structured population dynamics models which presents similarities with the method we develop, presented below.

At Dassault Systèmes, we designed an agile modeling framework and developed a model generator, able to generate simulatable models conforming to assumptions and population structure given as input [18].

The framework's formalism is an extension of the McKendrick equation, where the mortality term $\mu(a)\phi(t, a)$ is replaced with generic exchanges terms, such as

$$\int_{\mathbb{D}} \kappa(t, \phi(t), a, x, y)\phi(t, a, y) \, dy,$$

where $x$ (resp. $y$) represents the origin (resp. destination) state, which belongs to the state domain $\mathbb{D}$.
A conservative form of exchange terms, written as

$$ \int_D (k(t, \phi(t), a, x, y)\phi(t, a, y) - k(t, \phi(t), a, y, x)\phi(t, a, x)) \, dy, $$

is also studied, to model exchanges between the compartments of our structured and generic population.

To simulate structured models, we operate a semi-discretization on the age derivative, and get a dynamical system of ODE simulatable by Dassault Systèmes's solvers embedded in the Dymola Software [19]. The input format of Dymola being models written in the Modelica language [20], the modeler generates Modelica files by coupling an object-oriented implementation of the structured population formalism with a templating approach, using the Python library Jinja2 [21].

In our previous work, we applied this framework to the simulation of a structured predator-prey model, exploring various and interesting phenomena in terms of asymptotic behaviors. The issue of such population system's stability comes naturally in the discussion, wondering under which assumptions a population system converges to a fixed point, or ends up with a cyclic trajectory, or even presents a chaotic evolution.

The issue of stability and asymptotic behavior is one of the most common when dealing with the analysis of dynamical systems of ODEs. Linear stability refers to the ability of a system to come back to an equilibrium point when it undergoes small perturbations [22]. The usual approach consists in writing the jacobian matrix of the differential system and examining its eigenvalues. The sign of the real parts of eigenvalues gives insight on the equilibrium point's nature. See [23-25] for further perspectives on linear and nonlinear stability analysis. Beyond linear and local stability, many definitions of stability for dynamical systems have been given, such as the stability in the sense of Lyapunov [26], Poincaré [27, 28] or Lagrange [29, 30]. Lyapunov's method addresses stability problems both near and far from the equilibrium. It consists in finding a positive definite function, called a Lyapunov function, whose derivative's sign will give insight on the equilibrium point asymptotic stability. To go further into nonlinear stability analysis and Lyapunov method, see [30, 31]. We also cite [32] as an exhaustive reference on ordinary differential equations.

To tackle this issue of stability and asymptotic behavior, we got interested in computing the eigenvalues of classic demographic models. We observed that the eigenvalues were distributed in circular structures. Our goal is to better understand these structures. The fact that these structures appear even in a very simple problem such as the semi-discretized McKendrick equation supports the intuition that they are due to the discretization process.

This article is divided in three sections. In section 2, we derive from the well-known predator-prey model an age-structured, semi-discretized formulation, and carry numerical observations on the distribution of jacobian's eigenvalues at equilibrium points. We have been particularly surprised by the shape of such eigenvalues. They seem to form circular figures whose radius seems to be unbounded with the number of age classes, and tangent to the y-axis. However, studying this problem appeared
difficult, due to multiple nonlinearities present in the Lotka-Volterra equations. Thus, we tried to see if a similar phenomenon could be observed in a simpler problem. This is why, in section 3, we study the McKendrick-von Foerster equation, semi-discretized in age, and prove a theoretical result on the distribution of the system's eigenvalues. We show that, in some sense, and under some regularity assumptions, the set of eigenvalues is dense into the unit circle of the complex plane. We then discuss the role of the regularity assumptions Interestingly, this result also applies to the Leslie matrix. In section 4, we recall some classic results on the McKendrick equation, and prove that the eigenvectors form a complete basis. We observe numerically that the first eigenvalues of the discrete problem converge to the ones of the continuous problem.

## 2. Age-structured predator prey model

### 2.1. Predator-prey model and age-structuration

#### 2.1.1. The predator-prey model

The predator-prey model and the Lotka-Volterra equations are among the most known models in population dynamics. For the sake of exhaustiveness, let's recall the classic Lotka-Volterra equations [33, 34]:

$$ \begin{cases} \frac{d}{dt}x(t) = \alpha x(t) - \beta y(t)x(t) \\ \frac{d}{dt}y(t) = \delta x(t)y(t) - \gamma y(t) \end{cases} \qquad (2.1) $$

where:

* x is the number of preys,

* y is the number of predators,

* α is the natality rate for preys,

* β describes the effect of the presence of predators on the prey mortality rate,

* γ is the mortality rate for predators,

* δ describes the effect of the presence of preys on the predator natality rate.

It is well known that for positive parameters α, β, γ, δ, the system (2.1) has a periodic behavior (see Figure 1).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Prey-predator model</figcaption></figure>

We refer the reader to [35–37] for detailed analyses on the Lotka-Volterra equations.

### 2.1.2. Age-structuration of the predator-prey model

Age-structured predator-prey models have been widely studied in the literature on population dynamics. We may cite [38–40] for recent breakthroughs on this topic.

A classic approach to reach a further level of realism consists in introducing a structuration on the age of individuals. We introduce a second, positive variable *a* standing for the age, a transport term as a partial age derivative and adapt the shapes of rates linked to predation. We get the system (2.2a)-(2.2f):

$$ \partial_t x(t, a) + \partial_a x(t, a) + \left( \int_0^{+\infty} \beta(a, \tilde{a}) y(t, \tilde{a}) \, d\tilde{a} \right) x(t, a) = 0, \quad (2.2a) $$

$$ \partial_t y(t, a) + \partial_a y(t, a) + \gamma(a)y(t, a) = 0, \quad (2.2b) $$

$$ x(t, 0) = \int_{0}^{+\infty} \alpha(a)x(t, a) \, da, \quad (2.2c) $$

$$ y(t, 0) = \int_{0}^{+\infty} \int_{0}^{+\infty} \delta(a, \tilde{a}) x(t, \tilde{a}) y(t, a) \, d\tilde{a} \, da, \quad (2.2d) $$

$$ x(0, a) = x_0(a), \quad (2.2e) $$

$$ y(0, a) = y_0(a), \quad (2.2f) $$

where we added initial conditions $x_0$ and $y_0$. The reproduction terms now appear in boundary conditions (2.2c)-(2.2d), as they affect only the populations at age zero. The rates $\alpha, \gamma$ (resp. $\beta, \delta$) are to be chosen with compact support in $\mathbb{R}_+^*$ (resp. $\mathbb{R}_+^* \times \mathbb{R}_+^*$) to ensure existence and uniqueness of the solution of this system.
There exists a wide diversity of behaviors which can be observed in a such age-structured system. Contrary to Equation 2.1, which had parameters in $\mathbb{R}^4$, the space of parameters in Equation 2.2a-Equation 2.2f has infinite dimension. It is interesting noticing that the resulting behavior is not necessarily periodic; it may also, for example, converge to a fixed point (see Figure 2).

<figure><img src="image_5.png" /><figcaption>Figure 2. Age-structured prey-predator model</figcaption></figure>

This figure has been computed using a semi-discretization on age and Dassault Systèmes ODE solver, Dymola. In the next subsection, we explain the step of semi-discretization.

## 2.2. Semi-discretization of the age-structured predator-prey model

Many approaches exist to achieve a numerical simulation of Equation 2.2a-Equation 2.2f. The approach we took on consists in two steps :

1. derive from Equation 2.2a-Equation 2.2f an age semi-discretized formulation,

2. solve numerically the semi-discretized formulation thanks to ODE solvers.

We first truncate the age interval to [0, A], with A > 0 representing a maximal age. Let consider a discretization of [0, A] into N intervals $[a_{i-\frac{1}{2}}, a_{i+\frac{1}{2}}]$ for i in 0, ..., N - 1, with $a_{-\frac{1}{2}} = 0$ and $a_{N-\frac{1}{2}} = A$. The discretization will be regular, with $\Delta a_N = a_{i+\frac{1}{2}} - a_{i-\frac{1}{2}}$ independent of i.

The choice of a numerical scheme is not the heart of our work. What follows could be adapted with any numerical scheme. For the sake of simplicity, we choose an upwind scheme, where, for i in 0, ..., N - 1,

$$ egin{cases} x(t, a_{i+\frac{1}{2}}) & \text{is approximated by} \\ y(t, a_{i+\frac{1}{2}}) & \text{is approximated by} \end{cases} \frac{1}{\Delta a_N} \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} x(t, a) \, da \\ \frac{1}{\Delta a_N} \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} y(t, a) \, da \qquad (2.3) $$

We also define the aggregated population variables

$$ X_i(t) = \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} x(t, a) \, da, $$
$$Y_i(t) = \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} y(t, a) \, da.$$

Applying the semi-discretization (2.3) to the age-structured continuous system (2.2a)-(2.2f) leads to the dynamical system (2.4a)-(2.4d).

$$d_i X_i + \frac{X_i - X_{i-1}}{\Delta a_N} + \left( \sum_{k=0}^{N-1} \beta_{ik} Y_k \right) X_i = 0, \quad 1 \le i \le N-1, \quad (2.4a)$$

$$d_i X_0 + \frac{1}{\Delta a_N} X_0 - \sum_{k=0}^{N-1} \alpha_k X_k + \left( \sum_{k=0}^{N-1} \beta_{0k} Y_k \right) X_0 = 0, \quad (2.4b)$$

$$d_i Y_i + \frac{Y_i - Y_{i-1}}{\Delta a_N} + \gamma_i Y_i = 0, \quad 1 \le i \le N-1, \quad (2.4c)$$

$$d_i Y_0 + \frac{1}{\Delta a_N} Y_0 - \sum_{j=0}^{N-1} \left( \sum_{k=0}^{N-1} \delta_{jk} X_k \right) Y_j + \gamma_0 Y_0 = 0, \quad (2.4d)$$

where the boundary conditions (2.2c)-(2.2d) have been adapted into special terms in Equation 2.4b-Equation 2.4d, as input population flows into the zero-aged class.

The discrete parameters $(\alpha_i)_{0 \le i \le N-1}$, $(\beta_i)_{0 \le i \le N-1}$, $(\gamma_i)_{0 \le i \le N-1}$, $(\delta_i)_{0 \le i \le N-1}$ are derived from the continuous ones :

$$\left\{
\begin{aligned}
\alpha_i &= \frac{1}{\Delta a_N} \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} \alpha(a) \, da \\
\beta_{ij} &= \frac{1}{\Delta a_N^2} \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} \int_{a_{j-\frac{1}{2}}}^{a_{j+\frac{1}{2}}} \beta(a, \tilde{a}) \, d\tilde{a} \, da \\
\gamma_i &= \frac{1}{\Delta a_N} \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} \gamma(a) \, da \\
\delta_{ij} &= \frac{1}{\Delta a_N^2} \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} \int_{a_{j-\frac{1}{2}}}^{a_{j+\frac{1}{2}}} \delta(a, \tilde{a}) \, d\tilde{a} \, da
\end{aligned}
\right.$$

The choice of functions $\alpha, \beta, \gamma, \delta$ follows usual hypotheses, which may depend on the type of ecosystem one wants to model.

The limitation of our analysis to age semi-discretized models may seem restrictive, but our goal was not to study in a systematic way fully discrete numerical schemes. In modeling and engineering literature, problems often consist in ODEs with a discrete subdivision of the population into age classes. A classical example is the World3 model [41], where the population is divided into four age classes. We can notice that in this literature, the authors are not always aware of the connection between their age-structured equations and semi-discretized McKendrick formulation.

## 2.3. Numerical observations

A natural step, to study the asymptotic behavior of the system (2.4a)-(2.4d), consists in studying the eigenvalues of the jacobian at an equilibrium point. To reach such an equilibrium, we may numer-
ically run the simulation until the system seems to reach an equilibrium point, and study the jacobian
at this final point. We ensure that the system has converged by numerically testing that enough points
are sufficiently close before stopping the computation.

We are interested in the sign of the real part of the first eigenvalue, which characterizes the stability
properties of the system. As one can easily compute, thanks to linear algebra tools, all the eigenval-
ues of the system (2.4a)-(2.4d), we decided to give a look at the whole spectrum, and how its shape
depends on the number of age groups. We get the Figure 3.

<figure><img src="image_6.png" /><figcaption>(a) N = 100</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) N = 200</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) N = 300</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) N = 400</figcaption></figure>

**Figure 3.** Eigenvalues of the age-structured semi-discretized predator-prey model, at the end of the simulation, for various value of number of age groups *N*.
The eigenvalues distribution observed on Figure 3 is somewhat interesting, with circles tangent to the y-axis with a center which seems to translate along the x-axis, as N increases; we also observe inner structures, which go from points to small circles as N increases.

### 3. Semi-discretized McKendrick equation

The numerical observations on Figure 3 are of interest, but the predator-prey model and the age-structured Lotka-Volterra equations are very complex to study. To simplify the problem at hand, we will study a McKendrick equation. Contrary to the predator-prey model and its nonlinearities due to predation, the McKendrick equation is linear. As a result, the infinitesimal generator in the discrete age-structured formulation is linear, so studying the stability through jacobian eigenvalues is equivalent to studying the eigenvalues of the operator itself.

Numerically, we may formulate the following conjecture (see Figure 4) : asymptotically, the eigenvalues set in the complex plane seems to describe a circle, with center located at $(-\frac{1}{\Delta a_N}, 0)$, with radius $\frac{1}{\Delta a_N}$.
<figure><img src="image_10.png" /><figcaption>(a) N = 100</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) N = 200</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) N = 300</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) N = 400</figcaption></figure>

<caption>**Figure 4.** Eigenvalues of the semi-discretized McKendrick equation, for various value of number of age groups N.</caption>

### 3.1. Semi-discretization

Let us recall the well-known McKendrick equation :

$$ egin{cases} 
\partial_t \phi(t, a) + \partial_a \phi(t, a) + \mu(a) \phi(t, a) = 0, \quad \forall (t, a) \in ]0, T] \times \mathbb{R}_+^*, \\ 
\phi(t, 0) = \int_0^{+\infty} b(a) \phi(t, a) da, \quad \forall t \in ]0, T], \\ 
\phi(0, a) = \phi_0(a) \quad \forall a \in \mathbb{R}_+^*, 
\end{cases} 
\quad (3.1) $$
where $\phi$ is the population density, $\mu$ is the death rate, $b$ the birth rate, $a$ the age of individuals. The semi-discretization of Equation 3.1 with an upwind-scheme is straight-forward:

$$ \left\{ \begin{array}{l} d_t \psi_i + \frac{\psi_i - \psi_{i-1}}{\Delta a_N} + \mu_i \psi_i = 0 \quad 1 \le i \le N-1 \\ d_t \psi_0 + \frac{1}{\Delta a_N} \psi_0 - \sum_{k=0}^{N-1} b_k \psi_k + \mu_0 \psi_0 = 0 \end{array} \right. \qquad (3.2) $$

where:

* $N$ is the number of age groups,

* $\psi_i(t) = \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} \phi(t, a) da$, is a discrete population variable, between ages $a_{i-\frac{1}{2}}$ and $a_{i+\frac{1}{2}}$, for $i$ in $0, \dots, N-1$, following the notations of subsection 2.2,

* $b_i = \frac{1}{\Delta a_N} \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} b(a) da$,

* $\mu_i = \frac{1}{\Delta a_N} \int_{a_{i-\frac{1}{2}}}^{a_{i+\frac{1}{2}}} \mu(a) da$.

## 3.2. Eigenvalues problem

The eigenvalues problem associated to Equation 3.2 may be written as the following:

$$ \left\{ \begin{array}{l} - \left( \frac{\psi_i - \psi_{i-1}}{\Delta a_N} + \mu_i \psi_i \right) = \lambda \psi_i \quad 1 \le i \le N-1 \\ \Delta a_N \psi_0 - \sum_{k=0}^{N-1} b_k \psi_k = \lambda \psi_0 \end{array} \right. \qquad (3.3a) $$

$$ \left\{ \begin{array}{l} - \left( \frac{1}{\Delta a_N} \psi_0 - \sum_{k=0}^{N-1} b_k \psi_k + \mu_0 \psi_0 \right) = \lambda \psi_0 \\ \Delta a_N \psi_0 - \sum_{k=0}^{N-1} b_k \psi_k = \lambda \psi_0 \end{array} \right. \qquad (3.3b) $$

where $\lambda$ is a complex eigenvalues and $(\psi_i)_{0 \le i \le N-1}$ is an eigenvector in $\mathbb{C}^N$.
Some steps of computation give, for $i$ in $0, \dots, N-1$,

$$ \psi_i = \frac{\psi_0}{\prod_{k=1}^{i} (1 + \Delta a_N (\mu_k + \lambda))} $$

We put back this expression into Equation 3.3b, which gives after simplification

$$ -1 + \Delta a_N \sum_{i=0}^{N-1} \frac{b_i}{\prod_{k=0}^{i} (1 + \Delta a_N (\mu_k + \lambda))} = 0 $$

Multiplying this equation by $\prod_{k=0}^{N-1} (1 + \Delta a_N (\mu_k + \lambda))$ gives the following characteristic equation:
$$ -\prod_{k=0}^{N-1} (1 + \Delta a_N(\mu_k + \lambda)) + \Delta a_N \sum_{i=0}^{N-1} b_i \prod_{k=i+1}^{N-1} (1 + \Delta a_N(\mu_k + \lambda)) = 0 \quad (3.4) $$

From now, for the sake of realism, we will assume an hypothesis on fertile age groups: for any given number of age groups $N$ set in the semi-discrete model, there exists two integers $n$ and $m$ such that for $i < n - m$ or $i > n$, $b_i = 0$. $n$ is the index of the older fertile age group, $m$ is the width of the fertile age groups range. Because the fertility age range in a population is supposed to be independent from the age structure, we will assume the existence of reals $\alpha$ and $\beta$, strictly between 0 and 1, such that $n \sim \beta N$, $m \sim \alpha N$. From now, we will suppose that the problem may be expressed in terms of $n$, and will denote $\Delta a_N$ by $\Delta a_n$. We will denote by $C > 0$ a real such that for any $n \in \mathbb{N}^*$, $\Delta a_n Bn \le C$. We also achieve a rescaling of the eigenvalue, by setting $z = 1 + \Delta a_n \lambda$. The characteristic equation (3.4) rewrites, under these assumptions,

$$ \prod_{k=1}^{n} (z + \Delta a_n \mu_k) = \Delta a_n \sum_{i=n-m}^{n} b_i \prod_{k=i+1}^{n} (z + \Delta a_n \mu_k) \quad (E_n) $$

### 3.3. Asymptotic distribution of eigenvalues

We now present our main result on the circular distribution of eigenvalues for the semi-discretized McKendrick equation.

To achieve that, we assume an additional and strong hypothesis: the birth and death rates are piecewise-constant. This assumption is, to be sure, restrictive, but adapted to the shape of data on birth rates and death rates we get from official statistical institutes. In a perspective of numerical simulation and model calibration, this hypothesis is not much restrictive.

**Theorem 3.1** (Asymptotic properties of eigenvalues of Equation $E_n$). We assume piecewise-constant birth and death rates. Then we get the following results:

1. Let $(z_n)_{n \in \mathbb{N}}$ be a sequence of solutions of Equation $E_n$. Then

$$ \lim_{n \to +\infty} |z_n| = 1 $$

2. For any $z_0 = e^{i\theta_0}$ in $\mathbb{U}$, there exists a sequence $(z_n)_{n \in \mathbb{N}}$ of solutions to Equation $E_n$ converging to $z_0$.

This theorem shows a result of density of the solutions to Equation $E_n$ in the unit circle. We provide below a sketch of demonstration.

*Proof.* As birth and death rates $(b_i)_{1 \le i \le N}$ and $(\mu_i)_{1 \le i \le M}$ are assumed piecewise-constant, they are bounded by reals denoted respectively by $B$ and $M$.

1. • Let $R = \limsup_{n \to +\infty} |z_n|$. We prove by contradiction that $R \le 1$. If $R > 1$, up to a subsequence, we may suppose that for any $n \in \mathbb{N}^*$, $|z_n| \ge \frac{R+1}{2}$. We denote by $k_n \in [1, n]$ the index such that
$|z_n + \Delta a_n \mu_{k_n}| = \min_{1 \le k \le n} |z_n + \Delta a_n \mu_k|$. We have the following inequalities:

$$|z_n + \Delta a_n \mu_{k_n}|^n \le \prod_{k=1}^{n} |z_n + \Delta a_n \mu_k| \le \Delta a_n \sum_{i=n-m}^{n} b_i \prod_{k=i+1}^{n} |z_n + \Delta a_n \mu_k|,$$ $$|z_n + \Delta a_n \mu_{k_n}|^n \le \Delta a_n Bn (|z_n| + \Delta a_n M)^m,$$

$$|z_n| - \Delta a_n \mu_{k_n} \le |z_n + \Delta a_n \mu_{k_n}| \le (\Delta a_n Bn)^{\frac{1}{n}} (|z_n| + \Delta a_n M)^{\frac{m}{n}},$$

$$|z_n| \le C_n^{\frac{1}{n}} (|z_n| + \Delta a_n M)^{\frac{m}{n}} + \Delta a_n \mu_{k_n},$$

$$|z_n|^{1-\frac{m}{n}} \le C_n^{\frac{1}{n}} \left( 1 + \frac{\Delta a_n M}{|z_n|} \right)^{\frac{m}{n}} + \frac{\Delta a_n \mu_{k_n}}{|z_n|^{\frac{m}{n}}},$$

$$\left(\frac{R+1}{2}\right)^{1-\frac{m}{n}} \le C_n^{\frac{1}{n}} \left(1+\frac{\Delta a_n M}{|z_n|}\right)^{\frac{m}{n}} + \frac{\Delta a_n \mu_{k_n}}{|z_n|^{\frac{m}{n}}}.$$

Taking the limit as $n$ goes to $+\infty$, we get

$$\left(\frac{R+1}{2}\right)^{1-\alpha} \le 1,$$

which is absurd if $R > 1$, as $\alpha < 1$.

* Let $r = \liminf_{n \to +\infty} |z_n|$. We prove by contradiction that $r \ge 1$. If not, we can extract a subsequence, still denoted by $(z_n)_{n \in \mathbb{N}}$, such that $|z_n|$, which converges to a $L < 1$. Under the piecewise-constant hypothesis, the following computation holds on Equation $E_n$:

$$\prod_{k=1}^{n} (z_n + \Delta a_n \mu_k) = \Delta a_n \sum_{j=0}^{K-1} C_j \sum_{i=jp_n}^{(j+1)p_n-1} \left( \prod_{l=0}^{j-1} (z_n + \Delta a_n M_l)^{p_n} \right) \left( z_n + \Delta a_n M_j \right)^{i-jp_n},$$

where $K$ is the number of subintervals on which $b$ and $\mu$ are constant, $p_n$ the number of age groups on each subinterval, and $C_j$ (resp. $M_j$) is the common value of $b_{n-i}$ (resp. $\mu_{n-i}$) for $i$ between $jp_n$ and $(j+1)p_n - 1$,

$$\prod_{j=0}^{K-1} \prod_{k=jp_n}^{(j+1)p_n-1} (z_n + \Delta a_n \mu_{n-k}) = \Delta a_n \sum_{j=0}^{K-1} C_j \left( \prod_{l=0}^{j-1} (z_n + \Delta a_n M_l)^{p_n} \right) \frac{(z_n + \Delta a_n M_j)^{p_n} - 1}{z_n + \Delta a_n M_j - 1},$$

$$\prod_{j=0}^{K-1} (z_n + \Delta a_n M_j)^{p_n} = \Delta a_n \sum_{j=0}^{K-1} C_j \left( \prod_{l=0}^{j-1} (z_n + \Delta a_n M_l)^{p_n} \right) \frac{(z_n + \Delta a_n M_j)^{p_n} - 1}{z_n + \Delta a_n M_j - 1}.$$

We know that

$$\Delta a_n \underset{n \to +\infty}{\sim} \frac{A\beta}{n}$$
and that for $j \neq 0$, $\prod_{l=0}^{j-1} (z_n + \Delta a_n M_l)^{p_n}$ converges to 0 as well as $(z_n + \Delta a_n M_j)^{p_n}$. Then we may extract a subsequence, still denoted by $(z_n)_{n \in \mathbb{N}}$, and have a real $C_1$ such that

$$ \prod_{j=0}^{K-1} (z_n + \Delta a_n M_j)^{p_n} = \frac{C_1}{n} + o\left(\frac{1}{n}\right) $$

Taking the logarithm of each member gives

$$ \sum_{j=0}^{K-1} \ln|z_n + \Delta a_n M_j| = \frac{\ln(C_1) - \ln(n)}{p_n} + o(1). $$

As we have, for any $j \in [0, K-1]$,

$$ |z_n + \Delta a_n M_j| \le |z_n| + \Delta a_n M_j \xrightarrow{n \to +\infty} L < 1, $$

we will have for sufficiently large $n$,

$$ \sum_{j=0}^{K-1} \ln|z_n + \Delta a_n M_j| \le K \ln(L) < 0, $$

and, because $p_n$, the number of age groups on each subinterval, is proportional to the total number of age group $n$, we also have

$$ \lim_{n \to +\infty} \frac{\ln(C) - \ln(n)}{p_n} \xrightarrow{n \to +\infty} 0, $$

we get a contradiction.

2. We denote the average mortality rate $\frac{1}{n} \sum_{k=1}^{n} \mu_k$ by $\bar{\mu}$ and rewrite Equation $E_n$ as a fixed point equation:

$$ \left(z + \frac{\bar{\mu}}{n}\right)^n = \left(z + \frac{\bar{\mu}}{n}\right)^n - \prod_{k=1}^{n} \left(z + \frac{\mu_k}{n}\right) + \frac{\gamma}{n} \sum_{i=0}^{m} b_{n-i} \prod_{k=n-i+1}^{n} \left(z + \frac{\mu_k}{n}\right), \quad (3.5) $$

where $\gamma = A\beta$.

We prove the existence of a solution to Equation 3.5 (and consequently, to Equation $E_n$) in the angular sector $\mathcal{D}_{n,\theta_0}$ defined by

$$ \mathcal{D}_{n,\theta_0} = \left\{ z \in \mathbb{C}, |z| \in [0, 2], \arg(z) \in \left[\theta_0 - \frac{\pi}{n}, \theta_0 + \frac{\pi}{n}\right] \right\}. $$

In order to apply the Brouwer fixed point theorem*, we draw a bijection $\xi_n$ between $\mathcal{D}_{n,\theta_0}$ and $\overline{\mathcal{D}}(0, 2)$,

<footnote>*see e.g. [42] for details and demonstration</footnote>
We rewrite Equation 3.5 as

$$ \left(z + \frac{\bar{\mu}}{n}\right)^n = \varphi_n(z) \qquad (3.6) $$

and denote by $\chi_n$ a monovaluation of the inverse of the map $z \mapsto z^n$. We then rewrite Equation 3.6 as

$$ \tilde{z} = \xi_n \left( \chi_n \circ \varphi_n \circ \xi_n^{-1}(\tilde{z}) - \frac{\bar{\mu}}{n} \right) \qquad (3.7) $$

where $\tilde{z} = \xi_n(z)$ belongs to $\bar{\mathcal{D}}(0, 2)$. The right-hand side of Equation 3.7 is denoted by $f_n(\tilde{z})$ so that a stronger form of Equation 3.5 is equivalent to

$$ f_n(\tilde{z}) = \tilde{z}, \qquad (3.8) $$

One can easily show that $f_n$ is continuous. It remains to prove that $\bar{\mathcal{D}}(0, 2)$ is stable by $f_n$. We have, on the last term of $\varphi_n(z)$, for $z$ in $\bar{\mathcal{D}}(0, 2)$, the following estimation:

$$ \left| \frac{C}{n} \sum_{i=0}^{m} b_{n-i} \prod_{k=n-i+1}^{n} \left(z + \frac{\mu_k}{n}\right) \right| \le \frac{C}{n} B \sum_{i=0}^{m} \left(2 + \frac{M}{n}\right)^i \le \frac{CB}{1 + \frac{M}{n}} \frac{1}{n} \left(2 + \frac{M}{n}\right)^{m+1}, $$

$$ \left| \frac{C}{n} \sum_{i=0}^{m} b_{n-i} \prod_{k=n-i+1}^{n} (z + \mu_k) \right| = \left(2 + \frac{M}{n}\right)^{m+1} o\left(\frac{1}{n}\right). $$

The two first terms of $\varphi_n(z)$ are rewritten as

$$ \left(z + \frac{\bar{\mu}}{n}\right)^n - \prod_{k=1}^{n} \left(z + \frac{\mu_k}{n}\right) = \left(z + \frac{\bar{\mu}}{n}\right)^n \left(1 - \frac{\prod_{k=1}^{n} \left(z + \frac{\mu_k}{n}\right)}{\left(z + \frac{\bar{\mu}}{n}\right)^n}\right), $$

and after some computation, the following estimation holds:

$$ \frac{\prod_{k=1}^{n} \left(z + \frac{\mu_k}{n}\right)}{\left(z + \frac{\bar{\mu}}{n}\right)^n} = 1 + o\left(\frac{1}{n}\right), $$

so that

$$ |\varphi_n(z)| \le O\left(\frac{1}{n}\right) \left( \left(2 + \frac{\bar{\mu}}{n}\right)^n + \left(2 + \frac{M}{n}\right)^{m+1} \right). $$

One can then prove the existence of a $C_1 > 0$ such that, for $n$ large enough,

$$ |\varphi_n(z)| \le \frac{2C1}{n} \left(2 + \frac{M}{n}\right)^n, $$
$$|\chi_n \circ \varphi_n(z)| \le (2C_1)^{\frac{1}{n}} \left( 2 + \frac{M - 2\ln(n)}{n} + o\left(\frac{\ln(n)}{n}\right) \right),$$

$$|\xi_n(\chi_n \circ \varphi_n(z) - \frac{\bar{\mu}}{n})| \le 2 + \frac{M2\ln(2C_1) + \bar{\mu} - 2\ln(n)}{n} + o\left(\frac{\ln(n)}{n}\right),$$

and then conclude that $f_n(z) \in \bar{D}(0, 2)$.

By applying Brouwer's fixed point theorem, we then get a solution of Equation 3.8, which provides by applying $\xi_n$ backwards a solution of Equation 3.5 which belongs to $D_{n,\theta_0}$, i.e. which has argument converging to $\theta_0$ as $n$ goes to infinity.

The conjugation of the two steps shows that the constructed sequence $(z_n)_{n \in \mathbb{N}}$ converges to $z_0$ as $n$ goes to infinity.

□

We can verify that the piecewise regularity assumption of Theorem 3.1 is also crucial here. When taking random uniform identically distributed values of natality or mortality between 0 and 1, one can construct counter-examples that give us eigenvalues inside the unit disk. And this, even with a large number of age groups (see Figure 5). Although stochastic data can be prevalent in biological applications, they have to be excluded from this framework, as they would invalidate Theorem 3.1. A contrario, we are convinced that the result holds for natality and mortality rates with bounded variation, but this question stays open.
<figure><img src="image_14.png" /><figcaption>**Figure 5.** Computation of eigenvalues of Equation 3.3a-Equation 3.3b with random natality and mortality values, with N = 1000 age groups. The target circle is represented for reference.</figcaption></figure>

### 3.4. Relation with the Leslie matrix

In this subsection, we are interested in the dynamic of a population *P* divided into *N* age groups. The dynamic is ruled by the following recurrence equation:

$$P(t + \Delta t) = LP(t) \qquad (3.9)$$
where $\Delta t$ represents a time step and $P(t)$ is a vector of $\mathbb{R}^N$ with, for age group $i$, $P(t)_i$ the population of the $i^{th}$ age group. In the literature, the matrix $L$ is called the Leslie matrix [10]. It is a nonnegative matrix given by

$$L = \begin{pmatrix}
\Delta t b_0 & \cdots & \cdots & \cdots & \Delta t b_{N-1} \\
1 - \Delta t \mu_0 & 0 & \cdots & \cdots & 0 \\
0 & \ddots & \ddots & & \vdots \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & 0 & 1 - \Delta t \mu_{N-2} & 0
\end{pmatrix}$$

We first recall a classic result on the spectrum of the Leslie Matrix, which is an application of the Perron-Frobenius theorem.

**Theorem 3.2.** *The Leslie Matrix has one real eigenvalue which is equal to its spectral radius.*

Obviously there is a strong link between the Leslie matrix and the matrix of the dynamic system (3.2) resulting from the semi-discretization of the McKendrick equation, which will be denoted by $M$. Taking $\Delta a_N = \Delta t$ gives

$$M = \begin{pmatrix}
\frac{1}{\Delta t} + \mu_0 - b_0 & -b_1 & \cdots & \cdots & -b_{N-1} \\
-\frac{1}{\Delta t} & \frac{1}{\Delta t} + \mu_1 & 0 & \cdots & 0 \\
0 & \ddots & \ddots & \ddots & \vdots \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & 0 & -\frac{1}{\Delta t} & \frac{1}{\Delta t} + \mu_{N-1}
\end{pmatrix}$$

We denote by $J$ the following matrix:

$$J = \begin{pmatrix}
\mu_0 & 0 & \cdots & \cdots & 0 \\
-\mu_0 & \mu_1 & 0 & \cdots & 0 \\
0 & \ddots & \ddots & \ddots & \vdots \\
\vdots & \ddots & \ddots & \ddots & 0 \\
0 & \cdots & 0 & -\mu_{N-2} & \mu_{N-1}
\end{pmatrix}$$

Then one immediately sees that

$$\frac{P(t + \Delta t) - P(t)}{\Delta t} = \frac{L - I_N}{\Delta t} P(t) = (-M + J)P(t)$$

so that

$$L = I_N - \Delta t M + \Delta t J$$

As $\Delta t J = O(\Delta t)$, the Leslie Matrix $L$ is a perturbation of matrix $I_N - \Delta t M$, and we have, as $\Delta t$ goes to 0, $L$ close to $I_N - \Delta t M$.

It follows from Theorem 3.1 that $I_N - \Delta t M$ is a matrix whose eigenvalues are close to the unit circle when $N$ is high enough. Consequently, the eigenvalues of $L$ will also be close to the unit circle (as the map $\lambda \mapsto 1 - \Delta t \lambda$ is precisely the rescaling operation applied to the eigenvalues of $M$).
## 4. McKendrick equation spectrum

Let us recall the McKendrick equation introduced in the previous section:

$$ egin{cases} 
\partial_t \phi(t, a) + \partial_a \phi(t, a) + \mu(a)\phi(t, a) = 0, \quad \forall (t, a) \in ]0, T] \times \mathbb{R}_+^*, \\ 
\phi(t, 0) = \int_0^{+\infty} b(a)\phi(t, a) \, da, \quad \forall t \in ]0, T], \\ 
\phi(0, a) = \phi_0(a) \quad \forall a \in \mathbb{R}_+^*, 
\end{cases} 
\quad (3.1) $$

and the assumptions presented in the introduction: $b$ is zero out of a fertility interval $[a_0, a_1]$, and for any $t \ge 0$, $\phi(t, \cdot)$ is zero out of an age interval $[0, A]$, where $A > 0$ is the maximum age. The boundary term may then be rewritten $\int_{a_0}^{a_1} b(a)\phi(t, a) \, da$.

For what follows, we will also assume that the mortality coefficient $\mu(t, a)$ and the natality coefficient $b(t, a)$ are nonnegative $L^\infty(]0, A]$ functions.

It is well known that the problem (3.1) is well posed. Moreover, one can exhibit an analytical solution.

**Theorem 4.1 (Analytical solution of McKendrick equation).** There exists one unique solution to the McKendrick equation. Moreover, if the infertility interval is such that $0 < a_0 < a_1 < A < +\infty$, then the solution is given by the following analytical formula, for all $t < a_0$:

$$ \phi(t, a) = 1_{a>t} \psi_I(t, a) + 1_{a \le t} \psi_B(t, a) $$

with

$$ \psi_I(t, a) = \phi_0(a - t) \exp \left( \int_0^t -\mu(u, a - t + u) \, du \right) $$

and

$$ \psi_B(t, a) = \left( \int_{a_0}^{a_1} \psi_I(t - a, \tilde{a}) \, da \right) \exp \left( - \int_{t-a}^t \mu(u, a - t + u) \, du \right) $$

By recurrence, the formula can be extended on each time intervals of the form $[ia_0, (i + 1)a_0[$.

In what follows, the parameters $\mu$ and $b$ will be independent of the time. Let us introduce the following spectral problem$^\dagger$ in $L^2(]0, A]$ associated with problem (3.1), where the infinite bound in the natality term may be replaced by $A$:

$$ egin{cases} 
-d_a\phi(a) - \mu(a)\phi(a) = \lambda\phi(a) \quad \forall a \in ]0, A[, \\ 
\phi(0) = \int_0^A b(a)\phi(a) \, da. 
\end{cases} 
\quad (4.1) $$

This problem is not self-adjoint, and as often in this kind of situation, it is useful to introduce its adjoint spectral problem

$$ egin{cases} 
d_a\phi^*(a) - \mu(a)\phi^*(a) + b(a)\phi^*(0) = \lambda^*\phi^*(a) \quad \forall a \in ]0, A[, \\ 
\phi^*(A) = 0. 
\end{cases} 
\quad (4.2) $$

<footnote>$^\dagger$Introducing this spectral problem and the following one on $L^2(]0, A]$ or $L^2([0, A])$ would be equivalent.</footnote>
The eigenvectors of the problems (4.1) and (4.2), if they exist, can be computed explicitly. Indeed, if $(\lambda, \phi)$ is an eigenpair solution of Equation 4.1, then

$$ \phi(a) = \phi(0)\exp\left(-\int_{0}^{a} (\lambda + \mu(\tilde{a})) d\tilde{a}\right). $$

Moreover, the eigenvalue $\lambda$ is solution of the following characteristic equation, which is obtained from the boundary condition:

$$ 1 = \int_{0}^{A} b(a)\exp\left(-\int_{0}^{a} (\lambda + \mu(\tilde{a})) d\tilde{a}\right) da. \quad (4.3) $$

If $(\lambda^*, \phi^*)$ is an eigenpair solution of Equation 4.2, then

$$ \phi^*(a) = \phi^*(0) \int_a^A b(\tilde{a}) \exp\left(-\int_a^{\tilde{a}} (\lambda^* + \mu(\tilde{a})) d\tilde{a}\right) d\tilde{a} $$

and the eigenvalue satisfies the same characteristic Equation 4.3 as the direct problem. We shall see later more complete results on the eigenpairs of problems (4.1) and (4.2). Les us say for now that there exists a unique real eigenvalue to Equation 4.1, often called the first or leading eigenvalue, which is associated with a positive eigenvector.

The reason why the spectral problem is important is that the global dynamic of the population will be dictated by the sign of the real part of the first eigenvalue of problem (4.1).

**Theorem 4.2.** Let $(\lambda_1, \phi_1)$ denote the first eigenpair of problem (4.1). Then for any $\varepsilon > 0$, there exists a constant $c$ and a time $T$ such that for any $t > T$, we have

$$ \|\phi(t, \cdot) - \langle \phi_1^*, \phi_I \rangle \exp(\lambda_1 t) \phi_1\|_{L^2([0, A[)} \le c \exp((\lambda_1 - \varepsilon)t). $$

For a proof of this proposition, the reader will usefully read [43].

## 4.1. Compactness of the McKendrick operator

We study the compact nature of the McKendrick operator $T$, defined in the Theorem 4.3 below.

**Theorem 4.3.** The operator $T$, defined by

$$ T: \left( \begin{array}{cc} L^2([0, A[) & \rightarrow L^2([0, A[) \\ f & \mapsto \phi \end{array} \right), $$

where $\phi$ is solution of

$$ \left\{ \begin{array}{l} d_a \phi(a) + \mu(a) \phi(a) = f(a), \\ \phi(0) = \int_0^A b(a) \phi(a) da, \end{array} \right. \quad (4.4) $$

is compact. Its eigenvalues are the inverse of the eigenvalues of Equation 4.1. Moreover, $T$ is positive.

*Proof.* The analytical solution of problem (4.4) is given by the following formula:

$$ \phi(a) = \left( \phi(0) + \int_{0}^{a} \exp\left(\int_{0}^{s} \mu(v) dv\right) f(s) ds \right) \exp\left(-\int_{0}^{a} \mu(s) ds\right), $$
with

$$ \phi(0) \left( 1 - \int_{0}^{A} b(a) \exp \left( - \int_{0}^{a} \mu(s) \, ds \right) da \right) = \int_{0}^{A} b(a) \int_{0}^{a} \exp \left( \int_{a}^{s} \mu(v) \, dv \right) f(s) \, ds \, da. $$

From this the following a priori estimate can be clearly obtained:

$$ \begin{cases} |\phi(0)| \le C \|f\|_2, \\ \|\phi\|_2 \le C \|f\|_2, \\ \|d_a \phi\|_2 \le C \|f\|_2. \end{cases} $$

This ensures the compactness of $T$. In fact, it is possible that the factor $1 - \int_0^A b(a) \exp(-\int_0^a \mu(s) \, ds) \, da$ is equal to 0. In this situation, the right-hand side term must also be equal to 0. Notice that 0 is then the real eigenvalue of problem (4.1), and that the condition satisfied by $f$ precisely corresponds to the Fredholm alternative. The solution of Equation 4.4 will be defined up to a function proportional to the eigenvector associated with the null eigenvalue. In the $L^2$ space quotiented by the first eigenspace, the a priori estimate will still hold. Then the derivative will also be bounded, and the compactness ensured. The positivity is obvious considering the analytical formula. $\square$

The compactness of operator $T$ allows using the following classic functional analysis theorem which gives us information on the eigenvalues and eigenvectors of problem (4.1) (see for instance [44]).

**Theorem 4.4.** Let $T$ be a compact operator. Then its spectrum is at most a countable set, with no accumulation point different from 0. Each nonzero value $\lambda$ of the spectrum is an eigenvalue of $T$ of finite multiplicity. The conjugate value $\bar{\lambda}$ is an eigenvalue of the adjoint operator $T^*$ with the same multiplicity.

## 4.2. Spectral properties

Finally, let us write a summary of the results that are known on the problems (4.1) and (4.2). A demonstration of these results can again be read in [43].

**Theorem 4.5 (Spectral properties McKendrick equation).** [leftmargin=1em]

* The spectral problems (4.1) and (4.2) have one and only one real eigenvalue, solution of the characteristic equation (4.3).

* All the other eigenvalues arise as pairs of complex and conjugate numbers. Their real part is greater than the real part of the unique real eigenvalue.

* There exists at most a countable number of eigenpairs.

* If there is an infinity of eigenvalues, their module tends to infinity.

* If we denote by $(\lambda_k, \phi_k)$ and $(\bar{\lambda}_k, \phi_k^*)$ the eigenpairs of Equation 4.1 and Equation 4.2, then these eigenvectors can be normalized so that they are orthogonal:

$$ \langle \phi_l^*, \phi_k \rangle = \delta_{kl} $$

* In every vertical strip of finite width, there exists at most a finite number of eigenvalues.
Of course, one question remains unanswered. Do the eigenvectors of problem (4.1) form a complete basis? As in the heat equation for instance, we would like to know if we can write for any $f \in L^2([0, A])$:

$$ egin{cases} f = \sum_{k=0}^{+\infty} \langle \phi_k^*, f \rangle \phi_k, \\ T f = \sum_{k=0}^{+\infty} \frac{1}{\lambda_k} \langle \phi_k^*, f \rangle \phi_k. \end{cases} $$

We will show that this result is true.

**Theorem 4.6.** The eigenvectors of the spectral problem (4.1) form a complete basis of the space $L^2([0, A])$.

*Proof.* The operator $T$ is a compact operator. Consequently, its spectrum is at most countable and the following decomposition holds for any meromorphic function $f$ (see [45] for instance):

$$ f(T) = \sum_{\lambda \in \Lambda_T} \frac{1}{2i\pi} \oint_{C_\lambda} f(z) (zI - T)^{-1} dz, $$

where the sum runs over the values in the spectrum $\Lambda_T$ of $T$, and $C_\lambda$ is a contour encircling the isolated eigenvalue $\lambda$. For $f = I$, this yields

$$ T = \sum_{\lambda \in \Lambda_T} \lambda \frac{1}{2i\pi} \oint_{C_\lambda} (zI - T)^{-1} dz + \frac{1}{2i\pi} \oint_{C_\lambda} (z - \lambda) (zI - T)^{-1} dz \quad (4.5) $$

Let us introduce the following operators

$$ T_{\lambda,0} = \frac{1}{2i\pi} \oint_{C_\lambda} (zI - T)^{-1} dz, $$

$$ T_{\lambda,1} = \frac{1}{2i\pi} \oint_{C_\lambda} (z - \lambda) (zI - T)^{-1} dz. $$

The operator $T_{\lambda,0}$ is a projection operator on the eigenspace corresponding to the eigenvalue $\lambda$, and $T_{\lambda,1}$ a nilpotent operator. The decomposition formula (4.5) is the so-called Dunford decomposition. We know explicitly the formula for the eigenvectors corresponding to a given eigenvalue. Obviously, any eigenspace is of dimension 1. Consequently, if the nilpotent term is in fact null, the Dunford decomposition implies that there is an infinite number of eigenvalues, and that their eigenspaces form a complete basis allowing to decompose the identity and to diagonalize the operator $T$. It remains to see if the nilpotent term is indeed null. Remembering that the eigenvalues of $T$ are the inverse of the eigenvalues introduced in Equation 4.1, we want to know if we can find a solution $f$ to the equation

$$ (T - \frac{1}{\lambda_k} I) f = \phi_k, $$

where $T\phi_k = \frac{1}{\lambda_k}\phi_k$. Then the Fredholm alternative implies that $\langle \phi_k^*, \phi_k \rangle = 0$ which is not true since its value is equal to 1. Thus, the nilpotent space is null, which concludes the proof. □
To end this section, let us give two numerical examples. In the first example, the natality and mortality rates are constant. It is easy then to compute the eigenvectors once the eigenvalues are known. Nevertheless, the eigenvalues must be calculated numerically. We will denote this set of eigenvalues by "continuous spectrum". The spectrum displays a very classic cusp shape as can be seen on Figure 6.

<figure><img src="image_15.png" /><figcaption>**Figure 6.** Complex eigenvalues of the problem (4.1), with constant natality and mortality rates.</figcaption></figure>

Plotting on the same figure both discrete and continuous spectrums shows an interesting link between the two problems. At a large scale, it seems that the continuous spectrum is tangent to the circle constituting the discrete spectrum (Figure 7a). At a small scale, we observe a correspondence in the shapes of the two spectrums. We highlight the excellent convergence of the first eigenvalue of the discrete spectrum towards continuous spectrum's one. The following eigenvalues have an acceptable convergence error. Indeed, the convergence of the discrete spectrum to a circle seems to be nonuniform, as the discrete spectrum presents the characteristic cusp shape (Figure 7b). This observation paves the way to a deeper study of the links between discrete and continuous spectrums. The conjecture that the continuous spectrum can be obtained from the discrete one, taking the limit on the number of age groups, seems to be reasonable.
<figure><img src="image_16.png" /><figcaption>Figure 7. Continuous and discrete spectrums of the McKendrick equation</figcaption></figure>

In the second example, the mortality is constant, but the natality is concentrated in a Dirac operator $b(a) = b\delta_{a=a_0}$. Then the eigenpairs are given by

$$ egin{cases} 
\lambda_k = \mu - \frac{\log(b)}{a_0} + i \frac{2k\pi}{a_0}, \\ 
\phi_k(a) = \exp\left(\left(\mu - \frac{\log(b)}{a_0} + i \frac{2k\pi}{a_0}\right)a\right), 
\end{cases} $$

so that we could say that the eigenvectors form a Fourier basis. The evolution problem has periodic solutions, up to a multiplicative growing or decaying exponential term. Interestingly, this situation is not at all a pure mathematical curiosity, as such reproduction strategies are employed by certain species, like some bamboo species that have a cycle of decades.

## Discussion

In this paper, we present theoretical and numerical studies on a semi-discretized McKendrick equation eigenvalues. This work aims at enlightening a specific phenomenon occurring with the distribution of the eigenvalues of the McKendrick semi-discretized equation. Although we managed to demonstrate the numerical conjectures, the theoretical results do not exhaust the issue. One could seek a deeper explanation on why such figures arise when the number of age groups goes to infinity. Moreover, what influence does the numerical scheme or the equation itself have on the discrete spectrum? How could one link the spectra of discrete and continuous McKendrick equation? What extension could one draw with the predator-prey model, or even other population models?

We acknowledge that studying the simple linear McKendrick equations does not exhaust the discussion on nonlinear systems more generally. We started by studying the predator-prey model and got
interesting and, to say the least, surprising numerical conjectures. We did not seek rigorous analysis of such conjectures, and focused on a more accessible problem to develop theoretical methods. The question of extending the proof of the circular repartition of eigenvalues to more complex and nonlinear population models stays open and could be studied in future works.

The above-mentioned theoretical result enlightens the importance of choosing wisely the number of age classes in or semi-discrete model. Can one trust a model with a few age classes? Is there a risk that the model exhibits wrong stability properties in this situation? Instead of working all along the study with a fixed number of age classes, we propose formal and theoretical tools to apprehend the impact of such a modeling choice. This issue is well-known from PDE experts regularly dealing with mesh convergence.

In this paper, we focused only on the upwind scheme to discretize the age derivative. Other schemes would likely lead to comparable results. Such study could be the topic of further research. This circular or disk shape of spectrum seems to be ubiquitous, particularly in the study of random polynomial roots [46, 47].

## Future work

As presented in the previous discussion, many interesting topics arise from our study. The generalization of such theoretical results to an age-structured Lotka-Volterra model, and in a larger way, to this type of equations (i.e. transport equations with exchange terms) is one of these topics. To achieve that, one could linearize the age-structured Lotka-Volterra model at an equilibrium point, to get the spectrum of the continuous problem and compare it with the discrete one. Other numerical schemes could produce similar theoretical and numerical results. Understanding better the computation and characterization of equilibrium points, and handling the nonlinearities of the predator-prey model which were eclipsed in the McKendrick equation, represents a crucial issue. It would also be a topic of interest to deepen the interpretation and analysis of the eigenvalues distribution; especially, explain why circles appear with the discrete formulation and not with the continuous one, and how one can properly get the second spectrum from the first by taking the limit on the number of age groups.

## Conclusion

In this work, we presented numerical and theoretical considerations on spectral properties of population dynamics models. We took simple and academic population dynamics models with an age structure, starting from predator-prey models extended within a McKendrick formalism, then undergoing a semi-discretization process. We proved a density result on the discrete formulation's eigenvalues, as the number of age groups goes to infinity. We showed the property that the eigenfunctions of McKendrick spectral problem form a complete basis of the space $L^2(]0, A[)$.
These considerations enlighten the stability analysis of such population models, where the eigenvalues distribution play a fundamental role.
# Acknowledgments

This work has been carried out under a collaboration between Sorbonne Université and Dassault Systèmes.

# Conflict of interest

* This work was funded by Dassault Systèmes and the Association Nationale de la Recherche et de la Technologie (ANRT).

* The authors certify that they do not have any conflict of interest.

* The authors have no financial or non-financial interests to disclose.

* The authors provide their full consent for publication.

# References

1. A.G McKendrick, Applications of Mathematics to Medical Problems, *Proceedings of the Edinburgh Mathematical Society*, **44** (1925), 98–130, https://doi.org/10.1017/S0013091500034428.

2. W.O Kermack and A.G McKendrick, A Contribution to the Mathematical Theory of Epidemics, *Proceedings of the Royal Society of London.*, **115** (1927), 700–721, https://doi.org/10.1098/rspa.1927.0118.

3. F Hoppenstaedt, *Mathematical Theories of Populations: Demographics, Genetics and Epidemics*, no. 20 in CBMS -NSF Regional Conference Series in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1975, https://doi.org/10.1137/1.9781611970487.

4. G.F Webb, *Theory of nonlinear age-dependent population dynamics*, no. 89 in Monographs and textbooks in pure and applied mathematics, M. Dekker, New York, 1985.

5. F Sharpe and A Lotka, A problem in age-distribution, *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science*, **21** (1911), 435–438.

6. W.R Thompson, On the Reproduction of Organisms with overlapping Generations, *Bulletin of Entomological Research*, **22** (1931), 147–172, https://doi.org/10.1017/S0007485300029837.

7. L.C Cole, The Population consequences of life history phenomena, *The Quarterly Review of Biology*, **29** (1954), 103–137.

8. H von Foerster, Some Remarks on Changing Populations, in *The Kinetics of Cellular Proliferation* (ed. J.F Stohlman). Grune and Stratton, New York, 1959, 382–407.

9. E Cannan, The Probability of a Cessation of the Growth of Population in England and Wales During the Next Century, *The Economic Journal*, **5** (1895), 505, https://doi.org/10.2307/2956626.

10. P.H Leslie, On the use of matrices in certain population mathematics, *Biometrika*, **33** (1945), 183–212, https://doi.org/10.1093/biomet/33.3.183.

11. B.L Keyfitz and N Keyfitz, The McKendrick partial differential equation and its uses in epidemiology and population study, *Mathematical and Computer Modelling*, **26**, https://doi.org/10.1016/S0895-7177(97)00165-9.
12. M.E Gurtin and R.C Maccamy, Non-linear age-dependent population dynamics, *Archive for Rational Mechanics and Analysis*, **54** (1974), 281–300, https://doi.org/10.1007/BF00250793.

13. M Iannelli, *Mathematical theory of age-structured population dynamics*, no. 7 in Applied mathematics monographs, Giardini editori e stampatori, Pisa, 1995.

14. O Diekmann and P Getto, Boundedness, global existence and continuous dependence for nonlinear dynamical systems describing physiologically structured populations, *Journal of Differential Equations*, **215** (2005), 268–319, https://doi.org/http://doi.org/10.1016/j.jde.2004.10.025.

15. B Perthame, *Transport Equations in Biology*, Frontiers in Mathematics, Birkhäuser Basel, Basel, 2007, https://doi.org/10.1007/978-3-7643-7842-4.

16. O Diekmann, M Gyllenberg, J.A.J Metz and H.R Thieme, The "cumulative" formulation of (physiologically) structured population models, in *Evolution equations, control theory, and biomathematics: proceedings of the Han-sur-Lesse conference* (eds. P Clément and G Lumer). no. 155 in Lecture notes in pure and applied mathematics, Marcel Dekker, New York, 1994.

17. O Diekmann, M Gyllenberg, J.A.J Metz and H.R Thieme, On the formulation and analysis of general deterministic structured population models I. Linear Theory, *Journal of Mathematical Biology*, **36** (1998), 349–388, https://doi.org/10.1007/s002850050104.

18. L Attias, V Siess and S Labbé, An agile modeling framework for population dynamics, *Mathematics and Computers in Simulation*, **234** (2025), 113–134, https://doi.org/10.1016/j.matcom.2025.02.013.

19. D Brück, H Elmqvist, H Olsson and S Mattsson, Dymola for Multi-Engineering Modeling and Simulation, in *Proceedings of the 2nd International Modelica Conference* (ed. M Otter), 2002, 55–1–55–8.

20. P Fritzson and V Engelson, Modelica — A unified object-oriented language for system modeling and simulation, in *ECOOP’98 — Object-Oriented Programming* (ed. E Jul). Springer Berlin Heidelberg, Berlin, Heidelberg, 1998, 67–90.

21. A Ronacher, *Jinja2 Documentation*, Technical report, 2014, URL http://jinja.pocoo.org.

22. D.G Magiros, Stability concepts of dynamical systems, *Information and Control*, **9** (1966), 531–548, https://doi.org/10.1016/S0019-9958(66)80014-1.

23. R Bellman, *Stability theory of differential equations*, Dover publications, inc. edition, Dover Publications, Mineola, N.Y., 1953.

24. J.P LaSalle, *The Stability of Dynamical Systems*, no. 25 in CBMS -NSF Regional Conference Series in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1976, https://doi.org/10.1137/1.9781611970432.

25. J.P Demailly, *Analyse numérique et équations différentielles*, Collection Grenoble Sciences, EDP Sciences, Les Ulis, France, 2006.

26. A.M Lyapunov, *The general problem of the stability of motion*, Taylor & Francis, London ; Washington, DC, 1992.

27. D.G Magiros, On a Class of Precessional Phenomena and Their Stability in the Sense of Liapunov, Poincaré and Lagrange, in *Selected Papers of Demetrios G. Magiros* (ed. S.G Tzafestas). Springer Netherlands, Dordrecht, 1985, 432–439, https://doi.org/10.1007/978-94-009-5368-0_36.
28. D.G Magiros, On stability definitions of dynamical systems, *Proceedings of the National Academy of Sciences*, **53** (1965), 1288-1294, https://doi.org/10.1073/pnas.53.6.1288.

29. E.P Gyftopoulos, Lagrange Stability and Liapunov's Direct Method, in *Proceedings of Symposium on Reactor Kinetics and Control*, 1963.

30. N.P Bhatia and G.P Szegö, *Stability theory of dynamical systems*, Classics in mathematics, Springer, Berlin ; New York, 1970.

31. M.R Roussel, Stability Analysis for ODEs, in *Nonlinear Dynamics: A hands-on introductory survey*. 2053-2571, Morgan & Claypool Publishers, San Rafael, CA, USA, 2019, 3–1 to 3–12.

32. V.I Arnol'd, *Ordinary differential equations*, Springer, Berlin ; New York, 1992.

33. A.J Lotka, Elements of Physical Biology, *Science Progress in the Twentieth Century* (1919-1933), **21** (1926), 341–343.

34. V Volterra, Fluctuations in the Abundance of a Species considered Mathematically, *Nature*, **118** (1926), 558–560, https://doi.org/10.1038/118558a0.

35. P.J Wangersky, Lotka-Volterra Population Models, *Annual Review of Ecology and Systematics*, **9** (1978), 189–218, https://doi.org/10.1146/annurev.es.09.110178.001201.

36. Y Takeuchi, *Global Dynamical Properties of Lotka-Volterra Systems*, World Scientific, Shizuoka University, Japan, 1996, https://doi.org/10.1142/2942.

37. R Cherniha and V Davydovych, Construction and application of exact solutions of the diffusive Lotka–Volterra system: A review and new results, *Communications in Nonlinear Science and Numerical Simulation*, **113**, https://doi.org/10.1016/j.cnsns.2022.106579.

38. Y Lu, K.A Pawelek and S Liu, A stage-structured predator-prey model with preda-tion over juvenile prey, *Applied Mathematics and Computation*, **297** (2017), 115–130, https://doi.org/10.1016/j.amc.2016.10.035.

39. X Zhang and Z Liu, Hopf bifurcation analysis in a predator-prey model with predator-age structure and predator-prey reaction time delay, *Applied Mathematical Modelling*, **91** (2021), 530–548, https://doi.org/10.1016/j.apm.2020.08.054.

40. S Bentout, S Djilali and A Atangana, Bifurcation analysis of an age-structured prey–predator model with infection developed in prey, *Mathematical Methods in the Applied Sciences*, **45** (2022), 1189–1208, https://doi.org/10.1002/mma.7846.

41. D.H Meadows, D.L Meadows, J Randers and W.W Behrens III, *The Limits to Growth: A Report for the Club of Rome's Project on the Predicament of Mankind*, Universe Books, New York, 1972, URL https://www.library.dartmouth.edu/digital/digital-collections/limits-growth.

42. W.M Boothby, On Two Classical Theorems of Algebraic Topology, *The American Mathematical Monthly*, **78** (1971), 237–249, https://doi.org/10.1080/00029890.1971.11992736.

43. J Song, Y Jing Yu, Z Xian Wang, S.J Hu, Z.X Zhao, J.Q Liu, X Feng De and G.T Zhu, Spectral properties of population operator and asymptotic behaviour of population semigroup, *Acta Mathematica Scientia*, **2** (1982), 139–148, https://doi.org/10.1016/S0252-9602(18)30629-5.

44. T Kato, *Perturbation Theory for Linear Operators*, vol. 132 of Classics in Mathematics, Springer Berlin Heidelberg, Berlin, Heidelberg, 1995, https://doi.org/10.1007/978-3-642-66282-9.
45. P.M Riechers and J.P Crutchfield, Beyond the spectral theorem: Spectrally decomposing arbitrary functions of nondiagonalizable operators, *AIP Advances*, **8**, https://doi.org/10.1063/1.5040705.

46. P Erdos and P Turan, On the Distribution of Roots of Polynomials, *The Annals of Mathematics*, **51** (1950), 105, https://doi.org/10.2307/1969500.

47. G.A Mezincescu, D Bessis, J.D Fournier, G Mantica and F.D Aaron, Distribution of roots of random real generalized polynomials, *Journal of Statistical Physics*, **86** (1997), 675–705, https://doi.org/10.1007/BF02199115.

<figure><img src="image_17.png" /></figure>

AIMS Press

©2025 Laurent Attias, Stéphane Labbé, Vincent Siess, licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0)