- Email: [email protected]
Mathematical Modeling Applied to Sustainable Management of Marine Resources
Available online at www.sciencedirect.com
ScienceDirect Procedia Engineering 194 (2017) 337 – 344
10th International Conference on Marine Technology, MARTEC 2016
Mathematical Modeling Applied to Sustainable Management of Marine Resources Md. Haider Ali Biswasa,∗, Md Rajib Hossaina , Mitun Kumar Mondala a Mathematics
Discipline, Science Engineering and Technology School, Khulna University, Khulna-9208, Bangladesh
Abstract In this work, we formulate and study a nonlinear mathematical model of fishery management to understand the dynamics of a fishery resource system in an aquatic environment that consists of two zones; one is free fishing zone and another is reserve zone where fishing is strictly prohibited. We have analyzed the model by finding the existence of equilibrium points: biological and bionomic, dynamical behavior of equilibrium points and also derived the conditions of stability and instability of the system. The behavior of this dynamical model of marine fishery management has been illustrated by the numerical simulations to establish the presented analytical results. c 2017 2017The The Authors. Published by Elsevier Ltd.is an open access article under the CC BY-NC-ND license © Authors. Published by Elsevier Ltd. This (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the organizing committee of the 10th International Conference on Marine Technology. Peer-review under responsibility of the organizing committee of the 10th International Conference on Marine Technology. Keywords: mathematical model; natural resources; marine resources; fishery models; non-linear differential equations; numerical analysis
1. Introduction In recent years, environmental issues including air and water pollution, climate change, overexploitation of marine ecosystems, exhaustion of fossil resources and conservation of biodiversity are receiving major attention from the public, stakeholders and scholars from the local to the planetary scales. It is now well recognized that human activities yield major ecological and environmental stresses with irreversible loss of species, destruction of habitat or climate catastrophes as the most dramatic examples of their effects. From this point of view, it is of great concern and environmental challenges today to preserve, conserve and manage the renewable natural resources in the marine and coastal areas for the sustainable development in Bangladesh. Sustainable development emphasizes the need to organize and control the dynamics and the complex interactions between man, production activities, and natural resources in order to promote their coexistence and their common evolution. Efficient and sustainable management of natural resources and the control measure of such systems mainly depends on essential understanding the mechanisms of their evolution over time where mathematical modeling in terms of nonlinear differential equations (ODEs) can play a significant role. As a result mathematical modeling is broadly used to discuss different types of real phenomena which lead to design better prediction, prevention, management and control strategies. For example, Biswas et al. [3] studied the potential impacts of Global Climate Change in Bangladesh. ∗
Corresponding author. Tel.: +88 01711948396; fax: +88 041 731244. E-mail address: [email protected]
1877-7058 © 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the organizing committee of the 10th International Conference on Marine Technology.
doi:10.1016/j.proeng.2017.08.154
338
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
Mondal and Biswas [15] developed a mathematical model to describe the transmission of Nipah virus between bats and human. Neilan and Lenhart [17] showed the application optimal control strategy in disease modelling. Biswas and Haque [5] discussed necessary of nonlinear dynamical system to control the infectious disease. We also refer readers to [1, 4, 6, 8, 9, 16] and the references within for knowing the more applications of mathematical modeling. During the last decades, many researchers presented their mathematical models as to optimizing fish production because of, mathematical modeling is the art of translating problems from an application area into traceable mathematical formulation with theoretical and numerical analysis. From of them, Leung and Wang [13] presented a mathematical model for commercial fishing, Clark et al. [7] studied the effects of irreversibility of capital investment upon optimal exploitation policies for renewable resources stocks, Kitabatake [11] developed a dynamic model for fishery resources with predator-prey relationship based on observational data for Lake Kasumigaura in Japan, Mesterton-Gibbons [14] also investigated an optimal policy to maximize the present value from the combined harvest of predator and prey, Zhang et al. [19] proposed and analyzed a model to study the optimal harvesting policy of a stage structured problem and derived necessary and sufficient condition for the coexistence and extinction of species. Recently, Song and Chen [18] discussed the optimal harvesting policy and stability for a two-species competitive system and derived conditions for the existence of a globally asymptotically stable positive equilibrium and a threshold of harvesting for the mature population, Dubey et al. [10] proposed a dynamic model for a single-species fishery which depends partially on a logistically growing resource. From the literature discussed above and to the best of our views, we have proposed a mathematical model to optimize the production of fish with the help of system of nonlinear differential equation. We have discussed the equilibrium points, dynamical behavior of those points and obtained the conditions of stability and instability of the proposed model. Finally, we have discussed the numerical simulations with graphical representation of the proposed mathematical model. 2. Model formulation We have considered a fishery resource consisting of two zones: a free fishing zone and a reserved area zone. We study a prey-predator system in a two patch environment, one accessible to both prey and predators and the other one being a refuge for the prey. Each patch is supposed to be homogenous. The prey refuge constitutes a reserve area of prey and no fishing is permitted in the reserved zone while the unreserved zone area is an open-access fishery zone. We suppose that the prey migrate between the two patches randomly. It is assumed that, in each zone growth of fish population follows logistic model. Keeping these in view, the model becomes dx x = rx 1 − (1) − σ1 x + σ2 y − qEx dt K dy Y = sy 1 − (2) + σ1 x − σ2 y dt L Here, x(t) be the biomass densities of the same population inside the unreserved area and y(t) be the biomass densities of the same population inside the reserved area, respectively at a time t. Let E be the total effort applied for harvesting the fish population in the unreserved area and the value of qE is called fishing mortality. Again r and s are the intrinsic growth rate of fish subpopulation inside the unreserved and reserved areas respectively. Let the fish subpopulation of the unreserved area migrate into reserved area at a rate σ1 and the fish subpopulation of the reserved area migrate into unreserved area at a rate σ2 . The crying capacities of fish species in unreserved and reserved area K and L respectively, q is the catch ability coefficient of fish species in the unreserved area. If there is no migration of fish population from reserved area to unreserved area (i.e.σ2 = 0) and r − σ1 − qE < 0, then ∂x ∂t < 0. Similarly, if there is no migration of fish population from unreserved area to reserved area (i.e.σ1 = 0) and s − σ2 < 0, then ∂y ∂t < 0. 3. Model analysis The nonlinear system in equations (1)-(2) has qualitatively analyzed in this section so as to find the equilibrium points and the dynamical behavior of that points as well as test the stability at positivity analysis equilibrium points.
339
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
3.1. Positivity analysis Now we will prove that all the variable in the system model equations (1) and (2) are positive. Lemma 1. If x(t) ≥ 0, y(t) ≥ 0 then the solution of x(t) and y(t) of the model system of equations (1) and (2) are positive. Proof: To prove this lemma, we get from the equation (1),
dx dt
= rx(1 − kx ) − σσ1 x + σ2 y − qex
dx e(σ1 +qE)tdt = e(σ1 +qE)t Now, dx dt + σ1 x + qEx ≥ 0 ⇒ dt + (σ1 + qE)x ≥ 0 and the integrating factor, IF = (σ1 +qE)t (σ1 +qE)t (σ1 +qE)t Multiplying by IF on both sides, we get, e e +e (σ1 + qE)x ≥ 0 ⇒ e(σ1 +qE)t ≥ c1 ⇒ x ≥ −(σ1 +qE)t , c1 e where, c1 is an integer constant. For the value of constant c1 we applying the initial conditions, when t = 0, then x(t) = x(0), hence x(0) ≥ c1 . Now putting the value of c1 in equation we get, x(t) ≥ x(0)e−(σ1 +qE)t . Since x(t) ≥ 0 and ∂y dy Y σ1 +qE ≥ 0, hence x(t) ≥ 0 for all t ≥ 0. Again from the equation (2), we get, ∂t = sy(1− L )+σ1 x−σ2 y ⇒ ∂t +σ2 y ≥ 0 and the integrating factor IF = eσ2 dt = eσ2 t . d σ2 t σ2 t σ2 t −σ2 t Now, multiplying eσ2 t on both sides, we have, eσ2 t dy , dt + e σ2 y ≥ 0 ⇒ dt e y ≥ 0 ⇒ e y ≥ c2 ⇒ y ≥ c2 e where c2 is an integer constant. For the value of c2 , we applying initial conditions, when t = 0, then y(t) = y(0), hence y(0) ≥ c2 . Now putting the value of c2 in equation we get, y(t) ≥ y(0)e−σ2 t , since y(t) ≥ 0 and σ2 ≥ 0. Hence x(t) ≥ 0 for all t ≥ 0. Therefore it is proved that x(t) and y(t) are positive for all t ≥ 0. 3.2. Existence of equilibrium Now we obtain the steady states of the system by equating the derivatives on the left hand sides to zero and solving the resulting algebraic equations. This gives two possible steady states, namely p(0, 0) and p(x∗ , y∗ ). Here (x∗ , y∗ ) are the positive solutions of the following algebraic equations x (3) rx 1 − − σ1 x + σ2 y − qEx = 0 K Y sy 1 − (4) + σ1 x − σ2 y = 0 L Eliminating y from equation (3) and (4), we get, x kσ1 x + qEkx − krx + rx2 σ2 y = σ1 x + qEx − rx 1 − ⇒y= k σ2 k
(5)
Putting the value y in equation (4) kσ1 x + qEkx − krx + rx2 kσ1 x + qEkx − krx + rx2 kσ1 x + qEkx − krx + rx2 s 1− x − σ + σ =0 1 2 σ2 k σ2 k σ2 k L s(r − σ1 − qE)2 (s − σ2 )r sr2 (s − σ2 ) 3 2 x + 2sr(r − σ1 − qE)x − − (r − σ1 − qE) − σ1 = 0 x+ ⇒ Lσ2 σ2 Lσ1 2 k2 Lσ2 2 ⇒ ax3 + bx2 + cx + d = 0 where a =
sr2 s(r − σ1 − qE)2 (s − σ2 )r (s − σ2 ) , b = 2sr(r − σ − qE), c = − and d = (r − σ1 − qE) − σ1 1 Lσ2 σ2 Lσ1 2 k2 Lσ2 2
The above equation has a unique positive solution x = x∗ if the following inequalities hold s(r − σ1 − qE)2 (s − σ2 )r , (s − σ1 )(r − σ1 − qE) < σ1 σ2 < Lσ2 k Substituting x∗ , y∗ from (3) and to be positive, we must have, x∗ > (k/r)(r − σ1 − qE)
(6)
340
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
3.3. Dynamical behavior of equilibria The dynamical behavior of equilibria can be studied by computing vibrational matrices corresponding to each equilibrium. In equation (1), we note that the trivial equilibrium p0 is unstable. Using the Routh-Hurwitz criteria, it is easy to check that all eigenvalues of the vibrational matrix corresponding to p∗ have negative real parts, and hence p∗ is locally asymptotically stable in the plane. This implies that we can find a small circle with center p∗ such that any solution (x(t), y(t)) of system (2), which is inside the circle at some time t = t1 will remain inside the circle for all t ≥ t1 and will tend to (x∗ , y∗ ) as t → ∞. In the following lemma we show that all solutions of model (2) are positive and uniformly bounded. Lemma 2. The set φ = {(x, y)R2 2 : w = x + y ≤ ϑμ } is a region of attraction for all solutions initiating in the K L interior of the positive quadrant, where ϑ is a positive constant and μ = 4r (r − ϑ − eq)2 + 4s (s + ϑ)2 . Proof: Let, w(t) = x(t) + y(t) and ϑ > 0 be a constant. Then we have, dw K rx2 sy2 L + (r − ϑ − eq)2 x2 − + (s + ϑ)y ≤ (r − ϑ − eq)2 + (s + ϑ)2 = μ =− dt K L 4r 4s By the theory of differential inequality [2], we have, 0 < w(t), x(t), y(t) ≤ ϑμ (1 − e−ϑt ) + w(0), x(0), y(0), e−ϑt , when t → ∞ and 0 < w < ϑμ . In the following theorems we will discuss global stability of equilibria of p∗ . Theorem 1. The equilibrium point p∗ is globally asymptotically stable. Proof: Let us consider the follwing definite function about p∗ x y V = x − x∗ ln ∗ + K1 y − y∗ − y∗ ln ∗ x y ∗
y∗ σ2 s r ∗ 2 ∗ 2 2 − xσ∗ xy xy − yx∗ < 0. This shows Differentiating V with respect to time t we get, dV dt = − K (x − x ) − x∗ σ1 L y − y that dV dt is negative and hence by Liapunov’s [12], it follows that the positive equilibrium is globally asymptotically stable. 3.4. Stability of steady states Let the two function u = u(x, y) = dx dt and v = v(x, y) = and (2). Then the Jacobian of the system isJ(x,y)
⎧ ∂u ⎪ ⎪ ⎨ ∂x =⎪ ∂v ⎪ ⎩ ∂x
⎫
∂u ⎪ ⎪ ∂y ⎬ ∂v ⎪ ⎪ ⎭ ∂y
=
r−
dy dt
2rx k
from the steady state (x∗ , y∗ ) = (0, 0) of the equation (1) − σ1 − qE σ2 σ1 s − 2y L − σ2
For the steady states (x, y) = (0, 0), we get, J(0,0) =
r − σ1 − qE σ2 σ1 s − σ2
The characteristics equation of the matrix with eigenvalue ω is, r − σ1 − qE − ω 0 |J − ωI| = 0 s − σ2 − ω Now, r − σ1 − qE − ω = 0 or, s − σ2 − ω = 0 ∴ ω1 = r − σ1 − qE, ω2 = s − σ2 Here the eigenvalue ω of the matrix “J” determine the stability of the states. Depending on ω following stability conditions: 1. If the eigenvalue ω < 0, then the steady state is stable, 2. If the eigenvalue ω > 0, then the system is unstable.
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
4. Numerical results and discussion In this section, we have numerically solved our proposed models and discussed their stability and instability. In Sections, we have discussed the numerical simulation of model (1), (2) respectively. We have already discussed existence of equilibrium, positivity analysis and stability of steady states. The model is simulated using the parameter values the total effort applied for harvesting the fish population in the reserve area E = 0.2. The carrying capacities of fish species in unreserved and reserved area respectively K = 110 and L = 200. The intrinsic growth rate of fish subpopulation are respectively r = 0.03 and s = 0.04. The catch ability fish coefficient q = 0.001. The fish subpopulation of unreserved and reserved area are respectively σ1 = 0.05 and σ1 = 0.05 and qE denote the fishing mortality rate of the model (1). Then the result is displayed in the Fig. 1. The displayed graph depicts the biomass densities of reserved area unreserved area where the population increased and decreased within 90 days. First of all the Fig. 1 has shown the population of unreserved area and reserved area are same on the 9 days. In 90 days the population of biomass densities of unreserved area increases above 100. At the same time the population of biomass densities of reserved area increases above 100. Again, Fig. 2 the population of reserved area increases above 900. At the same time the population of unreserved area increases above 400.and unreserved area are same 10 days. If we input the value of the total effort applied for harvesting the fish population and fish catch ability coefficient increases then the population decreases. If the population of the total effort applied for harvesting the fish population and fish catch ability coefficient decreases so the population decreases. Similarly by changing the values of parameter then the graph is changed. The Fig. 3 and Fig. 4 are changing the parameter values then the graph is changed. If we are changing the parameter values of intrinsic growth rate of unreserved and reserved area respectively and are increasing then the population are increasing.
Fig. 1. The graph of fishery model of the population with time (90 Days)
Fig. 2. Variation of population with time (90 Days) for E = 0.002, K = 1000, L = 200, r = 0.03, s = 0.04 and q = 0.001
341
342
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
Fig. 3. Variation of population with time (90 Days) for E = 0.002, K = 1000, L = 200, r = 0.299, s = 0.004 and q = 0.00001
Fig. 4. Variation of population with time (90 Days) for E = 0.002, K = 1000, L = 200, r = 0.03, s = 0.4 and q = 0.00001
Fig. 5. Variation of population with time (90 Days) for E = 4, K = 1000, L = 200, r = 0.03, s = 0.004 and q = 0.1
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
Fig. 6. Variation of population with time (90 Days) for E = 2.8, K = 1000, L = 200, r = 0.3, s = 0.006 and q = 0.2
If we are using the parameter values of carrying capacities of fish spices in the unreserved and reserved area fish production respectively K and L area increasing then the graph is increasing and by changing the parameter values K and L are increasing then the fish production is increasing. Similarly, the Fig. 5 and Fig. 6 below the X axis. So the population decreases. The maximum fish production is depend on crying capacities intrinsic growth rate of the unreserved and reserved area and harvesting rate or fishing mortality. In this Figure we observed that, the harvesting rate of fishing mortality is increasing then the graph is decreasing. Again the fishing mortality is decreasing then the graph is increasing. In this figure we observed that, by changing the values of parameter the characteristics figure is changed. So the model is fishery management model it can be find out of minimum and maximum fish production of reserved and unreserved area. 5. Conclusion In this mathematical modeling the system has been assumed that the two ecosystems; one free fishing zone and the other reserved zone. In this paper, we analysed the fish production by using many mathematical models. It has been found that the total user’s cost of harvest per unit of effort equals to the discounted value of the future marginal profit of the effort at the steady-state level. It has also been noted that if the discounted rate increases, then the economic rent decreases and even may tend to zero if the discounted rate tend to infinity. The possibility that the predator species in this model be itself harvested for commercial purposes, or that its non-use value be an important element in the overall economic value of the marine ecosystem have not been considered in this paper. Their inclusion in the model would surely lead to characterize another set of trade-offs which are bound to be at the centre of marine reserve policies. We have derived existence of equilibrium and globally asymptotically stable and unstable. Finally, we have been numerical simulations of fish production with graphical representation with different values of parameter. References [1] L. Berezansky, L. Idels, M. Kipnis, Mathematical model of marine protected areas, IMA J. Appl. Math. 76 (2011) 312-325. [2] G. Birkhoff, G.S. Rota, Ordinary Differential Equations, Ginn, Boston, 1882. [3] M.H.A. Biswas, T. Rahman, N. Haque, Modeling the Potential Impacts of Global Climate Change in Bangladesh: An Optimal Control Approach, J. Fundam. Appl. Sci. 8 (2016) 1-19. [4] M.H.A. Biswas, On the Evolution of AIDS/HIV Treatment: An Optimal Control Approach, Curr. HIV Res. 12 (2014) 1-12. [5] M.H.A. Biswas, M.M. Haque, Nonlinear Dynamical Systems in Modeling and Control of Infectious Diseases, Springer Pro. Math. & Stat. Vol. 164 (2016). [6] M.H.A. Biswas, L.T. Paiva, M.D.R. de Pinho, A SEIR model for control of infectious diseases with constraints, Math. Bios. and Engi. 11 (2014) 761-784. [7] C.W. Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resources, Wiley, New York, 1976. [8] T. Das, R.N. Mukherjee, K.S. Chaudhuri, Bioeconomic harvesting of a prey-predator fishery, J. Bio. Dyn. 3 (2009) 447-462. [9] M. De Lara, L. Doyen, Sustainable Management of Natural Resources: Mathematical Models and Methods, Springer, Verlag Berlin, 2008. [10] B. Dubey , P. Chandra, P. Sinha, A Model for Fishery Resource with Reserve Area, Nonl. Ana.: Real Worlld Appl. 4 (2003) 625-637.
343
344 [11] [12] [13] [14] [15] [16] [17] [18] [19]
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344 Y. Kitabatake, A Dynamic Predator-Prey Model for Fishery Resources: a Case of Lake Kasumigaura, Environ, Plan. A. 14 (1982) 225-235. J. La Salle, S. Lefschetz, Stability by Liapunovs Direct Method with Applications, Academic Press, New York, London, 1961. A. Leung, A. Wang, Analysis of models for commercial fishing: mathematical and economical aspects, Econometrica. 44 (1976) 295-303. M. Mesterton-Gibbons, On the Optimal Policy for Combined Harvesting of Predator Prey, Nat. Res. Model. 3 (1996) 235-244. M.K. Mondal, M.H.A. Biswas, Modeling the Transmission Dynamics of the Nipah Virus Infection, Proc. of 1st Int. Con. on Math. and Its Appl., Mathematics Discipline, Khulna University, Khulna, Bangladesh, 23rd December 2015, pp. 22-26. J.D. Murray, Mathematical Biology: Biomathematics, Second Edition, Springer, New York, 1989. R.M. Neilan, S. Lenhart, An Introduction to Optimal Control with an Application in Disease Modeling, DIMACS Series in Disc. Math. 75 (2010) 67-81. X. Song, L. Chen, Optimal Harvesting and Stability for a Two-species Competitive System with Stage Structure, Math. Biosci. 170 (2001) 173-186. X. Zhang, L. Chen, A.U. Neumann, The Stage Structured Predator-prey Model and Optimal Harvesting Policy, Math. Biosci. 168 (2000) 201-210.
ScienceDirect Procedia Engineering 194 (2017) 337 – 344
10th International Conference on Marine Technology, MARTEC 2016
Mathematical Modeling Applied to Sustainable Management of Marine Resources Md. Haider Ali Biswasa,∗, Md Rajib Hossaina , Mitun Kumar Mondala a Mathematics
Discipline, Science Engineering and Technology School, Khulna University, Khulna-9208, Bangladesh
Abstract In this work, we formulate and study a nonlinear mathematical model of fishery management to understand the dynamics of a fishery resource system in an aquatic environment that consists of two zones; one is free fishing zone and another is reserve zone where fishing is strictly prohibited. We have analyzed the model by finding the existence of equilibrium points: biological and bionomic, dynamical behavior of equilibrium points and also derived the conditions of stability and instability of the system. The behavior of this dynamical model of marine fishery management has been illustrated by the numerical simulations to establish the presented analytical results. c 2017 2017The The Authors. Published by Elsevier Ltd.is an open access article under the CC BY-NC-ND license © Authors. Published by Elsevier Ltd. This (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the organizing committee of the 10th International Conference on Marine Technology. Peer-review under responsibility of the organizing committee of the 10th International Conference on Marine Technology. Keywords: mathematical model; natural resources; marine resources; fishery models; non-linear differential equations; numerical analysis
1. Introduction In recent years, environmental issues including air and water pollution, climate change, overexploitation of marine ecosystems, exhaustion of fossil resources and conservation of biodiversity are receiving major attention from the public, stakeholders and scholars from the local to the planetary scales. It is now well recognized that human activities yield major ecological and environmental stresses with irreversible loss of species, destruction of habitat or climate catastrophes as the most dramatic examples of their effects. From this point of view, it is of great concern and environmental challenges today to preserve, conserve and manage the renewable natural resources in the marine and coastal areas for the sustainable development in Bangladesh. Sustainable development emphasizes the need to organize and control the dynamics and the complex interactions between man, production activities, and natural resources in order to promote their coexistence and their common evolution. Efficient and sustainable management of natural resources and the control measure of such systems mainly depends on essential understanding the mechanisms of their evolution over time where mathematical modeling in terms of nonlinear differential equations (ODEs) can play a significant role. As a result mathematical modeling is broadly used to discuss different types of real phenomena which lead to design better prediction, prevention, management and control strategies. For example, Biswas et al. [3] studied the potential impacts of Global Climate Change in Bangladesh. ∗
Corresponding author. Tel.: +88 01711948396; fax: +88 041 731244. E-mail address: [email protected]
1877-7058 © 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the organizing committee of the 10th International Conference on Marine Technology.
doi:10.1016/j.proeng.2017.08.154
338
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
Mondal and Biswas [15] developed a mathematical model to describe the transmission of Nipah virus between bats and human. Neilan and Lenhart [17] showed the application optimal control strategy in disease modelling. Biswas and Haque [5] discussed necessary of nonlinear dynamical system to control the infectious disease. We also refer readers to [1, 4, 6, 8, 9, 16] and the references within for knowing the more applications of mathematical modeling. During the last decades, many researchers presented their mathematical models as to optimizing fish production because of, mathematical modeling is the art of translating problems from an application area into traceable mathematical formulation with theoretical and numerical analysis. From of them, Leung and Wang [13] presented a mathematical model for commercial fishing, Clark et al. [7] studied the effects of irreversibility of capital investment upon optimal exploitation policies for renewable resources stocks, Kitabatake [11] developed a dynamic model for fishery resources with predator-prey relationship based on observational data for Lake Kasumigaura in Japan, Mesterton-Gibbons [14] also investigated an optimal policy to maximize the present value from the combined harvest of predator and prey, Zhang et al. [19] proposed and analyzed a model to study the optimal harvesting policy of a stage structured problem and derived necessary and sufficient condition for the coexistence and extinction of species. Recently, Song and Chen [18] discussed the optimal harvesting policy and stability for a two-species competitive system and derived conditions for the existence of a globally asymptotically stable positive equilibrium and a threshold of harvesting for the mature population, Dubey et al. [10] proposed a dynamic model for a single-species fishery which depends partially on a logistically growing resource. From the literature discussed above and to the best of our views, we have proposed a mathematical model to optimize the production of fish with the help of system of nonlinear differential equation. We have discussed the equilibrium points, dynamical behavior of those points and obtained the conditions of stability and instability of the proposed model. Finally, we have discussed the numerical simulations with graphical representation of the proposed mathematical model. 2. Model formulation We have considered a fishery resource consisting of two zones: a free fishing zone and a reserved area zone. We study a prey-predator system in a two patch environment, one accessible to both prey and predators and the other one being a refuge for the prey. Each patch is supposed to be homogenous. The prey refuge constitutes a reserve area of prey and no fishing is permitted in the reserved zone while the unreserved zone area is an open-access fishery zone. We suppose that the prey migrate between the two patches randomly. It is assumed that, in each zone growth of fish population follows logistic model. Keeping these in view, the model becomes dx x = rx 1 − (1) − σ1 x + σ2 y − qEx dt K dy Y = sy 1 − (2) + σ1 x − σ2 y dt L Here, x(t) be the biomass densities of the same population inside the unreserved area and y(t) be the biomass densities of the same population inside the reserved area, respectively at a time t. Let E be the total effort applied for harvesting the fish population in the unreserved area and the value of qE is called fishing mortality. Again r and s are the intrinsic growth rate of fish subpopulation inside the unreserved and reserved areas respectively. Let the fish subpopulation of the unreserved area migrate into reserved area at a rate σ1 and the fish subpopulation of the reserved area migrate into unreserved area at a rate σ2 . The crying capacities of fish species in unreserved and reserved area K and L respectively, q is the catch ability coefficient of fish species in the unreserved area. If there is no migration of fish population from reserved area to unreserved area (i.e.σ2 = 0) and r − σ1 − qE < 0, then ∂x ∂t < 0. Similarly, if there is no migration of fish population from unreserved area to reserved area (i.e.σ1 = 0) and s − σ2 < 0, then ∂y ∂t < 0. 3. Model analysis The nonlinear system in equations (1)-(2) has qualitatively analyzed in this section so as to find the equilibrium points and the dynamical behavior of that points as well as test the stability at positivity analysis equilibrium points.
339
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
3.1. Positivity analysis Now we will prove that all the variable in the system model equations (1) and (2) are positive. Lemma 1. If x(t) ≥ 0, y(t) ≥ 0 then the solution of x(t) and y(t) of the model system of equations (1) and (2) are positive. Proof: To prove this lemma, we get from the equation (1),
dx dt
= rx(1 − kx ) − σσ1 x + σ2 y − qex
dx e(σ1 +qE)tdt = e(σ1 +qE)t Now, dx dt + σ1 x + qEx ≥ 0 ⇒ dt + (σ1 + qE)x ≥ 0 and the integrating factor, IF = (σ1 +qE)t (σ1 +qE)t (σ1 +qE)t Multiplying by IF on both sides, we get, e e +e (σ1 + qE)x ≥ 0 ⇒ e(σ1 +qE)t ≥ c1 ⇒ x ≥ −(σ1 +qE)t , c1 e where, c1 is an integer constant. For the value of constant c1 we applying the initial conditions, when t = 0, then x(t) = x(0), hence x(0) ≥ c1 . Now putting the value of c1 in equation we get, x(t) ≥ x(0)e−(σ1 +qE)t . Since x(t) ≥ 0 and ∂y dy Y σ1 +qE ≥ 0, hence x(t) ≥ 0 for all t ≥ 0. Again from the equation (2), we get, ∂t = sy(1− L )+σ1 x−σ2 y ⇒ ∂t +σ2 y ≥ 0 and the integrating factor IF = eσ2 dt = eσ2 t . d σ2 t σ2 t σ2 t −σ2 t Now, multiplying eσ2 t on both sides, we have, eσ2 t dy , dt + e σ2 y ≥ 0 ⇒ dt e y ≥ 0 ⇒ e y ≥ c2 ⇒ y ≥ c2 e where c2 is an integer constant. For the value of c2 , we applying initial conditions, when t = 0, then y(t) = y(0), hence y(0) ≥ c2 . Now putting the value of c2 in equation we get, y(t) ≥ y(0)e−σ2 t , since y(t) ≥ 0 and σ2 ≥ 0. Hence x(t) ≥ 0 for all t ≥ 0. Therefore it is proved that x(t) and y(t) are positive for all t ≥ 0. 3.2. Existence of equilibrium Now we obtain the steady states of the system by equating the derivatives on the left hand sides to zero and solving the resulting algebraic equations. This gives two possible steady states, namely p(0, 0) and p(x∗ , y∗ ). Here (x∗ , y∗ ) are the positive solutions of the following algebraic equations x (3) rx 1 − − σ1 x + σ2 y − qEx = 0 K Y sy 1 − (4) + σ1 x − σ2 y = 0 L Eliminating y from equation (3) and (4), we get, x kσ1 x + qEkx − krx + rx2 σ2 y = σ1 x + qEx − rx 1 − ⇒y= k σ2 k
(5)
Putting the value y in equation (4) kσ1 x + qEkx − krx + rx2 kσ1 x + qEkx − krx + rx2 kσ1 x + qEkx − krx + rx2 s 1− x − σ + σ =0 1 2 σ2 k σ2 k σ2 k L s(r − σ1 − qE)2 (s − σ2 )r sr2 (s − σ2 ) 3 2 x + 2sr(r − σ1 − qE)x − − (r − σ1 − qE) − σ1 = 0 x+ ⇒ Lσ2 σ2 Lσ1 2 k2 Lσ2 2 ⇒ ax3 + bx2 + cx + d = 0 where a =
sr2 s(r − σ1 − qE)2 (s − σ2 )r (s − σ2 ) , b = 2sr(r − σ − qE), c = − and d = (r − σ1 − qE) − σ1 1 Lσ2 σ2 Lσ1 2 k2 Lσ2 2
The above equation has a unique positive solution x = x∗ if the following inequalities hold s(r − σ1 − qE)2 (s − σ2 )r , (s − σ1 )(r − σ1 − qE) < σ1 σ2 < Lσ2 k Substituting x∗ , y∗ from (3) and to be positive, we must have, x∗ > (k/r)(r − σ1 − qE)
(6)
340
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
3.3. Dynamical behavior of equilibria The dynamical behavior of equilibria can be studied by computing vibrational matrices corresponding to each equilibrium. In equation (1), we note that the trivial equilibrium p0 is unstable. Using the Routh-Hurwitz criteria, it is easy to check that all eigenvalues of the vibrational matrix corresponding to p∗ have negative real parts, and hence p∗ is locally asymptotically stable in the plane. This implies that we can find a small circle with center p∗ such that any solution (x(t), y(t)) of system (2), which is inside the circle at some time t = t1 will remain inside the circle for all t ≥ t1 and will tend to (x∗ , y∗ ) as t → ∞. In the following lemma we show that all solutions of model (2) are positive and uniformly bounded. Lemma 2. The set φ = {(x, y)R2 2 : w = x + y ≤ ϑμ } is a region of attraction for all solutions initiating in the K L interior of the positive quadrant, where ϑ is a positive constant and μ = 4r (r − ϑ − eq)2 + 4s (s + ϑ)2 . Proof: Let, w(t) = x(t) + y(t) and ϑ > 0 be a constant. Then we have, dw K rx2 sy2 L + (r − ϑ − eq)2 x2 − + (s + ϑ)y ≤ (r − ϑ − eq)2 + (s + ϑ)2 = μ =− dt K L 4r 4s By the theory of differential inequality [2], we have, 0 < w(t), x(t), y(t) ≤ ϑμ (1 − e−ϑt ) + w(0), x(0), y(0), e−ϑt , when t → ∞ and 0 < w < ϑμ . In the following theorems we will discuss global stability of equilibria of p∗ . Theorem 1. The equilibrium point p∗ is globally asymptotically stable. Proof: Let us consider the follwing definite function about p∗ x y V = x − x∗ ln ∗ + K1 y − y∗ − y∗ ln ∗ x y ∗
y∗ σ2 s r ∗ 2 ∗ 2 2 − xσ∗ xy xy − yx∗ < 0. This shows Differentiating V with respect to time t we get, dV dt = − K (x − x ) − x∗ σ1 L y − y that dV dt is negative and hence by Liapunov’s [12], it follows that the positive equilibrium is globally asymptotically stable. 3.4. Stability of steady states Let the two function u = u(x, y) = dx dt and v = v(x, y) = and (2). Then the Jacobian of the system isJ(x,y)
⎧ ∂u ⎪ ⎪ ⎨ ∂x =⎪ ∂v ⎪ ⎩ ∂x
⎫
∂u ⎪ ⎪ ∂y ⎬ ∂v ⎪ ⎪ ⎭ ∂y
=
r−
dy dt
2rx k
from the steady state (x∗ , y∗ ) = (0, 0) of the equation (1) − σ1 − qE σ2 σ1 s − 2y L − σ2
For the steady states (x, y) = (0, 0), we get, J(0,0) =
r − σ1 − qE σ2 σ1 s − σ2
The characteristics equation of the matrix with eigenvalue ω is, r − σ1 − qE − ω 0 |J − ωI| = 0 s − σ2 − ω Now, r − σ1 − qE − ω = 0 or, s − σ2 − ω = 0 ∴ ω1 = r − σ1 − qE, ω2 = s − σ2 Here the eigenvalue ω of the matrix “J” determine the stability of the states. Depending on ω following stability conditions: 1. If the eigenvalue ω < 0, then the steady state is stable, 2. If the eigenvalue ω > 0, then the system is unstable.
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
4. Numerical results and discussion In this section, we have numerically solved our proposed models and discussed their stability and instability. In Sections, we have discussed the numerical simulation of model (1), (2) respectively. We have already discussed existence of equilibrium, positivity analysis and stability of steady states. The model is simulated using the parameter values the total effort applied for harvesting the fish population in the reserve area E = 0.2. The carrying capacities of fish species in unreserved and reserved area respectively K = 110 and L = 200. The intrinsic growth rate of fish subpopulation are respectively r = 0.03 and s = 0.04. The catch ability fish coefficient q = 0.001. The fish subpopulation of unreserved and reserved area are respectively σ1 = 0.05 and σ1 = 0.05 and qE denote the fishing mortality rate of the model (1). Then the result is displayed in the Fig. 1. The displayed graph depicts the biomass densities of reserved area unreserved area where the population increased and decreased within 90 days. First of all the Fig. 1 has shown the population of unreserved area and reserved area are same on the 9 days. In 90 days the population of biomass densities of unreserved area increases above 100. At the same time the population of biomass densities of reserved area increases above 100. Again, Fig. 2 the population of reserved area increases above 900. At the same time the population of unreserved area increases above 400.and unreserved area are same 10 days. If we input the value of the total effort applied for harvesting the fish population and fish catch ability coefficient increases then the population decreases. If the population of the total effort applied for harvesting the fish population and fish catch ability coefficient decreases so the population decreases. Similarly by changing the values of parameter then the graph is changed. The Fig. 3 and Fig. 4 are changing the parameter values then the graph is changed. If we are changing the parameter values of intrinsic growth rate of unreserved and reserved area respectively and are increasing then the population are increasing.
Fig. 1. The graph of fishery model of the population with time (90 Days)
Fig. 2. Variation of population with time (90 Days) for E = 0.002, K = 1000, L = 200, r = 0.03, s = 0.04 and q = 0.001
341
342
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
Fig. 3. Variation of population with time (90 Days) for E = 0.002, K = 1000, L = 200, r = 0.299, s = 0.004 and q = 0.00001
Fig. 4. Variation of population with time (90 Days) for E = 0.002, K = 1000, L = 200, r = 0.03, s = 0.4 and q = 0.00001
Fig. 5. Variation of population with time (90 Days) for E = 4, K = 1000, L = 200, r = 0.03, s = 0.004 and q = 0.1
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344
Fig. 6. Variation of population with time (90 Days) for E = 2.8, K = 1000, L = 200, r = 0.3, s = 0.006 and q = 0.2
If we are using the parameter values of carrying capacities of fish spices in the unreserved and reserved area fish production respectively K and L area increasing then the graph is increasing and by changing the parameter values K and L are increasing then the fish production is increasing. Similarly, the Fig. 5 and Fig. 6 below the X axis. So the population decreases. The maximum fish production is depend on crying capacities intrinsic growth rate of the unreserved and reserved area and harvesting rate or fishing mortality. In this Figure we observed that, the harvesting rate of fishing mortality is increasing then the graph is decreasing. Again the fishing mortality is decreasing then the graph is increasing. In this figure we observed that, by changing the values of parameter the characteristics figure is changed. So the model is fishery management model it can be find out of minimum and maximum fish production of reserved and unreserved area. 5. Conclusion In this mathematical modeling the system has been assumed that the two ecosystems; one free fishing zone and the other reserved zone. In this paper, we analysed the fish production by using many mathematical models. It has been found that the total user’s cost of harvest per unit of effort equals to the discounted value of the future marginal profit of the effort at the steady-state level. It has also been noted that if the discounted rate increases, then the economic rent decreases and even may tend to zero if the discounted rate tend to infinity. The possibility that the predator species in this model be itself harvested for commercial purposes, or that its non-use value be an important element in the overall economic value of the marine ecosystem have not been considered in this paper. Their inclusion in the model would surely lead to characterize another set of trade-offs which are bound to be at the centre of marine reserve policies. We have derived existence of equilibrium and globally asymptotically stable and unstable. Finally, we have been numerical simulations of fish production with graphical representation with different values of parameter. References [1] L. Berezansky, L. Idels, M. Kipnis, Mathematical model of marine protected areas, IMA J. Appl. Math. 76 (2011) 312-325. [2] G. Birkhoff, G.S. Rota, Ordinary Differential Equations, Ginn, Boston, 1882. [3] M.H.A. Biswas, T. Rahman, N. Haque, Modeling the Potential Impacts of Global Climate Change in Bangladesh: An Optimal Control Approach, J. Fundam. Appl. Sci. 8 (2016) 1-19. [4] M.H.A. Biswas, On the Evolution of AIDS/HIV Treatment: An Optimal Control Approach, Curr. HIV Res. 12 (2014) 1-12. [5] M.H.A. Biswas, M.M. Haque, Nonlinear Dynamical Systems in Modeling and Control of Infectious Diseases, Springer Pro. Math. & Stat. Vol. 164 (2016). [6] M.H.A. Biswas, L.T. Paiva, M.D.R. de Pinho, A SEIR model for control of infectious diseases with constraints, Math. Bios. and Engi. 11 (2014) 761-784. [7] C.W. Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resources, Wiley, New York, 1976. [8] T. Das, R.N. Mukherjee, K.S. Chaudhuri, Bioeconomic harvesting of a prey-predator fishery, J. Bio. Dyn. 3 (2009) 447-462. [9] M. De Lara, L. Doyen, Sustainable Management of Natural Resources: Mathematical Models and Methods, Springer, Verlag Berlin, 2008. [10] B. Dubey , P. Chandra, P. Sinha, A Model for Fishery Resource with Reserve Area, Nonl. Ana.: Real Worlld Appl. 4 (2003) 625-637.
343
344 [11] [12] [13] [14] [15] [16] [17] [18] [19]
Haider Ali Biswas et al. / Procedia Engineering 194 (2017) 337 – 344 Y. Kitabatake, A Dynamic Predator-Prey Model for Fishery Resources: a Case of Lake Kasumigaura, Environ, Plan. A. 14 (1982) 225-235. J. La Salle, S. Lefschetz, Stability by Liapunovs Direct Method with Applications, Academic Press, New York, London, 1961. A. Leung, A. Wang, Analysis of models for commercial fishing: mathematical and economical aspects, Econometrica. 44 (1976) 295-303. M. Mesterton-Gibbons, On the Optimal Policy for Combined Harvesting of Predator Prey, Nat. Res. Model. 3 (1996) 235-244. M.K. Mondal, M.H.A. Biswas, Modeling the Transmission Dynamics of the Nipah Virus Infection, Proc. of 1st Int. Con. on Math. and Its Appl., Mathematics Discipline, Khulna University, Khulna, Bangladesh, 23rd December 2015, pp. 22-26. J.D. Murray, Mathematical Biology: Biomathematics, Second Edition, Springer, New York, 1989. R.M. Neilan, S. Lenhart, An Introduction to Optimal Control with an Application in Disease Modeling, DIMACS Series in Disc. Math. 75 (2010) 67-81. X. Song, L. Chen, Optimal Harvesting and Stability for a Two-species Competitive System with Stage Structure, Math. Biosci. 170 (2001) 173-186. X. Zhang, L. Chen, A.U. Neumann, The Stage Structured Predator-prey Model and Optimal Harvesting Policy, Math. Biosci. 168 (2000) 201-210.