A cubic nonlinear population growth model for single species: theory, an explicit–implicit solution algorithm and applications

In this paper, we extend existing population growth models and propose a model based on a nonlinear cubic differential equation that reveals itself as a special subclass of Abel differential equations of first kind. We first summarize properties of the time-continuous problem formulation. We state the boundedness, global existence, and uniqueness of solutions for all times. Proofs of these properties are thoroughly given in the Appendix to this paper. Subsequently, we develop an explicit–implicit time-discrete numerical solution algorithm for our time-continuous population growth model and show that many properties of the time-continuous case transfer to our numerical explicit–implicit time-discrete solution scheme. We provide numerical examples to illustrate different behaviors of our proposed model. Furthermore, we compare our explicit–implicit discretization scheme to the classical Eulerian discretization. The latter violates the nonnegativity constraints on population sizes, whereas we prove and illustrate that our explicit–implicit discretization algorithm preserves this constraint. Finally, we describe a parameter estimation approach to apply our algorithm to two different real-world data sets.


Motivation
Differential equations play a vital role in many sciences. Especially, time development of populations within different time frames is of big interest to fields such as biology [13,31,37] and, more specifically, epidemiology [1-4, 9-11, 40, 45, 46] and population dynamics [7, 10, 14-16, 24, 26, 30, 48]. The temporal change of population sizes is not only of scientific interest, but different possibilities of this evolution impact the life of all humans.
Regarding Verhulst's famous logistic model in 1838 [43], we shortly summarize different models often discussed in the literature [20].
1) Let us begin with a simple linear ordinary differential equation for the total population size N 1 for all t ≥ 0. Here N 1 (t) is the first derivative of the population size N 1 with respect to time. Separation of variables under the assumption A > 0 leads to the exponentially increasing solution for all t ≥ 0. For A > 0, we obtain This is why other models should be sought for population dynamics.
2) A quadratic model could be proposed as an alternative. This model reads for the total population size N 2 for all t ≥ 0. The solution is given by for all t < 1 A·N 2,0 . This model includes one blowup at t b = 1 A·N 2,0 , that is, lim t→t b N 2 (t) = +∞.
Since it does not seem realistic that population sizes infinitely increase and since there are only limited materials on Earth, we expect that this model needs further modification.
3) Verhulst proposed his logistic model in 1838, which reads for the total population size N 3 for all t ≥ 0. Mathematica yields N 3 (t) = A · N 3,0 · e A·t A -B · N 3,0 + B · N 3,0 · e A·t as the solution which can be also obtained by a partial fraction decomposition. We obtain lim t→∞ N 3

(t) =
A B for the limiting behavior, where the result corresponds to the carrying capacity. This model has been successfully applied in many situations. However, it does not account for phenomena such as oscillations due to competition between different species [10]. 4) Another possibility might be the application of delay-differential equations such as Hutchinson model N 4 (t) = A · N 4 (t) -B · N 4 (t) · N 4 (tτ ) with time delay τ > 0 and its variants as proposed in [49]. Although this model accounts for various different dynamics depending on the time delay, we might have to deal with the time delay τ > 0 for numerical simulations. 5) Abel differential equations of first kind N 5 (t) = a 1 (t) · N 5 (t) + a 2 (t) · N 2 5 (t) + a 3 (t) · N 3 5 (t) ( 6 ) have been widely investigated [5,8,21]. Even some closed-form solutions for special classes of Abel differential equations of first kind are known [17,28,29]. These equations have been used in many different areas of science [41]. However, these types of differential equations are rarely applied to population dynamics of single species. This explains our motivation for this paper, where we expand Verhulst's logistic model by two modifications. First, we multiply the quadratic term with a time-dependent tuning function and extend Verhulst's logistic model by an additional cubic term, which leads us to a specific class of Abel differential equations of first kind. Our model reads under appropriate conditions to be specified.

Contributions and outline
Based on the aforementioned disadvantages of the presented models, we propose a nonautonomous nonlinear first-order differential equation, which uses an expansion of Verhulst's model by a cubic term and introduces multiplication of the quadratic term by a time-dependent tuning function. Our contributions can be summarized as follows: 1) We state our time-continuous problem formulation in Sect. 2.2; 2) Afterward, we analyze this time-continuous problem formulation and show the nonnegativity and boundedness of possible solutions, global uniqueness in time, and global existence in time based on an inductive usage of Banach's fixed-point theorem in Sect. 2. The proofs can be found in the Appendix; 3) As our last contribution regarding our time-continuous problem formulation, we show a stability bound with respect to initial conditions and data. Hence, we conclude that our problem formulation is well posed on compact time intervals [0, T] for T > 0. The proof of this statement is thoroughly described in the Appendix; 4) As our first major contribution, we introduce an explicit-implicit numerical solution algorithm in Sect. 3 and prove that many properties of the time-continuous case transfer to its time-discrete variant. Additionally, we thoroughly prove that it converges linearly toward the solution of our time-continuous problem; 5) Finally, as our second and final major contribution, we first provide one numerical example to compare the classical explicit Eulerian discretization scheme to our proposed explicit-implicit numerical solution algorithm. In addition, we give further numerical examples for different behaviors of solutions in Sect. 4. Afterward, we describe a parameter estimation approach and conclusively present two examples on real-world data with user-chosen time-dependent functions to illustrate usefulness of our proposed model with respect to real-world applications in Sect. 5. In Sect. 6, we sum up our results and look out on possible future research directions.

Time-continuous problem formulation
In this section, we first sum up some mathematical background material for our analysis. After that, we describe our time-continuous problem formulation for population dynamics in one species. Conclusively, we analyze certain properties of our proposed model in detail. Detailed proofs of these statements can be found in the Appendix.

Mathematical background material
To especially state the global existence and global uniqueness of the solution of our nonautonomous first-order nonlinear population model, we need to introduce some theoretical background material regarding nonlinear ordinary differential equations. This subsection heavily relies on the structure of [45], which is based on [34,35,39]. Let us first recall the Lipschitz continuity of a function on Euclidean spaces.
for all x, y ∈ S. Here, · denotes a suitable norm on the corresponding Euclidean space.
Let U ⊂ R m be open, and let F : U − → R n . We call F locally Lipschitz continuous if for every point x 0 ∈ U, there exists a neighborhood V of x 0 such that the restriction of F to V is Lipschitz continuous on V .
We consider the initial-value problem where z(t) = (x 1 (t), . . . , x n (t)) denotes our solution vector. Our vectorial function is represented by G(t, z(t)) = (g 1 (t, z(t)), . . . , g n (t, z(t))), and z 0 ∈ R n are given initial conditions. To conclude the global existence, we can apply the following theorem, which is a consequence of Grönwall's lemma.
for all z(t) ∈ R n , then the solution of the initial value problem (8) exists for all t ∈ [0, ∞). Moreover, for every finite T ≥ 0, we have Additionally, we need Banach's fixed point theorem to derive the global uniqueness.
. Let T : X − → X be a strict contraction, that is, there exists a constant K ∈ [0, 1) such that (Tx, Ty) ≤ K · (x, y) for all x, y ∈ X. Then the map T has a unique fixed point.
Finally, we need a modification of Grönwall's lemma for continuous dependence on perturbations of initial values and data. This might be also known to readers by stability analysis.
for all t ∈ I, then for all t ∈ I.

Formulation
First, we state the following assumptions on the model: is a bounded function, that is, there exist positive constants m f and M f such that m f ≤ f (t) ≤ M f for all t ≥ 0. Hence our nonlinear initial value problem for our population model reads for t ≥ 0. In the rest of this section, we want to prove some important properties of our model.

Boundedness
As one important consequence, we obtain the boundedness and positivity of solutions to the initial value problem (13).
Lemma 2.5 Solutions of our model (13) are bounded, that is, there exist constants N min and N max such that N min ≤ N(t) ≤ N max for all t ≥ 0. In particular, solutions to positive initial values remain positive for all t ≥ 0.

Global existence
In addition to the boundedness, we conclude the existence of solutions to the initial value problem (13) for all t ≥ 0.

Theorem 2.6
Solutions to our population growth model (13) exist for all t ≥ 0.

Global uniqueness
As already stated, our proof of global uniqueness in time is heavily based on Banach's fixed-point theorem, which has proven a valuable tool in different mathematical problem settings [44,45,50].

Theorem 2.7
The initial value problem (13) possesses a unique solution for all t ≥ 0.

Continuous dependence on initial conditions and data
Let us first provide a simple stability bound, where only perturbations of initial conditions are considered. In the following, · ∞ denotes the maximum norm on a finite-dimensional Euclidean space.
for all t ∈ [0, T], where Now, we want to extend this stability by additionally assuming perturbations with respect to our prescribed data.
is a time-dependent coefficient, and

Explicit-implicit time-discrete problem formulation
In this section, we propose an explicit-implicit time-discrete solution algorithm for (13). We show unique solvability and demonstrate that many properties of the time-continuous case transfer to our time-discrete scheme.

Time-discrete problem formulation
with n+1 = t n+1t n for all n ∈ {1, . . . , M -1} and given initial population size N 0 . Here f n+1 is an abbreviation for f n+1 = f (t n+1 ). The term of explicit-implicit solution algorithm refers to the right-hand-side of (16). Whereas the first two summands are treated explicitly, the last summand is a mixture of N n and N n+1 . Hence, the proposed solution algorithm lies somewhere between fully explicit and fully implicit numerical solution schemes. Thus, these algorithms are called explicit-implicit solution algorithms [18,45].

Solvability
Now, we can prove that our time-discrete model (16) is uniquely solvable. (16) is uniquely solvable.

Theorem 3.1 The time-discrete explicit-implicit population growth model
Proof We first observe from (16) that Rearranging yields This implies We conclude unique solvability for all n ∈ {1, . . . , M -1}. (16) is bounded, that is, there exist constants N num,min and N num,max such that

Theorem 3.2 The unique solution of our time-discrete population growth model
for all n ∈ {1, . . . , M}.
Proof 1) We first show that N num,min > 0. By (17), we observe by induction that Hence we obtain N num,min > 0. This means that our time-discrete solution is nonnegative for all n ∈ {1, . . . , M}.
Additionally, we see that our solution sequence is monotonically decreasing if and only if This yields Square addition as in the time-continuous case results in From (20) and the nonnegativity, we see that our time-discrete solution is monotonically decreasing if and only if as in the time-continuous case. Since the function f is bounded below by zero, this implies for our lower bound.
2) Similarly, we can conclude for our upper bound. This provides both bounds, and the proof is finished.
The result of Theorem 3.2 implies that our time-discrete population growth model (16) is unconditionally stable and preserves the nonnegativity of population sizes N . Proof Since this proof is relatively technical, we briefly describe our strategy. We start with consideration of local errors between time-continuous and time-discrete solutions. Afterward, we have to take into account that errors propagate in time. Finally, we need to accumulate these errors.

Convergence
1) For our investigation of local errors, we assume that (t p , N p ) = (t p , N(t p )) and consider the time interval [t p , t p+1 ] for arbitrary p ∈ {1, . . . , M -1}. Here, we only consider one time step and denote the time-discrete solution at t p+1 by N p+1 . By (17), we have Application of the triangle inequality yields after some manipulations. By the mean value theorem of calculus, there are τ 1 , These equations yield where we set We finally obtain for local errors.
2) In reality, (t p , N p ) does not lie on the time-continuous solution exactly. For that reason, we must investigate how the procedural error N p -N(t p ) propagates to the (p + 1)th time step. These estimates are carried out in this and the following steps.
We have to consider |N p+1 -N p+1 |. We observe that By the triangle inequality, we obtain and p+1 ≤ 1, we can conclude that Let us define This yields for error propagation in time from the pth time step to the (p + 1)th time step.
3.1) First, we want to prove inductively that for all p ∈ {0, . . . , M -1}. Let p = 0. Inequality (26) is obviously fulfilled. Let p = 1 to understand the concept. By application of the triangle inequality and inequalities (24) and (25) we observe that We now assume that We now want to show (26). We obtain and this finishes our inductive proof. 3.2) We define := max p∈N p . We infer that Table 1 Algorithmic summary of our explicit-implicit time-discretion population growth model If the initial conditions of our time-continuous and time-discrete models coincide and → 0, then the time-discrete solution converges linearly toward the time-continuous solution of the time-continuous population growth model. Finally, our statement is proven.

Algorithmic summary
We summarize our algorithmic approach in Table 1.

Numerical examples
Here, we present four examples. In our first example, we compare the classical explicit Eulerian discretization to our proposed explicit-implicit numerical solution algorithm.
Three synthetic examples show that different behaviors are possible.

Example 1: comparison of time discretization algorithms
First, we compare the classical explicit Eulerian time discretization algorithm to our explicit-implicit numerical solution algorithm • Function f (t) ≡ 1. In Fig. 1, we can clearly see that the classical Eulerian time discretization scheme becomes unstable and we obtain negative population sizes, whereas our proposed explicitimplicit numerical solution algorithm remains stable and provides positive population sizes for all times. It is well known in numerical analysis that explicit time-stepping methods are always conditionally stable [22,23,25]. Hence, we have to carefully select possible time-discretization schemes with respect to equations we want to solve.

Example 2: monotonically increasing solution
We present an example that is monotonically increasing. We consider the following inputs: •  Fig. 2, we see that we can get results that look similar to monotonically increasing solutions from Verhulst's logistic model.

Example 3: solution with changing sign in derivative
We give an example where the derivative of our population model changes its sign in time.
For that purpose, we propose the following input parameters: •   We observe from Fig. 3 that it is possible to obtain solutions that first monotonically increase and later monotonically decrease in time.

Example 4: oscillating solution
We want to construct an oscillating example where the following parameters are taken into account: •

Parameter estimation and real-world applications
For real-world applications, we need to present a parameter estimation approach. We first describe our method and apply it to two real-world data sets from Human World and Japanese populations. From our observations, we especially derive our time-dependent coefficient functions we later apply to both data sets.

Parameter estimation
We consider the proposed explicit-implicit time discretization method where we set with respect to the nonnegativity constraints on the sought parameters, where · 2 is the Euclidean norm. We apply the function lsqnonneg of GNU Octave to solve this nonnegative linear least-squares problem [19]. For further details on numerical optimization, we refer the interested readers to [33].

Example 5: human world population
The data are taken from the United Nations [42] for the human total population size including both sexes. We investigate three different settings where parameters and given functions are user-chosen or, in one case, are estimated. For more sophisticated parameter estimation techniques, which are outside the scope of our paper, we refer the interested readers to further works [1,11,12,47], as these techniques might be an interesting extension of this paper.
The following parameters are applied in our first user-chosen setting: In our third setting, we apply the following parameters, where A 3 and C 3 were estimated by our least-squares approach, and B 3 is the arithmetic mean of the depicted data points in Fig. 5: • A 3 = 0.018758, B 3 = 0.000632, and C 3 = 0.000154; • Time step size = 1.0 year; • Initial population size N 1 = 2.536431 billion people; • Function f 3 (t) = 1. All chosen functions g j (t) = B j · f j (t) for j ∈ {1, 2, 3} are shown in Fig. 5.
The results of all three settings are depicted in Fig. 6. We clearly observe that it is possible to construct fits with different monotonicity assumption regarding future predictions.
The errors between known real-world data and our portrayed settings are portrayed in Fig. 7. All graphs well describe the given data. However, as observed in Fig. 6, all three predictions provide different long-time behaviors.

Example 6: Japanese population
Data for our last example are again taken from the United Nations [42] for human total population size including both sexes in Japan. The following inputs are considered for this example, where a simplified function f is user-chosen by the results of our least-squares estimation approach as depicted in Fig. 8: •  As already mentioned, our estimation results for our time-dependent coefficient function can be seen in Fig. 8.
The results of this setting are portrayed in Fig. 9, which depicts a rapidly decaying population size in Japan. In contrast to the traditional Verhulst model, which only accounts for monotone behavior, our model is able to capture different intervals of monotonicity in our given data.

Conclusion and outlook
In this work, we proposed an extended population growth model based on a nonlinear first-order differential equation with quadratic and cubic terms of population size. At the beginning, we proved the nonnegativity and boundedness globally in time. After that, we established the global existence and uniqueness of the solution of our model in time. Subsequently, we discretized our model and provided an explicit-implicit solution algorithm. We showed that many properties of our time-continuous model transfer to our time-discrete variant. Finally, we applied our results to four synthetic and two real-world examples to demonstrate usefulness and variability of our presented model. As a conclusion, we saw that explicit time-stepping methods suffer from time restrictions and may violate the nonnegativity constraints on populations size, whereas our proposed explicitimplicit numerical solution algorithm preserves the nonnegativity.
Our model seems to be flexible. This flexibility might be extended in future work by multiplication of all terms on the right-hand side by time-dependent tuning functions. Compare, for example, [27,32,36,38]. Additionally, parameter estimation for models in mathematical biology could be regarded as further work because its manual choosing is Figure 9 Results for our prediction of the Japanese population size. Here t = 0 corresponds to the year 1950, and t = 69 corresponds to the year 2019. This is the time frame for our data points a delicate task [1,11,12,47]. However, more sophisticated techniques for solving inverse problems are out of the scope of this paper.
As a final extension, our model can be modified by application of fractional differential operators [6].

Proof of Lemma 2.5
First, we see that for all t ≥ 0. Pick arbitrary t 0 ∈ [0, ∞) such that Then Hence, we choose N min := min{N 0 , A C } > 0 as a lower bound for the total population size. Second, we further notice that Examining N (t) = 0 from (32) yields for arbitrary t ≥ 0 where N (t) = 0. Thus we can choose the constant N max := max{N 0 , } as an upper bound for the total population size due to the fact that Finally, the two previous steps of our proof imply that the solutions of our population growth model (13) stay bounded and positive.

Proof of Theorem 2.6
We apply Theorem 2.2. Applying the maximum norm and triangle inequality and using the boundedness of the function f and the solution by Lemma 2.5, we obtain and Theorem 2.2 proves our statement.

Proof of Theorem 2.7
Let t ∈ [0, τ ] be arbitrary for a chosen τ > 0. Let us assume that both functions N : [0, ∞) − → [0, ∞) and N : [0, ∞) − → [0, ∞) solve our population growth model (13). We see that By the boundedness of the solution from Lemma 2.5, we obtain for all t ∈ [0, τ ]. This yields the uniqueness on the time interval [0, τ ] by Banach's fixedpoint theorem. Inductively, we can apply a similar argument as in our first step to obtain the uniqueness on time intervals [k · τ , (k + 1) · τ ] for all k ∈ N 0 . Here t 0 = k · τ is our new initial time point. This completes the proof of global uniqueness in time.

Proof of Theorem 2.8
We first observe that B · M f · N a (τ ) + N b (τ ) · N a (τ ) -N a (τ ) dτ for all t ∈ [0, T]. Defining α 1 and α 2 as stated in our claim of Theorem 2.8, we obtain and hence, as both α 1 and α 2 are constant in time, (14) holds after application of Theorem 2.4, which proves our statement.
Again, we observe that N a (t) = N a,0 + t 0 A a · N a (τ ) + B a · f a (τ ) · N 2 a (τ ) -C a · N 3 a (τ ) dτ and We will estimate all terms I, II, and III.
By the inequality of our first step, we conclude that By multiple application of the inequality of our first step, we infer that ≤ B a · max{N a,max ; N b,max } 2 · t · f a (t)f b (t) ∞ + max{M f ,a ; M f ,b } · max{N a,max ; N b,max } 2 · t · |B a -B b | + t 0 2 · B a · max{M f ,a ; M f ,b } · max{N a,max ; N b,max } · N a (τ ) -N b (τ ) dτ .