Numerical integration of fully nonlinear size-structured population models

Applied Numerical Mathematics 50 (2004) 291–327 www.elsevier.com/locate/apnum

Numerical integration of fully nonlinear size-structured population models ✩ O. Angulo a,∗ , J.C. López-Marcos b a Departamento de Matemática Aplicada, Universidad de Valladolid, C/ Fco. Mendizabal, 1, 47014 Valladolid, Spain b Departamento de Matemática Aplicada, Universidad de Valladolid, Valladolid, Spain

Available online 6 March 2004

Abstract Characteristics’ schemes of second order are introduced for the numerical solution of size-structured population models. The models are fully nonlinear: all the vital functions (growth, mortality and fertility rates) depend on the total population. The schemes are completely analysed: consistency, stability, existence and convergence are established. Some numerical experiments are also reported in order to show numerically the results proved in our analysis and to study the efficiency of a selection nodes technique.  2004 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Fully nonlinear size-structured population models; Characteristics’ methods; Selection grid points; Convergence

1. Introduction In this paper we study numerical schemes based on integration along characteristic curves for nonlinear models describing the evolution of a size-structured population; i.e., a population where individuals are distinguished by an internal variable which is formally named “size”, but in reality it can have the physical meaning of length, weight, mass, maturity, etc. In particular we consider models that consist of a nonlinear partial differential equation with nonlocal terms (the population balance law)

 
 (1.1) ut + g x, Ig (t), t u x = −µ x, Iµ (t), t u, 0 < x < 1, t > 0, ✩

Supported in part by projects VA002/01,VA063/04 (Consejéria de Educatión y Culfura de la Junta de Castilla y León) and project DGESIC-DGES BFM2002-01250. * Corresponding author. E-mail addresses: [email protected] (O. Angulo), [email protected] (J.C. López-Marcos). 0168-9274/$30.00  2004 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2004.01.007

292

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

a nonlinear and nonlocal boundary condition that represents the birth law 
g 0, Ig (t), t u(0, t) =

1


 α x, Iα (t), t u(x, t) dx,

t > 0,

(1.2)

0

and an initial condition u(x, 0) = u0 (x),

0
x
1,

(1.3)

where Iµ (t), Iα (t) and Ig (t) are defined by 1 Iµ (t) =

γµ (x)u(x, t) dx,

t  0,

(1.4)

γα (x)u(x, t) dx,

t  0,

(1.5)

γg (x)u(x, t) dx,

t  0.

(1.6)

0

1 Iα (t) = 0

1 Ig (t) = 0

The independent variables x and t represent size and time, respectively. The maximum size individuals may obtain is normalized to one. The dependent variable u(x, t) is the size-specific density of individuals with size x at time t. We assume that the size of any individual varies according to the next ordinary differential equation
 dx = g x, Ig (t), t . (1.7) dt The nonnegative functions g, α and µ represent the growth, the fertility and the mortality rates, respectively. These are usually called the vital functions and define the life history of an individual. Note that all the vital functions (g, µ and α) depend on the size x (the structuring internal variable) and on the time t. The explicit time dependence can reflect the influence of some environmental changes on the vital functions or a seasonal behaviour of the population. These functions also depend on the total amount of individuals in the population by means of the weighted total population sizes Ig (t), Iα (t) and Iµ (t), which represent the way of weighting the size distribution density in order to model the different influence of the individuals of different sizes on the life conditions. More details on physiologically structured models can be found in [21,27,15,12]. To the best of the authors’ knowledge, the mathematical properties of solutions of the problem (1.1)– (1.6) have not been established in its more general form yet, but there are some results for particular situations. For the linear (γg ≡ γα ≡ γµ ≡ 0) autonomous problem, Ito et al. [16] have proved existence, uniqueness and regularity of the solutions. For a class of growth rate functions g(x), the mathematical properties of the solutions for nonlinear problems can be obtained, by means of the Murphy’s trick [23, 21], from the properties of the solutions of nonlinear age-structured models [27,21,15,12]. Also, when the growth rate function has the form g(x, t) (γg ≡ 0), the existence and uniqueness of solutions for nonlinear problems with γα ≡ γµ have been proved in [18] and their continuous dependence in [17].

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

293

For the fully nonlinear autonomous problem, the existence and uniqueness have been established in [11] when γg ≡ γα ≡ γµ . In addition, the analysis carried out by these authors could be used to our general problem [25]. We note that the mathematical analysis of the solutions of the problem (1.1)–(1.6) in its general setting is beyond the scope of this paper and we will concentrate on the numerical aspects of the problem. Thus, we assume existence and uniqueness of sufficiently smooth solutions under suitable hypotheses. In the following, we carry out a brief review of the numerical methods proposed for the solution of models such as (1.1)–(1.6). De Roos [24] introduced a semidiscrete method, the “escalator boxcar train”, in terms of momenta of the original density function. This method did not provide a direct approximation to the density function and its convergence has not been considered yet. For the autonomous linear problem (the growth rate depends only on size), Ito et al. [16] proposed a method based on the integration along the characteristic curves. They introduced the idea of a time-invariant spatial grid, called “natural grid” and demonstrated second order accuracy for initial conditions with compact support contained in [0, 1). Angulo and López-Marcos [8] considered a wide class of schemes of arbitrary order based on Runge–Kutta methods and integration along the characteristic curves. They proved convergence with the same order as the Runge–Kutta method employed. Also, they made a deep analysis of the natural grid pointing out the problem of the accumulation of grid nodes close to the maximum individual size. For the autonomous nonlinear problem, Huyer [14] introduced a numerical method based on the integration along the characteristics and on the use of the natural grid. Huyer proposed her scheme for a size-structured model in which classes of individuals with different growth rates were considered and demonstrated first order accuracy for initial conditions with compact support contained in [0, 1). Angulo and LópezMarcos [9] proposed a characteristic curves scheme based on the natural grid and on the use of an analytical representation of the theoretical solution that was introduced first by Abia for an age-structured model [2,3]. The analysis that they developed showed second order of convergence. The first attempt to obtain the numerical solution of the fully nonlinear problem (all the vital functions depend on the total population), was made by Ackleh and Ito [7] that proved convergence for a numerical method based on an upwind scheme. In [4], Ackleh et al. analysed the same upwind scheme for the solution of a coupled system of nonlinear size-structured populations. On the other hand, Ackleh and Deng [5,6] introduced an analytical technique based on an iterative process which gives monotone approximations by solving suitable linear problems. They proved the convergence of the approximations in different settings, autonomous and nonautonomous (the growth rate depend also on time), but using in all cases linear boundary conditions. Later, Angulo and López-Marcos [10] developed the first convergence analysis of a second order finite difference method, the box scheme, for the nonautonomous and nonlinear setting. Finally, the first explicit method of third order was proposed by Kostova [19] for a nonlinear autonomous model under a condition for the growth function that avoids the accumulation of grid nodes but that makes the model to be less realistic. A more detailed revision can be found in the work by Abia et al. [1]. In the present paper, we introduce second order numerical methods for the fully nonlinear model and we also carry out their convergence analysis. They are the highest order methods analysed for the fully nonlinear setting. Also, we prove the convergence of a particular adaptive moving mesh method. Throughout the paper we assume the following regularity conditions on the data functions and the solution of the problem (1.1)–(1.6): (H1) u ∈ C 2 ([0, 1] × [0, T ]), u(x, t)  0, x ∈ [0, 1), u(1, t) = 0, t  0. (H2) γµ , γα , γg ∈ C 2 ([0, 1]).

294

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

(H3) µ ∈ C 2 ([0, 1] × Dµ × [0, T ]), is nonnegative and Dµ is a compact neighbourhood of  1 γµ (x)u(x, t) dx, 0
t
T . 0

(H4) α ∈ C 2 ([0, 1] × Dα × [0, T ]), is nonnegative and Dα is a compact neighbourhood of  1 γα (x)u(x, t) dx, 0
t
T . 0

(H5) g ∈ C 3 ([0, 1] × Dg × [0, T ]), where Dg is a compact neighbourhood of  1 γg (x)u(x, t) dx, 0
t
T , 0

and g(x, z, t) > 0, x ∈ [0, 1) and g(1, z, t) = 0, t  0, z ∈ R. In addition, the characteristic curves x(t; t ∗ , x ∗ ) defined in (2.2) are continuous and differentiable with respect to the initial values (t ∗ , x ∗ ) ∈ [0, T ] × [0, 1]. Also, we assume that individuals cannot reach the maximum individual size. The above hypotheses can be based on three possible reasons. Firstly, biological assumptions such as the nonnegativity of the vital functions or, in (H5), to reflect that the size of any individual in the studied population [13,21,26] is strictly increasing during its lifetime and never reaches the maximum value. Secondly, the mathematical requirements to obtain existence and uniqueness of solutions for the problem (1.1)–(1.6), [16,11,7,18]. Lastly, the regularity properties needed in the numerical analysis to derive optimal rates of convergence [16,8–10,19]. The paper is organized as follows. Section 2 is devoted to introducing the numerical methods in detail. From Sections 3 to 6, we present the convergence analysis of the schemes introduced in Section 2. The third Section is devoted to introducing the discretization framework employed and some auxiliary results. In Section 4 the consistency analysis of the schemes is made and the fifth Section is devoted to study the stability properties of such methods. Next, in Section 6 we prove the convergence theorems. Finally, the last Section is devoted to presenting the numerical results.

2. Numerical integration along characteristics In this Section we describe methods that integrate numerically equations (1.1)–(1.6) along the characteristic curves. We begin rewriting the partial integro-differential equation in a more suitable form for its numerical treatment. So we define µ∗ (x, z1 , z2 , t) = µ(x, z1 , t) + gx (x, z2 , t), thus (1.1) has the form
 
ut (x, t) + g x, Ig (t), t ux (x, t) = −µ∗ x, Iµ (t), Ig (t), t u(x, t),

0 < x < 1, t > 0.

(2.1)

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

295

Next, we denote by x(t; t ∗ , x ∗ ) the characteristic curve of Eq. (2.1) that takes the value x ∗ at time t ∗ . This is the solution of the next initial value problem 
 

 x  t; t ∗ , x ∗ = g x t; t ∗ , x ∗ , Ig (t), t , t  t ∗ , 
(2.2) x t ∗; t ∗, x∗ = x∗. Then, we define the function 

 
w t; t ∗ , x ∗ = u x t; t ∗ , x ∗ , t ,

t  t ∗,

(2.3)

that satisfies the next initial value problem 



  d w t; t ∗ , x ∗ = −µ∗ x t; t ∗ , x ∗ , Iµ (t), Ig (t), t w t; t ∗ , x ∗ , dt 

w t ∗; t ∗, x∗ = u x∗, t ∗ ,

t  t ∗,

and therefore, it can be represented with the next formula   t 
∗ ∗


∗ ∗ w t; t , x = u x , t exp − µ∗ x τ ; t ∗ , x ∗ , Iµ (τ ), Ig (τ ), τ dτ .

(2.4)

(2.5)

t∗

(Recall that Ig (t) and Iµ (t) are defined by (1.4) and (1.6), respectively.) Now we introduce two numerical schemes that simultaneously integrate the coupled problems (2.2) and (2.5) together with the boundary condition (1.2). The first one increases the number of nodes at each time level because a new characteristic curve is considered in the integration. The second one avoids the increasing number of nodes by means of a selection of them at each time level. Both schemes are two-step methods. 2.1. Aggregation grid nodes method Let J and N be positive integers, we define the spatial and time discretization parameters as h = 1/J and k = T /N , respectively. The discrete time levels are tn = nk, 0
n
N , and the initial grid nodes are Xj0 = j h, 0
j
J . We suppose that the approximations to the theoretical solution in such nodes are known, Uj0 , 0
j
J . Thus, we denote     U0 = U00 , U10, . . . , UJ0−1 , UJ0 = 0 . X0 = X00 = 0, X10 , . . . , XJ0 −1 , XJ0 = 1 , We also suppose that at the first time level t1 ,   X1 = X01 = 0, X11 , . . . , XJ1 , XJ1 +1 = 1 ,

  U1 = U01 , U11, . . . , UJ1 , UJ1+1 = 0 ,

the grid nodes and the corresponding solution values, respectively, are known. Furthermore, Xj0 and Xj1+1 , 0
j
J − 1, are (numerically) in the same characteristic curve. We explain in our numerical experiments (Section 7), how we calculate X1 and U1 . We will obtain the numerical approximations at the time level t2 as follows. First we define the grid nodes X2 = {X02 = 0, X12 , . . . , XJ2 +1 , XJ2 +2 = 1} by numerical integration of (2.2) with the formulae 
1 1 1   3Q1 (X1 , γ 1g U1 ) − Q0 (X0 , γ 0g U0 ) k
k 2 1 , t1 + , (2.6) X1 = kg g 0, Q X , γ g U , t1 , 2 2 2

  (2.7) Xj2 = Xj0−2 + 2kg Xj1−1 , Q1 X1 , γ 1g U1 , t1 , 2
j
J + 1.

296

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

Then we calculate the corresponding approximations to the theoretical solution U2 = {U02 , U12 , . . . , UJ2+1 , UJ2+2 = 0}, by means of the following discretization of (2.5)  
  3Q1 (X1 , γ 1µ U1 ) − Q0 (X0 , γ 0µ U) k
, U12 = U01 exp −kµ∗ g 0, Q1 X1 , γ 1g U1 , t1 , 2 2 3Q1 (X1 , γ 1g U1 ) − Q0 (X0 , γ 0g U0 ) k , t1 + , (2.8) 2 2



  (2.9) Uj2 = Uj0−2 exp −2kµ∗ Xj1−1 , Q1 X1 , γ 1µ U1 , Q1 X1 , γ 1g U1 , t1 , 2
j
J + 1. Finally we derive the approximation U02 to u(0, t2 ) from a discrete version of the boundary condition (1.2) U02 =

Q2 (X2 , α(X2 , U2 )U2 ) . g(0, Q2 (X2 , γ 2g U2 ), t2 )

(2.10)

In Eqs. (2.6)–(2.10), we have used the next notation 





α X2 , U2 = α0 X2 , U2 , α1 X2 , U2 , . . . , αJ +1 X2 , U2 , αJ +2 X2 , U2 , where




  αj X2 , U2 = α Xj2 , Q2 X2 , γ 2α U2 , t2 , 0
j
J + 2,




 γ lφ = γφ X0l , γφ X1l , . . . , γφ XJl +l−1 , γφ XJl +l , l = 0, 1, 2, φ = g, µ, α;

and J +l
l l


 qjl Xl Vjl , Q X ,V = l

l = 0, 1, 2,

j =0

where qjl (Xl ), 0
j
J + l, l = 0, 1, 2, are the coefficients of composite quadrature rules of second order (at least) with nodes in Xl . Also, γ nφ Un , where φ can denote g, µ or α, n = 0, 1, 2; and α(X2 , U2 )U2 , represent the componentwise product of the corresponding vectors. We would like to note that we have employed the explicit midpoint rule to integrate (2.2), the midpoint quadrature rule to discretise (2.5) and extrapolation of the nonlocal terms in the equations corresponding to the computation of X12 and U12 . Also the value of U02 is defined with an implicit equation (2.10) which will be solved by means of an iterative procedure. Next, we describe the general time step, tn+2 , 0
n
N − 2. Now, we suppose that the numerical approximations at the previous time levels, tn and tn+1 , are known     Un = U0n , U1n , . . . , UJn+n−1 , UJn+n = 0 , Xn = X0n = 0, X1n , . . . , XJn +n−1 , XJn +n = 1 , and

  n+1 Xn+1 = X0n+1 = 0, X1n+1 , . . . , XJn+1 +n , XJ +n+1 = 1 ,   n+1 Un+1 = U0n+1 , U1n+1 , . . . , UJn+1 +n , UJ +n+1 = 0 .

We recall that Xjn and Xjn+1 +1 , 0
j
J + n − 1 are (numerically) in the same characteristic curve. First, we calculate the grid values at the time level tn+2   n+2 Xn+2 = X0n+2 = 0, X1n+2 , . . . , XJn+2 +n+1 , XJ +n+2 = 1 ,

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

297

by means of the numerical integration of (2.2) 
  k
n+2 n+1 , tn+1 , X1 = kg g 0, Qn+1 Xn+1 , γ n+1 g U 2

k , tn+1 + , 2 2

n+1 n+1 n+1   n+1 X ,γg U , tn+1 , 2
j
J + n + 1, Xjn+2 = Xjn−2 + 2kg Xjn+1 −1 , Q n+1 ) − Qn (Xn , γ ng Un ) 3Qn+1 (Xn+1 , γ n+1 g U

(2.11) (2.12)

and the approximations to the theoretical solution in these nodes at such time level   n+2 Un+2 = U0n+2 , U1n+2 , . . . , UJn+2 +n+1 , UJ +n+2 = 0 , using the discretization of (2.5)  

  n+2 n+1 ∗ k n+1 g 0, Qn+1 Xn+1 , γ n+1 , tn+1 , U1 = U0 exp −kµ g U 2 n+1 ) − Qn (Xn , γ nµ Un ) 3Qn+1 (Xn+1 , γ n+1 µ U

2 n+1 ) − Qn (Xn , γ ng Un ) 3Qn+1 (Xn+1 , γ n+1 g U

2




n+1 n+1 n+1  n+1 X ,γµ U , Ujn+2 = Ujn−2 exp −2kµ∗ Xjn+1 −1 , Q
n+1 n+1 n+1   n+1 X ,γg U , tn+1 , Q

, k , tn+1 + 2

,

2
j
J + n + 1.

(2.13)

(2.14)

We complete the equations of the time level tn+2 with the approximation U0n+2 to u(0, tn+2 ) by means of the discretization of the boundary condition (1.2) U0n+2 =

Qn+2 (Xn+2 , α(Xn+2 , Un+2 )Un+2 ) . n+2 ), t g(0, Qn+2 (Xn+2 , γ n+2 n+2 ) g U

(2.15)

We have used the same notation as in the integration at the time level t2 . Note that at every time level, we can use quadrature rules with different order (but at least second order), or based on diverse number of nodes or with nodes that have no relationship with the employed in rules of other time levels. (However they have to satisfy a suitable property, see Theorem 2). In our numerical experiments (see Section 7), we have used at each time level, the composite trapezoidal quadrature rule based on the whole grid nodes. The name aggregation grid nodes method is due to the fact that a new grid node is introduced at each time step. This node can be viewed like the numerical counterpart of the flux population from the minimum individual size (newborns). 2.2. Selection grid nodes method The following scheme is a variation of the previous one such that, by means of a selection of the grid nodes, the number of nodes does not increase at each time level. Thus, we try to lower the computational cost without lost of accuracy.

298

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

Let J and N be positive integers, we define the spatial and time discretization parameters as h = 1/J and k = T /N , respectively. The discrete time levels are tn = nk, 0
n
N , and the initial grid nodes are Xj0 = j h, 0
j
J . In the same way that for the aggregation grid nodes method, we assume that     U0 = U00 , U10, . . . , UJ0−1 , UJ0 = 0 , X0 = X00 = 0, X10 , . . . , XJ0 −1 , XJ0 = 1 , and

  X1 = X01 = 0, X11 , . . . , XJ1 , XJ1 +1 = 1 ,

  U1 = U01 , U11, . . . , UJ1 , UJ1+1 = 0 ,

are known and that Xj0 and Xj1+1 , 0
j
J − 1, are (numerically) in the same characteristic curve. Also,     U2 = U02 , U12 , . . . , UJ2+1 , UJ2+2 = 0 , X2 = X02 = 0, X12 , . . . , XJ2 +1 , XJ2 +2 = 1 , are defined by means of (2.6)–(2.10). Next we calculate Q2 (X2 , γ 2g U2 ) and Q2 (X2 , γ 2µ U2 ). We observe that at consecutive time levels, we are working with a different number of nodes because we have introduced a new node that fluxes through the boundary. So, at the time level t0 we have (J + 1) grid nodes, at t1 we have (J + 2) and at t2 we have (J + 3). Now, we eliminate the first grid node Xl2 that satisfies     2 X − X 2  = min X 2 − X 2 , (2.16) l+1

l−1

1
j
J +1

j +1

j −1

1 , the grid node in the same characteristic curve at t1 . Therefore we keep fixed and we also take out Xl−1 the number of nodes at the levels involved in the implementation of our two-step scheme: (J + 3) nodes for the time level reached in the integration and (J + 2) and (J + 1) for the previous ones. However, we do not recompute the approximations to the nonlocal terms at such time levels. Now we describe the general time step, tn+2 , 0
n
N − 2. We suppose that the approximations at the previous time levels tn and tn+1 , are known   n n   n n U0 , U1 , . . . , UJn−1 , UJn = 0 , X0 , X1 , . . . , XJn −1 , XJn = 1 ,

Qn (Xn , γ ng Un ), Qn (Xn , γ nµ Un ) and   n+1 X0 = 0, X1n+1 , . . . , XJn+1 , XJn+1 +1 = 1 ,



 U0n+1 , U1n+1 , . . . , UJn+1 , UJn+1 +1 = 0 ,

n+1 n+1 ), Qn+1 (Xn+1 , γ n+1 ). We recall that Xjn and Xjn+1 Qn+1 (Xn+1 , γ n+1 g U µ U +1 , 0
j
J − 1, are (numerically) in the same characteristic curve. In addition, the nodes considered at tn have two nodes less with respect to the moment when Xn was actually calculated, while the nodes used at tn+1 have only one node less than Xn+1 . First, we compute the grid values at the time level tn+2 ,   n+2 Xn+2 = X0n+2 = 0, X1n+2 , . . . , XJn+2 +1 , XJ +2 = 1 ,

by means of the numerical integration of (2.2) 
  k
n+2 n+1 , tn+1 , X1 = kg g 0, Qn+1 Xn+1 , γ n+1 g U 2 n+1 ) − Qn (Xn , γ ng Un ) 3Qn+1 (Xn+1 , γ n+1 g U

k , tn+1 + , 2 2

n+1 n+1 n+1   n+1 X ,γg U , tn+1 , 2
j
J + 1, Xjn+2 = Xjn−2 + 2kg Xjn+1 −1 , Q

(2.17) (2.18)

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

299

and the approximations to the theoretical solution in these nodes at such time level   n+2 Un+2 = U0n+2 , U1n+2 , . . . , UJn+2 +1 , UJ +2 = 0 , using the discretization of (2.5)  

  n+2 n+1 ∗ k n+1 g 0, Qn+1 Xn+1 , γ n+1 , tn+1 , U1 = U0 exp −kµ g U 2 n+1 ) − Qn (Xn , γ nµ Un ) 3Qn+1 (Xn+1 , γ n+1 µ U 2 n+1 ) − Qn (Xn , γ ng Un ) 3Qn+1 (Xn+1 , γ n+1 g U

2 Ujn+2

=

Ujn−2 exp






n+1 n+1 n+1  n+1 −2kµ Xjn+1 X ,γµ U , −1 , Q
n+1 n+1 n+1   n+1 X ,γg U , tn+1 , Q

, k , tn+1 + 2

, (2.19)

∗

2
j
J + 1.

(2.20)

We complete the equations at the time level tn+2 with the approximation U0n+2 to u(0, tn+2 ) using a discretization of the boundary condition (1.2) U0n+2 =

Qn+2 (Xn+2 , α(Xn+2 , Un+2 )Un+2 ) . n+2 ), t g(0, Qn+2 (Xn+2 , γ n+2 n+2 ) g U

(2.21)

n+2 n+2 ) and Qn+2 (Xn+2 , γ n+2 ). Note that for the time levels Now we calculate Qn+2 (Xn+2 , γ n+2 g U µ U tn , n  2, the quadrature rules use the same number of nodes (J + 3). This fact will play an important role in the convergence analysis. The notation used is similar to that employed in Section 2.1. Next, we eliminate the first grid node Xln+2 that satisfies     n+2 n+2  n+2  X min Xjn+2 (2.22) l+1 − Xl−1 = +1 − Xj −1 , 1
j
J +1

n+1 , the grid node in the same characteristic curve at the previous time level. and we take out Xl−1 This method can be viewed as a particular case of the aggregation grid nodes method because we can consider the selection procedure as the choice of the nodes employed in the quadrature rules at each time level tn , for n  2, which form a subgrid of the whole grid points. At the time levels tn , 0
n
2, the subgrid is the whole grid of the aggregation grid nodes method. When n  3, the subgrid is formed by the nodes obtained by integration along characteristics from the nodes selected at the previous time level together with the zero (minimum size). Note that 1 (maximum size) is always a node belonging to the subgrid because x = 1 is a characteristic curve (recall (H5)) and it cannot be removed with the criterion (2.22). From this perspective, we can still compute, hypothetically, the grid nodes (and the solution values at such nodes) that are not employed in the quadrature rules. The convergence analysis will be carried out in this statement. However, the nodes (and the solution values at them) that do not belong to the subgrid at the time level tn , 3
n
N , have no effect over the computation of the other nodes (and the solution values at them) at the subsequent time levels. Therefore, we only calculate the subgrid nodes (and the solution values at them). Hence, we save computational effort and, due to the accumulation of nodes close to the maximum individual size for the aggregation scheme, the selection strategy makes the accuracy not to diminish substantially.

300

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

3. Convergence analysis: Preliminaries In this Section we begin the analysis of numerical methods based on integration along characteristics that use quadrature rules with suitable properties to approximate the integral terms. Convergence will be obtained by means of consistency and nonlinear stability. First, we rewrite the numerical methods into the discretization framework developed by López-Marcos and Sanz-Serna [20]. We assume that the spatial discretization parameter, h, takes values in the set H = {h > 0: h = 1/J, J ∈ N}. Now, we suppose that the time step, k, satisfies k = rh, 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 space Ah =

N


 RJ +n−1 × RJ +n ,

n=0 J +n−1

is used for the approximations to the interior grid nodes and RJ +n for the approaches to where R the theoretical solution on them and on the left boundary node, at time level tn , 0
n
N . We also consider the space N
J −1 
J 
J +n−1  J J +1 N−1 ×R × R ×R × × RJ +n−1 , ×R R Bh = R n=2

where (RJ −1 × RJ ) × (RJ × RJ +1 ) is employed to compare with the initial approximations; RN−1 considers residuals that take place on the boundary node for every time step (except the first two);

theJ +n−1 (R × RJ +n−1 ), is used for the residuals which arise in the formulae that define the grid and N n=2 nodes and the solution values. We note that in the spaces, Ah and Bh , we consider neither the first and the last grid nodes nor the value of the solution in the last node because they are fixed values. Thus, both spaces have the same dimension. In order to measure the size of the errors, we define a∞ = max |aj |, 1
j
p

a = (a1 , a2 , . . . , ap ) ∈ Rp ,

J

+n−1    n V  = hVjn , 1

Vn ∈ RJ +n .

j =0

Now, we endow the spaces Ah and Bh with the following norms. If (y0 , V0 , . . . , yN , VN ) ∈ Ah , then       
0 0   y , V , . . . , yN , VN  = max max yn  , max Vn  . Ah ∞ ∞ 0
n
N

0
n
N

On the other hand, if (Y , Z , Y , Z , Z0 , Y , Z , . . . , YN , ZN ) ∈ Bh , then 
0 0 1 1   Y , Z , Y , Z , Z0 , Y2 , Z2 , . . . , YN , ZN  Bh 0

0

1

1

2

2

N N



 0  0  1  1    n               = Y ∞ + Z ∞ + Y ∞ + Z ∞ + Z0 ∞ + k Z ∞+ k Yn ∞ . n=2

For each h ∈ H , we define

n=2

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327


xh = x0 , x1 , x2 , . . . , xN ,

301


xn = x1n , . . . , xJn+n−1 ∈ RJ +n−1 ,

xj0 = j h, 1
j
J − 1,
 1
j
J + n − 1, 1
n
N; xjn = x tn ; tn−1 , xjn−1 −1 , ∗

xJn+n = ∗

(3.1) ∗

= 0, 1, n  0. Recall that x(t; t , x ) represents the theoretical solution of and we denote problem (2.2), t ∈ [0, T ], x ∈ [0, 1]. In addition, if u represents the theoretical solution of (1.1)–(1.6) we define

  un = un0 , un1 , . . . , unJ +n−1 ∈ RJ +n , uh = u0 , u1 , u2 , . . . , uN ,
 unj = u xjn , tn , 0
j
J + n, 0
n
N; x0n ∗

so u˜ h = (x0 , u0 , x1 , u1 , . . . , xN , uN ) ∈ Ah . (By hypothesis (H1), note that unJ +n = 0, n  0.) Proposition 1. Assume that hypotheses (H1), (H2) and (H5) hold, then the grid nodes defined in (3.1) satisfy that xjn+1 > xjn and there also exist two positive constants C1 and C2 such that C1 h
xjn+1 − xjn
C2 h,

0
j
J + n − 1, 0
n
N.

(3.2)

Proof. There are three different cases and we will use the mean-value theorem and the derivatives of the solutions of (2.2) with respect to the initial values. If n
j
J + n − 2, 1
n
N , xjn and xjn+1 are points belonging to characteristic curves with initial grid nodes as initial conditions at t = 0, then


 ∂ x(tn ; 0, ξ ) xjn+1 − xjn = x tn ; 0, xj0−n+1 − x tn ; 0, xj0−n = xj0−n+1 − xj0−n ∂x ∗   tn 
gx x(s; 0, ξ ), Ig (s), s ds , = h exp

(3.3)

0

ξ ∈ (xj0−n , xj0−n+1 ). < n
N , xjn and xjn+1

where If j as initial conditions, then

are points which belong to characteristic curves that have boundary nodes

∂ xjn+1 − xjn = x(tn ; tn−j −1 , 0) − x(tn ; tn−j , 0) = (tn−j −1 − tn−j ) ∗ x(tn ; ξ, 0) ∂t  tn   

gx x(s; ξ, 0), Ig (s), s ds , = kg 0, Ig (ξ ), ξ exp

(3.4)

ξ

where ξ ∈ (tn−j −1 , tn−j ). Finally, if n
N ,

 ∂ x(tn ; 0, ξ ) 1 − xJn+n−1 = x(tn ; 0, 1) − x tn ; 0, xJ0 −1 = 1 − xJ0 −1 ∂x ∗   tn 
gx x(s; 0, ξ ), Ig (s), s ds , = h exp 0

(3.5)

302

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

where ξ ∈ (xJ0 −1 , 1). Then we get (3.2) by means of (3.3)–(3.5) and the smoothness properties of the involved functions. 2 On the other hand, we denote by n




n

Q X ,V

n



=

J +n


 qjn Xn Vjn ,

j =0

the quadrature rule used at each time level. We have to note that the number of nodes considered at each time level is J + n + 1, that does not coincide with the number of nodes of Xn , 0
n
N . This is due to that the quadrature rules use the fixed values of the nodes X0n = x0n = 0, XJn +n = xJn+n = 1 and the prescribed value VJn+n = unJ +n = 0, 0
n
N . This notation is also valid in case that we consider quadrature rules whose nodes are, at each time level, a subgrid of Xn , 0
n
N . From now on, C will denote a positive constant, independent of h, k (k = rh), j (0
j
J + n) and n (0
n
N ); C possibly has different values in different places. We assume that the quadrature rules satisfy the next properties: (P1) |Iµ (tn ) − Qn (xn , γ nµ un )|
Ch2 , when h → 0, 0
n
N . (P2) |Iα (tn ) − Qn (xn , γ nα un )|
Ch2 , when h → 0, 0
n
N . (P3) |Ig (tn ) − Qn (xn , γ ng un )|
Ch2 , when h → 0, 0
n
N . 1  +n n n qj (x )α(xjn , Iα (tn ), tn )unj |
Ch2 , when h → 0, 0
n
N . (P4) | 0 α(x, Iα (tn ), tn )u(x, tn ) dx − Jj =0 n n (P5) |qj (x )|
qh, 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) Let R and p be positive constants with 1 < p < 2. The quadrature weights qjn are Lipschitz continuous functions on B∞ (xn , Rhp ), 0
j
J + n, 0
n
N . (P7) Let R and p be positive constants with 1 < p < 2. If yn , zn ∈ B∞ (xn , Rhp ) and Vn ∈ B∞ (un , Rhp ), then   J +n  

  

   qin yn − qin zn γφ zin Vin 
C yn − zn ∞ ,    i=0

when h → 0, 0
n
N , where φ = g, µ, α. (P8) Let R and p be positive constants with 1 < p < 2. If yn , zn ∈ B∞ (xn , Rhp ) and Vn ∈ B∞ (un , Rhp ), then   J +n  

  

   qin yn − qin zn αi zn , Vn Vin 
C yn − zn ∞ ,    i=0

when h → 0, 0
n
N . The properties (P1)–(P8) will be sufficient to carry out our convergence analysis. Therefore, our theorems hold for any quadrature rules that satisfy (P1)–(P8). The following result establishes that the composite trapezoidal rule used in our experiments satisfies these properties.

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

303

Theorem 2. Assume that the hypotheses (H1)–(H5) hold. If the quadrature rules are the composite trapezoidal quadrature on subgrids {xjnn }M(n) l=0 , 0
n
N with the property l

(SR) There exists a positive constant C such that, for h sufficiently small, xjnn − xjnn
Ch, 0
l
l+1

M(n) − 1, xjnn = 0, xjnn

M(n)

0

l

= 1, with {xjnn }M(n)−1 contained in xn , 0
n
N . l=1 l

Then, properties (P1)–(P8) hold. Proof. We note that the properties (P1)–(P6) can be easily derived under our assumptions by means of Proposition 1 and the properties of the composite trapezoidal quadrature rule. Recall that with our / {xjnn }M(n) notation, qjn (xn ) = 0, if xjn ∈ l=0 , 0
n
N . l On the other hand, the proofs of (P7) and (P8) are very similar. Therefore, we will establish (P7) and we will omit the proof of (P8). We define ∆n = yn − zn , 0
n
N . If we take into account that y0n = z0n = 0, yJn+n = zJn +n = 1 and n Vj n = VJn+n = 0, 0
n
N , then M(n)

J +n


n
n

 qj y − qjn zn γφ zjn Vjn j =0

=

 y nn − y nn j1

j0

2 +

=−

 y nn

jM(n)

n M(n)−1 zjnn − zjnn
 zjnn − zjnn
 n − yj n

 yjnl+1 l−1 l−1 1 0 n n − γφ zj n Vj n + − l+1 γφ zjnn Vjnn 0 0 l l 2 2 2 l=1

− yjnn

M(n)−1

2

−

zjnn

M(n)

− zjnn

M(n)−1

2


 γφ zjnn Vjnn M(n)

M(n)

M(n)−1
  1 n

n  n ∆j n γφ zj n Vj n − γφ zjnn Vjnn , l l+1 l+1 l−1 l−1 2 l=1

(3.6)

where φ = g, µ, α. Now, the hypotheses (H1)–(H5) and the property (SR), allow us to get, after some standard operational machinery, that
 2

n  n γφ zj n Vj n − γφ zjnn Vjnn l+1 l+1 l−1 l−1

n 


n 2
n 2 
2
n 2
3 γφ zj n − γφ xj n + 3 γφ zjnn − γφ xjnn Vj n Vj n l+1 l+1 l+1 l−1 l−1 l−1

n  n
n  n 2 + 3 γφ xj n Vj n − γφ xj n Vj n l+1 l+1 l−1 l−1

n 


n 2
n 2 
2
n 2
3 γφ zj n − γφ xj n + 3 γφ zjnn − γφ xjnn Vj n Vj n l+1 l+1 l+1 l−1 l−1 l−1

n 2
n


 2 2 2 + 9 γφ xj n + 9 γφ xjnn Vj n − unjn Vjnn − unjn l+1 l+1 l+1 l−1 l−1 l−1

n  n
2p
n  n 2  + 9 γφ xj n uj n − γφ xj n uj n
C h + h2 , 1
l
M(n) − 1. l+1

l+1

l−1

l−1

(3.7)

304

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

Thus, by means of the Cauchy–Schwarz inequality, (3.6), (3.7) and that M(n)
(C1 h)−1 , we obtain   J +n  



   qjn yn − qjn zn γφ zjn Vjn     j =0  M(n)−1  γ (znn )V nn − γ (znn )V nn  φ j φ j jl+1 jl−1  1 

l+1 l−1 n h∆j n =   l  2  l=1 h  M(n)−1 1/2  M(n)−1 1/2

1

 n
n  n 2 1
n 2 n γφ zj n Vj n − γφ zj n Vj n h ∆j n
l l+1 l+1 l−1 l−1 2 h l=1 l=1  M(n)−1 1/2

1

 n 
  2
C ∆ ∞ γφ zjnn Vjnn − γφ zjnn Vjnn l+1 l+1 l−1 l−1 h l=1  

  1/2 

C ∆n ∞ M(n) − 1 h2p−1 + h
C ∆n ∞ , as desired.

2

Let be R a positive constant and we denote by BAh (u˜ h , Rhp ) ⊂ Ah the open ball with center u˜ h and radius Rhp , 1 < p < 2. Next, we introduce the operator
 Φ h : BAh u˜ h , Rhp → Bh ,

 (3.8) Φ h y0 , V0 , . . . , yN , VN = Y0 , P0 , Y1 , P1 , P0 , Y2 , P2 , . . . , YN , PN , defined by the following equations: Y0 = y0 − X0 ∈ RJ −1 ,

(3.9)

P =V −U ∈R ,

(3.10)

Y =y −X ∈R ,

(3.11)

0

0

1

J

0

1

1

J

P =V −U ∈R 1

1

1

J +1

.

(3.12)

Vectors X0 , U0 , X1 and U1 represent approximations at t = 0 and t = k, respectively, to the initial grid nodes and to the theoretical solution at them. Also, P0n = V0n −

Y1n

Qn (yn , α(yn , Vn)Vn ) , g(0, Qn (yn , γ ng Vn ), tn )

(3.13)

 
  1 n k
n−1 y1 − kg g 0, Qn−1 yn−1 , γ n−1 = , tn−1 , g V k 2 n−1 n−2 ) − Qn−2 (yn−2 , γ n−2 ) 3Qn−1 (yn−1 , γ n−1 g V g V

2 Yjn =

k , tn−1 + 2


n−1 n−1
n−1 n−1 n−1   1 n yj − yjn−2 y ,γg V , tn−1 , −2 − 2kg yj −1 , Q 2k

 ,

2
j
J + n − 1,

(3.14) (3.15)

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

P1n =

Pjn =

  
  1 n k
n−1 V1 − V0n−1 exp −kµ∗ g 0, Qn−1 yn−1 , γ n−1 , tn−1 , g V k 2 n−1 n−1 n−1 n−1 n−2 ) 3Q (y , γ µ V ) − Qn−2 (yn−2 , γ n−2 µ V , 2  n−1 n−2 ) − Qn−2 (yn−2 , γ n−2 ) 3Qn−1 (yn−1 , γ n−1 k g V g V , tn−1 + , 2 2

(3.16)




 1 n ∗ n−1 n−1 n−1 n−1 Vj − Vjn−2 y , γ n−1 , µ V −2 exp −2kµ yj −1 , Q 2k
  n−1 , tn−1 , 2
j
J + n − 1, Qn−1 yn−1 , γ n−1 g V

2
n
N . Where, with the notation introduced in Section 2, 





α yn , Vn = α0 yn , Vn , α1 yn , Vn , . . . , αJ +n−1 yn , Vn , αJ +n yn , Vn ,


  αj yn , Vn = α yjn , Qn yn , γ nα Vn , tn , 0
j
J + n, 2
n
N.

305

(3.17) 2
n
N;

The next result shows that operator (3.8) is well defined. Proposition 3. Assume that hypotheses (H1)–(H5) hold and that the quadrature rules used in (3.13)– (3.17) satisfy the properties (P1)–(P8). If 

0 0 1 1 X , V , X , V , . . . , XN , VN ∈ BAh u˜ h , Rhp , where R is a fixed positive constant and 1 < p < 2, then, for h sufficiently small,
 Qn Xn , γ nφ Vn ∈ Dφ , 0
n
N,

(3.18)

and n+1 ) − Qn (Xn , γ nφ Vn ) 3Qn+1 (Xn+1 , γ n+1 φ V

2

∈ Dφ ,

0
n
N − 1, where φ = g, µ, α.

(3.19)

Proof. The definition of Qn , the hypotheses (H1)–(H5), the properties (P1)–(P6) and that Vn is bounded, allow us to obtain   n
n n n Q X , γ V − Iφ (tn ) φ 


Qn Xn , γ nφ Vn − Qn xn , γ nφ Vn   

 
 + Qn xn , γ nφ Vn − Qn xn , γ nφ un  + Qn xn , γ nφ un − Iφ (tn )   J +n   J +n 



n  n  




n 
n  n   n n n n n n qj x γφ Xj − γφ xj Vj  qj X − qj x γφ Xj Vj  + 
     j =0 j =0   J +n 


  n n n n n  qj x γφ xj Vj − uj  + o(1) +   j =0


CRhp (J + n + 1)γφ ∞ V∞ + CqRhp+1 (J + n + 1)V∞ + Rqhp+1 (J + n + 1)γφ ∞ + o(1),

0
n
N, h → 0.

(3.20)

306

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

Therefore (3.18) holds, for h sufficiently small. On the other hand, (3.19) is easily derived from (3.20) and the properties of the extrapolation procedure. 2 The operator Φ h clearly describes the aggregation grid nodes method, and
  Uh = X0 , U0 , X1 , U1 , . . . , XN , UN ∈ Ah , satisfies


 U h = 0 ∈ Bh , Φh 

(3.21)

if and only if it is a numerical solution of the scheme defined by (2.6)–(2.15). However, it is not clear how the selection grid nodes method can be written in terms of an operator defined as in (3.8). We will use the ideas presented at the end of the Section 2.2. First, we determine the indexes of the nodes used when they are considered as a subset, Sn , of the whole possible nodes, In , at each time level tn , 0
n
N . By definition,  In = {0, 1, 2, . . . , J + n}, 0
n
N, (3.22) S1 = I1 , S2 = I2 , Sn ⊂ In , n > 2. S0 = I0 , Let Xn ∈ B∞ (xn , Rhp ), n  2. For n > 2, if we have determined Si , 0
i
n − 1, and jln−1 is the first index belonging to Sn−1 such that   n−1   n−1  X n−1 − X n−1  X n−1 − X n−1 min (3.23) n−1 = n−1 , j j j j l+1

l−1

1
m
J +1

m+1

m−1

{j0n , j1n , . . . , jJn+2 }

(the selection criterion) then Sn =  n j1n = 1, j0 = 0,    n n−1  + 1, 2
m
l, jm = jm−1 n n−1  jm = jm + 1, l + 1
m
J + 1,    n jJ +2 = J + n.

is defined by

(3.24)

By means of (3.22)–(3.24) the selection procedure is completely determined. Next, we establish the following auxiliary result. Proposition 4. Let be Xn ∈ B∞ (xn , Rhp ), 0
n
N , where R is a fixed positive constant and 1 < p < 2, then     n x n − x nn   X nn − X nn  − 2Rhp , (3.25) jm+1 jm−1 jm+1 jm−1   n   n n n p X n − X n   x n − x n  − 2Rh , (3.26) j j j j m+1

m−1

m+1

m−1

1
m
J + 1. Proof. By the triangle inequalities,       n x n − x nn   X nn − X nn  − x nn − X nn + X nn − x nn  jm+1 jm−1 jm+1 jm−1 jm+1 jm+1 jm−1 jm−1    Xjnn − Xjnn  − 2Rhp , 1
m
J + 1. m+1

m−1

In a similar way (3.26) can be derived.

2

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

307

Now, we have to show that the subgrids defined by (3.22)–(3.24), satisfy the property (SR). This fact is a consequence of the following property of the indexes of the subgrid nodes. Proposition 5. Assume that the hypotheses (H1)–(H5) hold. Let be R and p positive constants with 1 < p < 2. If Xn ∈ B∞ (xn , Rhp ), 0
n
N , and K is a nonnegative integer with K > max {1, (2 + 4R)/C1 }. Then, for h sufficiently small, n − jmn
K, jm+1

0
m
J + 1, 0
n
N.

(3.27)

Proof. Suppose that (3.27) is not true. Thus, let n∗ be the first time level such that there exists l∗ , 1
l∗
J + 1, with jln∗∗+1 − jln∗∗ > K.

(3.28)

Clearly n∗ > 2. Now, if we consider the nodes xjnn∗∗ and xjnn∗∗ , the corresponding characteristic curves l∗

l∗ +1

= x(tn∗ −1 ; tn∗ , xjnn∗∗ ) and xjnn∗∗−1−1 = x(tn∗ −1 ; tn∗ , xjnn∗∗ ), cannot be consecutive in the nodes at tn∗ −1 , xjnn∗∗−1 −1 l∗ l∗ l∗ +1 l∗ +1 order of the subgrid at the time level tn∗ −1 due to the definition of the subgrids (3.22)–(3.24) and because and xjnn∗∗−1−1 correspond to n∗ is the first time level in which property (3.28) occurs. Therefore, xjnn∗∗−1 −1 l∗

l∗ +1

n∗ −1 n∗ −1 n∗ −1 x nn∗∗−1 −1 and x n∗ −1 , respectively, and the eliminated node at the time level tn∗ −1 is X n∗ −1 . Thus, X n∗ −1 is jl∗ −1

jl∗ +1

the first node that satisfies   n∗ −1  X n −1 − X nn∗ −1 min ∗ ∗ −1 = jl∗ +1

jl∗ −1

1
m
J +1

jl∗

 n∗ −1  X n −1 − X nn∗ −1  −1 . ∗ jm+1

(3.29)

∗ jm−1

−1 −1 + 1, jln∗∗+1 = jln∗∗+1 + 1 and Now using (3.24) we have that jln∗∗ = jln∗∗−1

 −1 −1 − jln∗∗−1 = jln∗∗+1 − 1 − jln∗∗ − 1 = jln∗∗+1 − jln∗∗ > K. jln∗∗+1

(3.30)

Then by means of (3.26), (3.29), (3.30) and Proposition 1, we obtain   n∗ −1   n∗ −1   n∗ −1 n∗ −1  n∗ −1  p p    X n −1 − X nn∗ −1 ∗ ∗ −1  X n∗ −1 − X n∗ −1  x n∗ −1 − x n∗ −1 − 2Rh > KC1 h − 2Rh , jm+1

jm−1

jl∗ +1

jl∗ −1

jl∗ +1

jl∗ −1

1
m
J + 1. By the above inequality and (3.25), we have   n∗ −1   n∗ −1 n∗ −1  p p   x n −1 − x nn∗ −1 ∗ ∗ −1  X n∗ −1 − X n∗ −1 − 2Rh > KC1 h − 4Rh , jm+1

jm−1

jm+1

jl∗

jm−1

1
m
J + 1.

Hence   [ J2 ]

 n∗ −1 
 J n∗ −1   + 1 KC1 h − 4Rhp > J h = 1, x n∗ −1 − x n∗ −1 > 1 j2m+2 j2m 2 m=0 which is a contradiction. Therefore the property (3.27) holds.

2

Thus, the property (SR) is immediate from Propositions 1 and 5.

308

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

Proposition 6. Assume that the hypotheses (H1)–(H5) hold. Let be R and p positive constants with 1 < p < 2. If Xn ∈ B∞ (xn , Rhp ), 0
n
N , and K is a nonnegative integer with K > max {1, (2 + 4R)/C1 }. Then, for h sufficiently small, the subgrids defined by the selection technique given by (3.22)–(3.24) satisfy the property (SR). Now, let Φ h be the operator defined with quadrature rules based on nodes which belong to subgrids defined by (3.22)–(3.24). Note that, as it was mentioned in Section 2, Φ h takes into account the whole possible nodes at each time level and their corresponding solution values. However, at each time level tn , n  3, the subgrid is formed by the zero node and, also, by the nodes which belong to the characteristic curves corresponding to the selected nodes at tn−1 . Hence, when a node is removed at a time level, the points belonging to its characteristic curve in the subsequent time levels cannot be subgrid nodes. Therefore, the nodes removed play no role in the numerical computation because those that carry out the numerical integration are determined by (3.22)–(3.24). If  Uh = (X0 , U0 , X1 , U1 , . . . , XN , UN ) ∈ BAh (u˜ h , Rhp ), satisfies
 (3.31) U h = 0 ∈ Bh , Φh  the nodes and the corresponding values of the solution at such nodes of  Uh , that belong to the subgrids defined in (3.22)–(3.24) are a numerical solution of the scheme defined by (2.17)–(2.22). On the other hand, the numerical solution of the scheme defined by (2.17)–(2.22) together with the values of the numerical integration along the characteristic curves that were removed by the selection criterion (ordered in the correct way) satisfy (3.31).

4. Consistency We define the local truncation error as lh = Φ h (u˜ h ) ∈ Bh , and we say that the discretization (3.8) is consistent if, as h → 0,   lim Φ h (u˜ h )Bh = lim lh Bh = 0. The next theorem establishes the consistency of the numerical scheme defined by Eqs. (3.9)– (3.17). Theorem 7. Assume that hypotheses (H1)–(H5) hold and that the quadrature rules considered satisfy properties (P1)–(P8). Then, as h → 0, the local truncation error satisfies,          
 Φ h (u˜ h ) = u0 − U0  + u1 − U1  + x0 − X0  + x1 − X1  + O h2 + k 2 . (4.1) Bh ∞ ∞ ∞ ∞ Proof. We denote Φ h (u˜ h ) = (Z0 , L0 , Z1 , L1 , L0 , Z2 , L2 , . . . , ZN , LN ). First, we set the bounds for Ln , 2
n
N . By means of (2.5) and (3.17), the regularity hypotheses (H1)–(H5), the properties of the quadrature rules (P1) and (P3), and the error bound of the mid-point quadrature rule, we have

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

309

   tn+2  n   n+2  

 1 ∗ n L 
u  exp − µ x τ ; tn , xj −2 , Iµ (τ ), Ig (τ ), τ dτ j 2k j −2  tn

 


    n+1 n+1 n+1 n+1 n+1 n+1 n+1 n+1 , Q , γ u , γ u x , Q x , t − exp −2kµ∗ xjn+1  n+1 µ g −1  C
2k

 tn+2  

  µ∗ x τ ; tn , xjn−2 , Iµ (τ ), Ig (τ ), τ dτ   tn

 

   − 2kµ∗ x tn+1 ; tn , xjn−2 , Iµ (tn+1 ), Ig (tn+1 ), tn+1   
n+1 
n+1 n+1
n+1 n+1 n+1   x ,γµ u , tn+1  + 2k µ xj −1 , Iµ (tn+1 ), tn+1 − µ xj −1 , Q 


  + 2k gx x n+1 , Ig (tn+1 ), tn+1 − gx x n+1 , Qn+1 xn+1 , γ n+1 un+1 , tn+1  j −1

j −1



g


 C 3 n+1  k + 2k Iµ (tn+1 ) − Qn+1 xn+1 , γ n+1 µ u 2k 
 n+1  + 2k Ig (tn+1 ) − Qn+1 xn+1 , γ n+1 g u 

C k 2 + h2 , 2
j
J + n + 1, 0
n
N − 2.




(4.2)

With respect to Ln1 , 2
n
N , we use (3.16) and the same arguments than for (4.2) together with the error bounds for the extrapolation considered and for the explicit Euler method to derive    tn+2  n+2  1  n+1  
L 
u  exp − µ∗ x(τ ; tn+1 , 0), Iµ (τ ), Ig (τ ), τ dτ 1 k 0  tn+1



∗

− exp −kµ




  k
n+1 g 0, Qn+1 xn+1 , γ n+1 , tn+1 , g u 2 n+1 ) − Qn (xn , γ nµ un ) 3Qn+1 (xn+1 , γ n+1 µ u 2 n+1 ) − Qn (xn , γ ng un ) 3Qn+1 (xn+1 , γ n+1 g u

2

, k , tn+1 + 2

 tn+2 
C  µ∗ x(τ ; tn+1 , 0), Iµ (τ ), Ig (τ ), τ dτ
 k  tn+1

     k k k k  , Ig tn+1 + , tn+1 + − kµ x tn+1 + ; tn+1 , 0 , Iµ tn+1 +  2 2 2 2       k k k  , tn+1 + + k µ x tn+1 + ; tn+1 , 0 , Iµ tn+1 + 2 2 2 ∗

   

310

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

  
  k
k k  n+1 , t − µ g 0, Qn+1 xn+1 , γ n+1 u + + , t , I t n+1 µ n+1 n+1 g 2 2 2      k

  k k n+1 , t u + + , t , I t + k µ g 0, Qn+1 xn+1 , γ n+1 n+1 µ n+1 n+1 g 2 2 2 
  k
n+1 , tn+1 , − µ g 0, Qn+1 xn+1 , γ n+1 g u 2  n+1 ) − Qn (xn , γ nµ un ) 3Qn+1 (xn+1 , γ n+1 k  µ u , tn+1 + 2 2       k k k , tn+1 + + k gx x tn+1 + ; tn+1 , 0 , Ig tn+1 + 2 2 2   
  k k  k
n+1 n+1 n+1 n+1 , tn+1 + x ,γg u , tn+1 , Ig tn+1 + − gx g 0, Q 2 2 2     
  k k k
n+1 n+1 n+1 n+1  , tn+1 + x ,γg u , tn+1 , Ig tn+1 + + k gx g 0, Q 2 2 2 
  k
n+1 , tn+1 , − gx g 0, Qn+1 xn+1 , γ n+1 g u 2  n+1 ) − Qn (xn , γ ng un ) 3Qn+1 (xn+1 , γ n+1 k  g u , tn+1 + 2 2        k k
2 
C k + x tn+1 + ; tn+1 , 0 − g 0, Ig (tn+1 ), tn+1  2 2 


  n+1 , tn+1  + k g 0, Ig (tn+1 ), tn+1 − g 0, Qn+1 xn+1 , γ n+1 g u 
 
 n+1  + Iµ (tn ) − Qn xn , γ nµ un  + Iµ (tn+1 ) − Qn+1 xn+1 , γ n+1 µ u  
 
 n+1 n+1 n+1 n+1  n n n n    x ,γg u + Ig (tn ) − Q x , γ g u + Ig (tn+1 ) − Q 

C k 2 + h2 , 0
n
N − 2.

(4.3)

Now in order to find an estimate for the boundary terms, due to the hypotheses (H1) and (H5), and the property (P3), we have that  

 
 g 0, Qn+2 xn+2 , γ n+2 un+2 , tn+2 − g 0, Ig (tn+2 ), tn+2 un+2 
Ch2 . (4.4) g 0 Also, the hypothesis (H4) and the properties (P2), (P4) and (P5) allow us to write   1 

n+2 n+2  n+2  
 n+2 n+2 x ,α x ,u u   α x, Iα (tn+2 ), tn+2 u(x, tn+2 ) dx − Q   0

  1 J

+n+2  


   qjn+2 xn+2 α xjn+2 , Iα (tn+2 ), tn+2 un+2
 α x, Iα (tn+2 ), tn+2 u(x, tn+2 ) dx − j    0

j =0

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

311

J +n+2  


n+2 n+2  n+2 

n+2  n+2
 n+2 n+2 n+2 n+2 + qj x α xj , Iα (tn+2 ), tn+2 uj − Q x ,α x ,u u    j =0


Ch . 2

(4.5)

Thus, by means of (3.13), the hypothesis (H5) and the inequalities (4.4), (4.5), we obtain  n+2 

 
  
L 
C g 0, Qn+2 xn+2 , γ n+2 un+2 , tn+2 un+2 − Qn+2 xn+2 , α xn+2 , un+2 un+2  g 0 0   

 

C g 0, Qn+2 xn+2 , γ n+2 un+2 , tn+2 − g 0, Ig (tn+2 ), tn+2 un+2  g

0

  1   
 

  +  α x, Iα (tn+2 ), tn+2 u(x, tn+2 ) dx − Qn+2 xn+2 , α xn+2 , un+2 un+2    0


Ch , 2

0
n
N − 2.

(4.6)

Next, analogous arguments than those used to derive (4.2), (4.3) lead us to establish the bound for the truncation errors produced by the grid nodes,  n+2  
Z 
C k 2 + h2 , 1
j
J + n + 1, 0
n
N − 2. (4.7) j We omit the proof because it is straightforward from the above arguments. Therefore, (4.1) follows from (4.2), (4.3), (4.6) and (4.7). 2

5. Stability Another notion that plays an important role in the analysis of the numerical method is the stability with h-dependent thresholds. For h ∈ H , let Rh be a real number (the stability threshold) with 0 < Rh < ∞; we say that the discretization (3.8) is stable for u˜ 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  h in that ball, Vh , W BAh (u˜ h , Rh ) is contained in the domain of Φ h , and, for all   
 
   h 
S Φ h   h . Vh − Φ h W Vh − W We begin with the next auxiliary result. Proposition 8. Assume that hypotheses (H1)–(H5) hold and that the quadrature rules considered satisfy properties (P1)–(P8). Let be yn , zn ∈ B∞ (xn , Rhp ) and Vn , Wn ∈ B∞ (un , Rhp ). Then, as h → 0,      n
n n n

 Q y , γ V − Qn zn , γ n Wn 
C Vn − Wn  + yn − zn  , 2
n
N, (5.1) φ φ 1 ∞ where φ could be g, µ and α. Proof. The triangle inequality, the hypothesis (H2), the properties (P5) and (P7) and that Wn ∞
C, yield

312

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

 n
n n n
 Q y , γ V − Qn zn , γ n Wn  φ φ 

 


Qn yn , γ nφ Vn − Wn  + Qn yn , γ nφ Wn − Qn zn , γ nφ Wn   J +n   J +n  J +n

 n  


n 

n 
n  n  


n
n
n 
n  n  n n n hVi − Wi  +  qi y γφ yi − γφ zi Wi  +  qi y − qi z γφ zi Wi 
C     i=0 i=0 i=0    
 n
C V − Wn 1 + yn − zn ∞ , 2
n
N, as desired.

2

Next, we introduce the theorem that establishes the stability of the discretization defined by the equations (3.9)–(3.17). Theorem 9. Assume that hypotheses (H1)–(H5) hold and that the quadrature rules considered satisfy properties (P1)–(P8). Then, the discretization is stable for u˜ h with Rh = Rhp , 1 < p < 2. Proof. We denote

 Φ h y0 , V0 , y1 , V1 , . . . , yN , VN = Y0 , P0 , P1 , Y1 , P0 , Y2 , P2 , . . . , YN , PN ,

 Φ h z0 , W0 , z1 , W1 , . . . , zN , WN = Z0 , R0 , Z1 , R1 , R0 , Z2 , R2 , . . . , ZN , RN , (y0 , V0 , y1 , V1 , . . . , yN , VN ), (z0 , W0 , z1 , W1 , . . . , zN , WN ) ∈ BAh (u˜ h , Rh ). Now, we set En = Vn − Wn ∈ RJ +n ,

∆n = yn − zn ∈ RJ +n+1 ,

0
n
N.

By means of (3.15), the hypothesis (H5) and (5.1) enable us to write    n+2   n  ∆ 
∆  + 2k Y n+2 − Z n+2  j −2 j j j   n+1
n+1 n+1 n+1 
  n+1    y ,γg V − Qn+1 zn+1 , γ n+1 + Ck ∆n+1 g W j −1 + Q        
∆nj−2  + Ck ∆n+1 ∞ + En+1 1 + Ck Yjn+2 − Zjn+2 , 2
j
J + n + 1, 0
n
N − 2.

(5.2)

Next, (3.14), the hypothesis (H5) and (5.1) allow us to obtain    n+2 
   
∆ 
k Y n+2 − Z n+2  + Ck k g 0, Qn+1 yn+1 , γ n+1 Vn+1 , tn+1 g 1 1 1
 
n+1 , tn+1  − g 0, Qn+1 zn+1 , γ n+1 g W 

 n+1 n+1  − Qn+1 zn+1 , γ n+1 + Qn+1 yn+1 , γ n+1 g V g W 

 + Qn yn , γ ng Vn − Qn zn , γ ng Wn           
k Y1n+2 − Z1n+2  + Ck ∆n+1 ∞ + En+1 1 + ∆n ∞ + En 1 

 n+1 n+1  + k Qn+1 yn+1 , γ n+1 − Qn+1 zn+1 , γ n+1 g V g W          
Ck ∆n+1  + En+1  + ∆n  + En  + k Y n+2 − Z n+2 , ∞

1

∞

Thus, when N  n > j  2, from (5.2), we have that

1

1

1

0
n
N − 2. (5.3)

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

313

• if j is even, j/2−1



     j/2−1   n En−1−2l  + ∆n−1−2l  Ck Y n−2l − Z n−2l , ∆ 
Ck j j −2l j −2l 1 ∞ l=0

(5.4)

l=0

• if j is odd, (5.3) yields (j −1)/2−1 (j −1)/2−1



       n   n−j +1  En−1−2l  + ∆n−1−2l  + Ck  + Ck Y n−2l − Z n−2l  ∆ 
∆ j 1 j −2l j −2l 1 ∞ l=0

l=0 (j −1)/2

       Y n−2l − Z n−2l 
Ck En−j −1 1 + ∆n−j −1 ∞ + Ck j −2l j −2l l=0 (j −1)/2

+ Ck

     En−1−2l  + ∆n−1−2l  . 1 ∞

(5.5)

l=0

Therefore, when N  n > j  1, by means of (5.4) and (5.5), we establish that  n−1  n−1 n



 m  m  m   n m ∆ 
C k E 1 + k ∆ ∞ + k Y − Z ∞ . j m=n−j −1

m=n−j −1

(5.6)

m=n−j +1

On the other hand, when J + n − 1  j  n  2, due to (5.2) it follows that • if n is even, n/2−1 n/2−1



       n  0  En−1−2l  + ∆n−1−2l  + Ck Y n−2l − Z n−2l , ∆ 
∆  + Ck j j −n j −2l j −2l 1 ∞ l=0

(5.7)

l=0

• if n is odd, (n−1)/2−1 (n−1)/2−1



        n  1 En−1−2l  + ∆n−1−2l  + Ck  + Ck Y n−2l − Z n−2l . ∆ 
∆ j j −n+1 j −2l j −2l 1 ∞ l=0

l=0

(5.8) Thus, when J + n − 1  j  n  2, (5.7) and (5.8) yield  n−1  n−1 n

 



 1      n   0 ∆ 
∆  + ∆  + C k Em 1 + k ∆m ∞ + k Ym − Zm ∞ . j ∞ ∞ m=1

m=1

Then, by means of (5.6) and (5.9), we can conclude that      n ∆ 
∆0  + ∆1  ∞ ∞ ∞  n−1  n−1 n

 



    +C k Em 1 + k ∆m ∞ + k Ym − Zm ∞ , m=0

m=0

On the other hand, from (3.17) we arrive at

m=2

(5.9)

m=2

2
n
N.

(5.10)

314

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

 n+2   n 


   E 
E  exp −2kµ∗ y n+1 , Qn+1 yn+1 , γ n+1 Vn+1 , Qn+1 yn+1 , γ n+1 Vn+1 , tn+1 j −2 µ g j j −1  



  n+1 n+1 n+1 n+1 + Wjn−2  exp −2kµ∗ yjn+1 y , γ n+1 , Qn+1 yn+1 , γ n+1 , tn+1 µ V g V −1 , Q

 n+1
n+1 n+1 n+1  
n+1 n+1 n+1 z , γ n+1 ,Q z ,γg W , tn+1  − exp −2kµ∗ zjn+1 µ W −1 , Q   (5.11) + 2k P n+2 − R n+2 , 2
j
J + n + 1, 0
n
N − 2. j

j

Now, by means of the hypotheses (H3) and (H5), we have


 
n+1 n+1 n+1 n+1 y , γ n+1 , Qn+1 yn+1 , γ n+1 , tn+1
1 + Ck, exp −2kµ∗ yjn+1 µ V g V −1 , Q 2
j
J + n + 1, 0
n
N − 2.

(5.12)

Thus, (5.1), (5.11), (5.12), the hypotheses (H3) and (H5), and that Wn ∞
C, enable us to write      n+2  E 
2k P n+2 − R n+2  + (1 + Ck)E n  j −2 j j j 

 n+1
n+1 n+1 n+1   n+1 n+1 n+1 y , γ n+1 ,Q y ,γg V , tn+1 + Ck Wjn−2 µ∗ yjn+1 µ V −1 , Q

 n+1
n+1 n+1 n+1   n+1 n+1 n+1 z , γ n+1 ,Q z ,γg W , tn+1  − µ∗ zjn+1 µ W −1 , Q         
(1 + Ck)E n  + Ck En+1  + ∆n+1  + Ck P n+2 − R n+2 , j −2

2
j
J + n + 1, 0
n
N − 2.

1

∞

j

j

(5.13)

Now, by means of (3.16), (5.1), a similar inequality to (5.12), the hypotheses (H3) and (H5), and that Wn+1 ∞ < C, we obtain    n+2   n+1 

  E 
E  exp −kµ∗ k g 0, Qn+1 yn+1 , γ n+1 Vn+1 , tn+1 , g 1 0 2 n+1 ) − Qn (yn , γ nµ Vn ) 3Qn+1 (yn+1 , γ n+1 µ V , 2  n+1 ) − Qn (yn , γ ng Vn ) 3Qn+1 (yn+1 , γ n+1 k g V , tn+1 + 2 2     
  k
n+1 , tn+1 , + W0n+1  exp −kµ∗ g 0, Qn+1 yn+1 , γ n+1 g V 2 n+1 ) − Qn (yn , γ nµ Vn ) 3Qn+1 (yn+1 , γ n+1 µ V , 2 n+1 ) − Qn (yn , γ ng Vn ) 3Qn+1 (yn+1 , γ n+1 k g V , tn+1 + 2 2  

  k n+1 , tn+1 , − exp −kµ∗ g 0, Qn+1 zn+1 , γ n+1 g W 2 n+1 ) − Qn (zn , γ nµ Wn ) 3Qn+1 (zn+1 , γ n+1 µ W , 2  n+1 ) − Qn (zn , γ ng Wn ) 3Qn+1 (zn+1 , γ n+1 k  g W , tn+1 + 2 2  
n+2 + k P1 − R1n+2

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

          
(1 + Ck)E0n+1  + Ck En 1 + ∆n ∞ + En+1 1 + ∆n+1 ∞   + k P n+2 − R n+2 , 0
n
N − 2. 1

1

315

(5.14)

Now from (3.13) and the hypothesis (H5) it follows that   n+2   n+2 E 
P − R0n+2  0 0   n+2 n+2  Q (y , α(yn+2 , Vn+2 )Vn+2 ) Qn+2 (zn+2 , α(zn+2 , Wn+2 )Wn+2 )    − + g(0, Qn+2 (yn+2 , γ g Vn+2 ), tn+2 ) g(0, Qn+2 (zn+2 , γ g Wn+2 ), tn+2 )   
P0n+2 − R0n+2  


 
  n+2 + C g 0, Qn+2 zn+2 , γ n+2 , tn+2 Qn+2 yn+2 , α yn+2 , Vn+2 Vn+2 g W
 
 

n+2 , tn+2 Qn+2 zn+2 , α zn+2 , Wn+2 Wn+2  − g 0, Qn+2 yn+2 , γ n+2 g V 

 
C Qn+2 yn+2 , α yn+2 , Vn+2 Vn+2  

 

  n+2 n+2 , tn+2 − g 0, Qn+2 yn+2 , γ n+2 , tn+2  × g 0, Qn+2 zn+2 , γ n+2 g W g V 

  n+2 , tn+2  + g 0, Qn+2 yn+2 , γ n+2 g V 


 
  × Qn+2 yn+2 , α yn+2 , Vn+2 Vn+2 − Qn+2 zn+2 , α zn+2 , Wn+2 Wn+2  , 0
n
N − 2.

(5.15)

Next, the hypothesis (H4), the property (P5), and that Vn ∞
C, enable us to obtain  n+2
n+2
n+2 n+2  n+2  
C, 0
n
N − 2. Q y ,α y ,V V

(5.16)

Also, by means of (5.1) and the hypothesis (H5) we have 

 

  g 0, Qn+2 zn+2 , γ n+2 Wn+2 , tn+2 − g 0, Qn+2 yn+2 , γ g Vn+2 , tn+2  g     
C En+2 1 + ∆n+2 ∞ , 0
n
N − 2.

(5.17)

Furthermore, the definition of αi , the hypothesis (H4) and (5.1) yield 
n+2 n+2 
 αi y , V − αi zn+2 , Wn+2    

  + C Qn+2 yn+2 , γ n+2 Vn+2 − Qn+2 zn+2 , γ n+2 Wn+2 
C ∆n+2 α α i       
C ∆n+2  + C En+2  + ∆n+2  , 0
n
N − 2.

(5.18)

i

1

∞

Next, by means of (5.18), the hypothesis (H4), the properties (P5) and (P8), and that W∞
C, we arrive at  n+2
n+2
n+2 n+2  n+2 

  Q y ,α y ,V V − Qn+2 zn+2 , α zn+2 , Wn+2 Wn+2    J +n+2  




   qin+2 yn+2 αi yn+2 , Vn+2 Vin+2 − Win+2 
   i=0   J +n+2 


n+2 n+2 

n+2  n+2   n+2 n+2 n+2 n+2 n+2 qi y αi y , V − qi z αi z , W Wi  +   i=0  J +n+2  
  


   qin+2 yn+2 − qin+2 zn+2 αi zn+2 , Wn+2 Win+2 
C En+2 1 +    i=0

316

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

 J +n+2  




n+2 n+2 
n+2  n+2   n+2 n+2 n+2 + qi y αi y , V − αi z , W Wi    i=0  n+2    n+2 
C E 1 + ∆ ∞ , 0
n
N − 2.

(5.19)

Therefore, we complete the derivation of the stability estimate for the boundary node taking into account (5.15)–(5.17) and (5.19), and the hypothesis (H5),        n+2 
 E 
C En+2  + ∆n+2  + P n+2 − R n+2 , 0
n
N − 2. (5.20) 0 0 0 1 ∞ Thus, when N  n > j  2, from (5.13), we obtain that • if j is even, j/2−1

     n n−2l  l  n−2l E 
(1 + Ck)j/2 E n−j  + Ck (1 + Ck) − R P j j −2l j −2l 0 l=0



     (1 + Ck)l En−1−2l 1 + ∆n−1−2l ∞ ,

j/2−1

+ Ck

(5.21)

l=0

• if j is odd, (5.14) yields (j −1)/2−1

     n n−2l  l  n−2l E 
(1 + Ck)(j −1)/2 E n−j +1  + Ck (1 + Ck) − R P j 1 j −2l j −2l l=0 (j −1)/2−1

+ Ck



     (1 + Ck)l En−1−2l 1 + ∆n−1−2l ∞

l=0

      
 n−j
(1 + Ck)(j +1)/2 E0  + Ck (1 + Ck)(j −1)/2 En−j −1 1 + ∆n−j −1 ∞  (j −1)/2

       n−2l  (1 + Ck)l En−1−2l 1 + ∆n−1−2l ∞ + (1 + Ck)l Pjn−2l . −2l − Rj −2l

(j −1)/2

+

l=0

l=0

(5.22) Therefore, when N  n > j  1, by means of (5.21) and (5.22), we establish that   n−1 n−1 n



 m  m  m     n m      E 
C E n−j  + k E 1+ k ∆ ∞+ k P −R ∞ . j 0 m=n−j −1

m=n−j −1

m=n−j +1

On the other hand, when J + n − 1  j  n  2, due to (5.13) it follows that • if n is even, n/2−1

     n n−2l  E 
(1 + Ck)n/2 E 0  + 2k (1 + Ck)l Pjn−2l j j −n −2l − Rj −2l l=0

(5.23)

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327



     (1 + Ck)l En−1−2l 1 + ∆n−1−2l ∞ ,

317

n/2−1

+ Ck

(5.24)

l=0

• if n is odd, (n−1)/2−1

     n n−2l   + 2k E 
(1 + Ck)(n−1)/2 E 1 (1 + Ck)l Pjn−2l j j −n+1 −2l − Rj −2l l=0



     (1 + Ck)l En−1−2l 1 + ∆n−1−2l ∞ .

(n−1)/2−1

+ Ck

(5.25)

l=0

Thus, when J + n − 1  j  n  2, we can conclude from (5.24) and (5.25) that   n−1 n−1 n



 1  m  m  m  0   n E 
C E  + E  + k E 1 + k ∆ ∞ + k P − Rm ∞ . j 1 1 m=1

m=1

(5.26)

m=2

Now, multiplying |Ejn | by h and summing in j , 0
j
J + n − 1, 2
n
N , from (5.20), (5.23) and (5.26) and that k = rh, we have n−1 +n−1  

 n J

 n  n  E   E  = hE n  + h E h + 0 j j 1 j =1

j =n

    
 
hP0n − R0n  + Ch En 1 + ∆n ∞  n−1 n−1



 n−j    h E0  + k Em 1 + +C j =1

+C

m=n−j −1

J

+n−1 j =n



m=n−j −1

m=1

m=1

n−1

n−1

  k Em 1 +

  + h∆n ∞ +

n−1

  k ∆m ∞ +

m=0

   
C E0 1 + E1 1 +



m=0

m=n−j +1

  k Pm − Rm ∞

  k ∆m ∞ +

m=2

 n−j  hE0 

n

    k Pm − Rm ∞ + hP0n − R0n 



m=2 n

m=0

+

n

j =1

m=0

n

  k ∆m ∞ +

  n−1 n−1 n



 0  1  m  m  m  m h E 1 + E 1 + k E 1 + k ∆ ∞ + k P − R ∞

     
C E0 1 + E1 1 + hEn 1 +



n−1

n

n−2

      
 n−j n−j  k Em 1 + h P0 − R0  + En−j 1 + ∆n−j ∞ j =1

  m  n   m n k P − R ∞ + hP0 − R0 

m=2

n n



        k Em 1 + k ∆m ∞
C E0 1 + E1 1 + m=0

m=0



318

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

 n n



 m  m   m m + k P − R ∞ + hP0 − R0  , m=2

Then

2
n
N.

m=2

 n n



     m    n   E 
C E0  + E1  + k E 1+ k ∆m ∞ 1 1 1 m=2

+

n

m=0

  m  m k P − R ∞ + P0 − R0 ∞ ,

2
n
N.

(5.27)

m=2

Thus, by means of the discrete Gronwall lemma,   n n



 0  1  m  m   n m E 
C E  + E  + k ∆ ∞ + k P − R ∞ + P0 − R0 ∞ , 1 1 1 m=0

2
n
N.

m=2

(5.28) Next, we substitute (5.28) into (5.10) to have   n  0  1     ∆ 
∆  + ∆  + C k E0  + k E1  ∞ ∞ ∞ 1 1 +

n−1

  m m



 0  1  l  l  l     k E 1 + E 1 + k ∆ ∞ + k P − R ∞ + P0 − R0 ∞

m=2

+

l=0

n−1 n



    k ∆m ∞ + k Ym − Zm ∞ m=0



l=2

m=2

 n

          k ∆m ∞
C ∆0 ∞ + ∆1 ∞ + E0 1 + E1 1 + m=2

 n n



 m    + k P − Rm ∞ + P0 − R0 ∞ + k Ym − Zm ∞ , m=2

2
n
N.

(5.29)

m=2

Again, by means of the discrete Gronwall Lemma, it follows that           n ∆ 
C ∆0  + ∆1  + E0  + E1  + P0 − R0 ∞ ∞

∞

∞

1

1

 n n



 m  m   m m + k P − R ∞ + k Y − Z ∞ , m=2

2
n
N.

m=2

Next, we substitute (5.30) in (5.28) to obtain            n E 
C ∆0  + ∆1  + E0  + E1  + P0 − R0 ∞ 1 ∞ ∞ 1 1

(5.30)

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

 n n



 m  m   m m + k P − R ∞ + k Y − Z ∞ , m=2

319

2
n
N.

(5.31)

m=2

And, finally we substitute (5.30) and (5.31) in (5.20), (5.23) and (5.26) to arrive at           n E 
C ∆0  + ∆1  + E0  + E1  + P0 − R0 ∞ ∞

∞

∞

1

1

 n n



 m  m   m m + k P − R ∞ + k Y − Z ∞ , m=2

2
n
N.

(5.32)

m=2

So, due to (5.30) and (5.32) we have 
0 0   ∆ , E , . . . , ∆ N , EN  A 
0 0 1 1 h  
C ∆ , E , ∆ , E , P 0 − R 0 , Y 2 − Z 2 , P 2 − R 2 , . . . , Y N − Z N , P N − R N  Bh .

2

6. Convergence The global discretization error is defined as Uh ∈ Ah . e˜ h = u˜ h −  We say that the discretization (3.8) is convergent if there exists h0 > 0 such that, for each h ∈ H with Uh for which, as h → 0, h
h0 , (3.21) has a solution      Uh Ah = lim e˜ h Ah = 0. lim u˜ h −  In our analysis, we shall use the following result of the general discretization framework introduced by López-Marcos and Sanz-Serna [20]. Theorem 10. Assume that (3.8) is consistent and stable with thresholds Rh . If Φ h is continuous in B(u˜ h , Rh ) and lh Bh = o(Rh ) as h → 0, then: (i) For h sufficiently small, the discrete equations (3.21) possess a unique solution in B(u˜ h , Rh ). (ii) As h → 0, the solutions converge and ˜eh Ah = O(lh Bh ). Finally, we enunciate the next theorem that establishes the convergence of the numerical method defined by Eqs. (3.9)–(3.17). Theorem 11. Assume that hypotheses (H1)–(H5) hold and that the quadrature rules considered satisfy properties (P1)–(P8). Then, for h sufficiently small, the numerical method defined by Eqs. (3.9)–(3.17) has a unique solution  Uh ∈ B(u˜ h , Rh ) and       
  Uh − u˜ h Ah
C x0 − X0 ∞ + u0 − U0 ∞ + x1 − X1 ∞  
 + u1 − U1 ∞ + O h2 + k 2 . (6.1)

320

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

The proof of Theorem 11 is immediately derived by means of the consistency (Theorem 7), the stability (Theorem 9) and Theorem 10. Next, we can establish an error bound for the numerical and the theoretical solution at the numerical values of the grid nodes. Theorem 12. Assume that hypotheses (H1)–(H5) hold and that the quadrature rules considered satisfy properties (P1)–(P8). For h sufficiently small, let be N
 RJ +n , ∈ u∗h = u0∗ , u1∗ , u2∗ , . . . , uN ∗ n=0

defined by





 un∗ = u X0n , tn , u X1n , tn , . . . , u XJn +n−1 , tn ∈ RJ +n ,

0
n
N , where Xjn , 0
j
J + n − 1, 0
n
N , are the grid nodes given by the scheme (3.9)– (3.17). Then,           
2 ll Un − un∗ ∞
C h2 + k 2 + x0 − X0 ∞ + u0 − U0 ∞ + x1 − X1 ∞ + u1 − U1 ∞ . (6.2) This theorem follows immediately from Theorem 11. In particular, if X0 = x0 and U0 = u0 and X1 and U1 are second order approximations to x1 and u1 , respectively, the proposed numerical schemes are second order accuracy. Remark. An interesting question is to consider how the numerical schemes analysed in this paper should be used to carry out the numerical integration of problems based on biological data. We restrict our attention to schemes that employ the composite trapezoidal quadrature rule to approximate the integral terms. First, we have to construct the data of the problem (the vital and the weight functions) in a such way that the properties (H2)–(H5) hold. Next, we consider an initial condition u0 ∈ C 2 [0, 1]. We have to take into account that three compatibility conditions between the initial and the boundary conditions are needed to have a twice continuously differentiable solution of the problem. If some of these conditions are not satisfied the solution has a singularity along the characteristic curve passing through point (0, 0), but it is twice continuously differentiable on the left and on the right of such curve. The kind of the singularity depends on the first compatibility condition that does not hold. The necessary condition for the continuity at the points (x(t; 0, 0), t), t  0 (recall (1.7)) is 
g 0, Ig (0), 0 u0 (0) =

1


 α x, Iα (0), 0 u0 (x) dx.

(6.3)

0

From a biological point of view, this condition implies that the birth law is also valid at t = 0. This condition seems reasonable in many situations and it has been used, for example, in the demographic model considered by Milner and Rabbiolo [22] and i n the size-structured model built by Sulsky [26] for a Gambussia affinis population. If (6.3) does not hold our convergence analysis cannot be applied, however convergence probably occurs. The second compatibility condition is needed for the continuity of the first derivatives at the points (x(t; 0, 0), t), t  0,

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

321




 
 gz 0, Ig (0), 0 Ig (0) + gt 0, Ig (0), 0 u0 (0) − g 0, Ig (0), 0 v(0) 1 =

 
  


αz x, Iα (0), 0 Iα (0) + αt x, Iα (0), 0 u0 (x) − α x, Iα (0), 0 v(x) dx,

(6.4)

0

where Iφ (0)

1 =

γφ (x)v(x) dx,

φ = g, α,

0

 


v(x) = g x, Ig (0), 0 u0 (x) x + µ x, Iµ (0), 0 u0 (x). If (6.3) holds but (6.4) does not, then the first derivatives have a (finite) jump at the points of the characteristic curve (x(t; 0, 0), t), t  0. When the local error is analysed, there is only a problem with the integral terms. For each time level tn , 0
n
N , if x(tn ; 0, 0) is a node of the quadrature rule then there is not a loss of accuracy. On the other hand, if x(tn ; 0, 0) is not a node of the composite rule then the trapezoidal rule on the subinterval that contains such point gives an approximation o(h) as h tends to zero. Thus, in any case, the local truncation error is lh Bh = o(h), as h tends to zero. Therefore, the nonlinear stability analysis has to be made with stability thresholds Rh such that Rh = o(h) and lh Bh = o(Rh ), as h tends to zero. Thus, at least first order of convergence can be established. On the other hand, if (6.3) and (6.4) hold the solution, in the worst case, has a jump on the second derivatives if the corresponding compatibility condition is not satisfied. But, in this situation, by means of similar arguments than those made for the previous one, we can derive that the local error is second order still and, hence, the second order of convergence is retained.

7. Numerical results We have carried out numerical experiments with the schemes defined in Section 2. In order to approximate the integral terms we have used the composite trapezoidal quadrature rule. We have obtained the initial conditions   X1 = X01 = 0, X11 , . . . , XJ1 , XJ1 +1 = 1 ,   U1 = U01 , U11 , . . . , UJ1 , UJ1+1 = 0 , by means of the well-known second order method given by the following equations: 
1/2 1/2 1/2 k 1/2 1 0 1/2 , Xj +1 = Xj + kg Xj +1 , Q X , γ g U , 2   
1/2 1/2 1/2 1/2
1/2 1/2 1/2 k 1/2 1 0 ∗ 1/2 , Uj +1 = Uj exp −kµ Xj +1 , Q X , γ µ U , Q X , γ g U , 2

(7.1) (7.2)

0
j
J − 1, and U01 =

Q1 (X1 , α(X1 , U1 )U1 ) , g(0, Q1 (X1 , γ 1g U1 ), t1 )

(7.3)

322

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

where we have used the notation of Section 2 and 1/2

1/2

X0 = 0,

1/2

XJ +1 = 1, UJ +1 = 0,

  k 1/2 Xj +1 = Xj0 + g Xj0 , Q0 X0 , γ 0g U0 , t0 , 2   k ∗
0 0
0 0 0 0
0 0 0  1/2 0 Uj +1 = Uj exp − µ Xj , Q X , γ µ U , Q X , γ g U , t0 , 2

(7.4) (7.5) (7.6)

0
j
J − 1, and 1/2

U0

=

Q1/2 (X1/2, α(X1/2, U1/2)U1/2 ) 1/2

g(0, Q1/2 (X1/2, γ g U1/2 ), k2 )

.

(7.7)

We have considered two different test problems that present meaningful nonlinearities (both from a mathematical and a biological point of view). The numerical integration for each numerical experiment was carried out on the time interval [0, 10]. Problem 1. The size-specific growth, fertility and mortality moduli are chosen as 

z
1 + exp −0.2 sin2 t exp 0.1 sin2 t , 2 1+z
 (1 + exp (−0.1 sin2 t))2 z , α(x, z, t) = 3x 2 1 − x 2 (z + 1)2 exp (−0.1 sin2 t)


g(x, z, t) = 0.1 1 − x 4

µ(x, z, t) = 0.2z

1 + x 2 + 2x 3 +

sin 2t 2

exp (−0.1 sin2 t)

.

The weight functions, φ = g, µ and α, are taken as  0
x
1/3,  c, 
3 2 γφ (x) = c(2 − 3x) 54x − 27x + 4 , 1/3 < x
2/3,  0, 2/3 < x
1, 1 , c= 21250 ln (5/3) − 21248 ln (4/3) − 71131/15 and we consider as the initial size-specific density the function u0 (x) =

1−x . 1+x

The problem (1.1)–(1.6) has the periodic solution u(x, t) given by u(x, t) =


 1−x exp −0.1 sin2 t . 1+x

Problem 2. We chose the data in this problem to obtain a solution which tends to a nontrivial equilibrium. Thus, the size-specific growth, fertility and mortality moduli are chosen as

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327


 g(x, z, t) = 0.225 1 − x 2

323

z t (1 + (0.16 + 0.22 exp (−0.225t 2 ))2 ) , 1 + z2 0.16 + 0.22 exp (−0.225t 2 )

z (1.16 + 0.22e−0.225t )2 α(x, z, t) = 0.225tx (1 − x) (1 + z)2 0.16 + 0.22e−0.225t 2 2

2

2

1 + e−0.225t , 61 − 88 log 2 + (38 log 2 − 79/3)e−0.225t 2 z . µ(x, z, t) = 1.35t 0.16 + 0.22e−0.225t 2 2

×

The weight functions, φ = g, µ and α, are taken as  0
x
1/3,  1, 
3 2 γφ (x) = (2 − 3x) 54x − 27x + 4 , 1/3 < x
2/3,  0, 2/3 < x
1, and we consider as the initial size-specific density the function u0 (x) =

(1 − x)2 1−x + . 4 (1 + x) (1 + x)3

Now, u(x, t) =

1 − x −0.225t 2 (1 − x)2 + e , 4 (1 + x) (1 + x)3

is the solution of (1.1)–(1.6). Since we know the exact solution of each problem, we can show numerically that our methods are second order accuracy by means of error tables. In addition, we can compare their computational efficiency. In Tables 1 and 2, we present the results obtained for the test problem 1 with the aggregation and the selection grid nodes method, respectively. Tables 3 and 4 show the results for the second test problem. In each entry in columns two to seven of Tables 1 and 3, the upper value represents the global error for the aggregation grid nodes method 
!  max u Xjn , tn − Ujn  . eh,k = max 0
n
N 0
j
J +n

In Tables 2 and 4, the upper value is the global error for the selection grid nodes scheme 
 
   eh,k = max max u Xj0 , t0 − Uj0 , max u Xj1 , t1 − Uj1 , 0
j
J

max

2
n
N

0
j
J +1


!!  max u Xjn , tn − Ujn  .

0
j
J +2

In all the tables, the lower number on the left is the CPU time measured in seconds and the lower number on the right is the order s of the method as computed from s=

log(e2h,2k /eh,k ) . log(2)

324

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

Table 1 Errors, CPU time (seconds) and order of aggregation grid nodes method for Problem 1 k\h

6.250E–2

6.250E–2

3.125E–2

3.161E–03 0.7

3.125E–2

8.594E–04 0.8

3.182E–03 2.5

1.563E–2

3.191E–03 8.9

7.813E–3

3.194E–03 32.6

3.906E–3 1.953E–3

2.00

7.991E–04 130.6

3.195E–03 513.1

2.00

7.984E–04 33.4

3.194E–03 129.0

2.00

7.959E–04 9.3

2.00

7.993E–04 516.1

3.277E–04 1.1

7.902E–04 2.8

1.563E–2

2.00

2.00

1.975E–04 10.1

2.00

1.990E–04 34.9

2.00

1.996E–04 133.6

2.00

1.998E–04 522.3

1.993E–04 1.5

2.148E–04 3.2

7.813E–3

2.00

2.00

5.367E–05 11.7

2.00

4.938E–05 38.0

2.00

4.975E–05 139.8

2.00

4.991E–05 534.7

1.675E–04 2.5

8.179E–05 4.1

3.906E–3

2.00

2.01

2.042E–05 15.0

2.00

1.341E–05 44.2

2.00

1.235E–05 152.1

2.00

1.244E–05 559.3

1.596E–04 4.3

4.962E–05 5.9

1.953E–3

2.00

4.168E–05 9.6

2.01

1.238E–05 21.6

2.00

5.101E–06 56.5

2.00

3.353E–06 176.7

2.00

3.086E–06 608.6

2.00

Table 2 Errors, CPU time (seconds) and order of selection grid nodes method for Problem 1 k\h

6.250E–2

3.125E–2

1.563E–2

7.813E–3

3.906E–3

1.953E–3

6.250E–2

3.363E–03

8.553E–04

2.729E–04

1.757E–04

1.621E–04

1.589E–04

0.3

0.1 3.125E–2

3.443E–03 0.5

0.3 1.563E–2

3.462E–03 3.476E–03 3.475E–03 3.478E–03 3.6

1.99 8.798E–04

3.4

1.8 1.953E–3

1.96 8.737E–04

1.7

0.9 3.906E–3

1.97 8.840E–04

0.9

0.5 7.813E–3

0.5 8.562E–04

1.99 8.797E–04

6.8

1.99

1.0

2.071E–04 1.0

2.05

2.065E–04 1.7

2.05

2.172E–04 3.3

2.03

2.180E–04 6.5

2.00

2.200E–04 13.1

2.00

1.9

6.770E–05 1.9

2.01

5.248E–05 3.4

1.98

5.250E–05 6.4

1.98

5.359E–05 12.8

2.02

5.455E–05 25.6

2.00

3.8

4.405E–05 3.7

2.00

1.668E–05 6.8

2.02

1.295E–05 12.7

2.02

1.267E–05 25.4

2.05

1.334E–05 50.8

2.01

4.026E–05 7.4

2.01

1.098E–05 13.4

2.00

4.160E–06 25.3

2.00

3.246E–06 50.6

2.00

3.166E–06 101.2

2.00

Each column and each row of the tables correspond with different values of the spatial and time discretization parameter, respectively. The results in the tables clearly confirm the expected second order of convergence for both methods. In Figs. 1 and 2, we show efficiency plots (corresponding to Problems 1 and 2, respectively) where we present the error (in the vertical axis) and the work, measured as the cpu-time spent in seconds (in the horizontal axis), in logarithmic scale. We present the results obtained when the ratio employed (r = hk ) is the most efficient for the aggregation grid nodes method among the values presented in each table.

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

325

Table 3 Errors, CPU time (seconds) and order of aggregation grid nodes method for Problem 2 k\h 1.563E–02 7.813E–03

6.250E–2 4.336E+13

6.184E+15

2.489E+14

8.3

9.1

1.791E+15 9.231E–01 97.7

1.953E–03

1.240E–02 320.8

9.766E–04

1.240E–02 1212.6

4.883E–04

1.563E–2

7.9 25.3 3.906E–03

3.125E–2

1.240E–02 4697.0

3.259E+20 25.9

–22.84 50.80 8.22 2.01 2.01

9.74 2.00 2.00

54.78

1.071E–03 337.6

9.74

1.925E–04 1245.6

7.705E–04 4728.1

–26.34

9.135E–01 106.6

7.704E–04 1228.6

3.086E–03 4710.7

68.27

2.48

1.925E–04 4759.9

3.906E–3 4.814E+14 13.7

2.120E+22 29.8

1.071E–03 328.7

3.086E–03 1219.6

–2.19

9.145E–01 101.5

3.086E–03 323.9

1.029E+14 10.6

2.831E+16 27.3

9.153E–01 98.9

7.813E–3

2.00

–17.18

9.133E–01 116.7

74.30

1.071E–03 355.7

9.74

4.809E–05 1278.9

4.48

4.812E–05 4823.7

3.926E+13 19.9

1.528E+19 35.1

1.953E–3

2.00

5.501E+18 45.5

–13.48

9.135E–01 137.1

63.86

1.071E–03 391.9

9.74

1.311E–05 1345.9

6.35

1.202E–05 4950.9

2.00

Table 4 Errors, CPU time (seconds) and order of selection grid nodes method for Problem 2 k\h

6.250E–2

3.125E–2

1.563E–2

7.813E–3

3.906E–3

1.953E–3

1.563E–02

4.677E–02

9.823E–03

2.407E–03

5.515E–02

8.785E–02

2.466E+14

0.5 7.813E–03

4.401E–02 0.9

3.906E–03

3.682E–02 1.8

1.953E–03

3.276E–02 3.5

9.766E–04

2.996E–02 6.5

4.883E–04

2.850E–02 11.3

0.9

1.7

9.774E–03 1.7

2.26

9.434E–03 3.3

2.22

8.064E–03 6.4

2.19

7.230E–03 11.2

2.18

6.816E–03 20.4

2.14

3.3

2.140E–03 3.0

2.20

2.034E–03 6.2

2.26

1.921E–03 11.5

2.30

1.731E–03 20.4

2.22

1.594E–03 39.1

2.18

6.3

5.772E–04 5.6

2.06

4.980E–04 10.9

2.10

4.225E–04 21.2

2.27

4.107E–04 39.8

2.23

3.818E–04 76.3

2.18

12.6

1.029E+00 10.9

-4.22

1.423E–04 21.5

2.02

1.191E–04 40.1

2.06

8.988E–05 78.4

2.23

9.617E–05 150.6

2.09

5.774E–04 21.5

7.25

3.099E–04 42.5

11.70

3.539E–05 76.4

2.01

2.895E–05 141.9

2.04

2.202E–05 282.8

2.03

We obtain r = 8 for the first problem, and r = 1 for the second one. Such figures show clearly the best behaviour of the selection grid nodes method. In both cases there is an improvement of the computational time spent in the calculations made to obtain the same global error. For the most expensive run, the CPU time cost is improved by about 45% for the first problem and by about 80% for the second. All the calculations were carried out on a Sun Sparc 11.

326

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

Fig. 1. Efficiency plot for Problem 1. Aggregation grid nodes method ( ) and selection grid nodes method ().

Fig. 2. Efficiency plot for Problem 2. Aggregation grid nodes method ( ) and selection grid nodes method ().

Acknowledgements The authors are grateful to Joan Saldaña for his comments about the analytical results related to the problem studied in the present paper.

O. Angulo, J.C. López-Marcos / Applied Numerical Mathematics 50 (2004) 291–327

327

References [1] L.M. Abia, O. Angulo, J.C. López-Marcos, Size-structured population dynamics models and their numerical solutions, Discrete Contin. Dyn. Syst. B 4 (3) (2004). [2] L.M. Abia, J.C. López-Marcos, Second order schemes for age-structured population equations, J. Biol. Syst. 5 (1997) 1–16. [3] L.M. Abia, J.C. López-Marcos, On the numerical integration of nonlocal terms for age-structured population models, Math. Biosci. 157 (1999) 147–167. [4] A.S. Ackleh, H.T. Banks, K. Deng, A finite difference approximation for a coupled system of nonlinear size-structured populations, Nonlinear Anal. 50 (2002) 727–748. [5] A.S. Ackleh, K. Deng, A monotone approximation for a nonlinear nonautonomous size-structured population model, Appl. Math. Comput. 108 (2000) 103–113. [6] A.S. Ackleh, K. Deng, Monotone scheme for nonlinear first-order hyperbolic initial-boundary value problems, Appl. Math. Lett. 13 (2000) 111–119. [7] A.S. Ackleh, K. Ito, An implicit finite difference scheme for the nonlinear size-structured population model, Numer. Funct. Anal. Optim. 8 (1997) 865–884. [8] O. Angulo, J.C. López-Marcos, Numerical schemes for size-structured population equations, Math. Biosci. 157 (1999) 169–188. [9] O. Angulo, J.C. López-Marcos, Numerical integration of nonlinear size-structured population equations, Ecol. Model. 133 (2000) 3–14. [10] O. Angulo, J.C. López-Marcos, Numerical integration of autonomous and nonautonomous nonlinear size-structured population models, Math. Biosci. 177–178 (2002) 39–71. [11] A. Calsina, J. Saldaña, A model of physiologically structured population dynamics with a nonlinear individual growth rate, J. Math. Biol. 33 (1995) 335–364. [12] J.M. Cushing, An Introduction to Structured Populations Dynamics, in: CMB-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, PA, 1998. [13] A.J. Fabens, Properties and fitting of the von Bertalanffy growth curve, Growth 29 (1965) 265–289. [14] W. Huyer, A size-structured population-model with dispersion, J. Math. Anal. Appl. 181 (1994) 716–754. [15] M. Iannelli, Mathematical Theory of Age-Structured Population Dynamics, in: Applied Mathematics Monographs, CNR, vol. 7, Giardini, Pisa, 1994. [16] K. Ito, F. Kappel, G. Peichl, A fully discretized approximation scheme for size-structured population models, SIAM J. Numer. Anal. 28 (1991) 923–954. [17] N. Kato, K. Sato, Continuous dependence results for a general model of size-dependent population dynamics, J. Math. Anal. Appl. 272 (2002) 200–222. [18] N. Kato, H. Torikata, Local existence for a general model of size-dependent population dynamics, Abstract Appl. Anal. 2 (1997) 207–226. [19] T. Kostova, An explicit third-order numerical method for size-structured population equations, Numer. Methods for Partial Differential Equations 19 (2003) 1–21. [20] J.C. López-Marcos, J.M. Sanz-Serna, Stability and convergence in numerical analysis III: Linear investigation of nonlinear stability, IMA J. Numer. Anal. 8 (1988) 71–84. [21] J.A.J. Metz, E.O. Dieckmann (Eds.), The Dynamics of Physiologically Structured Populations, Springer Lecture Notes in Biomath., vol. 68, Springer, Heildelberg, 1986. [22] F.A. Milner, G. Rabbiolo, Rapidly converging numerical algorithms for models of population dynamics, J. Math. Biol. 30 (1992) 733–753. [23] L.F. Murphy, A nonlinear growth mechanism in size-structured population dynamics, J. Theor. Biol. 104 (1983) 493–506. [24] A.M. de Roos, Numerical methods for structured population models: The escalator boxcar train, Numer. Methods Partial Differential Equations 4 (1988) 173–195. [25] J. Saldaña, Personal communication, 2004. [26] D. Sulsky, Numerical solution of structured population models II. Mass structure, J. Math. Biol. 32 (1994) 491–514. [27] G.F. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, Marcel Dekker, New York, 1985.