Original article Biomath 3 (2014), 1403241, 1–11

B f

Volume ░, Number ░, 20░░ 

BIOMATH

 ISSN 1314-684X

Editor–in–Chief: Roumen Anguelov  

B f

BIOMATH
h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum

Numerical Analysis of a Size-Structured
Population Model with a Dynamical Resource

Oscar Angulo∗, J.C. López-Marcos∗ and M.A. López-Marcos∗
∗Departamento de Matemática Aplicada

Universidad de Valladolid, Valladolid, Spain
Emails: oscar@mat.uva.es, lopezmar@mac.uva.es, malm@mac.uva.es

Received: 19 October 2013, accepted: 24 March 2014, published: 28 May 2014

Abstract—In this paper, we analyze the conver-
gence of a second-order numerical method for the
approximation of a size-structured population model
whose dependency on the environment is managed
by the evolution of a vital resource. Optimal con-
vergence rate is derived. Numerical experiments are
also reported to demonstrate the predicted accuracy
of the scheme. Also, it is applied to solve a problem
that describes the dynamics of a Daphnia magna
population, paying attention to the unstable case.

Keywords-structured population models; numeri-
cal methods; convergence; Daphnia magna

I. Introduction

Physiologically structured population models
are based on the use of one or more attributes that
structure the individuals in the population. Size is
one of the most natural and important attributes
structuring the population for many species: typi-
cal examples being fishes and trees. In such species
the ability of an individual to obtain the neces-
sary resources to survive and reproduce depends
strongly on its size. Structured population models
reflect the effect of physiological state of individ-
uals on the population dynamics. In addition, the
use of nonlinear structured population models al-

lows us to take into account the effect of competi-
tion for natural resources in the structured-specific
growth, mortality and fertility rates. We can find an
extensive study of physiologically structured pop-
ulation models, with analytical studies of aspects
such as existence and uniqueness, smoothness and
the asymptotic behaviour of solutions in [1], [2],
[3], [4], [5].

In this paper, we consider a size-structured pop-
ulation model nonlinearly coupled with an integro-
ordinary differential equation accounting for sub-
strate consumption and/or product formation. It
was introduced first in [6] for modeling a Daphnia
magna population. The model involves a nonlinear
hyperbolic partial differential equation

ut + (g(x, S (t), t) u)x = −µ(x, S (t), t) u, (1)

0 < x < xM (t), t > 0,

a nonlinear and nonlocal boundary condition
which reflects the reproduction process

g(0, S (t), t) u(0, t) =
∫ xM (t)

0
α(x, S (t), t) u(x, t) d x, (2)

t > 0,

and an initial size-distribution for the population

u(x, 0) = u0(x), 0 ≤ x ≤ xM (0). (3)

Citation: Oscar Angulo, J.C. López-Marcos, M.A. López-Marcos, Numerical Analysis of a
Size-Structured Population Model with a Dynamical Resource, Biomath 3 (2014), 1403241,
http://dx.doi.org/10.11145/j.biomath.2014.03.241

Page 1 of 11

http://www.biomathforum.org/biomath/index.php/biomath
http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

The influence of the environment on the life
history of the individuals is given by a function
S (t) which represents a physiological resource. Its
dynamics is managed by the next initial value
problem,

S ′(t) = f (S (t), I(t), t), t > 0,
(4)

S (0) = S 0,

which is coupled with (1)-(3). The evolution of
the resource also depends on the population, which
is performed by means of the nonlocal term I(t)
defined by

I(t) =
∫ xM (t)

0
γ(x, S (t), t) u(x, t) d x, t ≥ 0. (5)

Also, we consider that the maximum size xM (t)
that an individual could have at time t, changes
with time and its dynamics is described by

d
dt

xM (t) = g(xM (t), S (t), t), t > 0,
(6)

xM (0) = xM.

The independent variables x and t represent size
and time, respectively. The dependent variable
u(x, t) is the size-specific density of individuals
with size x at time t. In general, the size of
any individual varies according to the following
ordinary differential equation

d x
dt

= g(x, S (t), t). (7)

Functions g, α and µ represent the growth,
fertility and mortality rates, respectively. These are
usually called the vital functions and define the
life history of an individual. Functions α and µ
are nonnegative. Note that all the vital functions
(g, µ and α) depend on size x (the structuring
internal variable), on time t and on the value
of the resource at time t, which can reflect the
influence of the environmental changes on the
vital functions. Function f on the right-hand side
of (4) depends on the value of the resource at
time t, on the total amount of individuals in the
population by means of the weighted functional
I(t) (which represents the way of weighting the
size distribution density in order to model the

different influence of individuals of different sizes
on such dynamics) and on time t.

The paper is structured as follows. In the next
section we introduce some background on theoret-
ical and numerical results. In section III, we de-
scribe the numerical method, which is completely
analyzed in section IV. Finally, numerical results
that confirm the expected order of convergence
and show the biological example dynamics are
included in section V.

II. Preliminary results

Theoretical analysis (1)-(5) which includes ex-
istence, uniqueness and long-time behaviour, is
highly difficult. A theoretical study of this model
appeared first time in a H.Thieme’s [7] presenta-
tion and authors in [5] also pointed that, for an
analysis of such models, we had to perform as it
was made in [8]. However, only a simpler model,
in which a positive growth was employed, was
analyzed in [9].

The numerical solution to the model, due to
its obvious mathematical complexity, entails a
serious challenge. In [10], a simple modification
of the scheme succesfully employed in [11], for
the solution of general size-structured population
models, was considered. The Daphnia magna test
in section V was also introduced. Furthermore, the
theoretical steady state of the model for this test
was provided. This scheme was shown not suitable
for a long-time integration with this biological test.
In this work, the numerical method described in
section III was proposed. Some simulations was
performed to show values of the parameters which
led us to an asymptotically stable equilibrium and
to an asymptotically stable periodic situation, in
both cases solutions were bounded. The purpose
of the current work is to validate such numerical
method by means of an analysis of its convergence.

On the other hand, a modification of this nu-
merical procedure was introduced in [12] in order
to approximate singular asymptotic states for these
kinds of models. In this work, we showed stable
and unstable steady state solutions. The study of
these singular states is not the aim of present work.

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 2 of 11

http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

In our convergence analysis, we shall use the
general discretization framework introduced by
López-Marcos et al. [13]. We present it rather
tersely, the reader is referred to this paper and the
references therein for a more critical and detailed
treatment. Thus, once we introduced a fixed given
problem concerning a differential equation, we
denote u a theoretical solution to such a problem.
We denote Ũh a numerical approximation to u. The
subscript h reflects the dependency on a parameter
h (mesh size) which takes values in a set H of pos-
itive numbers with inf H = 0. The approximation
is reached by solving the discretized problem

Φh(Ũh) = 0 (8)

(this family of discrete problems with h is referred
as discretization), where, for each h ∈ H, the
mapping Φh is fixed with domain Dh ⊂ Ah and
taking values in Bh. Here, Ah and Bh are vector
spaces with the same finite dimension. We further
assume that, for each h ∈ H, we have chosen a
norm in both spaces and an element ũh ∈ Dh which
is a suitable discrete representation of u in Ah.

Furthermore, we introduce the global and local
discretization error, ẽh and lh, respectively, and the
consistency, stability and convergence properties
of a discretization. The following result is crucial
and uses a deep topological lemma due to Stet-
ter [14],

Theorem 1. Assume that (8) is consistent and
stable with thresholds Rh. If Φh is continuous in
B(ũh, Rh) and ‖lh‖Bh = o(Rh) as h → 0, then:

i) For h sufficiently small, the discrete equations
(8) possess a unique solution in B(ũh, Rh).

ii) As h → 0, the solutions converge and
‖ẽh‖Ah = O(‖lh‖Bh ).

III. The numericalmethod

The numerical method we employ to approx-
imate the solution to (1)-(5) is based on the
discretization of a representation of the solution
along the characteristic curves [10]. First of all
we rewrite the partial differential equation (1) in
a more suitable form for its numerical treatment.

So we define

µ∗(x, z, t) = µ(x, z, t) + gx(x, z, t).

Thus equation (1) has the form

ut + g(x, S (t), t) ux = −µ
∗(x, S (t), t) u, (9)

0 < x < xM (t), t > 0.

We denote by x(t; t∗, x∗) the characteristic curve of
equation (9) that takes the value x∗ at time t∗. Such
a characteristic curve is the solution to the initial
value problem

d
dt

x(t; t∗, x∗) = g(x(t; t∗, x∗), S (t), t), t ≥ t∗,
(10)

x(t∗; t∗, x∗) = x∗.

Now we consider the function that represents
the solution to (9) along the characteristic curves
u(x(t; t∗, x∗), t), t ≥ t∗, which satisfies the initial
value problem

d
dt

u(x(t; t∗, x∗), t)=

−µ∗ (x (t; t∗, x∗) , S (t), t) u(x(t; t∗, x∗), t), t ≥ t∗,

u(x(t∗; t∗, x∗), t∗) = u(x∗, t∗),

and, therefore, it can be represented in the follow-
ing integral form

u(x(t; t∗, x∗), t) = (11)

u(x∗, t∗)exp
{
−

∫ t
t∗
µ∗(x (τ; t∗, x∗) , S (τ),τ) dτ

}
, t ≥ t∗.

Given a constant step k > 0, we introduce the
discrete time levels tn = n k, n = 0, 1, 2, . . .. We
also take J a positive integer, as a parameter
related to the size variable which describes the
number of points in the uniform initial grid. The
diameter of such a mesh grid is h = xM/J, and
the initial grid nodes are X0j = j h, 0 ≤ j ≤ J.
In order to start the integration, we consider as
an approximation to the density at initial time
(t0), the grid restriction of the initial condition
in (3), U 0j = u0(X

0
j ), 0 ≤ j ≤ J. Also, we use S

0

in (4) as the initial value of the resource. Then,
the numerical method provides, at each discrete
time level, a mesh grid on the size interval in
Xn, the approximation to the density on such a

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 3 of 11

http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

mesh grid in Un and the approximation of the
value of the resource S n, from the approximations
we computed at the previous time level, by using
discretizations of equations (10), (11), (2), (4), (5).
Thus, for n = 0, 1, 2, . . ., the numerical solution at
time tn+1 = tn + k, is obtained from the known
values of the numerical solution at time tn as
follows,

Xn+10 = 0, (12)

Xn+1j+1 = X
n
j + k g(X

n+1/2
j+1 , S

n+1/2, tn+1/2), (13)

0 ≤ j ≤ J, (14)

S n+1 = S n + (15)

k f (S n+1/2,Q(Xn+1/2,γn+1/2.Un+1/2), tn+1/2), (16)
U n+1j+1 = U

n
j exp

{
−k µ∗(Xn+1/2j+1 , S

n+1/2, tn+1/2)
}
, (17)

0 ≤ j ≤ J, (18)

U n+10 =
Q(Xn+1,αn+1.Un+1)
g(Xn+10 , S

n+1, tn+1)
, (19)

where we have to compute approximations at time
level tn+1/2 = tn + k/2,

Xn+1/20 = 0,

Xn+1/2j+1 = X
n
j +

k
2

g(Xnj , S
n, tn), 0 ≤ j ≤ J,

S n+1/2 = S n +
k
2

f (S n,Q(Xn,γn.Un), tn),

U n+1/2j+1 = U
n
j exp

{
−

k
2
µ∗(Xnj , S

n, tn)
}
, 0 ≤ j ≤ J,

U n+1/20 =
Q(Xn+1/2,αn+1/2.Un+1/2)
g(Xn+1/20 , S

n+1/2, tn+1/2)
.

Note that the general step of the method increases
the number of grid points and also the dimension
of the vector with the numerical densities: at time
tn, we have J + 1 grid nodes in Xn and the (J + 1)-
dimensional vector Un, and at time tn+1 we obtain
J +2 grid nodes in Xn+1 and the (J +2)-dimensional
vector Un+1. In order to maintain the number of
grid points suitable to perform the next step, we
eliminate at time tn+1 the first grid node Xn+1l
which satisfies

|Xn+1l+1 − X
n+1
l−1 | = min1≤ j≤J+1

|Xn+1j+1 − X
n+1
j−1 |. (20)

We reproduce the same reduction in the corre-
sponding vector Un+1.

In the description of the method, we use the
following notation; vectors αp and γp contain the
evaluations of the functions α and γ in (2) and
(5), respectively, at the grid points in Xp, at the
resource value S p and at time t p. Products γp.Up
and αp.Up must be considered componentwise.
In order to approximate integrals over the inter-
val [0, xM (t p)], we use the composite trapezoidal
quadrature rule based on the grid points Xp =
[X p0 , X

p
1 , . . . , X

p
J ], that is

Q(Xp, Vp) =
J∑

j=1

X pj − X
p
j−1

2

(
V pj−1 + V

p
j

)
. (21)

Note that the method is implicit: all the expres-
sions provide explicit equations for the numerical
values at the highest time level, except those which
involve the numerical density U p0 at the first grid
point, but it is easy to implement the method in
an explicit form.

IV. Convergence Analysis

Below, we will analyze numerical methods
based on the integration along characteristics that
use a general quadrature rule with suitable proper-
ties to approximate the integral terms. The proofs
of every result are heavily laborious and they
would be included in a more technical work.

We assume the following regularity conditions
on the data functions and the solution to the
problem (1)-(5):

(H1) u ∈C2([0, xM (t)]× [0, T ]), u(x, t) ≥ 0, x ∈
[0, xM (t)], t ≥ 0.

(H2) S ∈C2([0, T ]), S (t) ≥ 0, t ≥ 0.
(H3) γ ∈C2([0, xM (t)]×D×[0, T ]), where D is a

compact neighbourhood of {S (t) , 0 ≤ t ≤ T}.

(H4) µ ∈C2([0, xM (t)]× D×[0, T ]), is nonneg-
ative and D is a compact neighbourhood of
{S (t) , 0 ≤ t ≤ T}.

(H5) α ∈C2([0, xM (t)]×D×[0, T ]), is nonneg-
ative and D is a compact neighbourhood of
{S (t) , 0 ≤ t ≤ T}.

(H6) f ∈ C2(D × DI × [0, T ]), is nonneg-
ative, D is a compact neighbourhood of

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 4 of 11

http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

{S (t) , 0 ≤ t ≤ T}, and DI is a compact neigh-
bourhood of{∫ xM (t)

0
γ(x, S (t)) u(x, t) d x, 0 ≤ t ≤ T

}
.

(H7) g ∈C3([0, xM (t)]×D×[0, T ]), where D is a
compact neighbourhood of {S (t) , 0 ≤ t ≤ T}
and g(0, s, t) ≥ C > 0, t ≥ 0, s ∈ R. In ad-
dition, the characteristic curves x(t; t∗, x∗) de-
fined in (7) are continuous and differentiable
with respect to the initial values (t∗, x∗) ∈
[0, T ] × [0, xM (t)].

The above hypotheses may be based on three pos-
sible reasons. First, biological assumptions such
as the nonnegativity of some of the vital functions
or, in (H7), to reflect that any individual in the
studied population could shrink [2]. Second, the
mathematical requirements to obtain the existence
and uniqueness of solutions for the problem (1)-
(5) [2]. Finally, the regularity properties needed in
the numerical analysis to derive optimal rates of
convergence [11].

We also assume that the spatial discretization
parameter, h, takes values in the set H = {h > 0 :
h = xM/J, J ∈N}. Now, we suppose that the time
step, k, satisfies k = r h, where r is an arbitrary and
positive constant, fixed throughout the analysis. In
addition, we set N = [T/k]. For each h ∈ H, we
define the spaces

Ah =

N∏
n=0

(
RJ+n ×RJ+n+1

)
×RN+1,

Bh =
(
RJ×RJ+1×R

)
×RN ×

N∏
n=1

(
RJ+n ×RJ+n

)
×RN.

Both spaces have the same dimension.
In order to measure the size of the errors, we

define ‖η‖∞ = max1≤ j≤p |η j|, η ∈ Rp, ‖Vn‖1 =∑J+n
j=0 h |V

n
j |, V

n ∈ RJ+n+1. Thus, we endow the
spaces Ah and Bh with the following norms. If(
y0, V0, . . . , yN, VN, a

)
∈Ah, then

‖
(
y0, V0, . . . , yN, VN, a

)
‖Ah =

max
(

max
0≤n≤N

‖yn‖∞, max
0≤n≤N

‖Vn‖∞,‖a‖∞
)
.

On the other hand, if(
Y0, Z0, A0, Z0, Y1, Z1, . . . , YN, ZN, A

)
∈Bh,

‖
(
Y0, Z0, A0, Z0, Y1, Z1, . . . , YN, ZN, A

)
‖Bh

=‖Y0‖∞ + ‖Z0‖∞ + |A0| + ‖Z0‖∞ +
N∑

n=1

k ‖Zn‖∞

+

N∑
n=1

k ‖Yn‖∞ +
N∑

n=1

k |An|.

Now, for each h ∈ H, we define xh =
(x0, x1, x2, . . . , xN ), xn = (xn1, . . . , x

n
J+n) ∈ R

J+n,
x0j = j h, 1 ≤ j ≤ J,

xnj = x(t
n; tn−1, xn−1j−1 ), 1 ≤ j ≤ J + n, 1 ≤ n ≤ N.

(22)
Also, xn+

1
2 = (x

n+ 12
1 , . . . , x

n+ 12
J+n+1) ∈R

J+n+1,

x
n+ 12
j = x(t

n+ 12 ; tn, xnj−1), 1 ≤ j ≤ J + n + 1, (23)

0 ≤ n ≤ N − 1. We denote xn0 = x
n+ 12
0 = 0, n ≥ 0.

In addition, if u represents the theoretical solution
to (1)-(5) we define uh = (u0, u1, u2, . . . , uN ), un =
(un0, u

n
1, . . . , u

n
J+n) ∈R

J+n+1,

unj = u(x
n
j, t

n), 0 ≤ j ≤ J + n, 0 ≤ n ≤ N, (24)

and un+
1
2 = (u

n+ 12
0 , u

n+ 12
1 , . . . , u

n+ 12
J+n+1) ∈R

J+n+2,

u
n+ 12
j = u(x

n+ 12
j , t

n), 0 ≤ j ≤ J + n + 1, 0 ≤ n ≤ N−1.
(25)

Finally, if S is the theoretical solution to (4) then
we define sh = (s0, s1, s2, . . . , sN ),

sn = S (tn), 0 ≤ n ≤ N, (26)

and
sn+

1
2 = S (tn+

1
2 ), 0 ≤ n ≤ N − 1. (27)

Therefore ũh = (x0, u0, x1, u1, . . . , xN, uN, sh) ∈
Ah.

Next, we introduce the discretization operator.
Let R be a positive constant and we denote by
BAh (ũh, R h

p) ⊂ Ah the open ball with center ũh
and radius R hp, 1 < p < 2,

Φh : BAh (ũh, R h
p) →Bh,

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 5 of 11

http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

Φh
(
y0, V0, . . . , yN, VN, a

)
=

(
Y0, P0, A0, P0, Y1, P1, . . . , YN, PN, A

)
, (28)

defined by the following equations:

Y0 = y0 − X0 ∈RJ, (29)

P0 = V0 − U0 ∈RJ+1, (30)

A0 = a0 − S 0 ∈R. (31)

Vectors X0, U0 and value S 0 represent approxi-
mations at t = 0, respectively, to the initial grid
nodes, to the theoretical solution at these points
and to the initial resource. Also,

Pn+10 = V
n+1
0 −

Q
(
yn+1,αn+1 · Vn+1

)
g
(
0, an+1, tn+1

) , (32)
Y n+1j+1 =

1
k

{
yn+1j+1

−ynj − k g(y
n+ 12 ,∗
j+1 , a

n+ 12 ,∗, tn+
1
2 )
}
, (33)

Pn+1j+1 =
1
k

{
V n+1j+1

(34)

− V nj exp
(
−k µ∗

(
y

n+ 12 ,∗
j+1 , a

n+ 12 ,∗, tn+
1
2

))}
, (35)

0 ≤ j ≤ J + n − 1,

An+1 =
1
k

{
an+1 − an

−k f
(
an+

1
2 ,∗,Q(yn+

1
2 ,∗,γn+

1
2 ,∗ · Vn+

1
2 ,∗), tn+

1
2

)}
, (36)

0 ≤ n ≤ N−1. Where, with the notation introduced
in Section III,

y
n+ 12 ,∗
j+1 = y

n
j +

k
2

g(ynj, a
n, tn), (37)

V
n+ 12 ,∗
j+1 = V

n
j exp

(
−

k
2
µ∗

(
ynj, a

n, tn
))
, (38)

0 ≤ j ≤ J + n − 1,

V
n+ 12 ,∗
0 =

Q(yn+
1
2 ,∗,αn+

1
2 ,∗ · Vn+

1
2 ,∗)

g(0, an+
1
2 ,∗, tn+

1
2 )

, (39)

an+
1
2 ,∗ = an +

k
2

f (an,Q(yn,γn · Vn), tn) , (40)

0 ≤ n ≤ N − 1. We denote by Q(X, V) =
M∑

l=0

ql(X) Vl the general quadrature rule employed

in (32)-(40). Note that, Φh takes into account
all the possible nodes and their corresponding
solution values at each time level, and it em-
ploys quadrature rules possibly based on a sub-
grid. If Ũh = (X0, U0, X1, U1, . . . , XN, UN, S) ∈
BAh (ũh, R h

p), satisfies

Φh(Ũh) = 0 ∈Bh, (41)

the nodes and the corresponding values of the
solution at such nodes of Ũh are a numerical
solution to the scheme defined by (13)-(19) when
the composite trapezoidal quadrature rule is given.
On the other hand, the numerical solution of the
scheme defined by (13)-(19) satisfies (41).

Henceforth, C will denote a positive constant,
independent of h, k (k = r h), j (0 ≤ j ≤ J + n)
and n (0 ≤ n ≤ N); C may have different
values in different places. Now, we assume that the
quadrature rules satisfy the following properties:

(P1) |I(tn) −Q (xn,γn · un)| ≤ C h2, when h →
0, 0 ≤ n ≤ N.

(P2)
∣∣∣∣I(tn−12 ) −Q(xn−12 ,γn−12 · un−12 )∣∣∣∣ ≤ C h2,

when h → 0, 1 ≤ n ≤ N.

(P3)

∣∣∣∣∣∣
∫ xM (tn )

0
α(x, S (tn), tn) u(x, tn) d x −Q (xn,αn · un)

∣∣∣∣∣∣
≤ C h2, when h → 0, 0 ≤ n ≤ N.

(P4)

∣∣∣∣∣∣∣∣
∫ xM (tn− 12 )

0
α(x, S (tn−

1
2 ), tn−

1
2 ) u(x, tn−

1
2 ) d x

−Q
(
xn−

1
2 ,αn−

1
2 · un−

1
2

)∣∣∣∣ ≤ C h2, when h → 0,
1 ≤ n ≤ N.

(P5) |q j(xn)| ≤ q h, where q is a positive
constant independent of h, k, j (0 ≤ j ≤ J + n)
and n (0 ≤ n ≤ N), for 0 ≤ j ≤ J + n,
0 ≤ n ≤ N.

(P6) |q j(xn−
1
2 )| ≤ q h, where q is a positive

constant independent of h, k, j (0 ≤ j ≤ J + n)
and n (0 ≤ n ≤ N), for 0 ≤ j ≤ J + n,
1 ≤ n ≤ N.

(P7) Let R and p be positive constants with
1 < p < 2. The quadrature weights q j are Lip-
schitz continuous functions on B∞(xn, R hp),

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 6 of 11

http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

0 ≤ j ≤ J + n, 1 ≤ n ≤ N and on
B∞(xn−

1
2 , R hp), 0 ≤ j ≤ J + n, 1 ≤ n ≤ N.

(P8) Let R and p be positive constants with
1 < p < 2. If yn, zn ∈ B∞(xn, R hp), Vn ∈
B∞(un, R hp) and an ∈ B∞(sn, R hp), then∣∣∣∣∣∣∣

J+n∑
i=0

(qi(yn) − qi(zn)) γ(zni , a
n, tn) V ni

∣∣∣∣∣∣∣ ≤
C‖yn − zn‖∞,

when h → 0, 0 ≤ n ≤ N.
(P9) Let R and p be positive constants with

1 < p < 2. If yn, zn ∈ B∞(xn, R hp), Vn ∈
B∞(un, R hp) and an ∈ B∞(sn, R hp), then∣∣∣∣∣∣∣
J+n∑
i=0

(qi(yn) − qi(zn)) α
(
zni , a

n, tn
)

V ni

∣∣∣∣∣∣∣ ≤
C‖yn − zn‖∞,

when h → 0, 0 ≤ n ≤ N.
(P10) Let R and p be positive constants with 1 <

p < 2. If yn−
1
2 , zn−

1
2 ∈ B∞(xn−

1
2 , R hp), Vn−

1
2 ∈

B∞(un−
1
2 , R hp), and an−

1
2 ∈ B∞(sn−

1
2 , R hp),

then∣∣∣∣∣∣∣
J+n+1∑

i=0

(
qi(yn−

1
2 ) − qi(zn−

1
2 )
)
γ
(
z

n−12
i , a

n−12 , tn−
1
2

)
V

n−12
i

∣∣∣∣∣∣∣
≤ C ‖yn−

1
2 − zn−

1
2‖∞,

when h → 0, 1 ≤ n ≤ N.
(P11) Let R and p be positive constants with 1 <

p < 2. If yn−
1
2 , zn−

1
2 ∈ B∞(xn−

1
2 , R hp), Vn−

1
2 ∈

B∞(un−
1
2 , R hp), and an−

1
2 ∈ B∞(sn−

1
2 , R hp),

then∣∣∣∣∣∣∣
J+n+1∑

i=0

(
qi(yn−

1
2 )− qi(zn−

1
2 )
)
α
(
z

n−12
i , a

n−12 , tn−
1
2

)
V

n−12
i

∣∣∣∣∣∣∣
≤ C‖yn−

1
2 − zn−

1
2‖∞,

when h → 0, 1 ≤ n ≤ N.

These are the enough properties the quadrature
rules have to satisfy to carry out our convergence
analysis. The following result establishes that the
composite trapezoidal rule used in our experiments
satisfies them.

Theorem 2. Assume that the hypotheses (H1)-
(H7) hold. If the quadrature rules are the compos-

ite trapezoidal quadrature on subgrids
{
xnjnl

}M(n)
l=0

,
0 ≤ n ≤ N with the property

(SR) There exists a positive constant C
such that, for h sufficiently small, xnjnl+1

− xnjnl
≤ C h,

0 ≤ l ≤ M(n) − 1, xnjn0
= 0, xnjnM(n)

= xJ+n, with{
xnjnl

}M(n)−1
l=1

contained in xn, 0 ≤ n ≤ N.

Then, properties (P1)-(P11) hold.

Now we introduce the following result over the
numerical values at the half-level time.

Proposition 1. Assume that the hypotheses (H1)-
(H7) hold and that the considered quadra-
ture rules satisfy properties (P1)-(P11). Let be
yn ∈ B∞(xn, R hp), Vn ∈ B∞(un, R hp) and
an ∈ B∞(sn, R hp). Then, as h → 0, yn+

1
2 ,∗ ∈

B∞(xn+
1
2 , R′ hp), an+

1
2 ,∗ ∈ B∞(sn+

1
2 , R′ hp), and

Vn+
1
2 ,∗ ∈ B∞(un+

1
2 , R′ hp) where xn+

1
2 , sn+

1
2 , and

un+
1
2 are defined by (23), (27) and (25), respec-

tively.

Now, we define the local discretization error as

lh = Φh(ũh) ∈Bh,

and we say that the discretization (28) is consis-
tent if, as h → 0,

lim‖Φh(ũh)‖Bh = lim‖lh‖Bh = 0.

The following theorem establishes the consis-
tency of the numerical scheme defined by equa-
tions (29)-(36).

Theorem 3. Assume that hypotheses (H1)-(H7)
hold and that the considered quadrature rules
satisfy properties (P1)-(P11). Then, as h → 0, the
local discretization error satisfies,

‖Φh(ũh)‖Bh = ‖u
0
− U0‖∞ + ‖x0 − X0‖∞

+ |s0 − S 0| + O(h2 + k2). (42)

Another notion that plays an important role in
the analysis of the numerical method is the sta-
bility with h-dependent thresholds. For h ∈ H, let
Rh be a real number ( the stability threshold) with

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 7 of 11

http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

0 < Rh < ∞: we say that the discretization (28)
is stable for ũh restricted to the thresholds Rh, if
there exist two positive constants h0 and S ( the
stability constant) such that, for any h ∈ H with
h ≤ h0, the open ball BAh (ũh, Rh) is contained in
the domain of Φh, and, for all Ṽh, W̃h in that ball,

‖Ṽh − W̃h‖ ≤ S ‖Φh(Ṽh) −Φh(W̃h)‖.

Below, we introduce the theorem that estab-
lishes the stability of the discretization defined by
equations (29)-(36).

Theorem 4. Assume that hypotheses (H1)-(H7)
hold and that the considered quadrature rules sat-
isfy properties (P1)-(P11). Then, the discretization
is stable for ũh with Rh = R hp, 1 < p < 2.

Finally, we define the global discretization error
as

ẽh = ũh − Ũh ∈Ah,

We say that the discretization (28) is convergent
if there exists h0 > 0 such that, for each h ∈ H
with h ≤ h0, (41) has a solution Ũh for which, as
h → 0,

lim‖ũh − Ũh‖Ah = lim‖ẽh‖Ah = 0.

We propose the following theorem which estab-
lishes the convergence of the numerical method
defined by equations (29)-(36).

Theorem 5. Assume that hypotheses (H1)-(H7)
hold and that the considered quadrature rules
satisfy properties (P1)-(P11). Then, for h suffi-
ciently small, the numerical method defined by
equations (29)-(36) has a unique solution Uh ∈
B(uh, Rh) and

‖Uh − uh‖Ah ≤ C
(
‖x0 − X0‖∞ + ‖u0 − U0‖∞

+|s0 − S 0| + O(h2 + k2)
)
. (43)

The proof of Theorem 5 is immediately derived
by means of the consistency (Theorem 3), the
stability (Theorem 4) and Theorem 1.

Next, we can establish an error bound for the
the numerical and the theoretical solution at the
numerical values of the grid nodes.

Theorem 6. Assume that hypotheses (H1)-(H7)
hold and that the considered quadrature rules
satisfy properties (P1)-(P11). For h sufficiently

small, let be u∗h = (u
0
∗, u

1
∗, u

2
∗, . . . , u

N
∗ ) ∈

N∏
n=0

RJ+n+1,

defined by

un∗ =
(
u(Xn0, t

n), u(Xn1, t
n), . . . , u(XnJ+n, t

n)
)
∈RJ+n+1,

0 ≤ n ≤ N, where Xnj , 0 ≤ j ≤ J + n, 0 ≤ n ≤ N,
are the grid nodes given by the scheme (29)-(36).
Then,

‖Un − un∗‖∞ ≤ C
(
‖x0 − X0‖∞ + ‖u0 − U0‖∞

+|s0 − S 0| + O(h2 + k2)
)
. (44)

This theorem follows immediately from The-
orem 5. Specifically, if X0 = x0, U0 = u0 and
S 0 = s0, the proposed numerical scheme is second-
order accurate.

At this moment, we have obtained convergence
of the numerical method (29)-(36) which does
not employ selection at each time level. Also, we
have proven the convergence of numerical methods
which employ a selection criterion, whenever the
positions, which are determined by the criterion
we have chosen, lead us to subgrids which satisfy
property (SR). For the criterion presented in this
paper, this property may be shown in two stages.
First, as proved in [11], it leads us to subgrids with
such a property when we applied it over nodes
which are in a neighbourhood of the theoretical
ones with radius R hp. In a second stage, it is
proven that the nodes, which in fact the numerical
method computes, are in such neighbourhoods. In
order to do this, it is enough to realize that such
nodes could be seen, up to each time level, as the
solutions obtained by a discrete operator which has
the form of that defined in (28).

V. Numerical results

We have carried out different numerical experi-
ments with the scheme defined in Section III. We
have considered a theoretical test problem with
meaningful nonlinearities (both from a mathemat-
ical and biological point of view). The numeri-
cal integration for the numerical experiment was

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 8 of 11

http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

carried out on the time interval [0, 10]. The size
interval was taken as [0, 1].

The size-specific growth, fertility
and mortality moduli are chosen as
g(x, z, t) = λ2

1+z
z

((
z

1+z

)2
− x2

)
+ xr1+z

(
29
30 −

z
k

)
,

α(x, z, t) = 32λ
1+

(
z

C( 2930 − zk )

)−29λ
30r

1+2
(

z
C( 2930 − zk )

)−29λ
30r

, µ(x, z, t) =

λ 1+zz

(
z

1+z + 2x
)
−

3r
1+z

(
29
30 −

z
k

)
. The weight

function is taken as γ(x, z, t) = x2 and, finally,
f (z, i, t) = rz

(
1 − zk

)
− rzi (1+z)

5

z5 (1+4e−λt ) . With this
choice of data functions, the problem (1)-(5) has
the following solution

u(x, t) =
(

S (t)
1 + S (t)

− x
)2
− e−λt


(

S (t)
1 + S (t)

)2
− x2

,
S (t) =

29
30

Ce29rt/30

1 + Ce29rt/30/k
,

r = 0.1, C = 24, k = 5, λ = 0.3. Since we
know the exact solution to the problem, we can
show numerically that our method is second-order
accurate by means of an error Table. In Table I,
each entry in columns two to five represents, at
the upper value, the global error

eh,k = max
{

max
0≤ j≤J

|u(X0j , t
0) − U 0j |, |S (t

0) − S 0|,

max
1≤n≤N

{
max
0≤ j≤J

|u(Xnj , t
n) − U nj |

}
, |S (tn) − S n|

}
and, at the lower number, the experimental order s
of the method as computed from s = log (e2h,2k/eh,k )log 2 .
Each column and each row of the Table corre-
spond to different values of the spatial and time
discretization parameter, respectively. The results
in the Table clearly confirm the expected second-
order convergence.

On the other hand, the numerical integration
of the model with an efficient method allows us
to consider a more realistic test problem. The
property of convergence in finite time interval is
important to carry out experiments in which the
long-time behaviour of the population is investi-
gated. Therefore, the numerical method has been
employed to describe the dynamics of a population

k\h 1.25e-2 6.25e-2 3.13e-2 1.56e-2
1.25e-2 1.67e-4 1.27e-4 1.22e-4 1.21e-4

6.25e-2 2.08e-4 4.18e-5 3.15e-5 3.03e-5
2.00 2.01 2.01

3.13e-2 2.02e-4 5.26e-5 1.05e-5 7.85e-6
1.98 2.00 2.00

1.56e-2 2.08e-4 5.05e-5 1.32e-5 2.62e-6
2.00 1.99 2.00

TABLE I
Error and experimental order of convergence.

of ectothermic invertebrates. This is the case of
the water flea, Daphnia magna. In this particular
case, the functions data are given by g(x, S, t) =
g( S1+S − x), µ(x, S, t) = µ, α(x, S, t) = α

S
1+S x

2,
f (S, i, t) = rS

(
1 − SK

)
− i S1+S , γ(x, S, t) = x

2,
and the values of the parameters are given by
g = 1, µ = 0.1, α = 0.75, r = 3, K = 8.3
and xM = 1 [6]. This set led us to unbounded
solutions [12]. Nevertheless, for the parameter
value g = 0.0075, the solution is bounded. As it
was pointed in [10], as the value of the parameter
K increases, the equilibrium state becomes unsta-
ble. We have performed a numerical experiment
with the discretization parameters k = 0.0625,
J = 4000 and the interesting value K = 9.64. In
this case, we observe the unstability of the equi-
librium and the solution evolving towards a cycled
situation (Figure 1). Taking into account that the
numerical solution is attracted to a limit cycle,
considering a sufficiently large time, the numerical
solution obtained after this long time integration
lies practically on such a cycle. In this way, the
numerical method provides an approximation to
the limit cycle. In Figure 2, the representation of
such a cycle in the tridimensional space defined by
the total population, the maximum individual size
and the dynamical resource, is drawn. From the
numerical results obtained for the total population,
the maximum size and the dynamical resource, we
can obtain an in depth analysis of these quantities
throughout a period of the limit cycle. For exam-
ple, we have estimated the period of the solution
by interpolation and it is about 64.6824.

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 9 of 11

http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

Fig. 1. Evolution of the numerical solution. Case of a
bounded unstable steady state.

46

48

50

52

54

0.76
0.77

0.78
0.79

0.8
0.81

0.82
0.83

0.84

2

3

4

5

6

maximum size
total population

re
s
o
u
rc

e

Fig. 2. Limit cycle appearing in the case of a bounded
unstable steady state.

VI. Conclusions

We have analyzed a second-order numerical
method for a problem that describes a population
with a possible shrinking size and with a depen-
dency on the environment managed by the evo-
lution of a vital resource. The second-order con-
vergence has been theoretically proven by means
of an argument of consistency and stability of the
scheme. We have reported numerical experiments
which demonstrate the predicted accuracy of the
scheme.

This knowledge leads us to employ this second-

order numerical method that, experimentally,
shows good stability properties for the study of
the behaviour of a real population. We have con-
sider the long time numerical approximation of a
particular model which describes the dynamics of
a Daphnia magna population. Moreover, when the
steady state is unstable we have used it to analyzed
the dynamics, appearing a limit cycle, and a good
approximation to the bifurcation process has been
obtained.

Acknowledgments

This work was supported in part by the Min-
isterio de Ciencia e Innovación (Spain), project
MTM2011-25238, and by the Junta de Castilla y
León (Spain), project VA191U13.

References

[1] G. F. Webb, Theory of Nonlinear Age-Dependent Popu-
lation Dynamics, Marcel Dekker, eds, New York, 1985.

[2] J. A. J. Metz and E. O. Dieckmann, editors, The Dynam-
ics of Physiologically Structured Populations, Springer
Lecture Notes in Biomathematics, 68, Springer, Heildel-
berg, 1986.

[3] M. Iannelli, Mathematical Theory of Age-Structured Pop-
ulation Dynamics, Applied Mathematics Monographs.
C.N.R., Giardini Editori e Stampatori, Pisa, 1995.

[4] J. M. Cushing. An Introduction to Structured Populations
Dynamics, CMB-NSF Regional Conference Series in
Applied Mathematics. SIAM, 1998.

[5] B. Perthame, Transport Equations in Biology, Birkhäuser
Verlag, Basel, 2007.

[6] A. M. de Roos, Numerical methods for structured popula-
tion models: The escalator boxcar train, Numer. Methods
Partial Differential Equations, 4 (1988), 173–195.

[7] H. Thieme,
http://math.la.asu.edu/∼dieter/workshop/schedule.html,
2003.

[8] A. Calsina and J. Saldaña. A model of physiologically
structured population dynamics with a nonlinear individ-
ual growth rate, J. Math. Biol., 33(4) (1995) 335–364.

[9] O. Diekmann, M. Gyllenberg, J.A.J. Metz, S. Nakaoka,
A.M. de Roos, Daphnia revisited: local stability and bi-
furcation theory for physiologically structured population
models explained by way of an example, J. Math. Biol.
61 (2010), 277-318.

[10] L. Abia, O. Angulo, J. C. López-Marcos and
M. A. López-Marcos, Long-Time Simulation of a Size-
Structured Population Model with a Dynamical Resource,
Math. Model. Nat. Phenom., 5 (2010), 1–21.

[11] O. Angulo and J.C. López-Marcos. Numerical integra-
tion of fully nonlinear size-structured population models,
App. Num. Math., 50 (2004) 291-327.

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 10 of 11

http://math.la.asu.edu/~dieter/workshop/schedule.html
http://dx.doi.org/10.11145/j.biomath.2014.03.241


O Angulo et al., Numerical Analysis of a Size-Structured Population Model...

[12] O. Angulo, J. C. López-Marcos and M. A. López-
Marcos, Numerical approximation of singular asymptotic
states for a size-structured population model with a dy-
namical resource, Math. Comput. Modelling, 54 (2011),
1693–1698.

[13] J. C. López-Marcos and J. M. Sanz-Serna, Stability
and convergence in numerical analysis III: Linear inves-
tigation of nonlinear stability, IMA J. Numer. Anal., 8
(1988), 71–84.

[14] H. Stetter, Analysis of discretization methods for ordi-
nary differential equations, Springer, Berlin, 1973.
http://dx.doi.org/10.1007/978-3-642-65471-8

[15] L.M. Abia, O. Angulo, J. C. López-Marcos and
M. A. López-Marcos, Numerical integration of a hier-
archically size-structured population model with contest
competition, J. Comp. Appl. Math., 258 (2014), 116–134.

Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 11 of 11

http://dx.doi.org/10.1007/978-3-642-65471-8
http://dx.doi.org/10.11145/j.biomath.2014.03.241

	Introduction
	Preliminary results
	The numerical method
	Convergence Analysis
	Numerical results
	Conclusions
	References