Fish populations dynamics with nonlinear stock-recruitment renewal conditions

Applied Mathematics and Computation 277 (2016) 101–110

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Fish populations dynamics with nonlinear stock-recruitment renewal conditions Gabriela Marinoschi a,∗, Angela Martiradonna b a b

Institute of Mathematical Statistics and Applied Mathematics, Calea 13 Septembrie 13, Bucharest 050711, Romania Dipartimento di Matematica, Università degli Studi di Bari Aldo Moro, Via E. Orabona 4, Bari 70125, Italy

a r t i c l e

i n f o

a b s t r a c t

MSC: 35K60 35K90 35Q92 47J35

The dynamics of a fish population with age-structure and space diffusion is studied under a renewal condition represented by various nonlocal nonlinear stock-recruitment functions, instead of the standard linear birth condition. This population dynamics model is approached as a Cauchy problem for an evolution equation with an unbounded nonlinear operator in a Hilbert space. The domain of the operator contains specific restrictions induced by the definition of the stock-recruitment function which make not possible the proof of the m-accretiveness property. Its lack is compensated by some other essential properties proved in the paper, which allow the proof of the existence and uniqueness of the solution. The semigroup formulation of the problem ensures the convergence of a time-difference scheme used for providing some numerical simulations which can give information about the stock, recruitment and fishing strategy.

Keywords: Evolution equations Accretive operators Nonlinear population dynamics Age-structure Stock-recruitment function

© 2016 Elsevier Inc. All rights reserved.

1. Introduction Models of fish population dynamics with age-structure and distribution in space allow to check the validity of the methods of assessment and fisheries management. Simulations based on them can be used to describe the current state of marine resources, to make projections on the stock in the future and to evaluate the best fishing strategies (see e.g., [15,21]). To describe the dynamics of a population exploited by fishery, it is common to use the concepts of stock and recruitment (see e.g., [7,16,18]). The fish stock is defined as a subpopulation of a species of fish which lives in specific geographic area exploited to fishery. Each stock is independent on other stocks of the same species: emigration and immigration events are ignored. All individuals of the stock are part of the same evolution process: they born, grow, reproduce and die with the same intrinsic parameters (see [8]). The number of newborns becoming part of the stock is called recruitment. The relation between the parental stock and the resulting recruitment is known as stock-recruitment function. Starting from initial discrete models, as those of Ricker and Beverton-Holt, continuous empirical curves, giving indication about some particular aspects of recruitment formation, were introduced in the literature. Following [16], we refer our study to some of the next models

R(S ) = b1 S e−b2 S (Ricker), R(S ) =

S S (Shepherd), (Beverton–Holt), R(S ) = b1 + b2 S b 1 + b2 S b3

R(S ) = b1 Sb3 e−b2 S (Saila–Lorda), R(S ) = b1 Sb3 (Cushing), where b1 , b2 , b3 > 0. ∗

Corresponding author. Tel.: +40 213182433. E-mail addresses: [email protected], [email protected] (G. Marinoschi), [email protected] (A. Martiradonna).

http://dx.doi.org/10.1016/j.amc.2015.12.041 0096-3003/© 2016 Elsevier Inc. All rights reserved.

102

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

Here, S denotes the spawning stock and R(S) the recruitment from S. These parameters may be estimated from given data by using various fitting techniques and some other mathematical tools (see e.g., [14]). In order to use a function providing a relationship consistent with reality, certain constraints have to be imposed. The functions should give nonnegative values for R, should not return a positive value for R for S = 0, and should not give infinite values of R, except possibly for infinite S. The stock-recruitments functions satisfying these conditions are called admissible (see [16]). Consequently, all previous functions should be defined for S ≥ 0. The Saila–Lorda model is a generalization of Ricker’s model, the third parameter b3 allowing more variety in shape of curve than the Ricker’s one, even if it was often found that the best-fitting equation is inadmissible. The Cushing equation use is limited by the fact that the function is unbounded as S increases, as it happens with the Shepherd model for b3 < 1. Moreover, the Cushing and Saila–Lorda models with b3 < 1 exhibit an infinite slope at zero S, so that all these cases are avoided in our study since they are not relevant in practice. In case of difficulty of parameters estimation of the previous stock-recruitment functions, two further models are often used in fish stock assessment (see [4]). These are the hockey-stick model

R(S ) = α min(S, S∗ ),

α , S∗ > 0,

where recruitment is assumed proportional to spawner and it becomes constant when the latter arises the given threshold level S∗ and the quadratic hockey-stick model

⎧ ⎪  ⎨α S, R (S ) = α S − ⎪ ⎩ ∗ αS ,

(S−(1−δ )S ) 4δ S∗

∗ 2



0 ≤ S ≤ ( 1 − δ )S ∗ ,

( 1 − δ )S ∗ < S ≤ ( 1 + δ )S ∗ S > ( 1 + δ )S ∗ ,

with α , S∗ > 0 and δ ∈ (0, 1), which allow a smoother transition between the two parts of the hockey-stick. Previous studies in population dynamics taking into account a nonlinear birth condition with R of particular forms are provided e.g., in [5,6,11,12,19,20], where the interest was generally in the stability of the stationary solutions and in the analysis of the stability region. In [17] the existence and uniqueness of the solution is proved using a parabolic regularization via the age variable combined with the Leray–Schauder fixed point theorem. In this paper we consider an age-structured model for fish population dynamics, in a bounded 3D domain, with the McKendrick equation for the population density (see [1,13]), but with the birth equation described by a nonlinear stockrecruitment function. Due to the previous comments on the stock-recruitment functions we shall consider as basic for our study the models of Ricker and Shepherd (with b3 ≥ 1) , the latter becoming Beverton–Holt for b3 = 1. Our theory can also support Saila–Lorda model with b3 ≥ 1, and the two hockey-stick models. We suppose that all vital parameters, fertility, mortality and diffusion term depend both on age and space. The aim of this work is to prove the well-posedness of the model by using semigroup and variational techniques. The model is treated as a Cauchy problem for an evolution equation with an unbounded nonlinear operator A in a Hilbert space. We note that the domain of A contains specific restrictions induced by the definition of the stock-recruitment function and so it cannot be shown that A is quasi m-accretive (as in previous papers in the literature, see e.g., [9,10]), in order to derive by the general theory of evolution equations with accretive operators in Hilbert spaces, that the Cauchy problem for the evolution equation has a solution. The way followed to provide the existence result is the proof of some other essential properties of A, in Proposition 3.1 and 3.2, these ensuring the solution existence and uniqueness given in Theorem 3.3. Theorem 4.1 is an additional general result stating the solution existence and uniqueness under weaker assumptions on the initial datum and indicating a technique for numerical simulations which are provided at the end. 2. Functional statement of the model Let  be an open bounded subset of RN , N ≤ 3 with a C1 boundary  = ∂ . We denote by p(t, a, x) the density of the population depending on the time t ∈ (0, T), age a ∈ (0, a+ ) and the space position x ∈ , with T, a+ < +∞. The functions μ0 (a) and μ(a, x) stand for the natural mortality and the extra mortality respectively, the latter being essentially due to fishery in this context. The function m(a, x) describes the sexual maturity function, K(a, x) is the diffusion coefficient and f(t, a, x) is a possible source term. For completing the model parameters we consider two types of stock-recruitment functions, corresponding to the Ricker and Shepherd functions. We propose the following model describing the dynamics of the fish population

⎧ pt + pa + μ0 (a ) p + μ(a, x ) p − ∇ · (K (a, x )∇ p) = f (t, a, x ) ⎪ ⎪ ⎪ ⎨ p(0, a, x ) = p0 (a, x )  +  a ⎪ p(t, 0, x ) = R 0 p(t, a, x )m(a, x )da ⎪ ⎪ ⎩ K (a, x )∇ p · ν = 0

in (0, T ) × (0, a+ ) × , in (0, a+ ) × , in (0, T ) × ,

(2.1)

on (0, T ) × (0, a+ ) ×  ,

where R stands for an admissible stock-recruitment function. Here, ν is the unit outward normal to  , pt and pa are the partial derivatives with respect to t and a, ∇ = ( ∂∂x , . . . , ∂ ∂x ) is the gradient operator and ∇ · is the divergence operator. 1

N

The vanishing flux on the boundary indicates that there is no exchange of population with the exterior of the domain .

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

103

We observe that if R is linear one gets the standard renewal condition used in population dynamics (see e.g., [1,13]). In all results we shall provide further we assume that:



 μ0 ∈ L1loc (0, a+ ), μ0 (a ) ≥ 0 a.e. a ∈ 0, a+ ,

a+ 0

μ0 (a )da = +∞,

(2.2)

μ ∈ L∞ ((0, a+ ) × ), 0 ≤ μ(a, x ) ≤ μ∞ a.e. (a, x ) ∈ (0, a+ ) × ,

(2.3)

m ∈ L∞ ((0, a+ ) × ),

0 ≤ m(a, x ) ≤ m∞ a.e. (a, x ) ∈ (0, a+ ) × ,

(2.4)

K ∈ L∞ ((0, a+ ) × ),

0 < K0 ≤ K (a, x ) ≤ K∞ a.e. (a, x ) ∈ (0, a+ ) × ,

(2.5)

there exists RM , CR > 0 such that 0 ≤ R(S ) ≤ RM ,

(2.6)

|R(S ) − R(S )| ≤ CR |S − S|,

(2.7)

where RM and CR are the upper bound and Lipschitz constant, respectively. We note that (2.6) and (2.7) are obeyed by Ricker and Shepherd functions. The function of Beverton–Holt is a particularization of Shepherd with b3 = 1. Moreover, Saila-Lorda with b3 ≥ 1, hockey-stick and quadratic hockey-stick models verify assumptions (2.6) and (2.7), so that they can be included in this study. For example, CR is equal to b1 (1 + b2 ) for the Ricker model,

2+b3 b1

for the Shepherd model,

b1 b3 e1−b3 b −1

b23

((b3 − 1 )b3 −1 + ebb33 ) for Saila–Lorda model with b3 > 1, α for the hockey-stick

δ model and α 1+ δ for the quadratic hockey-stick model. a Condition (2.2) guarantees that the survival probability π0 (a ) = e− 0 μ0 (σ )dσ is essentially bounded in (0, a+ ) and vanishes + at maximum age a . We perform the following changes of functions

p0 (a, x ) = p0 (a, x )/π0 (a ), f (t, a, x ) = f (t, a, x )/π0 (a ), p(t, a, x ) = p(t, a, x )/π0 (a ), (a, x ) = K (a, x ), (a, x ) = m(a, x )π0 (a ), μ (a, x ) = μ(a, x ), K m and obtain the equivalent system for the normalized functions indicated before by the symbol ∼, but without writing it, for convenience. The system reads

⎧ p + pa + μ(a, x ) p − ∇ · (K (a, x )∇ p) = f (t, a, x ) ⎪ ⎪ t ⎪ ⎨ p(0, a, x ) = p0 (a, x )  +  a ⎪ t, 0, x = R p t, a, x m ( a, x ) da p ( ) ( ) ⎪ 0 ⎪ ⎩ K (a, x )∇ p · ν = 0

in (0, T ) × (0, a+ ) × , in (0, a+ ) × ,

(2.8)

in (0, T ) × , on (0, T ) × (0, a+ ) ×  .

We denote

H = L2 (), V = H 1 (), V = (H 1 ()) , where H1 () is the Sobolev space with the standard norm and V is its dual, and we recall that V ⊂ H ⊂ V with compact injections. We also use the notation

a+ = (0, a+ ) × , H = L2 (0, a+ ; H ), V = L2 (0, a+ ; V ), V = L2 (0, a+ ; V ). We define the operator A : D(A ) ⊂ H → H, with the domain



D (A ) =

u ∈ H; u ∈ V, ua ∈ V , Au ∈ H, u(a, x ) ≥ 0 a.e. (a, x ) ∈

by

(Au, ϕ )H =



a+ 0

ua (a ), ϕ (a ) V ,V da +

 a+



 , u(0, x ) = R

μ(a, x )u ϕ dx da +

a+

 a+

0

a+

 u(a, x )m(a, x )da

K (a, x )∇ u · ∇ϕ dx da

 a.e. x ∈  ,

(2.9)

for any u ∈ D(A), ϕ ∈ V. Here, (·, · )H and ·, · V ,V represent the scalar product in H and the pairing between V and V, respectively. For simplicity, sometimes we shall not indicate the function arguments in the integrals. Given this functional framework we introduce the following Cauchy problem

 dp dt

(t ) + Ap(t ) = f (t )

p( 0 ) = p0 .

a.e. t ∈ (0, T )

(2.10)

We aim at proving that problem (2.10) has a strong solution p ∈ C ([0, T ]; H ), meaning that p(t) ∈ D(A) a.e. t ∈ (0, T), p is absolutely continuous on any compact subinterval of (0, T) and it satisfies (2.10). We also observe that a strong solution p of (2.10) is a solution of (2.8) in the sense of distributions and so we can focus on the study of (2.10).

104

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

Variants of model (2.1) for R linear, but involving a strong nonlinear dependence of μ and m on the solution have been studied e.g., in [9,10], for stratified one-dimensional domains , on the basis of the quasi m-accretiveness of the operator defined for those models. In this case we cannot prove this property for the operator (2.9), due to the nonnegativeness requirement in the domain of A. 3. Main results First, we give a definition. If E is a set we define (see [2], p. 30) the convex hull of the set E by



convE =

k 

σi ui ; ui ∈ E, 0 ≤ σi ≤ 1,

i=1

k 



σi = 1

i=1

and we observe that convE = convE. Indeed, the inclusion convE ⊂ convE is obvious. Let us take η ∈ convE. Then, η = k k n n n i=1 σi ui ; ui ∈ E. For each ui we have ui = lim ui with ui ∈ E. It follows that ηn := i=1 σi ui ∈ convE and so η = lim ηn ∈ n→∞

n→∞

convE. The next two propositions refer to the properties of A and the main property which replaces the lack of quasi maccretiveness of A is given in Proposition 3.2. These allow the existence result in Theorem 3.3. Proposition 3.1. A is quasi-accretive and closed in H. Proof. We have to show (see [2], p. 98) that there exists ω ∈ (0, ∞) such that

(Au − Av, u − v )H ≥ −ωu − vH for any u, v ∈ D(A ). We introduce F : H → H by



(F u )(x ) = R

a+ 0

(3.1)



m(a, x )u(a, x )da , for u ≥ 0.

(3.2)

Then, (Fu)(x) ≥ 0 a.e. x ∈ , and by a simple calculus we get that F is Lipschitz with the Lipschitz constant M having different values corresponding to each admissible model,

F u − F uH ≤ Mu − uH , u, u ∈ H, u, u ≥ 0.

√ More exactly, M = CR a+ m∞ , where CR is the Lipschitz constant of the function R. Then,

(Au − Av, u − v )H =



a+ 0

(ua − va )(a ), (u − v )(a ) V ,V da +



a+

μ(u − v )2 dx da +

(3.3)

 a+

K |∇ (u − v )|2 dx da

 a+ 1 1 u(a+ ) − v(a+ )2H − u(0 ) − v(0 )2H + K0 ∇ (u(a ) − v(a ))2H da 2 2 0 1  1 M2 + K0 u − v2H + K0 u − v2V ≥ u(a+ ) − v(a+ )2H − 2 2 ≥

1

≥−

2

M2 + K0



u − v2H ,

(3.4) (3.5)

where we used the relation u(0 ) − v(0 )H = F (u ) − F (v )H ≤ M2 u − v2H . This proves that A is quasi-accretive, or ωaccretive, with the constant ω = 12 M2 + K0 . Hence, for λ > ω the operator λI + A is accretive on H, where I is the identity on H. Now, we prove that the operator λI + A is closed for some λ > λ0 . We consider the sequences {un }n ≥ 1 ⊂ D(A), {ηn }n≥1 ⊂ H and u, η ∈ H such that 2

un → u,

2

ηn = λun + Aun → η strongly in H, as n → ∞.

(3.6)

We have to prove that u ∈ D(A) and λu + Au = η. First, let us observe that {un }n and {un, a }n are Cauchy sequences, in V and V , respectively. Indeed, by multiplying the equation

λ(un − um ) + Aun − Aum = ηn − ηm

(3.7)

by un − um , we get

λun − um 2H + (Aun − Aum , un − um )H ≤

λ 2

ηn − ηm 2H +

1 un − um 2H , 2λ

whence we obtain by (3.4)



 1 − ω un − um 2H + K0 un − um 2V ≤ ηn − ηm 2H . 2 2λ

λ

(3.8)

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

105

As {ηn }n is convergent in H it follows that {un }n is a Cauchy sequence in V for λ > λ0 = 2ω. Thus, it converges strongly in V to u. To prove that {un, a }n is a Cauchy sequence in V we write by (3.7)

un,a (a ) − um,a (a ), ψ V ,V = (ηn (a ) − ηm (a ), ψ )H − λ(un (a ) − um (a ), ψ )H −(μ(a )(un (a ) − um (a )), ψ )H − (K (a )∇ (un (a ) − um (a )), ∇ψ )H and compute



un,a − um,a 2V = ≤3



a+

 sup



ψ ∈V,ψV ≤1

0

    2  un,a (a ) − um,a (a ), ψ V ,V  da

(3.9)

2 un − um 2V . ηn − ηm 2H + (λ + μ∞ )2 un − um 2H + K∞

Thus, {un, a }n turns out to be Cauchy in V and it follows that un, a → ua strongly in V , as n → ∞. This proves that ua ∈ V , hence u ∈ C ([0, a+ ]; H ). Moreover, by (3.6), we get that at limit, u(a, x) ≥ 0 a.e. (a, x ) ∈ a+ . By Arzelà–Ascoli theorem we deduce that

un (a ) → u(a ) strongly in V , uniformly in a ∈ [0, a+ ],  a+ while, by the boundary condition un (0, x ) = R( 0 un (a, x )m(a, x )da ) we see that {un (0)}n is bounded in H and so, un (0) → ζ weakly in H. By the limit uniqueness we get that ζ = u(0 ). Actually, we still get using (3.3) that un (0, · ) = F (un )(· ) → F (u )(· ) = u(0, · ) strongly in H. Moreover, by passing to the limit in (λun + Aun , ψ )H = (ηn , ψ )H for any ψ ∈ V, on the basis of all previous convergences, we get (λu + Au, ψ )H = (η, ψ )H for any ψ ∈ V. Thus, we obtain that η = λu + Au, so Au ∈ H. In conclusion we have proved that u ∈ D(A) and A is closed, as claimed.  Proposition 3.2. The operator A has the property



convD(A ) ⊂

R(λI + A ), f or some λ0 > 0.

(3.10)

λ>λ0

Here I is the identity operator and R(λI + A ) is the range of the operator λI + A. Proof. We prove (3.10) in three steps. We begin with the characterization of the closure of D(A) proving that

D(A ) = {u ∈ H; u(a, x ) ≥ 0 a.e. (a, x ) ∈ a+ }.

(3.11)

Then we prove that

D (A ) ⊂



R(λI + A ), for some

λ0 > 0,

(3.12)

λ>λ0

and finally we show (3.10). Step 1. Let us consider the problem

du + A0 (a )u = 0, u(0 ) = u0 , da

(3.13)



where A0 (a): V → V , A0 (a )u =  K (a, x )∇ u · ∇ψ dx for any ψ ∈ V. If u0 ∈ H, it is known that this problem has a solution which we denote u(a, x; u0 ) that belongs to L2 (0, a+ ; V ), its derivative with respect to a is in L2 (0, a+ ; V ) and the linear application that maps u0 to this solution, u0 → u(a, x; u0 ) has the property

u(a, ·; u0 ) − u(a, ·; v0 )H ≤ u0 − v0 H . Moreover, the solution u is nonnegative if u0 is nonnegative. Let f ∈ H, f(a, x) ≥ 0 a.e. (a, x ) ∈ a+ , and consider a sequence fn ∈ C ∞ ([0, a+ ) × ), fn (a, x) ≥ 0 a.e. (a, x ) ∈ a+ , such that fn → f strongly in H. We also consider a function ρ ∈ C ∞ (R+ ) such that ρ (r ) = 0 for 0 ≤ r ≤ 1 and ρ (r ) = 1 for r > 2 (see e.g., [3,10]). Let us define the mapping n : V → H by

(n g)(a, x ) = ρ (na ) fn (a, x ) + (1 − ρ (na ))u(a, x; F g),

(3.14)

where u(a, x; Fg) is the solution to (3.13) corresponding to the initial condition u0 = F g, with (F g)(x ) =  a+ R( 0 m(a, x )g(a, x )da ), for g ∈ H, g ≥ 0 a.e. on a+ , see (3.2). We notice that

n g ∈ V = L2 (0, a+ ; V ), (n g)a ∈ V = L2 (0, a+ ; V ), n g ≥ 0 a.e. in a+ . Moreover, considering another solution u(a, x; Fh) to (3.13), for which the initial condition is Fh, we have

n g − n h2H = (1 − ρ )(u(·, ·; F g) − u(·, ·; F h ))2H

(3.15)

106

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

 ≤

2/n



0

 ≤



2/n 0

 

[(1 − ρ (na ))(u(a, x; F g) − u(a, x; F h ))] dx da 2

2 n

(u(a, x; F g) − u(a, x; F h ))2 dx da ≤ F g − F h2H ≤

2M2 g − h2H , n

(3.16)

by (3.3). In the above calculations we took into account that for na > 2 we have 1 − ρ (na ) = 0 and for na ≤ 2 we have 1 − ρ (na ) ≤ 1, since 0 ≤ ρ (na) ≤ 1. Consequently, for n sufficiently large the integral on [2/n, a+ ] vanishes. In fact we obtained that n is a contraction on H and so the equation n g = g has a unique solution, which we denote further by gn (since it depends on n),

gn (a, x ) = ρ (na ) fn (a, x ) + (1 − ρ (na ))u(a, x; F gn ). It can be seen that gn (0, x ) = u(0, x; F gn ) = (F gn )(x ) = R( of n gn we see that gn ∈ D(A). Next, we calculate

 a+ 0

m(a, x )gn (a, x )da ). Taking into account the previous properties

gn − fn H = (1 − ρ ) fn + (1 − ρ )u(·, ·; F gn )H  2/n  1/2  2/n  1/2 2 ≤ [(1 − ρ (na )) fn ]2 dx da + 1 − ρ ( na )) u ( a, x; F g ) dx da [( n ]   0 0   ≤

2  f n H + n

2 Mgn H → 0 as n → ∞. n

Therefore fn = gn a.e. on H and we have lim fn = lim gn = f in H, proving thus the density (3.11). n→∞

n→∞

Step 2. For (3.12) we fix g ∈ D(A ) and show that for any λ > λ0 > 0 there exists u ∈ D(A) such that

λu + Au = g. To this aim, for λ

(3.17)

> 0 let us define the operator Bλ (a): V → V ,

Bλ (a )v, ψ V ,V :=





(λ + μ(a, x ))v(x )ψ (x )dx +





K (a, x )∇v(x ) · ∇ψ (x )dx,

for any v, ψ ∈ V and consider the Cauchy problem

 du

( a ) + Bλ ( a ) u ( a ) = g ( a ) , u(0 ) = F u.

a.e. a ∈ (0, a+ )

da

(3.18)

We assert that this problem has a unique strong solution u ∈ C ([0, a+ ]; H ) ∩ W 1,2 ([0, a+ ]; V ) ∩ V. To prove this we proceed by the Banach fixed point theorem. We fix z ∈ H in

 dw

( a ) + Bλ ( a ) w ( a ) = g ( a ) , w(0 ) = F z,

a.e. a ∈ (0, a+ )

da

(3.19)

define φ : H → H by φ (z ) = w and show that it is a contraction on H. Problem (3.19) has a unique strong solution w ∈ C ([0, a+ ]; H ) ∩ W 1,2 ([0, a+ ]; V ) ∩ V. This follows by applying Theorem 4.17 in [2], p. 177, because Bλ (a) is linear, continuous from V in V and coercive, for λ > K0 ,

Bλ (a )vV ≤ max ((λ + μ∞ ), K∞ )vV , Bλ (a )v, v V ,V =

 

λv2 dx +

 

μ(a, x )v2 dx +



K (a, x )|∇v| dx ≥ (λ − K0 )v2H + K0 vV2 , 2



for any v ∈ V . Moreover, since Bλ depends on a and it is continuous, the function a → Bλ (a )w(a ) is measurable from [0, a+ ] to V . Then, writing (3.19) for w and w corresponding to z and z, making the difference, multiplying by w(a ) − w(a ) and integrating with respect to a we get

1 w(a ) − w(a )2H + (λ − K0 ) 2



0

a

w(s ) − w(s )2H ds + K0



0

a

1 2

w(s ) − w(s )V2 ds ≤ F z − F z2H ≤

M2 z − z2H 2

by (3.3), hence, for λ > K0 , we get

w(a ) − w(a )2H ≤ M2



a 0

z(s ) − z(s )2H ds.

(3.20)

We show that this implies that φ is a contraction in C ([0, a+ ]; H ) by using instead of the standard norm on this space the equivalent norm wb = supa∈[0,a+ ] (e−δ a w(a )H ), with some δ > 0. We write by (3.20)

e−2δ a w(a ) − w(a )H ≤ M2 z − zb e−2δ a 2

2



a 0

e2δ s ds ≤

M2 z − z2b (1 − e−2δa ) 2δ

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

107

and so we obtain

M2 z − z2b . 2δ

w − w2b ≤

Choosing now δ large enough we get that φ (z ) = w is a contraction, so that z = w is the unique fixed point, for λ > K0 . In conclusion we can replace z by w in (3.19) and show that (3.18) has a unique solution, as claimed. Thus, for λ > 2ω > K0 it follows that (3.12) is satisfied with λ0 = 2ω. Since g ∈ D(A ) we have g(a, x) ≥ 0 a.e. (a, x ) ∈ a+ and we show that u(a, x) the solution to (3.17) is nonnegative a.e. in a+ . To this end we multiply (3.17) by the negative part u− , noticing that u(0 ) = F (u ) ≥ 0 so that u− (0 ) = 0. Integrating with respect to a we get by a few calculations, using Stampacchia Lemma

    2 1 − u (a )2H + (λ + μ )(u− )2 dx + K ∇ u−  dx = − gu− dx ≤ 0. 2   

Thus, it follows that u− (a ) = 0, implying that u(a, x) ≥ 0 for all a ∈ [0, a+ ] a.e. x ∈ . Step 3. Finally, (3.10) is immediately obtained if we note that the nonnegative cone D(A ) is convex and so every g ∈ convD(A ) belongs to D(A ) which satisfies (3.12).  Theorem 3.3. Assume that f ∈ W 1,1 (0, T ; H ) and p0 ∈ D(A). Then, problem (2.10) has a unique solution p ∈ C ([0, T ]; H ) ∩ W 1,∞ (0, T ; H ) ∩ L∞ (0, T ; D(A )), p(t, a, x) ≥ 0 for any t ∈ [0, T ], a.e. (a, x ) ∈ a+ , which satisfies the estimate

 p(t ) + 2 H



t 0

 p( τ )  d τ ≤ 2 V

1 Kmin

   T 2 2 2  p 0 H +  f (τ )H dτ e(M +2K0 +1 )t

(3.21)

0

for any t ∈ [0, T], where Kmin = min (1, 2K0 ) and M is the Lipschitz constant for F. If p and p are solutions to (2.10) corresponding to p0 , p0 ∈ D(A ) and f, f ∈ W 1,1 (0, T ; H ), then

   T 2  p(t ) − p(t )2H ≤  p0 − p0 2H +  f (τ ) − f (τ )2H dτ e(M +2K0 +1 )t , ∀t ∈[0, T ].

(3.22)

0

Proof. The existence of the solution to (2.10) follows by Theorem 4.7 in [2], p. 146. The solution nonnegativeness is ensured since p(t) ∈ D(A) for any t ∈ [0, T ]. To prove estimate (3.21), we multiply (2.10) by p(t) and integrate over (0, t), obtaining

 p(t )2H −  p0 2H + 2



t 0

(Ap(τ ), p(τ ))H dτ ≤

Since A0 = 0, by (3.4) we have

 p(t )2H +



t 0

≤  p0 2H +

 p(τ )(a+ )2H dτ + 2K0



T

0

t 0

0



 p(t ) ≤   + p0 2H



T

0

 f ( τ )



2 Hd

T

 f (τ )2H dτ +



t 0

 p(τ )2H dτ .

 p(τ )2V dτ

 f (τ )2H dτ + (M2 + 2K0 + 1 )

which by Gronwall’s Lemma implies 2 H





τ e (M

2

 0

t

 p(τ )2H dτ ,

(3.23)

) , for any t ∈ [0, T ].

+2K0 +1 t

Substituting in (3.23) we get (3.21). To prove estimate (3.22) we multiply

d ( p(t ) − p(t )) + Ap(t ) − Ap(t ) = f (t ) − f (t ) dt by p(t ) − p(t ) and integrate over (0, t), getting by (3.4)

 p(t ) − p(t )2H ≤  p0 − p0 2H − 2 ≤  p0 − p0 2H + (M2 + 2K0 + 1 )



t

0



0

(Ap(τ ) − Ap(τ ), p(τ ) − p(τ ))H dτ + 2

t

 p(τ ) − p(τ )2H dτ +



T 0

 0

t

( f (τ ) − f (τ ), p(τ ) − p(τ ))H dτ

 f (τ ) − f (τ )2H dτ

for each t ∈ [0, T]. Then, (3.22) follows from the Gronwall Lemma and the uniqueness of the solution of system (2.10) is an obvious consequence of this estimate.  4. Numerical results and discussion In particular, if f ≡ 0 we can apply the Crandall–Liggett Theorem (see [2], p. 131) and state Theorem 4.1. Assume that f ≡ 0 and p0 ∈ D(A ). Then, problem (2.10) has a unique mild solution p ∈ C ([0, T ]; H ) given by

p(t ) = lim

n→∞



I+

t A n

−n

p0

108

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

Fig. 1. Graphics of the Ricker function R1 (left) and Shepherd function S1 (right).

uniformly in t on compact intervals. This theorem allows the numerical computation of the solution from the time finite-difference scheme

pi+1 + hApi+1 = pi for i = 0, . . . , n − 1, p0 = p0 (a, x ), where h = t/n for any t ≥ 0 and pi (a, x ) = p(ih, a, x ). Here, we present some numerical simulations performed on the basis of an algorithm implemented in Comsol Multiphysics v. 3.5a (FLN License 1025226), in order to put into evidence the behavior induced by different types of stockrecruitment functions. To this end we make the system (2.1) dimensionless by considering the characteristic values L for length, a+ for age and time, 1/a+ for the vital rates μ0 , μ and m, Kc for the diffusion coefficient and pc for the density and introducing the dimensionless variables and functions indicated by ∗ :

x = x∗ L, t = t ∗ a+ , a = a∗ a+ , p = p∗ pc , K = K ∗ Kc , m = m∗ /a+ ,

μ0 = μ∗0 /a+ , μ = μ∗ /a+ .

Performing these transformations we get the dimensionless system which has the same form as (2.1) but with K replaced by KK ∗ where the dimensionless parameter K = Kc a+ /L2 . The dimensionless domain is (a∗ , x∗ ) ∈ (0, 1) × (0, 1). The simulations are done with a space dependent initial distribution. According to Zhang et al. (see [22]) the diffusion rate of fish individuals is related to the body size, so that K is taken age-dependent with the same growing behavior as in [22]. The sexual maturity m is assumed to increase with age up to a maximum value and then to level off. The data considered below refer to the dimensionless functions, but for simplicity they are written without the superscript ∗ ,

K = 1,

μ0 ( a ) =

K (a, x ) = 10−4

1 , a+ − a

 a α 0.1



μ(a, x ) = 0, f = 0, m(a, x ) = 20

5a, (a, x ) ∈ (0, 0.2 ) × (0, 1 ) 1, otherwise,

, α = 1.41, p0 (a, x ) = 10−3 y.

The code was run for two cases referring to R: – R1: Ricker case R(S ) = Se−S (b1 = b2 = 1 ); – S1: Shepherd case R(S ) = S/(1 + S ) (b1 = b2 = b3 = 1 ). The graphics of these functions are drawn in Fig. 1. The results showing the contour plots of the solution (not divided by π 0 (a)) to the dimensionless system (2.1) at some times (t = 0.1, 0.5, 1), in a life period, are plotted in Fig. 2 for the R1 model, and in Fig. 3 for the S1 model, respectively. Comparing the graphics in Fig. 1 with the corresponding ones in Fig. 2 we observe that in the Shepherd case the solution has higher values due to the fact that R(S) in this case levels off at its maximum. In both cases (and especially in the Ricker case) the very young population can be viewed in isolated patches, while the adult and old individuals are spread in a larger spatial range.  a+ The graphics of the stock S(t, x ) = 0 m(a, x ) p(t, a, x )da, and recruitment p(t, 0, x ) = R(S(t, x )) represented with respect to x, at times t = 0, 0.2, 0.4, 0.6, 0.8, 1 (indicated in the legend), are plotted in Fig. 4 for the R1 model, and in Fig. 5 for the S1 model.

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

109

Fig. 2. Contour plots for Ricker model R1, R(S ) = SeS , at times t = 0.1, t = 0.5, t = 1.

Fig. 3. Contour plots for Shepherd model S1, R(S ) = S/(1 + S ), at times t = 0.1, t = 0.5, t = 1.

Fig. 4. The functions S(t, x) and p(t, 0, x ) = R(S (t, x )) for Ricker model R1.

The same behavior for the solution can be observed in Figs. 4 and 5. The stock in the Shepherd case (Fig. 5) exhibits higher values than in Ricker case, but in both cases it increases in time and space. However, while the recruitment is always increasing in time and space in the Shepherd model, it increases in the Ricker case only at the beginning (times t = 0, t = 0.2) and has a decreasing shape in time and space in the second part of the life period, remaining at lower values than in the Shepherd model. In this paper we proved the well-posedness of a fish-population model with nonlinear admissible stock-recruitment functions. The novelty is in the proof of the existence result, based on some essential properties of the operator defined in the associated evolution equation. The semigroup formulation of the problem ensures the convergence of a time-difference scheme useful for numerical computations which may offer information about the evaluation of the stock and on the strategy of fishing activity.

110

G. Marinoschi, A. Martiradonna / Applied Mathematics and Computation 277 (2016) 101–110

Fig. 5. The functions S(t, x) and p(t, 0, x ) = R(S (t, x )) for Shepherd model S1.

Acknowledgments G.M. acknowledges the supportof the Grant CNCS–UEFISCDI, project number PN-II-ID-PCE-2011-3-0045. A.M. is very grateful to Maria Teresa Spedicato (COISPA Tecnologia & Ricerca, Via dei trulli 18-20, 70126 Bari, Italy) for guide and support to research. The authors would like to thank to the anonymous reviewer for the helpful suggestions related to the numerical simulations. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

S. Anita, ¸ Analysis and Control of Age-Dependent Population Dynamics, Kluwer Academic Publishers, Dordrecht, 2000. V. Barbu, Nonlinear Differential Equations of Monotone Type in Banach Spaces, Springer, New York, 2010. V. Barbu, M. Iannelli, The semigroup approach to non-linear age-structured equations, Rend. Istit. Mat. Univ. Trieste Suppl. XXVIII (1997) 59–71. N.J. Barrowman, R.A. Myers, Still more spawner-recruitment curves: the hockey stick and its generalizations, Can. J. Fish. Aquat. Sci. 57 (2000) 665–676. F. Bekkal-Brikci, K. Boushaba, O. Arino, Nonlinear age structured model with cannibalism, Discret. Contin. Dyn. Syst. Ser. B 7 (2007) 201–218. J. Bhattacharyya, S. Pal, The role of space in stage-structured cannibalism with harvesting of an adult predator, Comput. Math. Appl. 66 (2013) 339–355. E.L. Cadima, Fish stock assessment manual, FAO Fish. Tech. Pap. 393 (2003) 161. FAO, Rome. S.X. Cadrin, L.A. Kerr, S. Mariani, Stock Identification Methods: Applications in Fisheries Science, second ed., Elsevier Academic Press, 2014. C. Cusulin, M. Iannelli, G. Marinoschi, Age-structured diffusion in a multi-layer environment, Nonlinear Anal. Real World Appl. 6 (2005) 207–223. C. Cusulin, M. Iannelli, G. Marinoschi, Convergence in a multi-layer population model with age-structure, Nonlinear Anal. Real World Appl. 8 (2007) 887–902. G. DiBlasio, Asymptotic behavior of an age-structured fish population, Comput. Math. Appl. 9 (1983) 377–381. G.D. Blasio, M. Iannelli, E. Sinestrari, J.M. Biol., Approach to equilibrium in age structured populations with an increasing recruitment process, J. Math. Biol. 13 1982 371–382. M. Iannelli, A. Pugliese, An Introduction to Mathematical Population Dynamics, Springer, 2014. S. Ion, D. Marinescu, S.G. Cruceanu, V. Iordache, A data porting tool for coupling models with different discretization needs, Environ. Model. Softw. 62 (2014) 240–252. G. Lembo, A. Alvaro, F. Fiorentino, S. Martino, M.-T. Spedicato, ALADYM: An age and length-based single species simulator for exploring alternative management strategies, Aquat. Living Resour. 22 (2009) 233–241. T.C. Iles, A review of stock-recruitment relationships with reference to flatfish populations, Neth. J. Sea Res. 32 (1994) 399–420. O. Nakoulima, A. Omrane, J. Velin, A nonlinear problem for age-structured population dynamics with space diffusion, topological methods in nonlinear analysis, J. Juliusz Schauder Cent. 17 (2001) 307–319. T.J. Quinn, R.B. Deriso, Quantitative Fish Dynamics, Oxford University Press, 1999. Y. Saito, C. Jung, Y. Lim, Effects of cannibalism on a basic stage structure, Appl. Math. Comput. 217 (2010) 2133–2141. F.J. Solis, Self-limitation, fishing and cannibalism, Appl. Math. Comput. 135 (2003) 39–48. M.-T. Spedicato, J.C. Poulard, C.Y. Politou, K. Radtke, G. Lembo, P. Petitgas, Using the ALADYM simulation model for exploring the effect of management scenarios on fish population metrics, Aquat. Living Resour. 23 (2010) 153–165. L. Zhang, U.H. Thygesen, M. Banerjee, Size-dependent diffusion promotes the emergence of spatiotemporal patterns, Phys. Rev. E 90 (2014) 012904.