Numerical study of size-structured population models: A case of <i>Gambussia affinis</i>

Oscar Angulo; Angel Durán; Juan Carlos López-Marcos

doi:10.1016/j.crvi.2004.11.007

Biological modelling / Biomodélisation

Numerical study of size-structured population models: A case of Gambussia affinis
[Étude numérique de modèles de population structurés en taille : le cas Gambussia affinis]

Présenté par : Pierre Auger

Oscar Angulo ^1,² ; Angel Durán ² ; Juan Carlos López-Marcos ²

¹ C/Fco. Mendizabal 1, 47014 Valladolid, Spain
² Departamento de Matemática Aplicada, Universidad de Valladolid, 47014 Valladolid, Spain

Comptes Rendus. Biologies, Volume 328 (2005) no. 4, pp. 387-402.

Résumés

Anglais
Français

We study, from a numerical point of view, some properties of a model which describes the evolution of a population of Gambussia affinis. Our model includes sufficiently smooth vital functions. First we select, among four numerical methods of second order, the most appropriate in terms of adaptation to the problem. The most efficient method also reveals new properties of the model for long times, such as the tendency to periodicity, obtained with different initial conditions. We also discuss some advantages and deficiencies of the model.

Nous étudions, d'un point de vue numérique, quelques propriétés d'un modèle qui décrit l'évolution d'une population de Gambussia affinis. Notre modèle inclut des fonctions vitales suffisamment lisses. Nous sélectionnons d'abord, parmi quatre méthodes numériques du deuxième ordre, celle qui est la plus appropriée en termes d'adaptation au problème. La méthode la plus efficace révèle aussi de nouvelles propriétés du modèle aux temps longs, telles qu'une tendance à la périodicité, obtenue avec différentes conditions initiales. Nous discutons aussi quelques avantages et déficiences du modèle.

Métadonnées

Reçu le : 2004-02-20
Accepté le : 2004-11-15
Publié le : 2005-02-01

PMID

DOI : 10.1016/j.crvi.2004.11.007

Keywords: Size-structured population, Numerical integration, Efficiency, Long time integration, Gambussia affinis
Mot clés : Population structurée en taille, Intégration numérique, Efficacité, Intégration aux temps longs, Gambussia affinis

Affiliations des auteurs :

Oscar Angulo ^{1, 2} ; Angel Durán ² ; Juan Carlos López-Marcos ²

¹ C/Fco. Mendizabal 1, 47014 Valladolid, Spain
² Departamento de Matemática Aplicada, Universidad de Valladolid, 47014 Valladolid, Spain

@article{CRBIOL_2005__328_4_387_0,
     author = {Oscar Angulo and Angel Dur\'an and Juan Carlos L\'opez-Marcos},
     title = {Numerical study of size-structured population models: {A} case of {\protect\emph{Gambussia} affinis}},
     journal = {Comptes Rendus. Biologies},
     pages = {387--402},
     publisher = {Elsevier},
     volume = {328},
     number = {4},
     year = {2005},
     doi = {10.1016/j.crvi.2004.11.007},
     language = {en},
}

TY  - JOUR
AU  - Oscar Angulo
AU  - Angel Durán
AU  - Juan Carlos López-Marcos
TI  - Numerical study of size-structured population models: A case of Gambussia affinis
JO  - Comptes Rendus. Biologies
PY  - 2005
SP  - 387
EP  - 402
VL  - 328
IS  - 4
PB  - Elsevier
DO  - 10.1016/j.crvi.2004.11.007
LA  - en
ID  - CRBIOL_2005__328_4_387_0
ER  -

%0 Journal Article
%A Oscar Angulo
%A Angel Durán
%A Juan Carlos López-Marcos
%T Numerical study of size-structured population models: A case of Gambussia affinis
%J Comptes Rendus. Biologies
%D 2005
%P 387-402
%V 328
%N 4
%I Elsevier
%R 10.1016/j.crvi.2004.11.007
%G en
%F CRBIOL_2005__328_4_387_0

Oscar Angulo; Angel Durán; Juan Carlos López-Marcos. Numerical study of size-structured population models: A case of Gambussia affinis. Comptes Rendus. Biologies, Volume 328 (2005) no. 4, pp. 387-402. doi : 10.1016/j.crvi.2004.11.007. https://comptes-rendus.academie-sciences.fr/biologies/articles/10.1016/j.crvi.2004.11.007/

Version originale du texte intégral

1 Introduction

The purpose of this paper is to analyse, from a numerical point of view, some aspects in the evolution of a population of Gambussia affinis or mosquitofish as an example of a size-structured population model. We show how to study the dynamics of a model from numerical techniques. First, it is appropriate to use different numerical integrators and compare the corresponding results. Such comparisons allow us to achieve some certainty regarding different questions in the dynamics of the population. From this, we can choose the numerical method which shows the best adaptation to the problem. With the help of the most efficient scheme we study the evolution of the population in order to evaluate the values and deficiencies of the model proposed. One of the relevant properties discovered is that, with appropriate vital functions, the population tends to be periodic and stable.

One of the ways to study the dynamics of species in biology or ecology is based on the so-called structured models. In these models, one or more internal variables of the individuals, like age, size, etc., allow us to show the effect of the physiological state of the individuals on the population dynamics. At the same time, these models incorporate vital functions like mortality, fertility, etc., generally in nonlinear form (with nonlocal terms). The way these vital functions depend on the structural variables and density of population determines the adaptation of the model to the real behaviour of the species [1–4].

In many cases, size is the individual parameter which has more influence on the evolution of the population. The ability of individuals to survive and reproduce depends on it. This gave rise to the consideration of the so-called size-structured mathematical models. The first models of this type were proposed in [5,6] and generalized the well-known age-structured models [7,8]. A general size-structured population model can be described by the following equations:

u_{t} + {(g (x, I_{g} (t), t) u)}_{x} = - μ (x, I_{μ} (t), t) u x_{\min} < x < x_{\max}, t > 0

(1.1)

g (x_{\min}, I_{g} (t), t) u (x_{\min}, t) = \int_{x_{\min}}^{x_{\max}} α (x, I_{α} (t), t) u (x, t) d x, t > 0

(1.2)

u (x, 0) = u_{0} (x), x_{\min} ⩽ x ⩽ x_{\max}

(1.3)

I_{ϕ} (t) = \int_{x_{\min}}^{x_{\max}} γ_{ϕ} (x) u (x, t) d x t ⩾ 0, ϕ = α, μ, g

(1.4)

The independent variables x, t represent, respectively, size and time, with x in an interval

[x_{\min}, x_{\max}]

. The dependent variable

u (x, t)

represents the density of individuals in the population with size x at time t. We assume that the size of any individual evolves according to the ODE

\frac{d x}{d t} = g (x, I_{g} (t), t)

The so-called vital functions of the population, which determine the life history of an individual, are given by the functions g, α and μ, representing, respectively, growth, fertility and mortality rates. All the functions in the model are nonnegative and depend on size and time. This time dependence implies a possible seasonal behaviour of the population. Finally, in order to consider the influence of individuals with different size on life conditions, the vital functions can also depend on the total number of individuals. This is indicated by the integrals (1.4), $I_{g} (t)$ , $I_{α} (t)$ and $I_{μ} (t)$ , which average the density of the population by using specified functions $γ_{g}$ , $γ_{α}$ and $γ_{μ}$ .

Eq. (1.1) is a PDE with nonlocal terms representing the balance law of the population. It describes the evolution of the number of individuals with a given size due to growth and mortality. Eq. (1.2) is a nonlocal boundary condition which represents the birth law of the population. Finally, the initial distribution of the density is given by Eq. (1.3).

It is also usually assumed that

u (x_{\max}, t) = 0, t ⩾ 0

(1.5)

a condition which shows the fact that individuals in the population cannot reach their maximum individual size. In some cases, this property can be deduced from the behaviour of the growth function g [9].

The paper is structured as follows. In Section 2, we study the evolution of a population of mosquitofish, from the general equations (1.1)–(1.4). We introduce significative hypotheses in the model, the equations and the selection of the vital functions with their biological properties. The type of these functions and the initial condition is determined both from a biological and mathematical point of view. In the latter case, one of the most noticeable properties is regularity. Smooth functions are initially chosen for numerical reasons of convergence of the methods used. This selection is different from the one considered in Ref. [10] and has a relevant influence on the numerical results and also in a biological sense. Section 3 is a numerical study of the model. In general, the analytical resolution of models like (1.1)–(1.4) is not possible; thus it is necessary to use numerical integration to obtain information about the behaviour of the population under consideration. This study is based on two questions. First, we analyse the efficiency of four numerical methods for solving the problem. This allows us to choose the schemes which, in our opinion, are the most appropriate for the model considered. Also, with the methods selected, we introduce a numerical analysis of the model for long times. This previous numerical study leads us to consider a final biological discussion of the model. Here, we first describe the behaviour of the population in one year. Finally, the analysis for long times shows new topics such as the periodicity in the evolution of the population. This tendency to periodicity seems to be a property of the model since it occurs for different initial conditions. Some other phenomena, such as overcrowding or extinction, may take place if the mortality rate, which includes the only nonlinearity in the model, is changed.

2 Gambussia affinis model

The biological model which our study focuses on analyses the dynamics of a population of the so-called Gambussia affinis or mosquitofish. Sulsky [10] analyses some numerical questions in the integration of a model similar to (1.1)–(1.4) by building vital functions from experimental data, including seasonal hypotheses. Functions used by Sulsky are, in general, $C^{1}$ -piecewise, with an initial condition and a mortality rate represented by discontinuous functions.

Sulsky introduced some hypotheses in the model in order to simplify the study. First, length is used in the model as the variable which represents individual size. This quantity is measured in milimetres, within a range between $x_{\min} = 9 mm$ and $x_{\max} = 63 mm$ . She also chose the day as a unit of measure for temporal evolution. Furthermore, she assumed that the model uses the same size distribution for males and females. Therefore, only the population of females is considered and $u (x, t)$ denotes the density of females with length x at time t. Finally, she only considered the influence of the number of individuals of different size in the mortality. This dependence is represented by the function

I_{μ} (t) = \int_{9}^{63} γ_{μ} (x) u (x, t) d x, t ⩾ 0

(2.1)

With these hypotheses, Sulsky's model for the mosquitofish is specified by the equations

u_{t} + {(g (x, t) u)}_{x} = - μ (x, I_{μ} (t), t) u 9 < x < 63, t > 0

(2.2)

g (9, t) u (9, t) = \int_{9}^{63} α (x, t) u (x, t) d x, t > 0

(2.3)

u (x, 0) = u_{0} (x), 9 ⩽ x ⩽ 63

(2.4)

The use of (2.1)–(2.4) to describe the dynamics of a population of Gambussia affinis requires the specification of the vital functions, fertility $(α)$ , growth $(g)$ and mortality $(μ)$ rates. Our proposal, based on the study of Sulsky [10], assumes, as she does, dependence of the vital functions with respect to seasons, with a periodic character. The most important difference in our study is the use of smoother data, functions with more regularity. This selection is motivated by two reasons. The first is of a numerical character: the numerical methods we use are of order two, hence it is theoretically necessary some regularity of data in order to make them work. We have carried out experiments with less regular data and the results do not change significatively although efficiency is lost. On the other hand, we think that the use of smoother functions is more realistic from a biological point of view and establishes new proposals in the model. Therefore, we have constructed $C^{2}$ -functions which are periodic in time (with annual periodicity) and that satisfy the compatibility conditions so that the solution to the problem is of a $C^{2}$ -type.

Now, we describe in detail our selection of the vital functions. We first consider the fertility rate which we suppose in the form $α (x, t) = α (x) T_{α} (t)$ . Here, we have used field data by Krumholtz [11], which relates the length of gravid females to the number of young they have. In order to fit Krumholtz's data, $(x_{i}, y_{i})$ , we choose the function $α (x)$ as a smooth spline which minimizes the functional

p \sum_{i} w_{i} {(y_{i} - α (x_{i}))}^{2} + (1 - p) \int_{x_{\min}}^{x_{\max}} {(α^{″} (x))}^{2} d x

Here we have used a MATLAB^© function called csaps and

p = 0.5

. This fitting is better than the one proposed by Sulsky [10], which does not take into account the decrease of the function to fit the last data (see Fig. 1).

On the other hand, it is necessary to include dependence on time, because fertility is affected by seasons. We choose the function $T_{α} (t)$ as

T_{α} (t) = {\begin{matrix} {(\frac{t}{30})}^{3} (1 - \frac{t - 30}{10} + \frac{{(t - 30)}^{2}}{150}) \\ for 0 ⩽ t ⩽ 30 \\ 1, for 30 ⩽ t ⩽ 90 \\ - {(\frac{t - 120}{30})}^{3} (1 + \frac{t - 90}{10} + \frac{{(t - 90)}^{2}}{150}) \\ for 90 ⩽ t ⩽ 120 \\ 0, for 120 ⩽ t ⩽ 365 \end{matrix}

with

T_{α} (t + 365 n) = T_{α} (t)

,

n = 0, 1, \dots

.

This function is periodic with period $T_{p} = 365$ days and includes the fact that during a year there is only a season of 120 days when the population is fertile. The breeding season begins in May $(t = n T_{p})$ (see Krumholtz data [11]). The birth rate increases regularly as a quintic polynomial during the first month until it reaches a maximum value for the next two months. Then it decreases to zero during the following month while no births occur during the rest of the year. The same behaviour takes place every year.

The growth rate is also assumed in the form $g (x, t) = g (x) T_{g} (t)$ . The function g is

g (x) = \frac{63}{80.2} (1 - \frac{x}{63}), for 9 ⩽ x ⩽ 63

This belongs to the type of growth curves proposed by von Bertalanffy [12] and takes into account that a newborn with a length of 9 mm takes about six weeks to reach an adult length of 31 mm. Furthermore, the function g is zero for the maximum length allowed (63 mm) and it is of the type which makes that (1.5) holds.

The growth also varies with the season. The function $T_{g} (t)$ is also $T_{p}$ -periodic with

T_{g} (t) = {\begin{matrix} 0.2 + 0.8 {(\frac{t}{30})}^{3} (1 - \frac{t - 30}{10} + \frac{{(t - 30)}^{2}}{150}) \\ for 0 ⩽ t ⩽ 30 \\ 1, for 30 ⩽ t ⩽ 90 \\ 0.2 - 0.8 {(\frac{t - 120}{30})}^{3} (1 + \frac{t - 90}{10} + \frac{{(t - 90)}^{2}}{150}) \\ for 90 ⩽ t ⩽ 120 \\ 0.2, for 120 ⩽ t ⩽ 365 \end{matrix}

T_{g} (t + 365 n) = T_{g} (t)

,

n = 0, 1, \dots

. Therefore, growth has a seasonal behaviour, the maximum being reached during the fertile months.

Now, we describe the mortality rate, which includes the nonlinearity of the model. Bostford et al. [13] estimate the instantaneous per capita mortality rate for females with a length larger than 30 mm. They consider the quantity

N (t) = \int_{x (t)}^{\infty} u (x, t) d x

as the number of individuals that survive from a certain size

x (t)

. They estimate the mortality rate for females larger than 30 mm as the slope of the natural logarithm of values of

N (t)

. Experiments estimate a mortality rate of 0.023 days⁻¹, suggesting that the population of adult survivors decreases as a Malthusian law.

Our proposal assumes a mortality rate of the form $μ (x, t, z) = μ (x, z) T_{μ} (t)$ , $z = I_{μ} (t)$ . The function $μ (x, z)$ is a regularization of the discontinuous function considered by Sulsky [10]

μ (x, z) = {\begin{matrix} 0.1 \exp (\frac{- C}{z}), for 9 ⩽ x ⩽ 31 \\ 0.1 \exp (\frac{- C}{z}) - (0.023 - 0.1 \exp (\frac{- C}{z})) \\ \times {(x - 31)}^{3} (1 - 3 (x - 32) (65 - 2 x)) \\ for 31 < x ⩽ 32 \\ 0.023, for 32 < x ⩽ 63 \end{matrix}

(2.5)

with

C = 2000

, where

I (t) = \int_{9}^{31} 2 u (x, t) d x

, which describes the number of immature fish in the population.

Mortality also varies seasonally. The function $T_{μ} (t)$ is

T_{μ} (t) = {\begin{matrix} 2 - {(\frac{t}{30})}^{3} (1 - \frac{t - 30}{10} + \frac{{(t - 30)}^{2}}{150}) \\ 0 ⩽ t ⩽ 30 \\ 1, 30 ⩽ t ⩽ 90 \\ 2 + {(\frac{t - 120}{30})}^{3} (1 + \frac{t - 90}{10} + \frac{{(t - 90)}^{2}}{150}) \\ 90 ⩽ t ⩽ 120 \\ 2, 120 ⩽ t ⩽ 365 \end{matrix}

We have that

T_{μ} (t + 365 n) = T_{μ} (t)

,

n = 0, 1, \dots

, and

T_{μ}

expresses the fact that a high mortality occurs during the nonfertile seasons.

The mortality function depends on the number of individuals by means of the quantity $I_{μ}$ described in (2.1). The weight function $γ_{μ}$ represents the way total population affects the mortality of the youngest. In our case, this dependence is exclusively due to the immature fish, showing the behaviour of juveniles with respect to predation and competition for limited resources. Explicitly,

γ_{μ} (x) = {\begin{matrix} 2, 9 ⩽ x ⩽ 30 \\ - 2 {(x - 31)}^{3} (1 + 3 (x - 30) (2 x - 59)) \\ 30 < x < 31 \\ 0, 31 ⩽ x ⩽ 63 \end{matrix}

3 Numerical approximation

As it was observed in Section 1, the analytical resolution of the Gambussia problem is not possible. Thus our study has to be carried out by using numerical integration and includes information of two types: the first one is quantitative. We can obtain results about the density of the population in a fixed time in terms of the size under certain vital conditions. The results can also be of a qualitative type, showing, for instance, a certain behaviour for long times.

This section has two parts. The first one establishes a competition between four numerical schemes, in terms of efficiency and adaptation to the biological problem. This allows us to select the integrator which in our opinion is the most appropriate for the case considered. Secondly, with this method, we discuss the results in a biological sense, describing the relevant properties of the model proposed.

The selection of the numerical method to be used for the experiments is a question to take into account. It is natural to consider those methods that, giving a good approximation to the solution to problems like (2.1)–(2.4), their computational efficiency is appropriate. Therefore, before obtaining conclusions about the population under study through numerical analysis, it is necessary to make a selection between possible numerical integrators, in terms of the various types of efficiency. Then, the analysis of the model from numerical techniques requires two stages: first, we have to consider different numerical integrators and reject those which give results with no biological sense or with little computational efficiency. A second stage consists of obtaining conclusions about the dynamics of the model, studying its values and deficiencies by using the best of the numerical methods employed. In the literature, we can find some studies in this way for different species, such as the monogonont rotifera whose sexual phase dynamics was studied in [14,15].

Numerical integration of size-structured models is still at a primary stage, especially in the analysis of the convergence of numerical schemes. Roughly speaking, there are basically two classes of numerical schemes for the integration of size-structured problems: the standard finite difference methods (Upwind, Lax–Wendroff, Box) [16] and those schemes based on integration along characteristic curves, where the PDE is transformed into a system of ODE, infinite dimensional, using the characteristic curves. Higher order methods are appropriate for obtaining the solution faster than the others. However, the regularity of the data in real problems is not as high as needed. So, second order schemes maintain a good compromise between the required smoothness and the efficiency of the schemes. The first study was written by de Roos [17], who introduced the escalator boxcar train, a semidiscrete scheme for integrals (momenta) of the density function. Ito et al. [18] introduced a method for the linear problem based on the integration along the characteristic curves by using an invariant grid called natural grid. They proved convergence of second order for initial conditions of compact support. Schemes based on the natural grid have been proposed by different authors, for the linear problem [19,20] and the nonlinear problem [21]. In [22], the approximation to the nonautonomous problem with a growth function depending on size has been analysed by using a second order scheme based on the box method. The analysis of the fully nonlinear problem was initiated by Ackleh et al. [23]. They introduced an implicit finite difference scheme of first order. The first methods of second-order were proposed by Angulo et al. [24]. These methods are based on the integration along the characteristic curves, with a dynamic construction of the grid to improve the efficiency. On the other hand, the diversity of the numerical methods used for the problem also generates works about efficiency [10,22]. One conclusion which can be deduced from them is that methods that integrate along the characteristic curves are probably more efficient. However, as it is revealed in [22], the choice of the numerical method probably depends on the particular problem to be solved. Finally, we can find an extensive revision in [26], focused on the physiologically-structured models, which extends the study made in [26] for the age-structured ones.

3.1 Numerical results

In our case, we have considered two difference methods of second order, the Lax–Wendroff scheme [10,22] and box method [22] and two characteristics methods, Aggregation Grid Nodes method (AGN) and Selection Grid Nodes method (SGN) [24] to study the problem (2.1)–(2.4). All the schemes are briefly described in Appendix A. The Lax–Wendroff method is an explicit scheme with an easy implementation on a uniform grid. The box method is an implicit scheme on a uniform grid. Therefore, a nonlinear system of equations at each time step has to be solved by means of an iterative technique [22]. Both AGN and SGN methods are explicit and implemented in a nonuniform spatial grid. In the case of the AGN scheme, the number of grid nodes grows at each time step, making the computational cost in the long-time numerical integration more expensive. The SGN method avoids this cost by means of a selection strategy applied to each time step. They are described and analysed in [24].

We assume an initial distribution of the form

u_{0} (x) = {\begin{matrix} 0, 9 ⩽ x ⩽ 34 \\ 5 {(1 + \frac{x - 38}{4})}^{3}, 34 ⩽ x ⩽ 38 \\ 5 + 15 \frac{x - 38}{4} + 15 {(\frac{x - 38}{4})}^{2} + 30 {(\frac{x - 38}{4})}^{3} \frac{x - 46}{4} \\ 38 ⩽ x ⩽ 42 \\ 5 {(2 - \frac{x - 38}{4})}^{3}, 42 ⩽ x ⩽ 46 \\ 0, 46 ⩽ x ⩽ 63 \end{matrix}

(3.1)

Note that the initial population only contains adults. This fact is congruent with the hypothesis that the initial time corresponds to the beginning of the fertile period of the population. It is also supposed that, if an individual which was born on the last day of the fertile period of the last year has survived, then its length is greater than 31 mm.

3.1.1 Efficiency study

First, we show some experiments in order to determine the adaptation of the methods to the problem. In Fig. 2 the evolution in a year of the total population

P (t) = \int_{9}^{63} u (x, t) d x

and the number of newborns is plotted. The figure is obtained with all the numerical methods employed. Nevertheless, we shall see that the behaviour of the schemes is quite different when we focus on how they approximate the density function.

Fig. 2
Evolution of total population and newborns with time (days).

For a final time of one year $(T = 365)$ we present, for each method, a figure which shows the value of the numerical solution as a function of the length at time $t = 365$ .

We begin with the results obtained by the Lax–Wendroff scheme. We should point out that this method must satisfy the Courant–Friedrichs–Levy (CFL) condition. A more complete discussion of this method can be found in general studies on the numerical solution to partial differential equations like, for example, the work of Morton and Mayers [16]. To the best of our knowledge, the convergence of such a scheme for the solution to size-structured population models has still to be carried out. Numerical experiments in the literature [22] suggest that the method is of second order. The most efficient results correspond to the values of $r = \frac{k}{h}$ , closer to the one determined by the CFL condition (in this case $r ⩽ 1.485$ ). Taking $r = 1$ in Fig. 3, we show the numerical solution with $N = 5840$ , $J = 864$ ( $k = 0.0625$ and $h = 0.0625$ ) and a computational cost of 8.49 CPU-time (the CPU-time will be always measured in seconds). We note that some spurious oscillations appear, giving, in some cases, a negative population that has no biological meaning and that could damage the long time integration. This behaviour is not acceptable for greater values of k and h, although these oscillations tend to disappear as the discretization parameters are smaller. The results are even worse for smaller values of r, as we have confirmed with different experiments.

Fig. 3
Population density with length. The Lax–Wendroff method, t=365, J=864, N=5840.

Next, we deal with the behaviour of the box method. In the form given in [22], the scheme is truly implicit and at each time step the nonlinear equations must be solved by some fixed-point iteration procedure. Again, the text of Morton and Mayers can be useful in finding more properties of this scheme and we should point out that the authors emphasize the difficulty of carrying out convergence analysis even in the simple cases [16, p. 110]. However, the numerical method is second-order accurate as h tends to zero assuming that the time step is of the form $k = r h$ , with r an arbitrary and positive constant fixed throughout the analysis. The stability and convergence analysis in the cases of an age-structured and a nonautonomous size-structured population model has been derived [25,26]. We should also point out that, contrary to the rest of the methods employed in this section, the box scheme does not make use of (1.5) in the case of size-structured population models. Therefore, this method can be used with more general models where this condition does not appear. Our experience with the box method leads us to use parameter values with a relation $r = \frac{k}{h}$ greater than the one considered for the Lax–Wendroff scheme. In Fig. 4, we show the numerical solution obtained with $N = 730$ , $J = 432$ ( $k = 0.5$ , $h = 0.125$ ) corresponding to $r = 4$ and a computational cost of 1.63 CPU-time. Here, we observe some oscillations in the final profile. This could provide wrong information about the solution. When using $r ⩽ 4$ , we need more CPU-time in order to make these oscillations disappear, since an appropriate number of nodes in the spatial grid is necessary. On the other hand, in our experiments with $r ⩾ 6$ , we need more CPU-time to obtain good profiles; in this case some oscillations, similar to those of the Lax–Wendroff method with less computational cost, appear.

Fig. 4
Population density with length. The box method, t=365, J=432, N=730.

Next, we have considered the numerical integration of the problem with the characteristics schemes introduced in [24]. We should note that, in the case of AGN method, the number of nodes grows at each time level because a new node is introduced at each time step. The SGN-scheme avoids the increasing number of nodes by using a selection of them at each time level. These two methods are the first two second-order schemes which solve the fully nonlinear model. The convergence of the schemes was carried out in [24]. In Fig. 5, we show the results obtained with both methods. We have used $N = 1460$ , $J = 54$ ( $k = 0.25$ , $h = 1$ ) with $r = 0.25$ , which implies a CPU-time of 1.24 for the first method and $N = 730$ , $J = 216$ ( $k = 0.5$ , $h = 0.25$ ) with $r = 2$ and 0.20 CPU-time for the second one. The values of r are presented as the most efficient in the exhaustive experimentation with different parameters. The profiles obtained in both cases are smooth and no oscillations appear. Besides, the computational cost is much smaller than the one needed for the finite difference methods. Therefore, Aggregation and Selection Grid Nodes methods are the most appropriate for the biological example, while the Lax–Wendroff scheme is the worst, since it needs much more CPU-time to obtain a profile of the solution without oscillations.

Fig. 5
Population density with length. AGN (J=54, N=1460) and SGN methods (J=216, N=730), t=365.

From a more quantitative point of view, the results are very similar. In Fig. 6, we present, in a logarithmic scale, an efficiency plot where we show the error as a function of the computational time (in seconds). In the case of the computation of the error, since we do not know the exact solution to this problem, we obtain a numerical approximation to the total population at time $t = 365$ with small enough discretization parameters, for each numerical scheme. These values are very closed among them and we can consider them as the theoretical one. Thus, for another values of the parameters h and k the corresponding error is obtained by comparing the approximation to the total population with the considered as the theoretical one. Next, we build, for each numerical scheme, a tableau with the errors corresponding to different values of the parameters k and h. Now, by means of the procedure employed in [22], we discovered the most efficient value of r (for each method) and, finally, we compare all the schemes in Fig. 6. These values of r have been introduced previously.

Fig. 6
Error vs. CPU-time. Lax–Wendroff method (∗), box method (□), AGN scheme (▵), SGN method (+).

Fig. 6 shows that the most efficient methods are the characteristics schemes, and the best one corresponds to the SGN method. The behaviour of the AGN method is better as less computational time is employed. It begins to be less efficient when the number of nodes grows with the corresponding computational cost increment. In this last case, the difference methods are close to the efficiency of the AGN scheme. Among the difference methods, the Lax–Wendroff scheme is more efficient than the box method. (We recall that this method is implicit.) However, the distance between them is not as large as we previously could suppose, because the box method provides more freedom to choose the parameter values (k and h) than the Lax–Wendroff scheme, since this scheme must satisfy the CFL condition.

Finally, we should note that our study is made with $C^{2}$ -piecewise data functions and with compatible initial and boundary data. If we were unable to reach such smoothness or the initial and boundary data were not compatible, the behaviour of the difference methods would be even worse, except for the box with $r = 1$ . On the other hand, the characteristics methods dynamics is good enough and they are shown to be as the most efficient schemes [15,25].

3.1.2 Long-time behaviour

We now study the evolution of the population of mosquitofish. The numerical integration for long times requires to control the spurious oscillations, in order to avoid inappropriate results. Thus we have analysed the integration of the biological example for a period of ten years $(T = 3650)$ by using the methods which integrate along the characteristic curves. As we have seen previously, they achieve a smooth solution profile in a more efficient way.

In Fig. 7, we show the evolution over ten years of the numerical approximation to the total population and the number of newborns obtained with the AGN and SGN methods. Here we have used $N = 14600$ , $J = 54$ , with a CPU-time of 12.4 in the case of the AGN method and $N = 3650$ , $J = 108$ and a CPU-time of 0.5 for the SGN method. This figure shows that the behaviour of mosquitofish population, with an initial distribution given by (3.1), tends to be periodic. This periodic evolution represents the influence of the seasons on the model and we have obtained this behaviour with different initial conditions, even with the one considered by Sulsky [10]. In fact, we have repeated the same experiment with the vital functions selected in [10] and the periodic behaviour has also been obtained, but with greater computational effort. This is one of the reasons why the regularity of vital functions is important.

Fig. 7
Evolution of total population and newborns for ten years. AGN and SGN methods.

Therefore, the periodic evolution of the numerical approximation seems to be a property of the model. The SGN method is, in this case, more efficient than the AGN method, since the first one achieves the periodic behaviour of the total population and the number of newborns with much less computational time.

The study made in the previous subsection, which includes an appropriate selection of the parameters k and h (that is, the parameter r), is necessary. If we employ the Lax–Wendroff method, the box scheme or the AGN and SGN methods with smaller values of r than those employed here, the solution to the problem could not be obtained for long values of T.

3.2 Biological discussion of the model

In this subsection, we study the model from a biological point of view. First we explain the vital functions involved in the problem and then we describe some biological conclusions from the numerical results.

The model employed for the solution to the Gambussia problem is structured: this means that the dynamics of the population are determined by the vital functions defined in Section 2. The fertility rate is considered as $α (x, t) = α (x) T_{α} (t)$ . Function $α (x)$ is built using the field data (see Fig. 1), where we assume that the fertile females have a size longer than 31 mm and the fertility increases with the size. The maximum is reached at 54 mm and then it declines to zero at the maximum size. In order to take into account the seasonal behaviour, this is multiplied by a periodic function $T_{α} (t)$ (with a period equal to 365 days) so that the fertile period are the months between May and August and the maximum fertility is reached in June and July (see Fig. 8). The growth rate is taken as $g (x, t) = g (x) T_{g} (t)$ (see Fig. 9) where $g (x)$ is chosen as a von Bertalanffy growth curve, which is a classical approach in ecology, and the parameters are obtained by the field data. The function $T_{g} (t)$ weights the growth in order to show that this is greater at the period of the year with more environmental facilities (June/July) and it maintains a minimum during the months between September and April (both inclusive) while May and August are considered transition months between the extreme situations. Finally, the mortality function is taken in the form of $μ (x, t) = μ (x) T_{μ} (t)$ (see Fig. 10). Its behaviour is also seasonal and it is supposed that the maximum mortality period spans the months between September and April and the minimum one corresponds to the months of June and July. According to the studies of Bostford [13], the mortality rate is constant for adults and the choice for the immature fish depends on the total immature population. This last property implies a life of competition among these individuals and demonstrates the influence of the dynamics of the population in the life conditions of immature fish. The mathematical expression of this fact is arbitrary but we will show numerically the existence of a periodic state which is asymptotically stable. This dynamics is coherent with the biological setting. We will also discuss later the role of the constant C in the population dynamics.

Fig. 9
(a) Shape of g(x) function. (b) Shape of $T_{g} (t)$ function.

Fig. 10
(a) Shape of μ(x,z) function. (b) Shape of $T_{μ} (t)$ function.

Another property of the vital functions is that the individuals do not reach the maximum value in finite time. Moreover, it is easy to show that the flux gu at the maximum size is null.

The previous discussion between the four numerical methods was necessary to establish the most appropriate one to deduce conclusions about the model. The numerical study in Section 3.1 shows that for this problem the most efficient scheme is the SGN method. Therefore, with this integrator, here we will make some remarks, in biological terms, about the model considered in this paper.

For a final time of a year $(T = 365)$ we present some figures which show the value of the numerical solution as a function of the length, at times $t = 30$ , $t = 60$ , $t = 90$ , $t = 120$ (end of the fertile period), $t = 150$ , $t = 200$ and $t = 365$ . We first analyse the results from Fig. 11. The choice made of the vital functions implies that at the end of the first month (May, $t = 30$ ), the greatest part of the population is immature (4053 immatures and 18 adults) and also the population is clearly divided between juveniles and adults. In the second month $(t = 60)$ , there is a fertile population which has been born in this season. This fact is maintained in the next two months ( $t = 90$ and $t = 120$ ), a situation which generates the greatest number of births in the season and the fact that the biggest part of the mature population is born in this season. At the end of the fertile period $(t = 120)$ we can show that the mature population of the previous season remains and it starts to diminish until it completely disappears over longer times. Also, in this period ( $t = 120$ , $t = 150$ ), we can show that the profile of the population, which could be considered as an archetype, begins to appear. We can say that the population at the end of the year is only composed of adults. (See Fig. 5.) Also, after the study of the monthly evolution of the newborns, we can affirm that these adults correspond to the live newborns of the last fertile month. The newborns in such a period do not reproduce in the year and they are the only ones with this property. Finally, the third fertile month is the most prolific.

Fig. 11
Evolution of population density with length. SGN method, J=216, N=730.

As has been pointed out before, the influence of the seasons is remarkable in Fig. 7. This shows that the total population and newborns tend to be periodic. This behaviour has been obtained with different initial conditions. Therefore, the existence of a periodic solution which is asymptotically stable is shown numerically. This means that the population stabilizes in a state with annual periodicity under the hypotheses considered and by using field data. This annual periodicity is a consequence of the dependence of the vital functions with respect to the seasons. The asymptotic state is determined by the life conditions imposed on the population. They have been deduced from real data where possible. Therefore, under simple hypotheses, and if there is no change in the vital conditions of the population, we have observed that the population reaches a state which allows a specific population with a specific size distribution to survive.

An open question concerns the capacity of the model to reproduce other behaviour of the population. This can depend on the vital functions selected. As an example, we focus on the mortality function. The variation in the parameter of the mortality rate in (2.5) can be interpreted as a change in the life conditions. This is observed in the model by the fact that the asymptotic state allows a larger or smaller number of individuals to survive.

We have already commented that if the constant C in (2.5) is bigger, maximum mortality is reached for greater size, which makes life conditions easier, with less competition and a greater number of resources for the immature population. This should imply growth in the population. In Fig. 12, we show the evolution of the population over ten years by using the SGN method, with a mortality rate of the form (2.5), where $C = 200, 000$ . The behaviour of the population again tends to be periodic. However, if we compare this with Fig. 7, we observe that the population grows in two orders of magnitude. This percentage shows a linear dependence of the population, with respect to the value of the constant C, on the mortality function (2.5).

Fig. 12
Evolution of total population and the newborns over ten years. SGN method. C=200,000.

Similarly, when the constant C is smaller, life conditions are harder. In Fig. 13, we show the evolution of the population of Gambussia affinis over ten years by using the SGN method with a mortality rate of the form (2.5), where $C = 20$ . We observe that when life conditions become harder, the total population of the asymptotic equilibrium does not manage to reach individual. This can be interpreted as the extinction of the population. However, the way this fact is shown mathematically, that is from a periodic function (which therefore does not go to zero) is a deficiency of the model. We will try to improve this in a future work.

Fig. 13
Evolution of total population and the newborns over ten years. SGN method. C=20.

Acknowledgments

The authors are very grateful to the referees for their suggestions to improve the paper.

This work has been supported by project BFM2002-01250 and project VA063/04 (Consejería de Educación y Cultura de la Junta de Castilla y León).

Appendix A Numerical methods

A.1 The Lax–Wendroff method [10]

Let J be a positive integer. Let the points of the grid in the size variable be $x_{j} = x_{\min} + j h$ , $0 ⩽ j ⩽ J$ , where $h = (x_{\max} - x_{\min}) / J$ is the grid diameter. We denote by k the time step, and the discrete time levels as $t_{n} = n k$ , $0 ⩽ n ⩽ N$ , $N = [T / k]$ . The sub-index j makes reference to the grid point $x_{j}$ and the super-index n to the time level $t_{n}$ . Finally, we denote by $U_{j}^{n}$ the numerical approximation to $u (x_{j}, t_{n})$ , $0 ⩽ j ⩽ J$ , $0 ⩽ n ⩽ N$ . We also consider that an approximation to the initial condition (1.3), $U^{0}$ , is given.

The Lax–Wendroff method is a two-stage scheme defined for each $n = 0, 1, \dots, N - 1$ . First we calculate the intermediate values

U_{j - 1 / 2}^{n + 1 / 2} = U_{j - 1 / 2}^{n} - \frac{k}{2 h} (g_{j}^{n} U_{j}^{n} - g_{j - 1}^{n} U_{j - 1}^{n}) - \frac{k}{2} μ_{j - 1 / 2}^{n} (U) U_{j - 1 / 2}^{n},

(A.1)

where

μ_{j - 1 / 2}^{n} (U) = μ (x_{j - 1 / 2}, Q_{h} (γ_{μ} U^{n}), t_{n}) x_{j - 1 / 2} = \frac{1}{2} (x_{j - 1} + x_{j}) U_{j - 1 / 2}^{n} = \frac{1}{2} (U_{j - 1}^{n} + U_{j}^{n}), j = 1, 2, \dots, J

and

g_{j}^{n} = g (x_{j}, t_{n}), j = 0, 1, \dots, J

The function

Q_{h} (U^{n})

,

n = 0, 1, \dots, N - 1

, denotes the composite trapezoidal quadrature rule [22]. In the second stage, we obtain the values

U_{j}^{n + 1}

,

j = 1, 2, \dots, J - 1

, as

U_{j}^{n + 1} = U_{j}^{n} - \frac{k}{h} (g_{j + 1 / 2}^{n + 1 / 2} U_{j + 1 / 2}^{n + 1 / 2} - g_{j - 1 / 2}^{n + 1 / 2} U_{j - 1 / 2}^{n + 1 / 2}) - k μ_{j}^{n + 1 / 2} (U) U_{j}^{n + 1 / 2}

U_{J}^{n + 1} = 0

where

U_{j}^{n + 1 / 2} = \frac{1}{2} (U_{j + 1 / 2}^{n + 1 / 2} + U_{j - 1 / 2}^{n + 1 / 2}) μ_{j}^{n + 1 / 2} (U) = μ (x_{j}, Q_{h} (γ_{μ} U^{n + 1 / 2}), t_{n} + \frac{k}{2}) j = 1, 2, \dots, J - 1

and

g_{j + 1 / 2}^{n + 1 / 2} = g (x_{j + 1 / 2}, t_{n} + k / 2), j = 0, 1, \dots, J - 1

Now,

Q_{h} (γ_{μ} U^{n + 1 / 2})

,

0 ⩽ n ⩽ N - 1

, are approximations to (2.1) by using the composite mid-point quadrature rule and the values

U_{j + 1 / 2}^{n + 1 / 2}

,

j = 0, \dots, J - 1

, defined in (A.1),

Q_{h} (γ_{μ} U^{n + 1 / 2}) = \sum_{j = 0}^{J - 1} h \cdot γ_{μ} (x_{j + 1 / 2}) U_{j + 1 / 2}^{n + 1 / 2}

. Finally, the approximation

U_{0}^{n}

to

u (x_{\min}, t_{n})

, for

1 ⩽ n ⩽ N

is calculated with the condition

g_{0}^{n} U_{0}^{n} = Q_{h} (α^{n} U^{n})

where

α^{n} = (α (x_{0}, t_{n}), α (x_{1}, t_{n}), \dots, α (x_{J}, t_{n}))

and

g_{0}^{n} = g (x_{\min}, t_{n})

,

1 ⩽ n ⩽ N

.

A.2 The Box method [22]

The parameters J, N, h and k are defined as in Appendix A.1.

We introduce the half integer grid points $x_{j - 1 / 2} = \frac{1}{2} (x_{j - 1} + x_{j})$ , $1 ⩽ j ⩽ J$ ; the mean value operators $U_{j}^{n + 1 / 2} : = \frac{1}{2} (U_{j}^{n + 1} + U_{j}^{n})$ , $U_{j - 1 / 2}^{n} = \frac{1}{2} (U_{j - 1}^{n} + U_{j}^{n})$ , $U_{j - 1 / 2}^{n + 1 / 2} = \frac{1}{2} (U_{j - 1}^{n + 1 / 2} + U_{j}^{n + 1 / 2})$ and the difference operator $D U_{j}^{n} = U_{j}^{n + 1} - U_{j}^{n}$ .

The box method is defined by

\frac{D U_{j}^{n} + D U_{j - 1}^{n}}{2 k} + \frac{g_{j}^{n + 1 / 2} U_{j}^{n + 1 / 2} - g_{j - 1}^{n + 1 / 2} U_{j - 1}^{n + 1 / 2}}{h} = - μ_{j - 1 / 2}^{n + 1 / 2} (U) U_{j - 1 / 2}^{n + 1 / 2}

(A.2)

g_{0}^{n + 1} U_{0}^{n + 1} = Q_{h} (α^{n + 1} U^{n + 1})

(A.3)

1 ⩽ j ⩽ J

,

0 ⩽ n ⩽ N - 1

, where

Q_{h}

represents the trapezoidal quadrature rule, and, for

0 ⩽ n ⩽ N - 1

:

g_{j}^{n + 1 / 2} = g (x_{j}, t_{n} + \frac{k}{2}), 0 ⩽ j ⩽ J g_{0}^{n + 1} = g (x_{\min}, t_{n + 1}) μ_{j - 1 / 2}^{n + 1 / 2} (U) = μ (x_{j - 1 / 2}, Q_{h} (γ_{μ} U^{n + 1 / 2}), t_{n} + \frac{k}{2}) 1 ⩽ j ⩽ J α^{n + 1} = (α (x_{0}, t_{n + 1}), α (x_{1}, t_{n + 1}), \dots, α (x_{J}, t_{n + 1}))

with

Q_{h} (γ_{μ} U^{n}) = \sum_{j = 1}^{J} h γ_{μ} (x_{j - 1 / 2}) U_{j - 1 / 2}^{n}, 0 ⩽ n ⩽ N

A.3 Aggregation Grid Nodes method (AGN) [24]

The parameters J, N, h and k are defined as in Appendix A.1.

The initial grid nodes are $X_{j}^{0} = x_{\min} + j h$ , $0 ⩽ j ⩽ J$ . We suppose that the approximations to the theoretical solution in such nodes are known, $U_{j}^{0}$ , $0 ⩽ j ⩽ J$ . We also suppose that at the first time level $t_{1}$ , the grid nodes, $X^{1}$ , and the corresponding solution values, $U^{1}$ , are known. Furthermore, $X_{j}^{0}$ and $X_{j + 1}^{1}$ , $0 ⩽ j ⩽ J - 1$ , are (numerically) in the same characteristic curve. Angulo and López-Marcos obtained the initial conditions by means of the well-known second-order method [24].

The numerical approximations at the time level $t_{n + 2}$ , $0 ⩽ n ⩽ N - 2$ are obtained as follows. We suppose that the numerical approximations at the previous time levels, $t_{n}$ and $t_{n + 1}$ , are known, $X^{n}$ , $U^{n}$ and $X^{n + 1}$ , $U^{n + 1}$ . Where $X_{j}^{n}$ and $X_{j + 1}^{n + 1}$ , $0 ⩽ j ⩽ J + n - 1$ belong (numerically) to the same characteristic curve. We introduce the notation

α (X^{n}) = (α (X_{0}^{n}, t_{n}), α (X_{1}^{n}, t_{n}), \dots, α (X_{J + n + 1}^{n}, t_{n})) γ_{μ}^{n} = (γ_{μ} (X_{0}^{n}), γ_{μ} (X_{1}^{n}), \dots, γ_{μ} (X_{J + n + 1}^{n}))

Also, $Q^{n} (X^{n}, V^{n}) = \sum_{j = 0}^{J + n + 1} q_{j}^{n} (X^{n}) V_{j}^{n}$ where $q_{j}^{n} (X^{n})$ , $0 ⩽ j ⩽ J + n + 1$ , are the coefficients of the composite quadrature rules of at least second order. This notation will be used throughout the subsection. First, the grid values at the time level $t_{n + 2}$ are calculated by

X_{0}^{n + 2} = x_{\min}, X_{J + n + 2}^{n + 2} = x_{\max}

(A.4)

X_{1}^{n + 2} = k g (\frac{k}{2} g (x_{\min}, t_{n + 1}), t_{n + 1} + \frac{k}{2})

(A.5)

X_{j}^{n + 2} = X_{j - 2}^{n} + 2 k g (X_{j - 1}^{n + 1}, t_{n + 1}) 2 ⩽ j ⩽ J + n + 1

(A.6)

and the approximations to the theoretical solution in these nodes at such time level using

U_{1}^{n + 2} = U_{0}^{n + 1} \exp (- k μ (\frac{k}{2} g (x_{\min}, t_{n + 1}), \frac{3 Q^{n + 1} (X^{n + 1}, γ_{μ}^{n + 1} U^{n + 1}) - Q^{n} (X^{n}, γ_{μ}^{n} U^{n})}{2}, t_{n + 1} + \frac{k}{2}))

(A.7)

U_{j}^{n + 2} = U_{j - 2}^{n} \exp (- 2 k μ (X_{j - 1}^{n + 1}, Q^{n + 1} (X^{n + 1}, γ_{μ}^{n + 1} U^{n + 1}), t_{n + 1})) 2 ⩽ j ⩽ J + n + 1

(A.8)

U_{J + n + 2}^{n + 2} = 0

(A.9)

The equations at the time level

t_{n + 2}

are completed with the approximation

U_{0}^{n + 2}

to

u (x_{\min}, t_{n + 2})

by means of the discretization of the boundary condition (1.2)

U_{0}^{n + 2} = \frac{Q^{n + 2} (X^{n + 2}, α (X^{n + 2}) U^{n + 2})}{g (x_{\min}, t_{n + 2})}

(A.10)

A.4 Selection Grid Nodes method (SGN) [24]

The following scheme considers a modification in the grid of the previous one so that, by using a selection of the grid nodes, the number of nodes does not increase at each time level. Thus, we try to reduce the computational cost without loss of accuracy.

The grid nodes and the numerical approximations at time $t_{2}$ , $X^{2}$ , $U^{2}$ , are defined by means of (A.4)–(A.10) for $n = 0$ . Next, we calculate $Q^{2} (X^{2}, γ_{μ}^{2} U^{2})$ .

At consecutive time levels, there are different numbers of nodes because a new node that fluxes through the boundary is introduced. So, at the time level $t_{0}$ , we have $(J + 1)$ grid nodes, at $t_{1}$ we have $(J + 2)$ and at $t_{2}$ we have $(J + 3)$ . Now, the first grid node $X_{l}^{2}$ that satisfies

| X_{l + 1}^{2} - X_{l - 1}^{2} | = \min_{1 ⩽ j ⩽ J + 1} | X_{j + 1}^{2} - X_{j - 1}^{2} |

is eliminated and, also

X_{l - 1}^{1}

, the grid node in the same characteristic curve at

t_{1}

, is taken out. The number of nodes at the levels involved in the implementation of our two-step scheme are kept fixed:

(J + 3)

nodes for the time level reached in the integration and

(J + 2)

and

(J + 1)

for the previous ones. However, the approximations to the nonlocal terms at such time levels are not recomputed.

Now, we suppose that the numerical approximations at time levels $t_{n}$ and $t_{n + 1}$ are known, and they are denoted by ${X_{0}^{n}, X_{1}^{n}, \dots, X_{J - 1}^{n}, X_{J}^{n} = x_{\max}}$ , ${U_{0}^{n}, U_{1}^{n}, \dots, U_{J - 1}^{n}, U_{J}^{n} = 0}$ , $Q^{n} (X^{n}, γ_{μ}^{n} U^{n})$ and ${X_{0}^{n + 1} = x_{\min}, X_{1}^{n + 1}, \dots, X_{J}^{n + 1}, X_{J + 1}^{n + 1} = x_{\max}}$ , ${U_{0}^{n + 1}, U_{1}^{n + 1}, \dots, U_{J}^{n + 1}, U_{J + 1}^{n + 1} = 0}$ , $Q^{n + 1} (X^{n + 1}, γ_{μ}^{n + 1} U^{n + 1})$ (note that $X_{j}^{n}$ and $X_{j + 1}^{n + 1}$ , $0 ⩽ j ⩽ J - 1$ , are, numerically, in the same characteristic curve). In addition, the grid considered at $t_{n}$ has lost two nodes with respect to the moment when $X^{n}$ was actually calculated, while the grid used at $t_{n + 1}$ has only one node less than $X^{n + 1}$ . The numerical grid nodes at the new time level $t_{n + 2}$ , are computed by means of (A.5) and (A.6), $2 ⩽ j ⩽ J + 1$ , and the approximations to the theoretical solution in these nodes are obtained using (A.7) and (A.8), $2 ⩽ j ⩽ J + 1$ . The equations at time level $t_{n + 2}$ are completed with the approximation $U_{0}^{n + 2}$ to $u (x_{\min}, t_{n + 2})$ using (A.10).

Now, we calculate $Q^{n + 2} (X^{n + 2}, γ_{μ}^{n + 2} U^{n + 2})$ . Note that, for the time levels $t_{n}$ , $n ⩾ 2$ , the quadrature rules always use the same number of nodes $(J + 3)$ . Finally, we eliminate the first grid node $X_{l}^{n + 2}$ that satisfies

| X_{l + 1}^{n + 2} - X_{l - 1}^{n + 2} | = \min_{1 ⩽ j ⩽ J + 1} | X_{j + 1}^{n + 2} - X_{j - 1}^{n + 2} |

and we take out

X_{l - 1}^{n + 1}

, the grid node in the same characteristic curve at the previous time level.

Bibliographie

[1] The Dynamics of Physiologically Structured Populations (J.A.J. Metz; E.O. Diekmann, eds.), Lecture Notes in Biomathematics, vol. 68, Springer, Heidelberg, 1986

[2] G.F. Webb Theory of Nonlinear Age-Dependent Population Dynamics, Marcel Dekker, New York, 1985

[3] M. Iannelli Mathematical Theory of Age-Structured Population Dynamics, Applied Mathematics Monographs, CNR, 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] J.W. Sinko; W. Streifer A new model for age-size structure of a population, Ecology, Volume 48 (1967) no. 6, pp. 910-918

[6] L.F. Murphy A nonlinear growth mechanism in size-structured population dynamics, J. Theor. Biol., Volume 104 (1983), pp. 493-506

[7] H. Von Förster Some remarks on changing populations (F. Stohlman, ed.), The Kinetics of Cellular Proliferation, Grune and Stratton, New York, 1959, pp. 382-407

[8] M.E. Gurtin; R.C. MacCamy Nonlinear age-dependent population dynamics, Arch. Ration. Mech. Anal., Volume 54 (1974), pp. 281-300

[9] A. Calsina; J. Saldaña A model of physiologically structured population dynamics with a nonlinear individual growth rate, J. Math. Biol., Volume 33 (1995) no. 4, pp. 335-364

[10] D. Sulsky Numerical solution of structured population models II. Mass structure, J. Math. Biol., Volume 32 (1994), pp. 491-514

[11] L.A. Krumholtz Reproduction in the western mosquitofish, Gambussia affinis (Baird & Girard) and its use in mosquito control, Ecol. Monogr., Volume 18 (1948), pp. 1-43

[12] A.J. Fabens Properties and fitting of the von Bertalanffy growth curve, Growth, Volume 29 (1965), pp. 265-289

[13] L.W. Bostford; B. Vondracek; T.C. Wainwright; A.L. Linden; R.G. Kope; D.E. Reed; J.J. Cech Population development of the mosquitofish, Gambussia affinis, in rice fields, Environ. Biol. Fish., Volume 20 (1987) no. 2, pp. 143-154

[14] O. Angulo; J.C. López-Marcos; M.A. López-Marcos A numerical simulation for the dynamics of the sexual phase of monogonont rotifera, C. R. Biologies, Volume 327 (2004), pp. 293-303

[15] O. Angulo, J.C. López-Marcos, M.A. López-Marcos, A numerical integrator for a model with discontinuous sink term: The dynamics of the sexual phase of monogonont rotifera, Nonlinear Anal. Real World Appl. (2005), in press

[16] K.W. Morton; D.F. Mayers Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, 1994

[17] A.M. de Roos Numerical methods for structured population models: The escalator boxcar train, Numer. Methods Partial Differ. Equations, Volume 4 (1988), pp. 173-195

[18] K. Ito; F. Kappel; G. Peichl A fully discretized approximation scheme for size-structured population models, SIAM J. Numer. Anal., Volume 28 (1991), pp. 923-954

[19] W. Huyer A size-structured population-model with dispersion, J. Math. Anal. Appl., Volume 181 (1994) no. 3, pp. 716-754

[20] O. Angulo; J.C. López-Marcos Numerical schemes for size-structured population equations, Math. Biosci., Volume 157 (1999), pp. 169-188

[21] O. Angulo; J.C. López-Marcos Numerical integration of nonlinear size-structured population equations, Ecol. Model., Volume 133 (2000), pp. 3-14

[22] O. Angulo; J.C. López-Marcos Numerical integration of autonomous and nonautonomous nonlinear size-structured population models, Math. Biosci., Volume 177–178 (2002), pp. 39-71

[23] A.S. Ackleh; K. Ito An implicit finite difference scheme for the nonlinear size-structured population model, Numer. Funct. Anal. Optim., Volume 8 (1997) no. 9–10, pp. 865-884

[24] O. Angulo; J.C. López-Marcos Numerical integration of fully nonlinear size-structured models, Appl. Numer. Math., Volume 50 (2004) no. 3–4, pp. 291-327

[25] L.M. Abia, O. Angulo, J.C. López-Marcos, Age-structured population dynamics models and their numerical solutions. Ecol. Model. (2005), in press

[26] L.M. Abia; O. Angulo; J.C. López-Marcos Size-structured population dynamics models and their numerical solutions, Discrete Contin. Dyn. Syst. Ser. B, Volume 4 (2004) no. 4, pp. 1203-1222

Commentaires - Politique