Research Article Received 15 April 2012
Published online in Wiley Online Library
(wileyonlinelibrary.com) DOI: 10.1002/mma.2675 MOS subject classification: 92D25; 35L02; 35K57; 65N99
Approximation of a population dynamics model by parabolic regularization‡ Mimmo Iannellia * † and Gabriela Marinoschib Communicated by P. Colli The basic linear model for describing an age structured population spreading in a limited habitat is considered with the purpose of investigating an approximation procedure based on parabolic regularization. In fact, a viscosity model is introduced by considering an appropriate approximating regularized parabolic problem and it is proved that the sequence of the approximating solutions tends to the solution to the original problem. The advantage of this approach is that it leads to the numerical solution of a parabolic problem that has more stable solutions than the hyperbolic-parabolic original problem and avoids the restrictions (compatibility conditions) needed to treat the latter. Moreover, for the solution of the approximating problem, it is possible to take advantage of established software packages dedicated to parabolic problems. Some examples of the approach are provided using COMSOL Multiphysics. Copyright © 2012 John Wiley & Sons, Ltd. Keywords: population dynamics; age structure; diffusion; numerical approximation
1. Introduction We consider the linear model of an age-structured population diffusing in a bounded habitat with no exchange of individuals with the outer environment @ @n @n @n C K.a, y/ C .a/n D g, in .0, T/ , (1) @t @a @y @y @n @n .t, a, 0/ D .t, a, L/ D 0, @y @y Z n.t, 0, y/ D
in .0, T/ .0, a /,
(2)
in .0, T/ .0, L/,
(3)
a
ˇ0 .a, y/n.t, a, y/da 0
n.0, a, y/ D n0 .a, y/
in
,
(4)
where t 0 represents time, a 2 Œ0, a age (we assume the maximum age a < C1), y 2 Œ0, L the space position, D .0, a / .0, L/. The state variable n.t, a, y/ is the age-space distribution of the number of individuals at time t, the function ˇ0 .a, y/ is the intrinsic fertility function, and .a/ is the intrinsic mortality rate of the population. Finally, the diffusion constant K.a, y/ drives the diffusion process and g.t, a, y/ represent a distributed in-flux. We note that the mortality rate .a/ is related to the survival probability ….a/ D e
Ra 0
./d
,
a Mathematics Department, Trento University, Via Sommarive 14, 38123 Trento, Italy b Institute of Mathematical Statistics and Applied Mathematics, Bucharest, Romania *Correspondence to: Mimmo Iannelli, Mathematics Department, Trento University, Via Sommarive 14, 38123 Trento, Italy. † E-mail:
[email protected] ‡ This work has been performed within the PRIN 2007 project ‘Mathematical Population Theory: methods, models, comparison with data’. Gabriela Marinoschi has been supported by the Grant CNCSIS PCCE-55/2008 ‘Sisteme diferentiale in analiza neliniara si aplicatii’.
Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI that must satisfy the condition ….a / D 0. Thus .a/ must be unbounded at a , and usually the assumption 2 L1loc .Œ0, a // is adopted. However, we perform the standard transformation p.t, a, y/ D
n.t, a, y/ ….a/
and obtain the following problem on p.t, a, y/ @p @p @ C @t @a @y
K.a, y/
@p @y
Df
in .0, T/ ,
@p @p .t, a, 0/ D .t, a, L/ D 0 @y @y Z
(5)
in .0, T/ .0, a /,
(6)
in .0, T/ .0, L/,
(7)
in ,
(8)
a
ˇ.a, y/p.t, a, y/da
p.t, 0, y/ D 0
p.0, a, y/ D p0 .a, y/
where, now, ˇ.a, y/ is the maternity function that summarizes the demographic process ˇ.a, y/ D ˇ0 .a, y/….a/,
(9)
and f .t, a, y/, p0 , respectively, are a normalized in-flux and a normalized initial datum, f .t, a, y/ D
g.t, a, y/ , ….a/
p0 .a, y/ D
n0 .a, y/ . ….a/
The transformed system (5)–(8), we will be concerned with, is equivalent to the standard form (1)-(4). In fact, once we have a solution for (5)–(8), we can compute n.t, a, y/ through the formula n.t, a, y/ D p.t, a, y/….a/.
(10)
Concerning the basic conditions on all the functions in (5)-(8), we assume K 2 C 1 ./, 0 < K0 K.a, y/ KC ,
(11)
ˇ 2 C 1 ./, 0 ˇ.a, y/ ˇC ,
(12)
f 2 C 1 .Œ0, T /.
(13)
System (1)–(4) constitutes the basic linear model for an age structured population (see for instance [1, 2]) spreading by a random mechanism. Its numerical treatment has been investigated by several authors adopting different strategies to approach the hyperbolic nature of the problem, due to the age structure. Early studies on this subject, focused on problems without diffusion, have designed several methods, mainly based on the methodology of characteristics. We may refer, for instance, to [3, 4] for a recent review, and to [5–7] for further recent work, presenting different approaches to linear and nonlinear problems, also dealing with the more general case of size structure. Problems including both age-structure and diffusion have been also considered from the numerical point of view ([8, 9] for very early approaches) and have recently attracted more attention in connection with the modeling of population mechanisms that are significantly depending on the interplay of both the age structure and the geographical spread [10–16]. Also, these methods are mainly based on integration along characteristics and are concerned with designing efficient ways to conjugate age and space discretization. The aim of this work is to develop a procedure for computing the solution to this problem by a non-standard technique. In this respect, we shall introduce a viscosity model by considering an appropriate approximating regularized parabolic problem and prove that the sequence of the approximating solutions tends to the solution to the original problem (5)–(8). The advantage of this method is that it leads to the numerical solution of a parabolic problem that has more stable solutions than the hyperbolic-parabolic original problem and avoids the restrictions (compatibility conditions) needed to treat the latter. Indeed, a condition of compatibility R a requiring p0 .a, y/ to verify (7) is hardly satisfied. For example, when p0 is a constant or a function of y then ˇ should be such that 0 ˇ.a/ da D 1, which introduces a restriction. We see that it is more difficult to find p0 depending also on a that satisfies (7). By introducing an approximating parabolic problem, there is no need to be constrained by such a compatibility condition. We shall prove the existence of a unique strong solution to the approximating problem indexed upon > 0 and then its convergence to a weak solution to the original problem (5)–(8), if p0 is sufficiently smooth. Also, we prove that if p0 2 L2 ./, then the mild solution to the approximating problem converges to a mild solution to the original problem. A stronger convergence in the latter case seems not possible to be obtained because of the fact that the approximating operator that is of elliptic type looses regularity by passing to the limit as ! 0. This is also put into evidence by the numerical results obtained with COMSOL Multiphysics [17]. Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI
2. Setting of the problem For simplicity, we denote V D H1 .0, L/, V 0 D .H1 .0, L//0 , where V is endowed with the standard norm, and the duality between these two spaces is denoted by h, iV 0 ,V . Then we consider V D L2 .0, a ; V/, V 0 D L2 .0, a ; V 0 /, with Z hf , iV 0 ,V D
a
0
hf .a/, .a/iV 0 ,V , for f 2 V 0 , 2 V.
(14)
This problem will be treated in the space E D L2 ./, within the same functional framework as in [18]. Namely, we define first A0 : D.A0 / V ! V 0 , by Z a @v ˇ.a, y/v.a, y/da, 2 V 0 , v.0, y/ D D.A0 / D v 2 V, @a 0 and hA0 v, iV 0 ,V D
@v , @a
Z C
V 0 ,V
K.a, y/vy y dyda,
8 2 V.
(15)
Next, we define the restriction of A0 to E A : D.A/ E ! E, D.A/ D fv 2 D.A0 /, A0 v 2 Eg ,
Av D A0 v
for any
(16)
v 2 D.A/.
Correspondingly, we have the Cauchy problem dp C Ap D f a.e. t 2 .0, T/ dt
(17)
p.0/ D p0 .
(18)
About this problem, we know that A is the infinitesimal generator of a semigroup etA such that tA e e˛t ,
(19)
with ˛ > 0, so that, for p0 2 E and f 2 C.Œ0, T; E/, (17)–(18) has a unique mild solution p.t/ D etA p0 C
Z
t
e.ts/A f .s/ds.
0
If moreover, p0 2 D.A/ and f 2 W 1,1 .0, T; L2 .//, then p 2 W 1,1 .Œ0, T; E/ \ L2 .0, T; D.A// and (5)–(8) are satisfied. Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI
3. The regularizing problem Let " be positive and consider the approximating parabolic problem @p" @2 p" @p" @ C " 2 @t @a @a @y
@p" K.a, y/ D f in .0, T/ , @y
@p" @p" .t, a, 0/ D .t, a, L/ D 0 in .0, T/ .0, a /, @y @y
@p" p" " @a
Z
(20)
(21)
a
.t, 0, y/ D
ˇ.a, y/p" .t, a, y/da in .0, T/ .0, L/,
(22)
0
@p" .t, a , y/ D 0 in .0, T/ .0, L/, @a
(23)
p" .0, a, y/ D p0 .a, y/ in .
(24)
To embed this problem in the space E, we consider the linear operator A" : D .A" / E ! E defined as A" v D
@2 v @v @ " 2 @a @y @a
@v K.a, y/ @y
(25)
on the domain D .A" / defined as the set of all the function v 2 H2 ./, satisfying the boundary conditions .i/ .ii/ .iii/
@v @v .a, 0/ D .a, L/ D 0 for a 2 Œ0, a , @y @y Z a @v ˇ.a, y/v.a, y/da for y 2 Œ0, L, v.0, y/ " .0, y/ D @a 0
(26)
v.a , y/ D 0 for y 2 Œ0, L.
In this way, we are led to the following abstract formulation of the approximating problem (20)–(24) dp" .t/ C A" p" .t/ D f .t/ a.e. t 2 .0, T/, dt
(27)
p" .0/ D p0 .
(28)
In order to prove well-posedness of this problem, first we note that Proposition 1 The operator A" is quasi-accretive and D.A" / D E. Proof By ., / and k k, we respectively denote the scalar product and the norm in E. Then, from the definition of the operator A" , (25) and (26), we have 2 2 Z L @v C " @v C v 2 .a , y/dy .A" v, v/ K0 @y @a 0 Z Z @ 2 1 ˇ.a, y/v.a, y/v.0, y/dyda v .a, y/dyda 2 @a Z L Z L 1 1 (29) v 2 .a , y/dy C v 2 .0, y/dy 2 0 2 0
Z 1 1 ˇ 2 a v 2 .a, y/ C v 2 .0, y/ dyda 2 C a 1 2 ˇC a kvk2 , 2 which shows that A" is quasi-accretive. Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI To prove that D.A" / is dense in E, we consider f 2 E and take the sequences ffn g, fn g such that fn 2 C01 .0, a / .0, L/ , fn ! f , in E, n 2 C 1 Œ0, a ,
n .0/ D 1,
n .a / D n0 .0/ D n0 .a / D 0,
n ! 0 in L2 .0, a /. Then, the sequence gn .a, y/ D cn .y/n .a/ C fn .a, y/,
(30)
with R a
cn .y/ D
ˇ.a, y/fn .a, y/da , R a 1 0 ˇ.a, y/n .a/da 0
satisfies all conditions (26) (notice that we have to take n sufficiently large in order to have the denominator different from 0). In fact, such conditions can be easily checked, noticing that cn .0/ D cn .L/ D cn0 .0/ D cn0 .L/ D 0. Finally, from (30) we also have gn ! f , in E. We note that, from the proof of the previous Proposition, it is clear that there exists a set D (containing functions in the form (30)) dense in E and such that D C 1 ./ ,
D D.A/ \ D .A / .
(31)
Now, we prove that A" is actually m-accretive. Namely, for g 2 E, we prove that the equation v C A" v D g,
(32)
has a solution v 2 D.A" /, for > 0 sufficiently large. To this aim, we consider the problem @2 v @v @ @v " 2 K.a, y/ Dg .i/ v C @a @a @y @y .ii/ .iii/ .iv/
@v @v .a, 0/ D .a, L/ D 0 @y @y Z a @v ˇ.a, y/!.a, y/da, v.0, y/ " .0, y/ D @a 0 @v .a , y/ D 0 @a
in ,
in .0, a /, (33) in .0, L/, in .0, L/,
where ! 2 L2 ./ is fixed. For this problem, resorting to the general theory for uniformly elliptic operators [19], we are able to state the following result. Proposition 2 Let ! 2 E. Then problem (33) has a unique weak solution v 2 H1 ./. If ! 2 H1 ./ then v 2 H2 ./. Then we have Theorem 1 The operator A" is quasi-m-accretive. Proof By Proposition 2, we can define a mapping T : E ! E, defining T ! as the solution of problem (33). Now, taking !, ! 2 H1 ./ and, respectively, denoting v D T ! and v D T !, we have @ @ @2 @ K.a, y/ .v v/ D 0. .v v/ C .v v/ " 2 .v v/ @a @a @y @y Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI Then, multiplying by .v v/, integrating by parts using (33) and following a similar procedure as in (29), we have kv vk2
2 a ˇC
2
k! !k2 .
2 a sufficiently large, the mapping T : E ! E is a contraction in the space E. As a conseThis inequality shows that, for > 0 D 12 ˇC quence, we have a fixed point in E that, thanks to Proposition 2, belongs to H2 ./. This latter statement shows that (32) has a solution for each g 2 E, and A" is quasi m-accretive.
In conclusion, the previous Theorem states that the operator A" is the infinitesimal generator of a strongly continuous semigroup etA" such that tA" (34) e e0 t . Thus, for any p0 2 E and f 2 C.Œ0, T; E/, tA"
p" .t/ D e
Z p0 C
t
e.ts/A" f .s/ds
(35)
0
is the unique mild solution to problem (20)-(24) and if p0 2 D.A" /, f 2 W 1,1 .0, T; L2 .//, then p" .t/ is such that
p" 2 C Œ0, T; L2 ./ \ L1 .0, T; D.A" // \ W 1,1 .0, T; L2 .//, and satisfies (20)–(24). Next section is devoted to prove the convergence of this solution to the solution of (17)–(18).
4. Convergence Here, we consider the problem of convergence of the (mild) solution of (27)–(28) to the (mild) solution of (17)–(18). Actually, it is enough to prove the weak convergence of the semigroups lim etA" p0 D etA p0 weakly in L2 .0, T; E/ .
(36)
"!0
where p0 2 E. In fact, from (35) and the uniform bound (34), we easily derive Z t Z t .ts/A" e f .s/ds D e.ts/A f .s/ds, weakly in L2 .0, T; E/. lim !0 0
0
Then we focus on (36). Actually we have Proposition 3 Let p0 2 E, then (36) holds. Proof Because the set D is dense in E (31), and the uniform bounds (19), (34) hold, it is enough to prove (36) for p0 2 D. Thus, because p0 2 D D .A /, from the general theory of m-accretive operators, we have that p .t/ D etA p0 belongs to C 1 .Œ0, T; E/ \ C .Œ0, T; D.A // (note that here p .t/ is the solution of (27)–(28) with f .t/ 0) and the following estimate holds dp e0 T kA p0 k M, .t/ (37) dt where @p0 @2 p0 @p0 @2 p0 C C C kKk 1 M D e0 T C ./ @a @a2 @y @y2 is independent of . Moreover, because p .t/ satisfies (20)–(24) (with f .t/ 0), multiplying (20) by p .t/ and integrating over , we obtain Z t Z t @p 2 @p 2 ds C K0 ds C kp0 k2 , (38) .s/ .s/ kp" .t/k2 C " @a @y 0 0 where we have integrated by parts, using the boundary conditions as in the estimate (29), and C is a constant independent of . Thus, we deduce the uniform boundedness, with respect to ", of the following sequences: fp" g in L2 .0, T; V/ , Copyright © 2012 John Wiley & Sons, Ltd.
(39) Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI
dp" , dt
p @p" in L2 .0, T; E/ . " @a
(40)
Moreover, because from (20) we have dp" @2 p" @p" C KC kp" .t/k .t/ " 2 .t/, .t/ V kkV , dt @a @a V 0 ,V for any 2 V, using (37), we obtain @ @p p .t/ " .t/ " @a @a
V0
M C KC p" .t/V ,
@p @ p" .t/ " .t/ , in the space L2 0, T; V 0 . @a @a From the previous estimates, we can select a subsequence (denoted in the same way as the sequence p" / satisfying the following weak convergence statements as " ! 0
so that, by (39), we also prove boundedness of the sequence
.i/
p" ! p
weakly
in
L2 .0, T; V/
.ii/
dp" dt
weakly
in
L2 .0, T; E/
weakly
in
weakly
in
L2 .0, T; E/ L2 0, T; V 0 .
.iii/ .iv/
! p @p " @a !
@ p" " @p !
@a @a
(41)
By (41, i) and (41, ii), we obtain dp dp ! weakly in L2 .0, T; E/ , dt dt and, by (41, iii) we have that "
@p" !0 @a
weakly in L2 .0, T; E/,
(42)
so that p "
@p" !p @a
weakly in L2 .0, T; E/.
(43)
Thus, by (41, iv), @p @ @p" p " ! @a @a @a
weakly in L2 0, T; V 0 .
(44)
Consequently, because @p @p p .t, 0, y/ D p .t, a, y/ @a @a Z a @ @p p .t, , y/d , C @ 0 @ we also have p" .t, 0, y/ "
@p" .t, 0, y/ ! p.t, 0, y/ weakly in L2 .0, T; V 0 /, @a
and, by (41, i) Z
Z
a
a
ˇ.a, y/p.t, a, y/da weakly in L2 .0, T; V/.
ˇ.a, y/p" .t, a, y/da ! 0
Thus, we conclude that p 2 L2 .0, T; V/,
0 @p @a
2 L2 .0, T; V 0 / and, because p" satisfies (22), going to the limit we have Z p.t, 0, y/ D
a
ˇ.a, y/p.a, y/da, 0
and in conclusion p.t/ 2 D.A0 /. Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI Now, for any 2 L2 .0, T; V/, we have Z lim
T
!0 0
hA p .t/, .t/iV 0 ,V dt
@ @p" p .t/ " .t/ , .t/ dt !0 @a @a 0 V 0 ,V ! Z TZ @p" @ K.a, y/ C .t/ .t/dydadt @y @y 0 Z TZ Z T @p @p @ .t/, .t/ dt C K.a, y/ .t/ .t/dydadt D @a @y @y 0 0 0 V ,V Z T D hA0 p.t/, .t/iV 0 ,V dt. Z
T
D lim
0
Thus, dp dp A p ! A0 p D 0 weakly in L2 .0, T; V 0 /, dt dt and, because
dp dt
2 L2 .0, T; E/, then A0 p.t/ 2 L2 .0, T; E/ so that p.t/ 2 D.A/,
dp .t/ C Ap.t/ D 0 a.e. t 2 .0, T/ dt
and p.t/ D etA p0 . Thus, (36) holds.
In conclusion, after the previous Proposition, we can state the following result. Theorem 2 Let p0 2 E and f 2 C.Œ0, T; E/, then the respective solutions p , p" 2 C.Œ0, T; E/ of problems (5)–(8) and (20)–(24) satisfy weak lim p" D p , in L2 .0, T; E/. "!0
5. Numerical tests The convergence result proved in the previous sections suggests the use of problem (20)–(24) to provide an approximation to problem (5)–(8). Concerning (20)–(24), we can either directly use a computational software for parabolic problems (we actually used COMSOL Multiphysics for some tests that we show further on) or perform an iterative procedure to avoid the implicit, non-local, nature of condition (22). In fact, such a condition, though supported by COMSOL Multiphysics may not be available in other packages or numerical methods that, instead, may be limited to Robin type conditions. Actually, taking any ! 2 E, we may consider an iterative procedure to build the sequence fpn g, setting p0 D ! and, for n > 0, taking pn as the solution of the problem @pn @2 pn @pn" @ C " " 2" @t @a @a @y
K.a, y/
@pn" @y
D f in .0, T/ ,
@pn" @pn .t, a, 0/ D " .t, a, L/ D 0 in .0, T/ .0, a /, @y @y Z a @pn pn" " " .t, 0, y/ D ˇ.a, y/pn1 " .t, a, y/da in .0, T/ .0, L/, @a 0 Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI @pn" .t, a , y/ D 0 in .0, T/ .0, L/, @a pn" .0, a, y/ D p0 .a, y/ in . In fact, proceeding with the argument used to derive (38), we obtain for n > 1 Z t 2 2 n n1 2 p .t/ pn1 p .s/ pn2 .t/ ˇC a .s/ ds, 0
which guarantees the convergence of the sequence in C.Œ0, T; E/. Concerning our numerical computations by COMSOL Multiphysics, we focus on particular problems for which the analytical solution can be determined and run some tests in order to compare the computed approximation to (20)–(24) with this analytical solution. Namely, we consider (5)–(8), in the domain D .0, 1/ .0, /, with the following particular maternity function ˇ.a/ D
4e2 a e .1 a/, 1 C e2
(45)
while, as far as the diffusion coefficient K.a, y/, the function f .t, a, y/ and the initial datum p0 .a, y/ are concerned, we consider the no-diffusion case K.a, y/ 0,
f .t, a, y/ 0,
p0 .a, y/ ea ,
(46)
and also the case 8 < K.a, y/ D .1 C a/.1 C y/, f .t, a, y/ D e.ta/ .1 C a/.sin y C cos y C y cos y/, : p0 .a, y/ D ea cos y.
a
b
c
d
(47)
Figure 1. Relative errors plotted against time in the case without diffusion.
Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI The two cases (46) and (47), respectively, correspond to the exact solutions ptest .t, a/ D e.ta/ ,
ptest .t, a, y/ D e.ta/ cos y.
(48)
Thus, we compute the solution p" .t, a, y/ to (20)–(24) corresponding to the two cases, and we calculate the relative difference against the analytical solution, that is, the relative error err D .p" .t, x, y/ ptest .t, a, y//=ptest .t, a, y/,
(49)
where the respective ptest are given by (48). The results relative to case (46) are shown in Figure 1 where (49) is plotted against time, at some points .a, y/ (with values a D 0, 0.25, 0.5, 0.75, 1 and y D 1) for different values of (values " D 0.1, 0.01, 0.001, 0.0001), respectively, corresponding to Figure 1a–d. In Table I, where the maximum relative error is shown, we can observe the high decrease of the relative errors with the decreasing of ". We also mention that, as decreases from 0.1 to 0.0001, the maximum of the absolute error decreases from 2 to 0.003. The results relative to case (47) are shown in Figure 2 where the variation (49) is displayed at the same points and for the same values of " as before. In Table I, the maximum relative error is shown for the different values of ". Similar tests for other values lead to comparable results, so that we can conclude that the computation of the solution by COMSOL Multiphysics is carried out with a good approximation.
Table I. Maximum relative error from the tests in Figure 1 and Figure 2, for different values of ".
Max. err (without diffusion) (%) Max. err (with diffusion) (%)
" D 101
" D 102
" D 103
" D 104
40 12
6 2
0.7 0.5
0.06 0.2
a
b
c
d
Figure 2. Relative errors plotted against time in the case with diffusion.
Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
M. IANNELLI AND G. MARINOSCHI
References 1. Webb GF. Theory of Nonlinear Age-Dependent Population Dynamics. Marcel Dekker: New York, 1985. 2. Iannelli M. Mathematical Theory of Age-Structured Population Dynamics, Applied Mathematics Monographs (C.N.R.) Giardini Editori e Stampatori: Pisa, Italy, 1994. 3. Abia LM, Angulo O, Lopez-Marcos JC. Size-structured population dynamics models and their numerical solutions. Discrete and Continuous Dynamical Systems B 2004; 4:1203–1222. 4. Abia LM, Angulo O, Lopez-Marcos JC. Age-structured population dynamics models and their numerical solutions. Ecological Modelling 2005; 188:112–136. 5. Kim M-Y. Discontinuous Galerkin methods for a model of population dynamics with unbounded mortality. SIAM Journal of Scientific Computing 2006; 27:1371–1393. 6. Kim M-Y, Selenge T. Age-time continuous Galerkin methods for a model of population dynamics. Journal of Computational and Applied Mathematics 2009; 223:659–671. 7. Angulo O, López-Marcos JC, López-Marcos MA, Milner FA. A numerical method for nonlinear age-structured population models with finite maximum age. Journal of Mathematical Analysis and Applications 2010; 361:150–160. 8. Lopez L, Trigiante D. A finite-difference scheme for a stiff problem arising in the numerical solution of a population dynamical model with spatial diffusion. Nonlinear Analysis, Theory, Methods and Applications 1985; 9:1–12. 9. Milner FA. A numerical method for a model of population dynamics with spatial diffusion. Computers and Mathematics with Applications 1990; 19:31–43. 10. Ayati B. A variable time step method for an age-dependent population model with non-linear diffusion. Siam Journal of Numerical Analysis 2000; 37:1571–1589. 11. Ayati BP, Dupont T. Galerkin methods in age and space for a population model with nonlinear diffusion. SIAM Journal on Numerical Analysis 2002; 40:1064–1076. 12. Deng Q, Hallam TG. Numerical approximations for an age structured model of a population dispersing in a spatially heterogeneous environment. Mathematical Medicine and Biology 2004; 21:247–268. 13. Kim M-Y. Global dynamics of approximate solutions to an age-structured epidemic model with diffusion. Advances in Computational Mathematics 2006; 25:451–474. 14. Kim M-Y. Discontinuous-continuous Galerkin methods for a structured model of a biological system. SIAM Journal of Scientific Computing 2008; 31:913–938. 15. Pelovska G. An improved explicit scheme for age-dependent population models with spatial diffusion. International Journal of Numerical Analysis and Modeling 2008; 5:466–490. 16. Cusulin C, Gerardo Giorda L. A numerical method for spatial diffusion in age-structured populations. Numerical Methods for Partial Differential Equations 2010; 26:253–273. 17. COMSOL Multiphysics v3.5a. Floating Network License 1025226. Comsol Sweden, 2009. 18. Cusulin C, Iannelli M, Marinoschi G. Age-structured diffusion in a multi-layer environment. Nonlinear Analysis: Real World Applications 2005; 6:207–233. 19. Grisvard P. Elliptic Problems in Nonsmooth Domains, Monographs and Studies in Mathematics 14, Pitman: Boston, 1985.
Copyright © 2012 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2012
