# Modelling and Simulation of Competition for Students’ Population with Holling Type II Response

Authors: Brian Odhiambo; Cyrus Ngari; Patrick Kimani; Peter Njori
Journal: Journal of Mathematics
 (2023)

---

## Abstract

The increase in the country’s population attracted the establishment of more schools, both public and private schools, to cater for the increasing number of students. However, there have been dynamics of students’ population both in public and private schools through transfer from one category of school to the other, through completion of the learning period, and through dropout due to unknown reasons which have subjected both the public and private schools to compete in order to maintain a good number of students. In this work, a modified Lotka–Volterra model of schools and nonenrolled entities population in the education system is studied. Private schools and nonenrolled entities play the role of a predator in public schools. Again, public schools and nonenrolled entities play the role of predators in private schools. Holling type II functional responses have been integrated in the analysis of the Lotka–Volterra model. The equilibrium points are established and their stability are determined using the Routh–Hurwitz criterion and eigenvalue method. Global stability has been done for the positive equilibrium point. Bifurcation is also done around the positive equilibrium point. Finally, a graphical illustration of various parameter is derived to show their effect on schools when they are varied. It is revealed that the increase in parametersθ2, θ3, and η3 greatly affects the schools population as they are the ones leading to predation in school. Therefore, proper strategies should be developed to focus on reducing the mentioned parameters to avoid leading schools’ population to extinct.


---

## Body

## 1. Introduction

Education stakeholders and strategists across countries work tirelessly to ensure that their respective country’s education system runs smoothly. Prey-predator models are usually used to study the interaction of animals and plants. This study applies the concept of a prey-predator model to study students’ population dynamics in schools. Prey-predator models can show different behaviors such as steady states, oscillations, and bifurcations depending on the parameters [1]. The study of the dynamics of prey-predator systems is one of the most dominant subjects in mathematical ecology due to its universal existence and importance [1]. The schools category population is considered in this study. The first category is the public schools which are owned and controlled by the government in their day-to-day operations. Teachers in such schools are government employees. The second category is private schools. United Nations Educational Scientific Cultural Organization (UNESCO) defines them as any school that is not operated by the government but is controlled and managed, whether for profit or not, by a private body like a community, foundation, faith-based organization, nongovernmental organizations (NGOs), private proprietor, or private enterprise [2]. This study also considers students who are out of school (nonenrolled) but are of school-going age as the third category and termed it as the nonenrolled category.The first two categories of schools face challenges like peer influence among students and parents causing transfer from one category to the other, population predation between categories, and competition for students [3]. These factors create population dynamics in the respective schools. The school population is an important parameter for economic development. Predicting the future population is a core factor for development planning and decision-making [4]. Quality of education has decreased in public schools owing to overcrowded classrooms attracting more private schools [5]. Private schools tend to concentrate in an area where public schools perform poorly [5]. In the same way Italian mathematician Volterra in the model of fishing argued that, when fishing was good, the number of fishermen increased. After sometime, fish declined and the number of fishermen also declined. After some time, the cycle repeated [6]. However, the development of private schools increases the competition in the schooling system [7]. This theoretically leads to efficiency gain in terms of both cost and quality as private and public schools compete to attract students.Due to prevailing situations in the country, where the country’s activities are geared towards achieving the Sustainable Development Goals (SDGs) which in addition is aimed at the realization of the vision 2030, the countries activities are thus focused on a better education system. This initiated the transfer of students to different categories of schools to acquire better education for self-reliance [3, 8]. Generally, the education sector shapes the nation through good leadership which can result to better development. For this matter, it is a secondary requirement to many people throughout the world. From the literature review, a range of social statistical techniques have been used to explain school population dynamics [2, 3, 5, 7–9]. But till now, very little research has been carried out on the interaction of public schools, private schools, and the nonenrolled entity population. Thus, it is important that approaches to modelling school population dynamics be considered based on past enrolment data incorporating the effects of competition and predation.In this paper, a modified prey-predator model for three categories, namely; public schools, private schools, and the nonenrolled entity has been formed. So the uniqueness in this paper is that both public and private schools act as a predator as well as a prey. Nonenrolled entity is purely a predator. In this work, our main goal is to study the dynamical behavior of these three categories of schools, which is very important for the economic development [4].Usually, Lotka–Volterra model comprises two equations, but this study presents the Lotka–Volterra model with three equations hence an improved one [10]. Holling extended the system of Lotka–Volterra and came up with the idea of functional response [11]. Thus, Holling type II functional response is used since our model adopts a logistic growth as proposed by [11]. Holling argued that when the prey is scattered in an area, the predator has to spend more of their time searching them as opposed to when they are concentrated in an area, where the predator will only have to spend most of their time handling the prey. Holling type II functional response with logistic growth is thus given by [1].Lotka–Volterra model can be used to show the competitive relationship between organisms, companies, markets, or organizations [12–19]. Thus, this study considers it as also much applicable in schools’ population competition. The main purpose of modelling population interactions is to understand what causes such fluctuations [6]. This is only possible when the parameters of the model are varied [1, 19]. This study is helpful to the country for future population prediction [20]. This can help the country to be prepared to cater for the future changes that are likely to happen in schools. However, parameter estimation in this model was challenging. This is due to difficulty in getting predation values as well as the nonenrolled entity values as this data is sensitive and lies with the government. These are the challenges faced in modelling.For the analysis of the model, the feasible region of the model is established and boundedness of the model is determined. The local stability of the equilibrium points determined with the help of a Wolfram Mathematica software is studied using the Routh–Hurwitz criterion [1, 21]. For the global stability, this study adopted the Lyapunov function. Finally, numerical simulation is done to verify the stabilities which cannot be determined analytically.


## 2. Model Formulation and Assumptions

In this work, the coexistence of the public schools, private schools, and nonenrolled entity is discussed. To formulate the model, we letxt, yt, and zt denotes the populations of public schools, private schools, and nonenrolled entity, respectively, at any given time t and they are subject to non-negative initial conditions x0=x0≥0, y0=y0≥0, and z0=z0≥0. The parameters γ and K represents the recruitment rate in the categories and the environmental carrying capacity, respectively.Figure1 shows a modified Lotka–Volterra model with Holling type II is formulated where public, private, and nonenrolled entities assume the role of prey and predator depending on the scenario. In this model, the predator (private schools) population grows by enrolment and transfer from public schools and again, the predator (nonenrolled) gains population as a result of decreased enrolment in schools and dropouts. The total recruitment into the said categories is denoted by π which is subdivided into the following categories πγ1 as the recruitment in public schools, πγ2 as the recruitment in the private schools, and (1−γ1−γ2) π for the nonenrolled entities. Since public and private schools consider their carrying capacity during enrolment, their model adopts a logistic growth [22], while the nonenrolled entities assume an exponential growth. The population dynamics in the categories is due to the predation effect known as Holling type II response [11], exit rate, and the competition as represented by η,μ, and θ, respectively. To develop an equation for a category, for example, public school, we take the enrolment which is carried out by considering the carrying capacity, then we subtract the exit, transfer to private, and drop out due to unknown factor, less predation from nonenrolled, and private schools, and then we add the population that is gained by the public when predating private schools. By this, we arrive at equation (1). The exit rate is a fraction of the study period which is a quarter for the secondary schools. For the nonenrolled entities, the study considers their exit rate as the period which a person cannot come back to school to study. Anybody who has not joined the two categories of schools is considered to be in the nonenrolled category even if they joined tertiary institutions.Figure 1
Model flow chart.Moreover, since a model is just a representation of the reality, it cannot capture all the features of the school setup. Therefore, this study considers the following assumptions in order to formulate the model:(i)
The study does not include special schools and group of schools in the country(ii)
Dropout rate due to other factors other the ones discussed in the model is insignificant(iii)
The rate of student retention in a class will be insignificant(iv)
Public and private schools cannot predate the nonenrolled entitiesFrom the above assumptions and illustrations, a mathematical model of the three categories has been developed in Figure1.


## 3. Model Equations

(1)dxdt=γ1πx1−xk1−μ+θ1+θ3x−η1xyB+x−η3xzB+x+η4xyB+y,(2)dydt=γ2πy1−yk2−μ+θ2y+θ1x+η1xyB+x−η2yzB+y−η4xyB+y,(3)dzdt=1−γ1−γ2πz+θ3x+θ2y+η3xzB+x+η2yzB+y−μz,with initial conditions x0≥0, y0≥0, and z0≥0.


## 4. Feasible Region of the Solution

Consider equation (1), by inspection method the term η4xy/B+y≥0 and γ1πx1−x/k1≥0. Thus, dx/dt≥−μ+θ1+θ3x−η1xy/B+x−η3xz/B+x. It follows, dx/x≥−μ+θ1+θ3−η1y/B+x−η3z/B+xdt. On integrating, lnx≥−μ+θ1+θ3−η1ys/B+xs−η3zs/B+xsds+c. xt≥x0e∫0t−μ+θ1+θ3+η1ys/B+xs+η3zs/B+xsds. Exponential of a negative is always positive.Consider equation (2), by inspection method, the term γ2πy1−y/k2,η1xy/B+x,  and θ1x is ≥0.Thus,yt≥y0e∫0t−μ+θ2+η2zs/B+ys+η4xs/B+ysds≥0, i.e., the exponential of a negative is positive.Similarly, for equation (3), we have zt≥z0e−μt≥0.


## 5. Boundedness of the Solution

For the system to be mathematically meaningful, it is necessary to show that its state variables are positive and bounded for all timet [1]. That is, the solution of the system with a positive initial value will remain positive for all time t≥0.(4)dNdt=γ1πx1−xk1−μ+θ1+θ3x−η1xyB+x−η3xzB+x+η4xyB+y+γ2πy1−yk2−μ+θ2y+η1xyB+x−η2yzB+y+θ1x−η4xyB+y+1−γ1−γ2πz+θ3x+θ2y+η3xB+x+η2yzB+y−μz.On simplifying the above,(5)dNdt=γN1−Nk4−μN,Nt=Cetγγ−μk4−etμ+etγCγ,t=0,C=N0−γk4+μk4+γN0.Thus,(6)Nt=etγγ−μk4N0etγγN0−etμ1−γ+μk4+γN0.


## 6. Steady States

To study the stability of our proposed model, the equilibrium points of the system need to be determined. The possible equilibrium to be considered are(7)−μ+π−2πxk1γ1−yBη1x+B2−zBη3x+B2+yη4y+B−θ1−θ3−xη1x+B−xyη4y+B2+xη4y+B−xη3x+ByBη1x+B2−yη4y+B+θ1−μ+π−2πyk2γ2+xη1x+B−zBη2y+B2−xBη4y+B2−θ2−yη2y+B−xzη3x+B2+zη3x+B+θ3−yzη2y+B2+zη2y+B+θ2−μ+πα+yη2y+B+xη3x+B.Case 1: No schools exist;x=0, y=0, z=0The trivial equilibrium points areE0=0,0,0Case 2: Only public schools exist;x≠0, y=0, z=0(8)γ1πx1−xk1−μ+θ1+θ3x=0.There are two equilibrium points sayE1=x1,x0 where(9)x0=0andx1=k1μ−πγ1+θ1+θ3πγ1.Case 3: Only private schools exist;x=0, y≠0, z=0.(10)γ2πy1−yk2−μ+θ2y=0.The equilibrium points areE2=y2,y0,wherey0=0 and y2=k2−μ+πγ2−θ2/πγ2.Case 4: No enrolment in schools;x=0, y=0, z≠0.(11)1−γ1−γ2πz−μz=0.The equilibrium points areE3=z0 and when α=μ,wherez0=0.Case 5: No private schools only;x≠0, y=0, z≠0.Three equilibrium points are obtainedE04=x04,z04, E14=x14,z14 and E24=x24,z24, where(12)x04=0andz04=0,x14=m−nqandz14=m−n−m−n/q+k1α1+k1β8qk1α5,x24=m+nqandz24=m+n−m−n/q+k1α1+k1β8qk1α5,wherem=α1α5+k1α1+β8φ3, n=k1α54k1α1α5−α1+β1φ3+α1α5+k1α1+β8φ32/k12α52, q=2α1φ3, α5=μ−α4, β5=μ−α4/θ3, β6=α1β5, β7=β1β5, β8=−β1+θ3, φ3xz=η3xz/B+x, α4=1−γ1−γ2π, β1=μ+θ1+θ3, α1=γ1π, and α2=γ2π.Case 6: Full enrolment;x≠0, y≠0, z=0.Three equilibrium points are obtained;E05=x05,y05, E15=x5,0, and E25=0,y5, where(13)x05=0andy05=0,x15=s15andy15=r4r2−r22+s15r3−s15r1+r5,x25=s25andy25=r4r2−r22+s25r3−s25r1+r5,wheres15=s1−s2/s3, s25=s1+s2/s3, s1=r3r42r5φ1−φ42+2α1α1−β1+r2r4−φ1+φ4ψ1, r1=α1ψ1, r2=α2−β2, r3=4α2ψ2, r4=1/2α2ψ2, r5=α1−β1+θ1, β1=μ+θ1+θ3, β2=μ+θ2, α1=γ1π, s2=r42φ1−φ42r32r42r52φ1−φ42−4r1r3α1−β1α1−β1+2r2r4−φ1+φ4+4r3r5α1α1−β1+r2r4−φ1+φ4ψ1+4r22α12ψ12, α2=γ2π, and s3=2r1r3r42φ1−φ42+2α12ψ12.Case 7: No public schools only;x=0, y≠0, z≠0.There are three equilibrium points, namely,E06y06,z06 and E16y16,z16, where,y06=0 and z06=0 which corresponds to the trivial equilibrium point.(14)y16=p1+−μ+α2p2p3andz16=p1+p2−μ+α2p3α1−β2+θ2−α1p1+p2−μ+α2ψ3p32μ−α2,y26=p1+p2μ−α2p3andz26=p1+p2μ−α2p3α1−β2+θ2−α1p1+p2μ−α2ψ3p32μ−α2,whereα1=γ2π, α2=1−γ1−γ2π, β2=μ+θ2, ψ2yz=η2yz/B+y, ψ3y=y/k2, p1=−β2+θ2ψ2+α1ψ2+μ−α2ψ3, p2=α1−β2+θ22ψ22/μ−α22+2α1α1−β2+θ2ψ2ψ3/μ−α2+α12ψ32, and p3=2α1ψ2ψ3.Case 8: Coexistence;x≠0, y≠0, z≠0.There are many equilibrium points areE07=x07,y07,z07, E17=x17,y17,z17, E27=x27,y27,z27, x07=0, y07=0, and z07=0. This corresponds to a trivial equilibrium point.(15)x17=0,y17=12−−μ+p1+α3φ2+α2−β2+θ2α2ψ2andz17=−β2+θ2φ2+α2φ2+−μ+α3ψ2+α2ψ2p12φ22,x27=p2p3−φ3−μ−p4+α3+p2p3φ3−p2θ3φ3p5,y27=0andz27=−μ−p4+α3+p2p3φ3−p2θ3φ3p5,x37=p2p3−φ3−μ−p4+α3+p2p3φ3+p2θ3φ3p5,y37=0andz37=−μ−p4+α3+p2p3φ3+p2θ3φ3p5,x47=p2hp3−m+nφ1−φ4h,y47=m+nhandz47=0,x57=p11φ6−p22+p66φ6/p77−p33+p44−p66φ6/p772p55,y57=p22+p66φ6/p77−p33+p44−p66φ6/p772p55andz57=p66p77,x67=p11φ6−p22+p66φ6/p77+p33+p44−p66φ6/p772p55,y67=p22+p66φ6/p77+p33+p44−p66φ6/p772p55andz67=p66p77,where β1=μ+θ1+θ3, β2=μ+θ2, p3=α1−β1, p11=μ−α3−θ2−θ3, p22=−α2+β2+θ1−p11φ3/φ6+p11φ5/φ6, p33=−4p11θ1−α2φ2−φ3+φ5/φ6, p44=α2−β2−θ1+p11φ3/φ6−p11φ5/φ6, p55=2−α2φ2−φ3+φ5, p66=−p112p552+2p11p22p55φ6−p222φ62+p33φ62+p442φ62, p77=2φ62−p11p55+p22φ6+p44φ6,p1=4α2−β2θ2φ22+−α2+β2+θ2φ2+α2μ−α3ψ22/α22ψ22, p2=1/α1ψ1, p5=2p2φ32, p4=4p22p3θ3φ32+−μ+α3+p2p3φ3−p2θ3φ32, γ1π=α1, γ2π=α2, 1−γ1−γ2π=α3, ψ1x=x/k1, ψ2y=y/k2, φ1xy=η1xy/B+x, φ2xy=η3yz/B+y, φ3xy=η3xz/B+x, φ4xy=η3xy/B+y, b1=b11+b12, b2=b21+b22+b23/b24, n2=1/2−n1α2φ2−φ32+φ3φ6, n3=−4α1θ1−β1θ1−zθ1φ5, n2=1/2−n1α2φ2−φ32+φ3φ6, n1=α1φ1, n3=−4α1θ1−β1θ1−θ1φ5, n4=−n1α2+n1β2−α1φ3+β1φ3+θ1φ3+α1φ6−β1φ6, b11=−n12n2θ2φ4+−n3n1α2φ2+φ32−φ3φ6n2φ3φ5−n1n2φ4+φ5α1−n2n4φ3−θ31+n2φ3φ3−2φ6+β1−1+n2φ3φ6, b12=β1−1+n2φ3φ6+n1−μ+α3+n2n4φ4+θ3φ3φ4+θ2φ5φ3−2φ6−β1φ4φ6, b21=n11/n124n12n2φ42+φ521+n2φ3φ3−2φ6+2n1n2φ4φ5−φ3+φ6α1−β1θ3+n2−n1θ2+θ3φ3−n4+β1φ6+−n3n1α2φ2+φ32−φ3φ6, b22=n12n2θ2φ4+−α1+β1φ5+φ5θ31+n2φ3φ3−2φ6−n2φ3−n4+β1φ6+−n3n1α2φ2+φ32−φ3φ6, b23=n1μ−α3+n2−n4φ4−θ3φ3φ4−θ2φ5φ3−2φ6+φ4β1φ6+−n3n1α2φ2+φ32−φ3φ62, b24=2n12n2φ42+φ521+n2φ3φ3−2φ6+2n1n2φ4φ5−φ3+φ6, m=α2−β2+p2p3φ1−p2θ1φ1−p2p3φ4+p2θ1φ4, h=2p2φ1−φ42+2α2ψ2, and n=α2−β2+p2p3−θ1φ1−φ42+4p2p3θ1p2φ1−φ42+α2ψ2.Given thatφ1, φ2, φ3, φ4, φ5, φ6, ≥0, and B>0.


## 7. Stability Analysis

### 7.1. Local Stability

The local stability of the system is determined by the nature of the eigenvalues of the variational matrix, that is, for negative eigenvalues imply a stable otherwise unstable. In the case where the matrix is large such that the nature of eigenvalues cannot be directly determined, we will apply Routh–Hurwitz criteria to determine stability.The system has eight solutions.E00,0,0,E1K1,0,0,E20,K2,0,E30,0,K3,E4x∗,y∗,0,E5x∗,0,z∗,E60,y∗,z∗,E7x∗,y∗,z∗.The variational matrix of the system is given by(16)b11−xη1x+B1−xyη4y+B42+xη4y+B4−xη3x+B3yB1η1y+B12−yη4y+B4b22−yη2y+B2xzη3x+B32+zη3x+B3+θ3−yzη2y+B22+zη2y+B2+θ2b33,where(17)b11=−μ+π−2πxk1γ1−yB1η1x+B12−zB3η3x+B32+yη4y+B4−θ1−θ3,b22=−μ+π−2πyk2γ2−xη1x+B1+zB2η2y+B22−xB4η4y+B42−θ2,andb33=−μ+πα+yη2y+B2+xη3x+B3.Theorem 1.
The trivial equilibrium point,E0, is stable if the eigenvalues of the variational matrix are all less than zero. Otherwise, the system will be unstable.Proof.
The variational matrix of (16) evaluated at E0 is given by(18)−μ+πγ1−θ1−θ300θ1−μ+πγ1−θ20θ3θ2−μ+πα.
By lettingρ=θ1+θ3 and α=1−γ1−γ2, eigenvalues are λ1=πα−μ,λ2=πγ2−θ2−μ,andλ3=πγ1−ρ−μ.
The equilibrium pointE00,0,0 is stable if and only if πγ2<θ2+μ, πα<μ, and πγ1<ρ+μ. Otherwise, the system will be unstable.Theorem 2.
Equilibrium point,E1, at the point x≠0, y=0, z=0, by Routh–Hurwitz criterion is locally asymptotically stable if A1,A3>0 and A1A2−A0A3>0. Otherwise, the system will be unstable.Proof.
x0=0 corresponds to a trivial equilibrium point.(19)x1=k1μ−πγ1+θ1+θ3πγ1.
The variational matrix is given by(20)a11a12−a13a21a220a31a32a33,where(21)a11=−μ+π−2πx1k1γ1−θ1−θ3,a12=−x1η1x1+B+x1η4B,a13=−x1η3x1+B,a21=θ1,a22=−μ+πγ2+x1η1x1+B−x1η4B−θ2,a31=θ3,a32=θ2anda33=−μ+πα+x1η3x1+B.
The characteristic equation is given by,p3+p2A1+pA2+A3=0, where A1=−a11+a22+a33, A2=a11a22+a22a33+a11a33−a12a21+a13a31 and A3=a13a21a32−a11a22a33+a12a21a33−a13a31a22.Theorem 3.
Equilibrium point,E2, at the point x=0, y≠0, z=0 is stable if and only if the eigenvalues of the variational matrix are less than zero.Proof.
The variational matrix is given by(22)a1100a21a220a31a32a33,where a11=−μ+πγ1−yη1/B+yη4/y+B−θ1−θ3, a21=yη1/B−yη4/y+B+θ1, a22=−μ+π−2πy/k2γ2−θ2, a23=−yη2/y+B, a31=θ3,a32=θ2, a33=−μ+πα+yη2/y+B.
The eigenvalues are−μ+πγ1−yη1/B+yη4/y+B−θ1−θ3, −μ+π−2πy/k2γ2−θ2 and −μ+πα+yη2/y+B.
The eigenvalues are stable if−μ+πγ1−yη1/B+yη4/y+B−θ1−θ3<0, −μ+π−2πy/k2γ2−θ2<0, and −μ+πα+yη2/y+B<0.Theorem 4.
Equilibrium point,E3, for point x=0, y=0, z≠0 is stable if and only if μ>πα, μ>πγ2+θ2 and μ>πγ1+ρ otherwise the system will be unstable.Proof.
The variational matrix is given by(23)−μ+πγ1−θ1−θ300θ1−μ+πγ1−θ20θ3θ2−μ+πα.
By choosingρ=θ1+θ3 and α=1−γ1−γ2, the eigenvalues are(24)λ1=πα−μ,λ2=πγ2−θ2−μ,λ3=πγ1−ρ−μ.
The eigenvalues will be negative if and only ifμ>πα, μ>πγ2+θ2 and μ>πγ1+ρ for it to be stable, otherwise unstable.Theorem 5.
Equilibrium point,E4, for case x≠0, y=0, z≠0 by the Routh–Hurwitz criterion is stable if D1,D3>0 and D1D2>D0D3 otherwise unstable.Proof.
The variational matrix of the equilibrium point atx≠0, y=0, z≠0 is given by(25)d11d12−d13d21d220d31d32d33,(26)d11=−μ+π−2πxk1γ1−zBη3x+B2−θ1−θ3,d12=−xη1x+B+xη4B,d13=xη3x+B,d21=θ1,d22=−μ+πγ2+xη1x+B−zη4B−xη4B−θ2,d31=−xzη3x+B2+zη3x+B+θ3,d32=zη2B+θ2,d33=−μ+πα+xη3x+B.
The characteristic equation is given byr3+r2D1+rD2+D3=0, where(27)D1=−d11+d22+d33,D2=d11d33+d22d33+d11d22−d12d21+d13d31,D3=d13d21d32−d11d22d33+d12d21d33−d13d31d22.Theorem 6.
Equilibrium point,E5, for the case x≠0, y≠0, z=0 is stable if C1,C3>0 and C1C2>C0C3, otherwise the system will be unstable.Proof.
The variational matrix for the equilibrium point atx≠0, y≠0, z=0 is given by(28)c11c12c13c21c22−c23c31c32c33,where(29)c11=−μ+p−2πxk1γ1−yBη1x+B2+yη4y+B−θ1−θ3,c12=−xη1x+B−xyη4y+B2+xη4y+B,c13=xη3x+B,c21=yBη1x+B2−yη4y+B+θ1,c22=−μ+π−2πyk2γ2+xη1x+B−xBη4y+B2−θ2,c23=yη2y+B,c31=θ3,c32=θ2,c33=−μ+πα+yη2y+B+xη3x+B.
The characteristic equation is given bym3+m2C1+mC2+C3=0, where C1=−c11+c22+c33, C2=c11c33+c22c33+c11c22−c12c21+c23c32+c13c31, and C3=c12c23c31−c11c22c33+c13c21c32+c12c21c33−c31c22c13.Theorem 7.
Equilibrium point,E6, for x=0, y≠0, z≠0, by Routh–Hurwitz is stable if G1,G3>0 and G1G2>G0G3, otherwise the system will be unstable.Proof.
The variational matrix for the equilibrium point atx=0, y≠0, z≠0 is given by(30)g1100g21g22−g23g31g32g33,where(31)g11=−μ+πγ1−yη1B−zη3B+yη4y+B−θ1−θ3,g22=−μ+π−2πyk2γ2−zBη2y+B2−θ2,g23=yη2y+B,g31=zη3B+θ3,g32=−yzη2y+B2+zη2y+B+θ2andg33=−μ+πα+yη2y+B.
The resulting characteristic equation is given bys3+s2G1+sG2+G3=0, where G1=−g11+g22+g33, G2=g11g33+g22g33+g11g22+g23g32, and G3=−g11g22g33+g11g23g32.Theorem 8.
Equilibrium point,E7, for the case x≠0, y≠0, z≠0 is stable if K1,K3>0 and K1K2>K3K0, otherwise the system will be unstable.Proof.
The variational matrix for the equilibrium point atx≠0, y≠0, z≠0 is given by(32)k11k12−k13k21k22−k23k31k32k33.where(33)k11=−μ+π−2πxK1γ1−yBη1x+B2−zBη3x+B2+yη4y+B−θ1−θ3,k12=−xη1x+B2−xyη4y+B2+xη4y+B,k13=xη3x+B,k21=−yBη4x+B2−yη4y+B+θ1,k31=−zxη3x+B2+zη3x+B+θ3,k23=yη2y+B,k22=−μ+π−2πyk2γ2+xη1x+B−zBη2y+B2−xBη4y+B2−θ2,k32=−zyη2y+B2+zη2y+B+θ2andk33=−μ+πα+yη2y+B+xη3x+B.
The resulting characteristic equation is given by(34)r3+r2K1+rK2+K3=0,where K1=−k11+k22+k33, K2=k11k33+k22k33+k11k22−k12k21+k23k32+k13k31, and K3=k12k23k31−k11k22k33+k13k21k32+k12k21k33−k31k22k13.


### 7.2. Global Stability

Linear stability analysis tells us how a system behaves near an equilibrium point. It does not however tell us anything about what happens farther away from equilibrium [23]. For higher dimensions, we consider using a Lyapunov function. This paper considers Lyapunov proposed by [1].Theorem 9.
Let(35)V=X−X∗22+δ1Y−Y∗22+δ2Z−Z∗22,be the Lyapunov function, where δ1,δ2>0 are to be taken properly such that the above function satisfies the Lyapunov conditions.V′E∗=0,  where E∗=X∗,Y∗,Z∗ and VX,Y,Z>0. The time derivative of V is dV/dt≤0 implies that the system is just stable and dV/dt<0 implies that E∗ is globally asymptotically stable.Proof.
(36)V=X−X∗22+δ1Y−Y∗22+δ2Z−Z∗22.
Getting a time derivative of (36),(37)dVdt=X−X∗dXdt+δ1Y−Y∗dYdt+δ2Z−Z∗dZdt,by substituting the values of dX/dt,dY/dt,dZ/dt results to(38)dVdt=X−X∗γ1πX1−Xk1−μ+θ1+θ3X−η1XYB+X−η3XZB+X+η4XYB+Y+δ1Y−Y∗γ2πY1−Yk2−μ+θ2Y+η1XYB+X−η2YZB+X+θ1X−η4XYB+Y+δ2Z−Z∗1−γ1−γ2πZ+θ3X+θ2Y+η3XZB+X+η2YZB+Y−μZ.
By factorizingX,Y,andZ,(39)dVdt=X−X∗γ1π1−Xk1−μ+θ1+θ3−η1YB+X−η3ZB+X+η4YB+YX+δ1Y−Y∗γ2π1−Yk2−μ+θ2+η1XB+X−η2ZB+X+θ1XY−η4XB+YY+δ2Z−Z∗1−γ1−γ2π+θ3XZ+θ2YZ+η3XB+X+η2YB+Y−μZ.
At special cases,X=X∗, Y=Y∗, and Z=Z∗ implying that X−X∗, Y−Y∗, and Z−Z∗ are also equilibrium points, then substituting results to(40)dVdt=X−X∗γ1π1−Xk1−μ+θ1+θ3−η1YB+X−η3ZB+X+η4YB+YX−X∗+δ1Y−Y∗γ2π1−Yk2−μ+θ2+η1XB+X−η2ZB+X+θ1XY−η4XB+YY−Y∗+δ2Z−Z∗1−γ1−γ2π+θ3XZ+θ2YZ+η3XB+X+η2YB+Y−μZ−Z∗.
Further simplifying,(41)dVdt=−X−X∗2γ1π−1+Xk1−μ+θ1+θ3+η1YB+X+η3ZB+X−η4YB+Y−δ1Y−Y∗2γ2π−1+Yk2+μ+θ2−η1XB+X−η2ZB+X−θ1XY+η4XB+Y−δ2Z−Z∗2−1+γ1+γ2π−θ3XZ−θ2YZ−η3XB+X−η2YB+Y+μ.
X−X∗2=0 is the same as X2−XX∗−XX∗+X∗2=0 implying XX−X∗=X∗X−X∗. Thus X=X∗. Similarly, Y=Y∗ and Z=Z∗. Since δ1>0 and δ2>0, then they can be chosen such that dV/dt≤0,  and that the equilibrium E∗ is globally asymptotically stable.


## 7.1. Local Stability

The local stability of the system is determined by the nature of the eigenvalues of the variational matrix, that is, for negative eigenvalues imply a stable otherwise unstable. In the case where the matrix is large such that the nature of eigenvalues cannot be directly determined, we will apply Routh–Hurwitz criteria to determine stability.The system has eight solutions.E00,0,0,E1K1,0,0,E20,K2,0,E30,0,K3,E4x∗,y∗,0,E5x∗,0,z∗,E60,y∗,z∗,E7x∗,y∗,z∗.The variational matrix of the system is given by(16)b11−xη1x+B1−xyη4y+B42+xη4y+B4−xη3x+B3yB1η1y+B12−yη4y+B4b22−yη2y+B2xzη3x+B32+zη3x+B3+θ3−yzη2y+B22+zη2y+B2+θ2b33,where(17)b11=−μ+π−2πxk1γ1−yB1η1x+B12−zB3η3x+B32+yη4y+B4−θ1−θ3,b22=−μ+π−2πyk2γ2−xη1x+B1+zB2η2y+B22−xB4η4y+B42−θ2,andb33=−μ+πα+yη2y+B2+xη3x+B3.Theorem 1.
The trivial equilibrium point,E0, is stable if the eigenvalues of the variational matrix are all less than zero. Otherwise, the system will be unstable.Proof.
The variational matrix of (16) evaluated at E0 is given by(18)−μ+πγ1−θ1−θ300θ1−μ+πγ1−θ20θ3θ2−μ+πα.
By lettingρ=θ1+θ3 and α=1−γ1−γ2, eigenvalues are λ1=πα−μ,λ2=πγ2−θ2−μ,andλ3=πγ1−ρ−μ.
The equilibrium pointE00,0,0 is stable if and only if πγ2<θ2+μ, πα<μ, and πγ1<ρ+μ. Otherwise, the system will be unstable.Theorem 2.
Equilibrium point,E1, at the point x≠0, y=0, z=0, by Routh–Hurwitz criterion is locally asymptotically stable if A1,A3>0 and A1A2−A0A3>0. Otherwise, the system will be unstable.Proof.
x0=0 corresponds to a trivial equilibrium point.(19)x1=k1μ−πγ1+θ1+θ3πγ1.
The variational matrix is given by(20)a11a12−a13a21a220a31a32a33,where(21)a11=−μ+π−2πx1k1γ1−θ1−θ3,a12=−x1η1x1+B+x1η4B,a13=−x1η3x1+B,a21=θ1,a22=−μ+πγ2+x1η1x1+B−x1η4B−θ2,a31=θ3,a32=θ2anda33=−μ+πα+x1η3x1+B.
The characteristic equation is given by,p3+p2A1+pA2+A3=0, where A1=−a11+a22+a33, A2=a11a22+a22a33+a11a33−a12a21+a13a31 and A3=a13a21a32−a11a22a33+a12a21a33−a13a31a22.Theorem 3.
Equilibrium point,E2, at the point x=0, y≠0, z=0 is stable if and only if the eigenvalues of the variational matrix are less than zero.Proof.
The variational matrix is given by(22)a1100a21a220a31a32a33,where a11=−μ+πγ1−yη1/B+yη4/y+B−θ1−θ3, a21=yη1/B−yη4/y+B+θ1, a22=−μ+π−2πy/k2γ2−θ2, a23=−yη2/y+B, a31=θ3,a32=θ2, a33=−μ+πα+yη2/y+B.
The eigenvalues are−μ+πγ1−yη1/B+yη4/y+B−θ1−θ3, −μ+π−2πy/k2γ2−θ2 and −μ+πα+yη2/y+B.
The eigenvalues are stable if−μ+πγ1−yη1/B+yη4/y+B−θ1−θ3<0, −μ+π−2πy/k2γ2−θ2<0, and −μ+πα+yη2/y+B<0.Theorem 4.
Equilibrium point,E3, for point x=0, y=0, z≠0 is stable if and only if μ>πα, μ>πγ2+θ2 and μ>πγ1+ρ otherwise the system will be unstable.Proof.
The variational matrix is given by(23)−μ+πγ1−θ1−θ300θ1−μ+πγ1−θ20θ3θ2−μ+πα.
By choosingρ=θ1+θ3 and α=1−γ1−γ2, the eigenvalues are(24)λ1=πα−μ,λ2=πγ2−θ2−μ,λ3=πγ1−ρ−μ.
The eigenvalues will be negative if and only ifμ>πα, μ>πγ2+θ2 and μ>πγ1+ρ for it to be stable, otherwise unstable.Theorem 5.
Equilibrium point,E4, for case x≠0, y=0, z≠0 by the Routh–Hurwitz criterion is stable if D1,D3>0 and D1D2>D0D3 otherwise unstable.Proof.
The variational matrix of the equilibrium point atx≠0, y=0, z≠0 is given by(25)d11d12−d13d21d220d31d32d33,(26)d11=−μ+π−2πxk1γ1−zBη3x+B2−θ1−θ3,d12=−xη1x+B+xη4B,d13=xη3x+B,d21=θ1,d22=−μ+πγ2+xη1x+B−zη4B−xη4B−θ2,d31=−xzη3x+B2+zη3x+B+θ3,d32=zη2B+θ2,d33=−μ+πα+xη3x+B.
The characteristic equation is given byr3+r2D1+rD2+D3=0, where(27)D1=−d11+d22+d33,D2=d11d33+d22d33+d11d22−d12d21+d13d31,D3=d13d21d32−d11d22d33+d12d21d33−d13d31d22.Theorem 6.
Equilibrium point,E5, for the case x≠0, y≠0, z=0 is stable if C1,C3>0 and C1C2>C0C3, otherwise the system will be unstable.Proof.
The variational matrix for the equilibrium point atx≠0, y≠0, z=0 is given by(28)c11c12c13c21c22−c23c31c32c33,where(29)c11=−μ+p−2πxk1γ1−yBη1x+B2+yη4y+B−θ1−θ3,c12=−xη1x+B−xyη4y+B2+xη4y+B,c13=xη3x+B,c21=yBη1x+B2−yη4y+B+θ1,c22=−μ+π−2πyk2γ2+xη1x+B−xBη4y+B2−θ2,c23=yη2y+B,c31=θ3,c32=θ2,c33=−μ+πα+yη2y+B+xη3x+B.
The characteristic equation is given bym3+m2C1+mC2+C3=0, where C1=−c11+c22+c33, C2=c11c33+c22c33+c11c22−c12c21+c23c32+c13c31, and C3=c12c23c31−c11c22c33+c13c21c32+c12c21c33−c31c22c13.Theorem 7.
Equilibrium point,E6, for x=0, y≠0, z≠0, by Routh–Hurwitz is stable if G1,G3>0 and G1G2>G0G3, otherwise the system will be unstable.Proof.
The variational matrix for the equilibrium point atx=0, y≠0, z≠0 is given by(30)g1100g21g22−g23g31g32g33,where(31)g11=−μ+πγ1−yη1B−zη3B+yη4y+B−θ1−θ3,g22=−μ+π−2πyk2γ2−zBη2y+B2−θ2,g23=yη2y+B,g31=zη3B+θ3,g32=−yzη2y+B2+zη2y+B+θ2andg33=−μ+πα+yη2y+B.
The resulting characteristic equation is given bys3+s2G1+sG2+G3=0, where G1=−g11+g22+g33, G2=g11g33+g22g33+g11g22+g23g32, and G3=−g11g22g33+g11g23g32.Theorem 8.
Equilibrium point,E7, for the case x≠0, y≠0, z≠0 is stable if K1,K3>0 and K1K2>K3K0, otherwise the system will be unstable.Proof.
The variational matrix for the equilibrium point atx≠0, y≠0, z≠0 is given by(32)k11k12−k13k21k22−k23k31k32k33.where(33)k11=−μ+π−2πxK1γ1−yBη1x+B2−zBη3x+B2+yη4y+B−θ1−θ3,k12=−xη1x+B2−xyη4y+B2+xη4y+B,k13=xη3x+B,k21=−yBη4x+B2−yη4y+B+θ1,k31=−zxη3x+B2+zη3x+B+θ3,k23=yη2y+B,k22=−μ+π−2πyk2γ2+xη1x+B−zBη2y+B2−xBη4y+B2−θ2,k32=−zyη2y+B2+zη2y+B+θ2andk33=−μ+πα+yη2y+B+xη3x+B.
The resulting characteristic equation is given by(34)r3+r2K1+rK2+K3=0,where K1=−k11+k22+k33, K2=k11k33+k22k33+k11k22−k12k21+k23k32+k13k31, and K3=k12k23k31−k11k22k33+k13k21k32+k12k21k33−k31k22k13.


## 7.2. Global Stability

Linear stability analysis tells us how a system behaves near an equilibrium point. It does not however tell us anything about what happens farther away from equilibrium [23]. For higher dimensions, we consider using a Lyapunov function. This paper considers Lyapunov proposed by [1].Theorem 9.
Let(35)V=X−X∗22+δ1Y−Y∗22+δ2Z−Z∗22,be the Lyapunov function, where δ1,δ2>0 are to be taken properly such that the above function satisfies the Lyapunov conditions.V′E∗=0,  where E∗=X∗,Y∗,Z∗ and VX,Y,Z>0. The time derivative of V is dV/dt≤0 implies that the system is just stable and dV/dt<0 implies that E∗ is globally asymptotically stable.Proof.
(36)V=X−X∗22+δ1Y−Y∗22+δ2Z−Z∗22.
Getting a time derivative of (36),(37)dVdt=X−X∗dXdt+δ1Y−Y∗dYdt+δ2Z−Z∗dZdt,by substituting the values of dX/dt,dY/dt,dZ/dt results to(38)dVdt=X−X∗γ1πX1−Xk1−μ+θ1+θ3X−η1XYB+X−η3XZB+X+η4XYB+Y+δ1Y−Y∗γ2πY1−Yk2−μ+θ2Y+η1XYB+X−η2YZB+X+θ1X−η4XYB+Y+δ2Z−Z∗1−γ1−γ2πZ+θ3X+θ2Y+η3XZB+X+η2YZB+Y−μZ.
By factorizingX,Y,andZ,(39)dVdt=X−X∗γ1π1−Xk1−μ+θ1+θ3−η1YB+X−η3ZB+X+η4YB+YX+δ1Y−Y∗γ2π1−Yk2−μ+θ2+η1XB+X−η2ZB+X+θ1XY−η4XB+YY+δ2Z−Z∗1−γ1−γ2π+θ3XZ+θ2YZ+η3XB+X+η2YB+Y−μZ.
At special cases,X=X∗, Y=Y∗, and Z=Z∗ implying that X−X∗, Y−Y∗, and Z−Z∗ are also equilibrium points, then substituting results to(40)dVdt=X−X∗γ1π1−Xk1−μ+θ1+θ3−η1YB+X−η3ZB+X+η4YB+YX−X∗+δ1Y−Y∗γ2π1−Yk2−μ+θ2+η1XB+X−η2ZB+X+θ1XY−η4XB+YY−Y∗+δ2Z−Z∗1−γ1−γ2π+θ3XZ+θ2YZ+η3XB+X+η2YB+Y−μZ−Z∗.
Further simplifying,(41)dVdt=−X−X∗2γ1π−1+Xk1−μ+θ1+θ3+η1YB+X+η3ZB+X−η4YB+Y−δ1Y−Y∗2γ2π−1+Yk2+μ+θ2−η1XB+X−η2ZB+X−θ1XY+η4XB+Y−δ2Z−Z∗2−1+γ1+γ2π−θ3XZ−θ2YZ−η3XB+X−η2YB+Y+μ.
X−X∗2=0 is the same as X2−XX∗−XX∗+X∗2=0 implying XX−X∗=X∗X−X∗. Thus X=X∗. Similarly, Y=Y∗ and Z=Z∗. Since δ1>0 and δ2>0, then they can be chosen such that dV/dt≤0,  and that the equilibrium E∗ is globally asymptotically stable.


## 8. Bifurcation Analysis

In the Lotka–Volterra system, many parameters are used to describe and formulate the system. But if parameters used in the model are changed, other types of dynamical behavior may occur and the critical parameter values at which such transitions happen are called bifurcation points [1]. This study addresses the transition factors of the prey-predator interactions. This transition with respect to small changes is called Hopf bifurcation. Hopf bifurcation occurs at the point where the system has nonhyperbolic equilibrium with purely imaginary eigenvalues [1].This study considersμ,θ as the bifurcation parameters and μ∗,θ∗ as the critical values of the bifurcation parameters.Theorem 10.
The positive equilibrium over real domainℝ is considered. Let Φ:0,∞⟶ℝ be the continuous differential function of μ. Then, Φμ=M1μM2μ−M3μ. Since μ∗ is the positive root of Φμ=0. Thus, Hopf bifurcation of the interior equilibrium EX∗,Y∗,Z∗ occurs at μ=μ∗ if and only if Φμ∗=0 and R2μ∗R4μ∗+R1μ∗R3μ∗≠0.Proof.
By the first condition above thatΦμ=0, the characteristic equation of the variational matrix in Theorem 8 can be written as r2+k2r+k1=0.
The roots of the equation areα1=−k1,α2=−ik2,α3=ik3.
There exist an interval ofμ∗+ξ,μ∗−ξ since Φμ∗ is continuous of all its roots.
The general form of complex roots is given by(42)θ1μ=ωμ+iξμandθ2μ=ωμ−iξμ.
We verify the transversality conditiondReθj/dμ≠0;j=1,2.
By substituting the general form of the complex roots into the characteristic equation and getting the derivative results to(43)ω3−ω2iξ−3ξ2ω+iξ3+σ1ω2−2iξωσ1−σ1ξ2+σ2ω−iξσ2+σ3.
By implicit differentiation of equation (43), we obtain(44)3ω2−3ξ2+2σ1ω+σ2ω′μ−6ωξ+2σ1ξξ′μ+σ1′ω−σ1′ξ2+σ2′ω+σ2′.
By lettingR1μ=3ω2−3ξ2+2σ1ω+σ2 , R2μ=6ωξ+2σ1ξ , R3=σ1′ω−σ1′ξ2+σ2′ω+σ2′, and R4=2σ1′ωξ+σ2′ξ.
Thus, we have(45)R1μω′μ−R2μξ′μ+R3μ=0,(46)R2μω′μ+R1μξ′μ+R4μ=0.
Forξ′μ at μ=μ∗, therefore dReθj/dμ=χ′μ∗.
By comparing the values ofξ′μ using equations (45) and (46).(47)R1μ∗+R3μ∗R2μ∗=−R2μ∗+R4μ∗R1μ∗.=−R2μ∗R4μ∗+R1μ∗R3μ∗/R12μ∗+R22μ∗≠0, only if R2μ∗R4μ∗+R1μ∗R3μ∗≠0. Thus, the transversality condition is satisfied and hence Hopf bifurcation occurs at μ=μ∗.


## 9. Numerical Simulation

In this section, the numerical analysis is done with the help of Wolfram Mathematica by considering the parameters values of the model cited in Table1. Due to the difficulty in parameter estimation, we relied heavily on previously published work related to this study. The initial values of the variables used to generate graphs in Figures 2–6 are x=3045227, y=217724 [24] in page 10 and since the data for the nonenrolled was unavailable, the study considered the school-going age and further extracted the data from the 2019 national census statistics data. From the numerical analysis, we have found that all the eight equilibrium points are unstable. Since our model lies in the real domain, we will not verify for imaginary equilibrium points.Table 1
Parameter and their description.
ParameterDescriptionValueSourceπSchools recruitment4,583,396γ1Enrolment rate in public schools0.6976[24]γ2Enrolment rate in private schools0.04988[24]k1Public schools carrying capacity3,546,153[24]k2Public schools carrying capacity297,114[24]θ1Transfer rate from public to private schools0.0009Estimatedθ2Transfer rate from private to nonenrolled0.0015[25]θ3Transfer rate from public to nonenrolled0.0015[25]μExit rate due to completion0.25Estimatedη1Private school predation rate on public schools0.0004Estimatedη2Nonenrolled entity predation rate on private schools0.0001[21]η3Nonenrolled entity predation rate on public schools0.0008[21]η4Public school predation rate on private schools0.0006EstimatedBSaturation constant1[26]Figure 2
Dynamics of public schools with respect toθ3 and η3.Figure 3
Dynamics of private school with respect toθ2.Figure 4
Dynamics of nonenrolled with respect toθ3 and η3.Figure 5
Dynamics of population (in millions) against each other at the same time interval. (a) Dynamics in public and private. (b) Dynamics in public and nonenrolled. (c) Dynamics in private and nonenrolled.
(a)(b)(c)Figure 6
Population projection in the long run.For Theorem1, the eigenvalues are 1.0884, −0.0315, and 2.77264; hence, the system is unstableTheorem2, the equilibrium point is −3.2503. Now, we have A1=−13.2227<0, A2=42.9696, and A3=1.2559 and by Routh–Hurwitz, the system is unstable at that equilibrium point.Theorem3, the equilibrium point −0.0425367, the eigenvalues are 2.7726, 0.0315, and 1.0888; hence, Theorem 3 is unstable since it has positive eigenvalues.Theorem4, the eigenvalues are 1.0884, −0.0315, and 2.77264, thus confirms that Theorem 4 is unstable.Theorem5, The equilibrium points are x14=3.2503 and z14=−0.1143 which are unstable since D1=1.7168, D2=−2.9637, and D3=−0.0999<0.x24=212.3727 and z24=930083.112; we have D1=188.24, D2=−205880, D3=224010.Theorem6, the equilibrium points are x15=−3.8359 and y15=−0.0155−4.783i. We have C1=−10.4076−9.6953i, C2=10.1555+100.899i, and C3=−0.005163−98.4232i.x25=10.3364 and y25=−0.0155−7.852i. We have C1=13.772−15.9126i, C2=−16.1924−219.157i, and C3=0.009540+257.578i. From the above-given theory, it confirms that the system is unstable for Theorem 6.Theorem7, the equilibrium points are y16=−7.763×10−13 and z16=1.97959×10−12. By Routh–Hurwitz, we have G1=−3.8296, G2=2.8963, and G3=0.09505. G3<0 implies the system is unstable.y26=−13255.0965 and z26=2.2493×109. We have G1=1.77256×106, G2=−4.8339×1010, and G3=5.26135×1010. Since G1G2−G3<0, the equilibrium points are unstable.Theorem8, the equilibrium points are x17=0,y17=−1656.91,z17=1867770.3. By Routh–Hurwitz criterion, we have K1=−964.05, K2=−3656941.28, and K3=3981178.6. Since K1<0, the equilibrium points are unstable.x27=3.25027,y27=0,z27=−0.004477. We have K1=1.08835, K2=−2.92308, and K3=−0.0987.x37=3.24852,y37=0,z37=7.98893. We have K1=1.71558, K2=−2.957, and K3=−0.1023. x47=3.25029,y47=0.04497,z47=0. We have K1=1.08897, K2=−2.85285, and K3=−0.2985.x57=−13291.8,y57=−0.001216,z57=1.1987×108. We have K1=76202.2, K2=−2.7233×108, and K3=2.9643×108.x67=1.8189×10−12,y67=−13291.8,z67=1.1987×108. We have K1=−10672.8, K2=−2.7233×108, and K3=2.9643×108. Hence, Theorem 8 is unstable.


## 10. Discussion of the Results

Just like disease, predation help in population regulation [15]. But in this study, predation from nonenrolled entity needs to be avoided at all cost to save the schools population by avoiding school dropouts from learners. In this study, the system of ODEs formed had eight equilibrium points. From the numerical simulation, all the eight equilibrium points were unstable. This implies that the system cannot tend to one of the states of having no schools’ population and no nonenrolled entity. In Figure 2, the parameter θ3 and η3 are varied from 0.0015 to 1.3 and 0.0008 to 1, respectively. In both cases as the value of the parameter increase the population of the public school and private school decreases. When the value of η3 is 1, the population of private school becomes extinct. This implies these parameters are the ones responsible for the decrease in schools’ population. In Figure 4 as θ3 and η3 increases the population of the nonenrolled entity increases. This is because θ3 and η3 are the parameters that lead to predation of the schools population as seen in Figure 2. These parameters further cause a decrease in the private population as seen in Figure 3. In Figure 5, public and private schools’ population are almost directly proportional to each other while for public and the nonenrolled, increases in the nonenrolled causes a decrease in the public population. This is due to predation of the nonenrolled as confirmed by Figure 4. The same situation happens in when private schools are plotted against the nonenrolled entity as seen in Figure 5. Figure 6 shows the variations of schools’ population in the long run which was plotted for six years, public and private schools will be directly proportional while for public and the nonenrolled, public school populations increase steadily faster but it reaches a time when the nonenrolled overwhelm the situation and thus public population remains constant. Figure 7 shows in the absence of the nonenrolled the behavior of public and private schools in the long run. Figure 8 show a bifurcation diagram with respect to μ. The bifurcation diagram confirms the theoretical results of Theorem 8. The system of schools and the nonenrolled entity tend to coexistence when μ>0.2 as in Theorem 8 illustrates the coexistence of the system.Figure 7
Public and private population in absence of nonenrolled entity.Figure 8
Bifurcation diagram forμ.


## 11. Conclusion

In this paper, we have studied a model of competition of students’ population growth. The first and second categories of the population are the public schools’ population and private schools’ population, respectively, which grows according to the logistic growth equation. The third category is the nonenrolled entity which grows exponentially. In this work, eight equilibrium points were determined and then studied for the local stability and global stability of the system. The local stability of the equilibrium points has been studied using the eigenvalue method and the Routh–Hurwitz criterion. A numerical simulation has been done to verify the local stability. From the numerical simulation, it was found that all the eight equilibrium points are unstable. However, the system was found to be globally asymptotically stable using the Lyapunov function. Finally, from the results, this study suggests that the government of Kenya through the ministry of education and any concerned body should consider implementing strategies that are aimed at reducing the parametersθ2, θ3, and η3. This will save the schools population from becoming to an extinct. This is important because education is the backbone of a country for the development purposes as the nation requires people who are wise and can help to solve the problems of the country. By employing strategies that mitigate these parameters, it will also help to prevent future occurrences of the same. Bifurcation was carried out for μ and which show coexistence when μ greater than 0.2.In this work, we have developed the Lotka–Volterra model for the system of ODEs to asses to effects of predation, school dropout, and competition among schools with Holling type II response. We used MATLAB software to generate graphs. However, concerning our future studies, we intend to do real data fitting to the model. Another direction that is also of interest is to do optimal control theory and stochastic modelling.



---

*Source: 4102242-2023-01-11.xml*