[datatype=bibtex] [fieldsource=doi,final] [fieldset=url,null] *A

A pseudospectral method for investigating the stability of linear population
models with two physiological structures

Alessia Andò111Area of Mathematics, Gran Sasso Science Institute, Viale F.
Crispi 7, 67100 L’Aquila, Italy444CDLab – Computational Dynamics Laboratory,
University of Udine<EMAIL_ADDRESS>Simone De
Reggi222Department of Mathematics, Computer Science and Physics, University of
Udine, Via delle Scienze 206, 33100 Udine, Italy444CDLab – Computational
Dynamics Laboratory, University of Udine,
<EMAIL_ADDRESS>
Davide Liessi222Department of Mathematics, Computer Science and Physics,
University of Udine, Via delle Scienze 206, 33100 Udine, Italy444CDLab –
Computational Dynamics Laboratory, University of Udine,
<EMAIL_ADDRESS>Francesca Scarabel333Department of
Mathematics, The University of Manchester, Oxford Rd, M13 9PL, Manchester,
United Kingdom444CDLab – Computational Dynamics Laboratory, University of
Udine, Italy555Joint UNIversities Pandemic and Epidemiological Research,
United<EMAIL_ADDRESS>
March 8, 2022

###### Abstract

The asymptotic stability of the null equilibrium of a linear population model
with two physiological structures formulated as a first-order hyperbolic PDE
is determined by the spectrum of its infinitesimal generator. We propose an
equivalent reformulation of the problem in the space of absolutely continuous
functions in the sense of Carathéodory, so that the domain of the
corresponding infinitesimal generator is defined by trivial boundary
conditions. Via bivariate collocation, we discretize the reformulated operator
as a finite-dimensional matrix, which can be used to approximate the spectrum
of the original infinitesimal generator. Finally, we provide test examples
illustrating the converging behavior of the approximated eigenvalues and
eigenfunctions, and its dependence on the regularity of the model
coefficients.

Keywords: bivariate collocation, infinitesimal generator, partial differential
equations, stability of equilibria, physiologically structured populations

2020 Mathematics Subject Classification: Primary: 37M99, 65L07, 65N25;
Secondary: 35L04, 37N25, 47D06, 92D25.

##  1 Introduction

Mathematical models are often used to describe the evolution of populations in
biology and epidemiology. An important class of models that has attracted
increased attention is that of structured population models, in which
individuals are characterized by one or more variables that describe the
i-state (i.e., the individual state) and determine the individual processes,
including for instance birth, growth and death. Examples of physiological
structures are age, size, spatial position or time since infection. We here
focus on a class of structured population models where the structuring
variables are continuous and the models are formulated as first-order
hyperbolic partial differential equations (PDEs) (see, e.g., [2, 23, 28, 32,
33] and references therein). In particular, we consider the case where
individuals are characterized by two different traits [31, 43].

In applications to population dynamics, the interest is often focused on the
long-term properties of the systems, for instance the existence of equilibrium
states and their stability. For linear models with two structures, it has been
proved that the stability of the zero equilibrium is determined by the
spectrum of the infinitesimal generator (IG) of the semigroup of solution
operators [31, 43]. Since the IG is an operator acting on an infinite-
dimensional space of functions, numerical techniques are required to obtain
finite-dimensional approximations of the operator and, in turn, of its
spectrum.

For the analysis of local stability of equilibria, pseudospectral methods have
been widely used both for delay equations [6, 7, 9, 14] and for PDE population
models with one structuring variable [4, 10, 34]. The main advantage of
pseudospectral methods is their typical spectral accuracy, by which the order
of convergence of the approximation error increases with the regularity of the
approximated function. In particular, the convergence is exponential for
analytic functions [42]. In the case of delay equations, this implies that the
spectrum of the IG is approximated with exponential order of convergence, as
the corresponding eigenfunctions are exponentials (see, e.g., [16, Proposition
3.4]). In the case of PDEs, the eigenfunctions are still exponential in time,
but the order of regularity with respect to the physiological variables
depends on the regularity of the model parameters, which therefore affects the
order of convergence of the approximation [43, Theorem 3.2]. A similar
behavior has been shown for the approximation of $R_{0}$ for structured
epidemic models [5, 8, 11], and in the approximation of the solution operators
[12, 13, 15].

For structured models with one single structuring variable, pseudospectral
methods have already been proposed to study the stability of the zero
equilibrium in [4]. In that paper the IG is approximated by combining
pseudospectral differentiation with the inversion of a (linear) algebraic
condition characterizing the domain of the operator. However, with more
structuring variables implementing this technique becomes substantially more
involved.

A different approach, which has been successfully employed in the context of
nonlinear PDEs with one physiological variable [34] and of renewal equations
[35, 36], consists in first reformulating the problem at hand via conjugation
with an integral operator, and then approximating the resulting transformed
operator via pseudospectral (also known as spectral collocation) techniques.
The advantages of this approach are mainly twofold: on the one hand, the
transformed operator acts on a space of absolutely continuous functions
(rather than the original space $L^{1}$), hence point evaluation, as well as
polynomial interpolation and collocation, are well defined; on the other hand,
the domain of the transformed operator is characterized by a trivial condition
(specifically, a zero boundary condition), which substantially simplifies the
numerical implementation. From a modeling point of view, the integrated state
has a clear interpretation as it represents the number of individuals whose
i-state is less than a given value.

Goal of this work is to introduce a numerical method for the stability
analysis of linear PDE population models with two structuring variables based
on this second approach. As far as we know, there is currently no other
numerical method available for this problem. We demonstrate the applicability
and computational efficiency of the new method with several numerical tests,
illustrating the convergence of the approximated eigenvalues to the exact ones
and supporting the conjecture of spectral accuracy.

The paper is organized as follows. In section 2 we introduce the prototype
model and the relevant solution operators and IG. Section 3 describes the
reformulation of the IG in terms of the integrated state. The reformulated IG
is discretized in section 4 and the resulting numerical method is applied to
some test problems in section 5. In section 6 we discuss the extension to
models with structuring variables evolving with nontrivial velocities and show
an example. Finally, we provide some concluding remarks in section 7. In
appendix A, for completeness, we apply the method to the case of one
structuring variable and illustrate it with a test model.

A MATLAB implementation of the method is available at
http://cdlab.uniud.it/software.

##  2 The prototype model

Let $x_{0},\bar{x},y_{0},\bar{y}\in\mathbb{R}$ such that $x_{0}<\bar{x}$ and
$y_{0}<\bar{y}$ and let
$\overline{R}\coloneqq[x_{0},\bar{x}]\times[y_{0},\bar{y}]$. We consider the
scalar first-order linear hyperbolic PDE

$\partial_{t}u(t,x,y)+\partial_{x}u(t,x,y)+\partial_{y}u(t,x,y)=-\mu(x,y)u(t,x,y),$
(2.1)

with boundary conditions

$\displaystyle u(t,x,y_{0})$
$\displaystyle=\iint_{\overline{R}}\alpha(x,\xi,\sigma)u(t,\xi,\sigma)\mathop{}\\!\mathrm{d}\xi\mathop{}\\!\mathrm{d}\sigma\eqqcolon
K_{\alpha}(u(t,\cdot,\cdot))(x),$ (2.2) $\displaystyle u(t,x_{0},y)$
$\displaystyle=\iint_{\overline{R}}\beta(y,\xi,\sigma)u(t,\xi,\sigma)\mathop{}\\!\mathrm{d}\xi\mathop{}\\!\mathrm{d}\sigma\eqqcolon
K_{\beta}(u(t,\cdot,\cdot))(y),$ (2.3)

where $u(t,x,y)$ is the density of the given population at time $t\geq 0$
depending on the two structuring variables $x\in[x_{0},\bar{x}]$ and
$y\in[y_{0},\bar{y}]$.

Following [31, Assumption 2.1], we assume that the model coefficients
$\mu\colon\overline{R}\to\mathbb{R}$,
$\alpha\colon[x_{0},\bar{x}]\times\overline{R}\to\mathbb{R}$ and
$\beta\colon[y_{0},\bar{y}]\times\overline{R}\to\mathbb{R}$ are nonnegative
$L^{1}$ functions, Lipschitz continuous on the interior of their domains, with
$\alpha(x,\cdot,\cdot)$ and $\beta(y,\cdot,\cdot)$ dominated by $L^{1}$
functions in $x$ and $y$, respectively, and $\mu$ bounded from above and
bounded away from $0$ (for a similar but more general approach, see [38,
section 4]). Observe that $\alpha(x,\cdot,\cdot)$ and $\beta(y,\cdot,\cdot)$
are essentially bounded, so the operators $K_{\alpha}$ and $K_{\beta}$ map
$L^{1}(\overline{R})$ to $L^{1}([x_{0},\bar{x}])$ and
$L^{1}([y_{0},\bar{y}])$, respectively, and are bounded.

Under these assumptions, for every $u_{0}\in L^{1}(\overline{R})$ the
initial–boundary value problem defined by (2.1)–(2.3) and
$u(0,\cdot,\cdot)=u_{0}$ admits a unique solution $u(t,\cdot,\cdot)\in
L^{1}(\overline{R})$ for $t\geq 0$. Moreover, the family of solution operators
$\\{T(t)\\}_{t\geq 0}$, defined by $T(t)u_{0}=u(t,\cdot,\cdot)$, forms a
strongly continuous semigroup of bounded linear operators in the Banach space
$L^{1}(\overline{R})$, see [43, Theorem 3.1] or [31, Theorem 2.3]. In
addition, $\\{T(t)\\}_{t\geq 0}$ is eventually compact [31, section 4].

The IG of the semigroup $\\{T(t)\\}_{t\geq 0}$ is the operator $A\colon
D(A)\to L^{1}(\overline{R})$ defined by

$A\phi=\lim_{t\to 0^{+}}\frac{1}{t}(T(t)\phi-\phi),$ (2.4)

where $D(A)\subset L^{1}(\overline{R})$ consists of the functions for which
the limit exists. Kang et al. [31, Remark 2.3] prove that the operator $A$
satisfies

$\displaystyle A\phi(x,y)$
$\displaystyle=-\partial_{x}\left[\phi(x,y)+\partial_{y}\int_{x_{0}}^{x}\phi(a,y)\mathop{}\\!\mathrm{d}a\right]-\mu(x,y)\phi(x,y)$
(2.5)
$\displaystyle=-\partial_{y}\left[\phi(x,y)+\partial_{x}\int_{y_{0}}^{y}\phi(x,a)\mathop{}\\!\mathrm{d}a\right]-\mu(x,y)\phi(x,y),$

for a.e. $x\in[x_{0},\bar{x}]$, a.e. $y\in[y_{0},\bar{y}]$, and every $\phi\in
D(A)$, and that $D(A)$ satisfies the inclusion

$\displaystyle D(A)\subset\biggl{\\{}$ $\displaystyle\phi\in
L^{1}(\overline{R})\;\bigg{|}$
$\displaystyle(x,y)\mapsto\int_{x_{0}}^{x}\phi(s,y)\mathop{}\\!\mathrm{d}s\text{
is absolutely continuous in $y$, for a.e. $x\in[x_{0},\bar{x}]$},$
$\displaystyle(x,y)\mapsto\left[\phi(x,y)+\partial_{y}\int_{x_{0}}^{x}\phi(s,y)\mathop{}\\!\mathrm{d}s\right]\text{
is absolutely continuous in $x$,}$
$\displaystyle\hphantom{(x,y)\mapsto\left[\phi(x,y)+\partial_{y}\int_{x_{0}}^{x}\phi(s,y)\mathop{}\\!\mathrm{d}s\right]{}}\text{
for a.e. $y\in[y_{0},\bar{y}]$},$ $\displaystyle\lim_{x\to
x_{0}^{+}}\left[\phi(x,y)+\partial_{y}\int_{x_{0}}^{x}\phi(s,y)\mathop{}\\!\mathrm{d}s\right]=K_{\beta}(\phi)(y)\text{
for a.e. $y\in[y_{0},\bar{y}]$},$
$\displaystyle(x,y)\mapsto\int_{y_{0}}^{y}\phi(x,s)\mathop{}\\!\mathrm{d}s\text{
is absolutely continuous in $x$, for a.e. $y\in[y_{0},\bar{y}]$},$
$\displaystyle(x,y)\mapsto\left[\phi(x,y)+\partial_{x}\int_{y_{0}}^{y}\phi(x,s)\mathop{}\\!\mathrm{d}s\right]\text{
is absolutely continuous in $y$,}$
$\displaystyle\hphantom{(x,y)\mapsto\left[\phi(x,y)+\partial_{x}\int_{y_{0}}^{y}\phi(x,s)\mathop{}\\!\mathrm{d}s\right]{}}\text{
for a.e. $x\in[x_{0},\bar{x}]$},$ $\displaystyle\lim_{y\to
y_{0}^{+}}\left[\phi(x,y)+\partial_{x}\int_{y_{0}}^{y}\phi(x,s)\mathop{}\\!\mathrm{d}s\right]=K_{\alpha}(\phi)(x)\text{
for a.e. $x\in[x_{0},\bar{x}]$},$
$\displaystyle\partial_{x}\left[\phi(x,y)+\partial_{y}\int_{x_{0}}^{x}\phi(s,y)\mathop{}\\!\mathrm{d}s\right]\in
L^{1}(\overline{R}),$
$\displaystyle\partial_{y}\left[\phi(x,y)+\partial_{x}\int_{y_{0}}^{y}\phi(x,s)\mathop{}\\!\mathrm{d}s\right]\in
L^{1}(\overline{R})\biggr{\\}},$

while [43, Remark 6.1] claims that equality holds. If $\phi\in D(A)$ is
sufficiently smooth, the action (2.5) of the operator $A$ simplifies and can
be expressed as

$A\phi=-(\partial_{x}\phi+\partial_{y}\phi+\mu\phi).$

The spectrum of the IG determines the stability of the zero
equilibrium.101010Note that since $A$ is an operator on a real Banach space,
in order to define and compute its spectrum the space and the operator need to
be complexified. For details see, e.g., [22, section III.7]. More precisely,
the latter is asymptotically stable if and only if the spectral abscissa of
$A$ is negative and it is unstable if the spectral abscissa is positive (see
[18, Theorem 9.5] and [24, Theorem VI.1.15]).

Observe that in the model defined by (2.1)–(2.3) the structuring variables
evolve with the same velocity as time. In sections 3 and 4 we present the
integral reformulation and the discretization restricting to this special
case, as in [31]. However, they can be applied to the case of nontrivial
velocities as described in section 6, where we also present an example.

##  3 Equivalent formulation in a space of absolutely continuous functions

To conveniently handle the boundary conditions in $D(A)$ from a numerical
point of view, inspired by the approach of [34] in the case of one structuring
variable, we argue in terms of the integrated state. In particular, we define
an isomorphism between $L^{1}(\overline{R})$ and a suitable space of functions
via integration. We then use this isomorphism and the semigroup
$\\{T(t)\\}_{t\geq 0}$ to construct an appropriate semigroup acting on a space
of functions with higher regularity.

With this goal in mind, we first recall the definition and some properties of
absolute continuity in the sense of Carathéodory. We refer the reader to [39]
for further details.

For $(x,y)\in\overline{R}$, define $R(x,y)\coloneqq[x_{0},x]\times[y_{0},y]$;
then $\overline{R}=R(\bar{x},\bar{y})$. A function $v$ defined on
$\overline{R}$ is _absolutely continuous in the sense of Carathéodory_ if and
only if there exist $e_{v}\in\mathbb{R}$, $f_{v}\in L^{1}([x_{0},\bar{x}])$,
$g_{v}\in L^{1}([y_{0},\bar{y}])$ and $h_{v}\in L^{1}(\overline{R})$ such that

$v(x,y)=e_{v}+\int_{x_{0}}^{x}f_{v}(a)\mathop{}\\!\mathrm{d}a+\int_{y_{0}}^{y}g_{v}(b)\mathop{}\\!\mathrm{d}b+\iint_{R(x,y)}h_{v}(a,b)\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b.$

Observe that the double integral in the last term is equal to the iterated
integral on the variables $a$ and $b$ in any order, thanks to Fubini’s
theorem. The space $AC(\overline{R})$ of absolutely continuous functions on
$\overline{R}$ in the sense of Carathéodory is a Banach space when equipped
with the norm $\lVert\cdot\rVert_{AC(\overline{R})}$ defined as

$\lVert v\rVert_{AC(\overline{R})}\coloneqq\lvert e_{v}\rvert+\lVert
f_{v}\rVert_{L^{1}([x_{0},\bar{x}])}+\lVert
g_{v}\rVert_{L^{1}([y_{0},\bar{y}])}+\lVert
h_{v}\rVert_{L^{1}(\overline{R})}.$

We consider a particular subspace of $AC(\overline{R})$, namely

$AC_{0}(\overline{R}):=\biggl{\\{}v\colon\overline{R}\to\mathbb{R}\;\bigg{|}\;v(x,y)=\iint_{R(x,y)}h_{v}(a,b)\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b\
\text{ for some }h_{v}\in L^{1}(\overline{R})\biggr{\\}}.$

Observe that $AC_{0}(\overline{R})$ is a Banach space, being a closed
subspace, and that $v(x_{0},y)=v(x,y_{0})=0$ for $v\in AC_{0}(\overline{R})$.
The operator $V\colon L^{1}(\overline{R})\to AC_{0}(\overline{R})$ defined by

$V\phi(x,y)=\iint_{R(x,y)}\phi(a,b)\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b$

defines an isomorphism between $L^{1}(\overline{R})$ and
$AC_{0}(\overline{R})$, with $V^{-1}\psi=\partial_{x}\partial_{y}\psi$ for all
$\psi\in AC_{0}(\overline{R})$. Observe that both $V$ and $V^{-1}$ are bounded
($\lVert V\rVert=\lVert V^{-1}\rVert=1$).

Note that, given a solution $u$ of (2.1)–(2.3),
$\iint_{R(x,y)}u(t,a,b)\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b$
represents the number of individuals whose structuring variables belong to
$R(x,y)$ at time $t$.

Returning now to (2.1)–(2.3), we define the family of operators
$\\{S(t)\\}_{t\geq 0}$ on $AC_{0}(\overline{R})$ as $S(t)\coloneqq
VT(t)V^{-1}$. Since $V$ and $V^{-1}$ are linear and bounded, recalling that
the composition of a compact and a bounded operator (in either order) is
compact, the operators $S(t)$ are in turn linear and bounded and form a family
with the same properties as $\\{T(t)\\}_{t\geq 0}$, namely they form a
strongly continuous and eventually compact semigroup on
$AC_{0}(\overline{R})$. Its IG is $B\colon D(B)\to AC_{0}(\overline{R})$, with
$B\coloneqq VAV^{-1}$ and $D(B)=VD(A)$.

With the aim of using $B$ to study the stability properties of (2.1)–(2.3), it
is important to understand the relation between the spectra of $A$ and $B$. By
[24, Corollary V.3.2(i)], since the semigroups $\\{T(t)\\}_{t\geq 0}$ and
$\\{S(t)\\}_{t\geq 0}$ are eventually compact, the spectra of both $A$ and $B$
are at most countable and consist only of eigenvalues (of finite algebraic
multiplicity). By [15, Proposition 4.1] $A$ and $B$ have the same nonzero
eigenvalues (with the same multiplicities) and corresponding eigenvectors.
Moreover, observe that $A$ is injective if and only if $B$ is, so $0$ is
either in both spectra or in none. Therefore, the spectra of $A$ and $B$
coincide. From [31, Proposition 3.1], noting that in that proof
$\chi=A\phi+\mu\phi$, for $\psi\in D(B)$ we can write

$\displaystyle B\psi(x,y)=$
$\displaystyle-\partial_{x}\psi(x,y)-\partial_{y}\psi(x,y)$ (3.1)
$\displaystyle+\int_{x_{0}}^{x}K_{\alpha}(\partial_{y}\partial_{x}\psi)(a)\mathop{}\\!\mathrm{d}a+\int_{y_{0}}^{y}K_{\beta}(\partial_{y}\partial_{x}\psi)(b)\mathop{}\\!\mathrm{d}b$
$\displaystyle-\iint_{R(x,y)}\mu(a,b)\partial_{y}\partial_{x}\psi(a,b)\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b.$

##  4 Pseudospectral discretization of the IG

In this section, we use pseudospectral methods with a tensorial approach to
obtain a finite-dimensional approximation of the operator $B$, whose spectrum
can be used to determine the stability properties of the system.

For $n$ and $m$ positive integers, let $\Pi_{n,m}^{0}$ be the space of
bivariate polynomials on $\overline{R}$ of degree at most $n$ in the first
variable and at most $m$ in the second variable, taking value $0$ at $x=x_{0}$
and $y=y_{0}$. These conditions are motivated by the fact that we use
polynomials in $\Pi_{n,m}^{0}$ as approximations of functions in
$AC_{0}(\overline{R})$. Let $\Theta_{X}=\\{x_{1},\ldots,x_{n}\\}$ be a mesh of
$n$ points in $(x_{0},\bar{x}]$, with $x_{0}<x_{1}<\dots<x_{n}=\bar{x}$, and
let $\Theta_{Y}=\\{y_{1},\ldots,y_{m}\\}$ be a mesh of $m$ points in
$(y_{0},\bar{y}]$, with $y_{0}<y_{1}<\dots<y_{m}=\bar{y}$. We approximate a
function $\psi\in AC_{0}(\overline{R})$ by a vector $\Psi\in\mathbb{R}^{nm}$
according to

$\psi(x_{i},y_{j})=\Psi_{i,j},\quad i=1,\ldots,n,\;j=1,\ldots,m,$

where the components of $\Psi$ are ordered according to the lexicographic
order of the double indices $(i,j)$.

Given $\Psi\in\mathbb{R}^{nm}$, let $\psi_{n,m}\in\Pi_{n,m}^{0}$ be the
polynomial interpolating $\Psi$ on $\Theta_{X}\times\Theta_{Y}$:

$\displaystyle\psi_{n,m}(x_{i},y_{j})=\Psi_{i,j},\quad
i=1,\ldots,n,\;j=1,\ldots,m.$

The finite-dimensional approximation of the operator $B$ is then
$B_{n,m}\colon\mathbb{R}^{nm}\to\mathbb{R}^{nm}$ defined as

$[B_{n,m}\Psi]_{i,j}\coloneqq(B\psi_{n,m})(x_{i},y_{j}),\quad
i=1,\ldots,n,\;j=1,\ldots,m.$ (4.1)

We can write more explicitly the entries of the matrix $B_{n,m}$ by using the
bivariate Lagrange representation of $\psi_{n,m}$, together with the explicit
action of the operator $B$ defined in (3.1). Let
$\\{\ell_{X,i}\\}_{i=0,\ldots,n}$ and $\\{\ell_{Y,j}\\}_{j=0,\ldots,m}$ be the
Lagrange bases of polynomials relevant to $\\{x_{0}\\}\cup\Theta_{X}$ and
$\\{y_{0}\\}\cup\Theta_{Y}$, i.e.,

$\displaystyle\ell_{X,i}(x)=\prod_{\begin{subarray}{c}k=0\\\ k\neq
i\end{subarray}}^{n}\frac{x-x_{k}}{x_{i}-x_{k}},\quad\ell_{Y,j}(y)=\prod_{\begin{subarray}{c}k=0\\\
k\neq j\end{subarray}}^{m}\frac{y-y_{k}}{y_{j}-y_{k}}.$

The polynomial $\psi_{n,m}$ can be written as

$\psi_{n,m}(x,y)=\sum_{i=1}^{n}\sum_{j=1}^{m}\ell_{X,i}(x)\ell_{Y,j}(y)\Psi_{i,j},\quad(x,y)\in\overline{R}.$

Note that indeed $\psi_{n,m}(x,y)=0$ for $x=x_{0}$ or $y=y_{0}$. Using (3.1)
and (4.1), we get

$\displaystyle{}_{k,l}=$
$\displaystyle-\sum_{i=1}^{n}\ell_{X,i}^{\prime}(x_{k})\Psi_{i,l}-\sum_{j=1}^{m}\ell_{Y,j}^{\prime}(y_{l})\Psi_{k,j}$
$\displaystyle+\int_{y_{0}}^{y_{l}}K_{\beta}\biggl{(}\sum_{i=1}^{n}\sum_{j=1}^{m}\ell_{X,i}^{\prime}\ell_{Y,j}^{\prime}\Psi_{i,j}\biggr{)}(b)\mathop{}\\!\mathrm{d}b$
$\displaystyle+\int_{x_{0}}^{x_{k}}K_{\alpha}\biggl{(}\sum_{i=1}^{n}\sum_{j=1}^{m}\ell_{X,i}^{\prime}\ell_{Y,j}^{\prime}\Psi_{i,j}\biggr{)}(a)\mathop{}\\!\mathrm{d}a$
$\displaystyle-\iint_{R(x_{k},y_{l})}\mu(a,b)\sum_{i=1}^{n}\sum_{j=1}^{m}\ell_{X,i}^{\prime}(a)\ell_{Y,j}^{\prime}(b)\Psi_{i,j}\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b.$

Using the linearity of $K_{\alpha}$ and $K_{\beta}$, it is easy to
characterize the entries of the matrix $B_{n,m}$. Let
$D_{X}\in\mathbb{R}^{n\times n}$ and $D_{Y}\in\mathbb{R}^{m\times m}$ be
defined as

$\displaystyle[D_{X}]_{i,j}$ $\displaystyle=\ell_{X,j}^{\prime}(x_{i}),$
$\displaystyle\quad i,j$ $\displaystyle=1,\ldots,n$
$\displaystyle[D_{Y}]_{i,j}$ $\displaystyle=\ell_{Y,j}^{\prime}(y_{i}),$
$\displaystyle\quad i,j$ $\displaystyle=1,\ldots,m.$

In other words, $D_{X}$ and $D_{Y}$ are the part of the differentiation
matrices associated with $\\{x_{0}\\}\cup\Theta_{X}$ and
$\\{y_{0}\\}\cup\Theta_{Y}$, respectively, deleting the first row and the
first column. The bivariate differentiation matrices in $x$ and $y$ are

$\mathbf{D}_{X}=D_{X}\otimes I_{m},\quad\mathbf{D}_{Y}=I_{n}\otimes D_{Y},$

where $\otimes$ denotes the Kronecker product. We can then write

$B_{n,m}=-\mathbf{D}_{X}-\mathbf{D}_{Y}+\mathbf{A}+\mathbf{B}-\mathbf{M},$

where $\mathbf{A},\mathbf{B},\mathbf{M}\in\mathbb{R}^{nm\times nm}$ are
defined by

$\displaystyle\mathbf{A}_{(k,l),(i,j)}$
$\displaystyle=\int_{x_{0}}^{x_{k}}K_{\alpha}(\ell_{X,i}^{\prime}\ell_{Y,j}^{\prime})(a)\mathop{}\\!\mathrm{d}a,$
(4.2) $\displaystyle\mathbf{B}_{(k,l),(i,j)}$
$\displaystyle=\int_{y_{0}}^{y_{l}}K_{\beta}(\ell_{X,i}^{\prime}\ell_{Y,j}^{\prime})(b)\mathop{}\\!\mathrm{d}b,$
(4.3) $\displaystyle\mathbf{M}_{(k,l),(i,j)}$
$\displaystyle=\iint_{R(x_{k},y_{l})}\mu(a,b)\ell_{X,i}^{\prime}(a)\ell_{Y,j}^{\prime}(b)\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b,$
(4.4)

for $k,i=1,\ldots,n$ and $l,j=1,\ldots,m$. Note that, if $\mu$ is constant,
the matrix $\mathbf{M}$ is diagonal with diagonal entries equal to $\mu$.

We finally note that, although the matrix $B_{n,m}$ is defined for any set of
nodes, the choice of the latter is critical to ensure the convergence of the
interpolating polynomials and, in turn, of the elements of the spectrum. In
the following numerical experiments, we choose the Chebyshev extremal points
in each interval $[x_{0},\bar{x}]$ and $[y_{0},\bar{y}]$.

In the univariate case, these nodes guarantee that the convergence rate of the
interpolating polynomial of degree $n$ is $O(n^{-k})$ if the interpolated
function is $C^{k}$ [42, Theorem 7.2], which implies that the order of
convergence is infinite if the function is smooth. Moreover, the convergence
rate is $O(c^{n})$ for some $0<c<1$ if the function is analytic [42, Theorem
8.2]. The two latter properties are often known as _spectral accuracy_ , see
[40, chapter 4] and [3, chapter 2]. Furthermore, observe that the relevant
differentiation matrices can be computed explicitly [40].

The classic result on the interpolation error being bounded by means of the
best uniform approximation error and the Lebesgue constant holds also in the
bivariate case. A multidimensional version of Jackson’s theorem on the best
uniform approximation error holds as well [37, Theorem 4.8]. Moreover, it is
easy to verify that the Lebesgue constant for the tensor-product Chebyshev
extremal nodes in $\overline{R}$ is the product of the univariate Lebesgue
constants in $[x_{0},\bar{x}]$ and $[y_{0},\bar{y}]$, hence it is $O(\log
n\log m)$. The tensor-product Chebyshev interpolation is thus near-optimal
also in the bivariate case.

Although a proof of convergence for the method is out of the scope of this
paper, we show that the order of convergence observed numerically for the
approximated eigenvalues and eigenvectors is consistent with the well-
established order of convergence of polynomial interpolation.

For implementation purposes, we observe that in general the integrals defining
$\mathbf{A}$, $\mathbf{B}$ and $\mathbf{M}$ in (4.2)–(4.4) cannot be computed
exactly. To approximate the integrals on $\overline{R}$ we use the
Clenshaw–Curtis cubature formula [19], which is based on Chebyshev extremal
points and is spectrally accurate [40, 41].

To compute the entries of $\mathbf{A}$ and $\mathbf{B}$, given a function $f$
defined on $[x_{0},\bar{x}]$ and $F\in\mathbb{R}^{n}$ such that
$F_{i}=f(x_{i})$, we consider the approximation

$\int_{x_{0}}^{x_{k}}f(a)\mathop{}\\!\mathrm{d}a\approx[D_{X}^{-1}F]_{k},\quad
k=1,\ldots,n,$

and similarly for functions defined on $[y_{0},\bar{y}]$ (see, e.g., [21]).
This approximation can be extended to the double integrals involved in the
entries of the matrix $\mathbf{M}$ for nonconstant $\mu$. More precisely,
given a function $\psi\in AC_{0}(\overline{R})$ and a vector
$\Psi\in\mathbb{R}^{nm}$ such that $\Psi_{k,l}=\psi(x_{k},y_{l})$, each
integral
$\int_{x_{0}}^{x_{k}}\int_{y_{0}}^{y_{l}}\psi(a,b)\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b$
can be approximated by the corresponding $(k,l)$-th entry of the vector
$\mathbf{D}_{X}^{-1}\mathbf{D}_{Y}^{-1}\Psi$.

##  5 Numerical experiments

In this section, we present several numerical experiments to investigate how
the spectrum of the finite-dimensional operator $B_{n,m}$ approximates the
spectrum of $B$, and in turn of $A$, of each problem at hand. For this
purpose, we select several parameter sets for which eigenvalues and
eigenfunctions of $A$ can be expressed analytically, and we study the
convergence of the approximated eigenvalues of $B_{n,m}$ to the analytic ones.
As for the eigenfunctions, we stress that, since $B_{n,m}$ represents an
approximation of the operator $B$, an eigenvector $\Psi$ of $B_{n,m}$ provides
an approximation $\psi_{n,m}$ of $V\phi$, where $\phi$ is an eigenfunction of
$A$; an approximation of $\phi$ is thus given by
$\partial_{x}\partial_{y}\psi_{n,m}$.

For each example we study the behavior for increasing $n=m$ of the absolute
error $\epsilon_{\lambda}$ on the known eigenvalue $\lambda$ and of the
absolute error $\epsilon_{\phi}$ in $L^{1}$ norm on the known eigenfunction
$\phi$, computed via Clenshaw–Curtis cubature.

In all examples we choose

$\alpha(x,\xi,\sigma)=\alpha_{1}(x)\gamma(\xi,\sigma),\quad\beta(y,\xi,\sigma)=\beta_{1}(y)\gamma(\xi,\sigma)$

in the boundary conditions (2.2)–(2.3), in order to simplify finding an
explicit eigenfunction.

We remark that the parameters are chosen in order to have an analytically
known eigenfunction with certain smoothness properties, without regard to any
specific biological interpretation.

To compute the spectrum of $B_{n,m}$ we use standard methods (namely MATLAB’s
eig function). Note that the approximated spectrum may contain spurious
eigenvalues (e.g., when $B$ has fewer eigenvalues than the dimension of
$B_{n,m}$); however, in our examples we only examine specific eigenvalues, so
that the possible spurious ones do not affect our analysis.

###  5.1 Analytic eigenfunctions

| Ex. 1.1 | Ex. 1.2 | Ex. 1.3 | Ex. 1.4
---|---|---|---|---
$\overline{R}$ | $[0,1]^{2}$ | $\left[\frac{\pi}{6},\frac{\pi}{2}\right]\times\left[\frac{\pi}{6},\frac{\pi}{4}\right]$ | $[0,2]\times[{-1},1]$ | $[0,2]\times[0,1]$
$\alpha_{1}(x)$ | $1$ | $\cos(x-\frac{\pi}{6})$ | $e^{x+1}$ | $e^{-x^{2}}$
$\beta_{1}(y)$ | $1$ | $\cos(\frac{\pi}{6}-y)$ | $e^{-y}$ | $e^{y}$
$\gamma(\xi,\sigma)$ | $1$ | $\left(\frac{\sqrt{2}-\sqrt{6}+2}{4}\right)^{-1}$ | $\frac{1}{4}e^{-\xi+\sigma}$ | $\left(\iint_{\overline{R}}e^{-a^{2}+b}\mathop{}\\!\mathrm{d}b\mathop{}\\!\mathrm{d}a\right)^{-1}$
$\mu(x,y)$ | $1$ | $1$ | $1$ | $2x+1$
$\lambda$ | $-1$ | $-1$ | $-1$ | $-2$
$\phi(x,y)$ | $1$ | $\cos(x-y)$ | $e^{x-y}$ | $e^{-x^{2}+y}$
Table 5.1: Parameters and resulting eigenvalue and eigenfunction in the
analytic cases.

We consider a first group of examples presenting an analytic eigenfunction.
The choices of the parameters and the resulting eigenvalue and eigenfunction
are listed in Table 5.1. Starting from Example 1.1, where all parameters are
constant, we gradually introduce nonconstant coefficients: $\alpha_{1}$ and
$\beta_{1}$ in Example 1.2, $\gamma$ in Example 1.3 and $\mu$ in Example 1.4.

Figure 5.1: Errors $\epsilon_{\lambda}$ and $\epsilon_{\phi}$ for Examples
1.1, 1.2, 1.3 and 1.4 defined in Table 5.1. Observe that for Example 1.1 with
$n=m=1$ the error $\epsilon_{\phi}$ is exactly $0$, hence it is not
represented in the logarithmic scale.

Considering Example 1.1, Figure 5.1 shows that the errors reach the machine
precision already for $n=m=2$. With $n=m=1$ the error $\epsilon_{\phi}$ on the
eigenfunction is exactly equal to $0$, which may be explained by the fact that
constant functions are interpolated exactly already by polynomials of degree
$0$. As $n=m$ increases, the errors increase, possibly due to numerical
instability.

Considering now Examples 1.2, 1.3 and 1.4, Figure 5.1 shows that both
$\epsilon_{\lambda}$ and $\epsilon_{\phi}$ decay with infinite order. Observe
that in Example 1.4 more nodes are needed to reach the error barrier than in
Examples 1.2 and 1.3, probably due to the approximation of the integrals
involving the nonconstant $\mu$.

###  5.2 Nonsmooth eigenfunctions

| Ex. 2.1 ($C^{2}$) | Ex. 2.2 ($C^{1}$) | Ex. 2.3 ($C^{0}$) | Ex. 2.4 (discontinuous)
---|---|---|---|---
$\overline{R}$ | $[0,1]\times[0,2]$ | $[0,1]\times[0,2]$ | $[0,1]\times[0,2]$ | $[0,1]\times[0,2]$
$\alpha_{1}(x)$ | $x^{2}\lvert x\rvert$ | $-x\lvert x\rvert$ | $\lvert x\rvert$ | $\mathbf{1}_{[0,+\infty[}(x)$
$\beta_{1}(y)$ | $y^{2}\lvert y\rvert$ | $y\lvert y\rvert$ | $\lvert y\rvert$ | $\mathbf{1}_{]-\infty,0]}(y)$
$\gamma(\xi,\sigma)$ | $\frac{5}{8}$ | $\frac{6}{7}$ | $\frac{3}{4}$ | $2$
$\mu(x,y)$ | $1$ | $1$ | $1$ | $1$
$\lambda$ | $-1$ | $-1$ | $-1$ | $-1$
$\phi(x,y)$ | $(x-y)^{2}\lvert x-y\rvert$ | $(x-y)\lvert x-y\rvert$ | $\lvert x-y\rvert$ | $\mathbf{1}_{[0,+\infty[}(x-y)$
Table 5.2: Parameters and resulting eigenvalue and eigenfunction in the
nonsmooth case. $\mathbf{1}_{A}$ is the indicator function of
$A\subset\mathbb{R}$. The regularity of the eigenfunction is shown in
parentheses in the titles.

For the second group of examples, we consider eigenfunctions which are not
smooth. The choices of the parameters and the resulting eigenvalue and
eigenfunction are listed in Table 5.2. Observe that the eigenfunctions have
the same regularity as the coefficients $\alpha_{1}$ and $\beta_{1}$, namely
$C^{2}$ for Example 2.1, $C^{1}$ for Example 2.2 and $C^{0}$ for Example 2.3.

We consider also Example 2.4 with discontinuous coefficients and
eigenfunction. In [31] it is assumed that $\alpha$ and $\beta$ are Lipschitz
continuous in order to prove that the semigroup $\\{T(t)\\}_{t\geq 0}$ is
eventually compact, but actually this hypothesis is used only in the interior
of the domains: thus, Example 2.4 still falls inside the scope of [31].
Observe also that the function $\mathbf{1}_{[0,+\infty[}(x-y)$, even if not
continuous, is in $D(A)$.

Figure 5.2: Errors $\epsilon_{\lambda}$ (left) and $\epsilon_{\phi}$ (right)
for Examples 2.1, 2.2, 2.3 and 2.4 defined in Table 5.2. The slopes of the
dashed gray lines are $-2$, $-3.5$, $-4.5$, $-5.5$ (left) and $-0.7$, $-1.7$,
$-2.7$, $-3.7$ (right), included to ease the interpretation of the plots.

Figure 5.2 suggests that the errors decay with finite order and that these
orders increase with the regularity of the eigenfunction. In particular,
focusing on Examples 2.1–2.3, we can observe that for both errors a loss of
one order of differentiability of the eigenfunction seems to correspond to a
loss of about one order of convergence (cfr. the dashed reference lines in
Figure 5.2). The convergence of $\epsilon_{\lambda}$ seems to be almost two
orders faster than that of $\epsilon_{\phi}$. To possibly explain this
difference, recall that we are actually collocating the eigenvalue problem for
$B$, which means that the eigenvalues are the same as $A$, but the
eigenfunctions correspond to integrals of the eigenfunctions of $A$, so the
comparison between the eigenfunctions involves differentiating the computed
ones.

As for Example 2.4, it appears that losing continuity of the eigenfunction
itself changes the order of convergence differently for $\epsilon_{\lambda}$
and $\epsilon_{\phi}$.

##  6 Structuring variables with nontrivial velocity

We have illustrated the method for systems in which both physiological
variables evolve at the same velocity as time. This should not be seen as too
restrictive, as systems with more general velocity terms [2, 23, 28, 32, 33]
can in some cases be reduced to (2.1)–(2.3) after a suitable scaling of
variables, so that similar theoretical results on the stability of the zero
solution hold [43]. In practice, however, the computation of the change of
variables, which in general is defined by the solution of an ODE system, may
be expensive, although necessary when the individual parameters (e.g., birth
and mortality rates) depend on the original (unscaled) variables. In this
case, directly approximating the original problem with nontrivial velocities
may be convenient from a computational point of view, as observed in [34].

In fact, the transformation via integration can be easily carried out for
problems of the form

$\displaystyle\partial_{t}u(t,x,y)+\partial_{x}(g_{X}(x)u(t,x,y))+\partial_{y}(g_{Y}(y)u(t,x,y))=-\mu(x,y)u(t,x,y),$
$\displaystyle\begin{aligned}
g_{Y}(y_{0})u(t,x,y_{0})&=K_{\alpha}(u(t,\cdot,\cdot))(x),\\\
g_{X}(x_{0})u(t,x_{0},y)&=K_{\beta}(u(t,\cdot,\cdot))(y),\end{aligned}$

where the positive functions $g_{X}(x)$ and $g_{Y}(y)$ describe the rates of
change of $x$ and $y$ in time. In this case, it is straightforward to verify
that, given the IG $A$, the operator $B=VAV^{-1}$ admits the representation

$\displaystyle B\psi(x,y)=$ $\displaystyle-
g_{X}(x)\partial_{x}\psi(x,y)-g_{Y}(y)\partial_{y}\psi(x,y)$
$\displaystyle+\int_{x_{0}}^{x}K_{\alpha}(\partial_{y}\partial_{x}\psi)(a)\mathop{}\\!\mathrm{d}a+\int_{y_{0}}^{y}K_{\beta}(\partial_{y}\partial_{x}\psi)(b)\mathop{}\\!\mathrm{d}b$
$\displaystyle-\iint_{R(x,y)}\mu(a,b)\partial_{y}\partial_{x}\psi(a,b)\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b,$

which is approximated by a matrix of the form

$B_{n,m}=-\mathbf{G}_{X}\mathbf{D}_{X}-\mathbf{G}_{Y}\mathbf{D}_{Y}+\mathbf{A}+\mathbf{B}-\mathbf{M},$

where $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{M}$ are the matrices defined in
section 4 and $\mathbf{G}_{X}=G_{X}\otimes I_{m}$,
$\mathbf{G}_{Y}=I_{n}\otimes G_{Y}$, with $G_{X}$ and $G_{Y}$ diagonal
matrices defined by $[G_{X}]_{i,i}=g_{X}(x_{i})$, $i=1,\dots,n$, and
$[G_{Y}]_{i,i}=g_{Y}(y_{j})$, $j=1,\dots,m$.

As an example let us consider

$\displaystyle\overline{R}\coloneqq\left[\frac{1}{2},\frac{3}{2}\right]\times\left[\frac{1}{2},2\right],\quad
g_{X}(x,y)\coloneqq x,\quad g_{Y}(x,y)\coloneqq\frac{y^{2}}{2},$
$\displaystyle\mu(x,y)\coloneqq y^{3}-2x^{2}-y+4,$
$\displaystyle\alpha(x,\xi,\sigma)\coloneqq\frac{e^{x^{2}-\frac{1}{4}}}{8c},\quad\beta(y,\xi,\sigma)\coloneqq\frac{e^{-y^{2}+\frac{1}{4}}}{2c},\quad
c=\frac{1}{4}\pi\iint_{\overline{R}}e^{a^{2}-b^{2}}\mathop{}\\!\mathrm{d}a\mathop{}\\!\mathrm{d}b,$

for which $\lambda=-5$ is an eigenvalue of the corresponding IG with
eigenfunction $\phi(x,y)=e^{x^{2}-y^{2}}$. We can observe in Figure 6.1 that
the errors computed by our method decay with infinite order even in this case.

Figure 6.1: Errors $\epsilon_{\lambda}$ and $\epsilon_{\phi}$ for the example
of section 6.

##  7 Concluding remarks

In this paper we proposed a numerical technique to analyze the stability of
the zero solution of linear population models with two structuring variables,
of the type considered in [31, 43].

Extensive numerical tests illustrate the convergence of the eigenvalues of the
finite-dimensional approximation to the true eigenvalues of the IG. The
numerical tests support the conjecture that the order of convergence of the
approximation depends on the regularity of the eigenfunction. A rigorous
theoretical proof of the convergence of the approximation is left to future
work.

Stability analysis requires not only that the eigenvalues of the IG are
approximated accurately, but also that no spurious eigenvalues are to the
right of the true spectral abscissa. In fact, in our examples we observe that
the approximation of the known eigenvalue is the numerical spectral abscissa,
suggesting that the method can be effectively used to study the stability.

Structured population models can also be formulated as renewal equations for
the population birth rate (or “recruitment function”) [1, 17, 20]. The renewal
equation formulation is particularly convenient from the theoretical point of
view as one can exploit the principle of linearized stability for nonlinear
equations, which can not be proved in general for the PDE formulation [1].
Results on the asymptotic behavior of solutions have also been recently
proved, under special assumptions, for renewal equations defined on a space of
measures, which makes it possible to consider a wider set of solutions
compared to PDEs [25, 26].

It would be interesting to apply pseudospectral methods in the framework of
renewal equations (admitting a state space of multivalued functions or even
measures), by extending the techniques developed for scalar renewal equations
[7, 12, 35]. However, when the structuring variables evolve with nontrivial
velocity, the renewal equation formulation requires to explicitly invert the
age-structure relation defined implicitly by an ODE system, which suffers from
the computational challenges highlighted in section 6. Hence, as explained
therein, directly tackling the PDE formulation may be computationally
convenient, as it bypasses the solution of the ODE system.

In this paper we restricted to structuring variables in bounded intervals, as
this allows to exploit the highly desirable convergence properties of
polynomial interpolation on Chebyshev nodes. However, unbounded domains are
common in the modeling literature (e.g., [2, 23, 30, 38]), for instance when
it is not easy to determine a suitable upper bound for a physiological
variable _a priori_ , or when the processes are naturally described by
probability distributions with unbounded support (e.g., exponential or Gamma
distributions).

In order to numerically treat these problems, truncating the domain can be
feasible sometimes, but the accuracy of the approximation would depend on both
the size of the truncated domain and the number of nodes in the domain. The
latter usually becomes very large because the choice of Chebyshev nodes does
not exploit the specific characteristics of the solutions, which usually
belong to exponentially weighted spaces [1]. For this reason, using
exponentially weighted interpolation and Laguerre-type nodes has proved
successful and more efficient than domain truncation in the case of delay
equations [27, 36]. It would be interesting to apply similar techniques [3] to
structured models in the PDE formulation with one or even two structuring
variables, although the latter brings in additional complications due to the
necessity to rely on multivariate interpolation.

##  Appendix A Models with one structuring variable

As recalled in the introduction, [4] already provides a pseudospectral method,
based on a different approach, to approximate the IG of models with one
structuring variable. In this appendix, for completeness, we adapt our
approach to this case, providing also an example.

We consider the scalar first-order linear hyperbolic PDE

$\partial_{t}u(t,x)+\partial_{x}u(t,x)=-\mu(x)u(t,x),$ (A.1)

with boundary condition

$u(t,x_{0})=\int_{x_{0}}^{\bar{x}}\beta(\sigma)u(t,\sigma)\mathop{}\\!\mathrm{d}\sigma\eqqcolon
K_{\beta}(u(t,\cdot)).$ (A.2)

As references on single structure models, see [29, 30, 44]; in particular, see
[30, section 1.2] for what concerns this appendix.

If $\mu,\beta\in L^{1}([x_{0},\bar{x}])$ are nonnegative and bounded, for
every $u_{0}\in L^{1}([x_{0},\bar{x}])$ the initial–boundary value problem
defined by (A.1)–(A.2) and $u(0,\cdot)=u_{0}$ admits a unique solution
$u(t,\cdot)\in L^{1}([x_{0},\bar{x}])$ for $t\geq 0$. Moreover, the family of
solution operators $\\{T(t)\\}_{t\geq 0}$, defined by $T(t)u_{0}=u(t,\cdot)$,
forms a strongly continuous and eventually compact semigroup of bounded linear
operators in the Banach space $L^{1}([x_{0},\bar{x}])$. Its IG $A\colon
D(A)\to L^{1}([x_{0},\bar{x}])$, defined as in (2.4), can be expressed as

$A\phi=-\phi^{\prime}-\mu\phi,$

and its domain $D(A)$ can be characterized as

$D(A)=\\{\phi\in AC([x_{0},\bar{x}])\;|\;\phi(x_{0})=K_{\beta}(\phi)\\}.$

Let us equip $AC([x_{0},\bar{x}])$ with the norm
$\lVert\cdot\rVert_{AC([x_{0},\bar{x}])}$ defined as $\lVert
f\rVert_{AC([x_{0},\bar{x}])}\coloneqq\lvert f(x_{0})\rvert+\lVert
f^{\prime}\rVert_{L^{1}([x_{0},\bar{x}])}$. Let $AC_{0}([x_{0},\bar{x}])$ be
the subspace of $AC([x_{0},\bar{x}])$ of functions that are null at $x_{0}$,
which is a Banach space, being a closed subspace. The operator $V\colon
L^{1}([x_{0},\bar{x}])\to AC_{0}([x_{0},\bar{x}])$ defined by
$V\phi(x)=\int_{x_{0}}^{x}\phi(\sigma)\mathop{}\\!\mathrm{d}\sigma$ defines an
isomorphism between $L^{1}([x_{0},\bar{x}])$ and $AC_{0}([x_{0},\bar{x}])$,
with $V^{-1}\psi=\psi^{\prime}$ for all $\psi\in AC_{0}([x_{0},\bar{x}])$.
Observe that both $V$ and $V^{-1}$ are bounded ($\lVert V\rVert=\lVert
V^{-1}\rVert=1$).

We define the family of operators $\\{S(t)\\}_{t\geq 0}$ on
$AC_{0}([x_{0},\bar{x}])$ as $S(t)\coloneqq VT(t)V^{-1}$. As in section 3, we
observe that they form a strongly continuous and eventually compact semigroup
on $AC_{0}([x_{0},\bar{x}])$ with IG $B\colon D(B)\to
AC_{0}([x_{0},\bar{x}])$, with $B\coloneqq VAV^{-1}$ and $D(B)=VD(A)$. We can
also derive the following expression for $B$, given $\psi\in D(B)$:

$B\psi=-\psi^{\prime}+K_{\beta}(\psi^{\prime})-V(\mu\psi^{\prime}).$

As in section 3, we can conclude that $A$ and $B$ have the same spectrum, at
most countable and consisting only of eigenvalues (of finite algebraic
multiplicity).

As an example, let us choose $[x_{0},\bar{x}]=[0,2]$, $\mu(x)\equiv 1$ and
$\beta(x)\coloneqq e^{-x}$. It can be shown that the only real eigenvalue of
the corresponding IG is the unique real solution of the equation

$\frac{1-e^{-2\lambda-4}}{\lambda+2}=1,$

which can be approximated to the machine precision with standard methods
(e.g., with MATLAB’s fzero we obtain $\lambda=-1.203187869979980$). The
relevant eigenfuction is

$\phi(x)=e^{-\int_{0}^{x}\mu(s)\mathop{}\\!\mathrm{d}s-\lambda
x}=e^{-(1+\lambda)x}.$

We can observe in Figure A.1 that the errors computed by our method decay with
infinite order.

Figure A.1: Errors $\epsilon_{\lambda}$ and $\epsilon_{\phi}$ for the example
of appendix A.

## Acknowledgements

We are grateful to Odo Diekmann for suggesting that we refer to the theory of
absolute continuity in the sense of Carathéodory.

The authors are members of INdAM Research group GNCS. Davide Liessi and
Francesca Scarabel are members of UMI Research group “Modellistica socio-
epidemiologica”. The work of Simone De Reggi and Davide Liessi was supported
by the Italian Ministry of University and Research (MUR) through the PRIN 2020
project (No. 2020JLWP23) “Integrated Mathematical Approaches to Socio-
Epidemiological Dynamics” (CUP: E15F21005420006). Francesca Scarabel is
supported by the UKRI through the JUNIPER modelling consortium (grant number
MR/V038613/1).

## References

  * [1] Carles Barril, Àngel Calsina, Odo Diekmann and József Z. Farkas  “On the formulation of size-structured consumer resource models (with special attention for the principle of linearised stability)”, 2021  arXiv:2111.09678 [math.AP]
  * [2] Fadia Bekkal Brikci, Jean Clairambault, Benjamin Ribba and Benoît Perthame  “An age-and-cyclin-structured cell population model for healthy and tumoral tissues”  In _Journal of Mathematical Biology_ 57.1, 2008, pp. 91–110  DOI: 10.1007/s00285-007-0147-x
  * [3] John P. Boyd  “Chebyshev and Fourier Spectral Methods”  Mineola, NY: Dover, 2001
  * [4] Dimitri Breda, Caterina Cusulin, Mimmo Iannelli, Stefano Maset and Rossana Vermiglio  “Stability analysis of age-structured population equations by pseudospectral differencing methods”  In _Journal of Mathematical Biology_ 54.5, 2007, pp. 701–720  DOI: 10.1007/s00285-006-0064-4
  * [5] Dimitri Breda, Simone De Reggi, Francesca Scarabel, Rossana Vermiglio and Jianhong Wu  “Bivariate collocation for computing ${R}_{0}$ in epidemic models with two structures”  In _Computers & Mathematics with Applications_, 2021  DOI: 10.1016/j.camwa.2021.10.026
  * [6] Dimitri Breda, Odo Diekmann, Mats Gyllenberg, Francesca Scarabel and Rossana Vermiglio  “Pseudospectral discretization of nonlinear delay equations: new prospects for numerical bifurcation analysis”  In _SIAM Journal on Applied Dynamical Systems_ 15.1, 2016, pp. 1–23  DOI: 10.1137/15M1040931
  * [7] Dimitri Breda, Odo Diekmann, Stefano Maset and Rossana Vermiglio  “A numerical approach for investigating the stability of equilibria for structured population models”  In _Journal of Biological Dynamics_ 7.sup1, 2013, pp. 4–20  DOI: 10.1080/17513758.2013.789562
  * [8] Dimitri Breda, Francesco Florian, Jordi Ripoll and Rossana Vermiglio  “Efficient numerical computation of the basic reproduction number for structured populations”  In _Journal of Computational and Applied Mathematics_ 384, 2021  DOI: 10.1016/j.cam.2020.113165
  * [9] Dimitri Breda, Philipp Getto, Julia Sánchez Sanz and Rossana Vermiglio  “Computing the eigenvalues of realistic Daphnia models by pseudospectral methods”  In _SIAM Journal on Scientific Computing_ 37.6, 2015, pp. A2607–A2629  DOI: 10.1137/15M1016710
  * [10] Dimitri Breda, Mimmo Iannelli, Stefano Maset and Rossana Vermiglio  “Stability Analysis of the Gurtin–MacCamy Model”  In _SIAM Journal on Numerical Analysis_ 46.2, 2008, pp. 980–995  DOI: 10.1137/070685658
  * [11] Dimitri Breda, Toshikazu Kuniya, Jordi Ripoll and Rossana Vermiglio  “Collocation of next-generation operators for computing the basic reproduction number of structured populations”  In _Journal of Scientific Computing_ 85, 2020  DOI: 10.1007/s10915-020-01339-1
  * [12] Dimitri Breda and Davide Liessi  “Approximation of eigenvalues of evolution operators for linear renewal equations”  In _SIAM Journal on Numerical Analysis_ 56.3, 2018, pp. 1456–1481  DOI: 10.1137/17M1140534
  * [13] Dimitri Breda and Davide Liessi  “Approximation of eigenvalues of evolution operators for linear coupled renewal and retarded functional differential equations”  In _Ricerche di Matematica_ 69.2, 2020, pp. 457–481  DOI: 10.1007/s11587-020-00513-9
  * [14] Dimitri Breda, Stefano Maset and Rossana Vermiglio  “Pseudospectral Differencing Methods for Characteristic Roots of Delay Differential Equations”  In _SIAM Journal on Scientific Computing_ 27.2, 2005, pp. 482–495  DOI: 10.1137/030601600
  * [15] Dimitri Breda, Stefano Maset and Rossana Vermiglio  “Approximation of eigenvalues of evolution operators for linear retarded functional differential equations”  In _SIAM Journal on Numerical Analysis_ 50.3, 2012, pp. 1456–1483  DOI: 10.1137/100815505
  * [16] Dimitri Breda, Stefano Maset and Rossana Vermiglio  “Stability of Linear Delay Differential Equations”, SpringerBriefs in Control, Automation and Robotics  New York: Springer, 2015  DOI: 10.1007/978-1-4939-2107-2
  * [17] Àngel Calsina, Odo Diekmann and József Z. Farkas  “Structured populations with distributed recruitment: from PDE to delay formulation”  In _Mathematical Methods in the Applied Sciences_ 39.18, 2016, pp. 5175–5191  DOI: 10.1002/mma.3898
  * [18] Ph…. Clément, H…. Heijmans, S. Angenent, C.. Duijn and B. Pagter  “One-Parameter Semigroups”, CWI Monographs 5  Netherlands: North-Holland Publishing Company, 1987
  * [19] C.. Clenshaw and A.. Curtis  “A method for numerical integration on an automatic computer”  In _Numerische Mathematik_ 2, 1960, pp. 197–205  DOI: 10.1007/BF01386223
  * [20] Odo Diekmann, Mats Gyllenberg, J… Metz, Shinji Nakaoka and André Marc Roos  “ _Daphnia_ revisited: local stability and bifurcation theory for physiologically structured population models explained by way of an example”  In _Journal of Mathematical Biology_ 61.2, 2010, pp. 277–318  DOI: 10.1007/s00285-009-0299-y
  * [21] Odo Diekmann, Francesca Scarabel and Rossana Vermiglio  “Pseudospectral discretization of delay differential equations in sun-star formulation: Results and conjectures”  In _Discrete and Continuous Dynamical Systems. Series S_ 13.9, 2020, pp. 2575–2602  DOI: 10.3934/dcdss.2020196
  * [22] Odo Diekmann, Stephan A. Gils, Sjoerd M. Verduyn Lunel and Hans-Otto Walther  “Delay Equations”, Applied Mathematical Sciences 110  New York: Springer, 1995  DOI: 10.1007/978-1-4612-4206-2
  * [23] Janet Dyson, Rosanna Villella-Bressan and Glenn Webb  “A nonlinear age and maturity structured model of population dynamics”  In _Journal of Mathematical Analysis and Applications_ 242.1, 2000, pp. 93–104  DOI: 10.1006/jmaa.1999.6656
  * [24] Klaus-Jochen Engel and Rainer Nagel  “One-Parameter Semigroups for Linear Evolution Equations”, Graduate Texts in Mathematics 194  New York: Springer, 2000  DOI: 10.1007/b97696
  * [25] Eugenia Franco, Odo Diekmann and Mats Gyllenberg  “Modelling physiologically structured populations: renewal equations and partial differential equations”, 2022  arXiv:2201.05323 [math.AP]
  * [26] Eugenia Franco, Mats Gyllenberg and Odo Diekmann  “One dimensional reduction of a renewal equation for a measure-valued function of time describing population dynamics”  In _Acta Applicandae Mathematicae_ 175, 2021  DOI: 10.1007/s10440-021-00440-3
  * [27] Mats Gyllenberg, Francesca Scarabel and Rossana Vermiglio  “Equations with infinite delay: Numerical bifurcation analysis via pseudospectral discretization”  In _Applied Mathematics and Computation_ 333, 2018, pp. 490–505  DOI: 10.1016/j.amc.2018.03.104
  * [28] Keith E. Howard  “A size and maturity structured model of cell dwarfism exhibiting chaotic behavior”  In _International Journal of Bifurcation and Chaos in Applied Sciences and Engineering_ 13.10, 2003, pp. 3001–3013  DOI: 10.1142/S0218127403008363
  * [29] Mimmo Iannelli  “Mathematical Theory of Age-Structured Population Dynamics”  Pisa: Giardini Editori e Stampatori, 1995
  * [30] Hisashi Inaba  “Age-Structured Population Dynamics in Demography and Epidemiology”  Singapore: Springer, 2017  DOI: 10.1007/978-981-10-0188-8
  * [31] Hao Kang, Xi Huo and Shigui Ruan  “On first-order hyperbolic partial differential equations with two internal variables modeling population dynamics of two physiological structures”  In _Annali di Matematica Pura ed Applicata_ 200, 2021, pp. 403–452  DOI: 10.1007/s10231-020-01001-5
  * [32] Pierre Magal and Shigui Ruan  “Structured Population Models in Biology and Epidemiology”, Lecture Notes in Mathematics 1936  Berlin, Heidelberg: Springer, 2008  DOI: 10.1007/978-3-540-78273-5
  * [33] “The Dynamics of Physiologically Structured Populations”, Lecture Notes in Biomathematics 68  Berlin, Heidelberg: Springer, 1986  DOI: 10.1007/978-3-662-13159-6
  * [34] Francesca Scarabel, Dimitri Breda, Odo Diekmann, Mats Gyllenberg and Rossana Vermiglio  “Numerical bifurcation analysis of physiologically structured population models via pseudospectral approximation”  In _Vietnam Journal of Mathematics_ 49, 2021, pp. 37–67  DOI: 10.1007/s10013-020-00421-3
  * [35] Francesca Scarabel, Odo Diekmann and Rossana Vermiglio  “Numerical bifurcation analysis of renewal equations via pseudospectral approximation”  In _Journal of Computational and Applied Mathematics_ 397, 2021  DOI: 10.1016/j.cam.2021.113611
  * [36] Francesca Scarabel and Rossana Vermiglio  “Equations with infinite delay: pseudospectral approximation of characteristic roots in an abstract framework”
  * [37] Martin H. Schultz  “${L}^{\infty}$-multivariate approximation theory”  In _SIAM Journal on Numerical Analysis_ 6.2, 1969, pp. 161–183  DOI: 10.1137/0706017
  * [38] Eugenio Sinestrari and Glenn F. Webb  “Nonlinear hyperbolic systems with nonlocal boundary conditions”  In _Journal of Mathematical Analysis and Applications_ 121.2, 1987, pp. 449–464  DOI: 10.1016/0022-247X(87)90255-1
  * [39] Jiří Šremr  “Absolutely continuous functions of two variables in the sense of Carathéodory”  In _Electronic Journal of Differential Equations_ 2010, 2010, pp. 1–11
  * [40] Lloyd Nicholas Trefethen  “Spectral Methods in MATLAB”, Software, Environments and Tools  Philadelphia: Society for IndustrialApplied Mathematics, 2000  DOI: 10.1137/1.9780898719598
  * [41] Lloyd Nicholas Trefethen  “Is Gauss Quadrature Better than Clenshaw–Curtis?”  In _SIAM Review_ 50.1, 2008, pp. 67–87  DOI: 10.1137/060659831
  * [42] Lloyd Nicholas Trefethen  “Approximation Theory and Approximation Practice”  Philadelphia: Society for IndustrialApplied Mathematics, 2013
  * [43] G.. Webb  “Dynamics of populations structured by internal variables”  In _Mathematische Zeitschrift_ 189, 1985, pp. 319–335  DOI: 10.1007/BF01164156
  * [44] Glenn F. Webb  “Theory of Nonlinear Age-Dependent Population Dynamics”, Monographs and Textbooks in Pure and Applied Mathematics 89  New York: Marcel Dekker, Inc., 1985