1 Introduction

1.1 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 [14, 911, 40, 45, 46] and population dynamics [7, 10, 1416, 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

$$ N_{1}^{\prime } ( t ) = A \cdot N_{1} ( t ), \qquad N_{1} ( 0 ) = N_{1, 0} > 0, $$
(1)

for the total population size \(N_{1}\) for all \(t \geq 0\). Here \(N_{1}^{\prime } ( 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

$$ N_{1} ( t ) = N_{1, 0} \cdot \exp ( A \cdot t ) $$
(2)

for all \(t \geq 0\). For \(A > 0\), we obtain

$$ \lim_{t \to \infty } N_{1} ( t ) = + \infty . $$

This is why other models should be sought for population dynamics.

2) A quadratic model could be proposed as an alternative. This model reads

$$ N_{2}^{\prime } ( t ) = A \cdot N_{2}^{2} ( t ), \qquad N_{2} ( 0 ) = N_{2, 0} > 0, $$
(3)

for the total population size \(N_{2}\) for all \(t \geq 0\). The solution is given by

$$ N_{2} ( t ) = \frac{N_{2, 0}}{1 - A \cdot N_{2, 0} \cdot t} $$
(4)

for all \(t < \frac{1}{A \cdot N_{2, 0}}\). This model includes one blowup at \(t_{b} = \frac{1}{A \cdot N_{2, 0}}\), that is,

$$ \lim_{t \to t_{b}} N_{2} ( t ) = + \infty . $$

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

$$ N_{3}^{\prime } ( t ) = A \cdot N_{3} ( t ) - B \cdot N_{3}^{2} ( t ), \qquad N_{3} ( 0 ) = N_{3, 0} > 0, $$
(5)

for the total population size \(N_{3}\) for all \(t \geq 0\). Mathematica yields

$$ N_{3} ( t ) = \frac{A \cdot N_{3, 0} \cdot \mathrm{e}^{A \cdot t}}{A - B \cdot N_{3, 0} + B \cdot N_{3, 0} \cdot \mathrm{e}^{A \cdot t}} $$

as the solution which can be also obtained by a partial fraction decomposition. We obtain

$$ \lim_{t \to \infty } N_{3} ( t ) = \frac{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}^{\prime } ( t ) = A \cdot N_{4} ( t ) - B \cdot N_{4} ( t ) \cdot N_{4} ( t - \tau ) $$

with time delay \(\tau > 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 \(\tau > 0\) for numerical simulations.

5) Abel differential equations of first kind

$$ N_{5}^{\prime } ( t ) = a_{1} ( t ) \cdot N_{5} ( t ) + a_{2} ( t ) \cdot N_{5}^{2} ( t ) + a_{3} ( t ) \cdot N_{5}^{3} ( 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

$$ N^{\prime } ( t ) = A \cdot N ( t ) + B \cdot f ( t ) \cdot N^{2} ( t ) - C \cdot N^{3} ( t ) $$

under appropriate conditions to be specified.

1.2 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. 1)

    We state our time-continuous problem formulation in Sect. 2.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. 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. 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. 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.

2 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.

2.1 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.

Definition 2.1

([39, Sect. 3.2])

Let \(m, n \in \mathbb{N}\). If \(S \subset \mathbb{R}^{m}\), a function \(\mathbf{F} : S \longrightarrow \mathbb{R}^{n}\) is called Lipschitz continuous on S if there exists a constant \(L \geq 0\) such that

$$ \bigl\lVert \mathbf{F} ( \mathbf{x} ) - \mathbf{F} ( \mathbf{y} ) \bigr\rVert _{\mathbb{R}^{n}} \leq L \cdot \lVert \mathbf{x} - \mathbf{y} \rVert _{\mathbb{R}^{m}} $$
(7)

for all \(\mathbf{x}, \mathbf{y} \in S\). Here, \(\lVert \cdot \rVert \) denotes a suitable norm on the corresponding Euclidean space.

Let \(U \subset \mathbb{R}^{m}\) be open, and let \(\mathbf{F} \colon U \longrightarrow \mathbb{R}^{n}\). We call F locally Lipschitz continuous if for every point \(\mathbf{x_{0}} \in U\), there exists a neighborhood V of \(\mathbf{x_{0}}\) such that the restriction of F to V is Lipschitz continuous on V.

We consider the initial-value problem

$$ \left \{ \begin{aligned} \mathbf{z}^{\prime } ( t ) & = \mathbf{G} \bigl( t, \mathbf{z} ( t ) \bigr), \\ \mathbf{z} ( 0 ) & = \mathbf{z_{0}}, \end{aligned} \right . $$
(8)

where \(\mathbf{z} ( t ) = ( x_{1} ( t ), \ldots , x_{n} ( t ) )\) denotes our solution vector. Our vectorial function is represented by \(\mathbf{G} ( t, \mathbf{z} ( t ) ) = ( g_{1} ( t, \mathbf{z} ( t ) ), \ldots , g_{n} ( t, \mathbf{z} ( t ) ) )\), and \(\mathbf{z_{0}} \in \mathbb{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.

Theorem 2.2

([39, Theorem 4.7.1])

If \(\mathbf{G} \colon [ 0, \infty ) \times \mathbb{R}^{n} \longrightarrow \mathbb{R}^{n}\) is locally Lipschitz continuous in both time and state variables and if there exist nonnegative real functions \(D \colon [ 0, \infty ) \longrightarrow [ 0, \infty )\) and \(K \colon [ 0, \infty ) \longrightarrow [ 0, \infty )\) such that

$$ \bigl\lVert \mathbf{G} \bigl( t, \mathbf{z} ( t ) \bigr) \bigr\rVert _{\mathbb{R}^{n}} \leq K ( t ) \cdot \bigl\lVert \mathbf{z} ( t ) \bigr\rVert _{\mathbb{R}^{n}} + D ( t ) $$
(9)

for all \(\mathbf{z} ( t ) \in \mathbb{R}^{n}\), then the solution of the initial value problem (8) exists for all \(t \in [ 0, \infty )\). Moreover, for every finite \(T \geq 0\), we have

$$ \bigl\lVert \mathbf{z} ( t ) \bigr\rVert _{\mathbb{R}^{n}} \leq \lVert \mathbf{z_{0}} \rVert _{\mathbb{R}^{n}} \cdot \exp \bigl( K_{\max} \cdot \vert t \vert \bigr) + \frac{D_{\max}}{K_{\max}} \cdot \bigl( \exp \bigl( K_{\max} \cdot \vert t \vert \bigr) - 1 \bigr) $$
(10)

for all \(t \in [ 0, T ]\), where

$$ D_{\max} = \max_{0 \leq s \leq T} \bigl\vert D ( s ) \bigr\vert \quad \textit{and} \quad K_{\max} = \max_{0 \leq s \leq T} \bigl\vert K ( s ) \bigr\vert . $$

Additionally, we need Banach’s fixed point theorem to derive the global uniqueness.

Theorem 2.3

([35, Theorem V.18])

Let \(( X, \varrho )\) be a complete metric space with metric \(\varrho \colon X \times X \longrightarrow [ 0, \infty )\). Let \(T \colon X \longrightarrow X\) be a strict contraction, that is, there exists a constant \(K \in [ 0, 1 )\) such that \(\varrho ( Tx, Ty ) \leq K \cdot \varrho ( x, y )\) for all \(x, y \in 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.

Theorem 2.4

([34, Theorem 1.3.1])

Let \(I := [ a, b ]\). Let \(u, f \colon I \longrightarrow [ 0, \infty )\) be continuous nonnegative functions. Let \(g \colon I \longrightarrow ( 0, \infty )\) be a continuous positive nondecreasing function. If

$$ u ( t ) \leq g ( t ) + \int _{a}^{t} f ( s ) \cdot u ( s ) \,\mathrm{d}s $$
(11)

for all \(t \in I\), then

$$ u ( t ) \leq g ( t ) \cdot \exp \biggl( \int _{a}^{t} f ( s ) \,\mathrm{d}s \biggr) $$
(12)

for all \(t \in I\).

2.2 Formulation

First, we state the following assumptions on the model:

  • \(A > 0\), \(B \geq 0\), \(C > 0\);

  • \(f \colon [ 0, \infty ) \longrightarrow ( 0, \infty )\) is a bounded function, that is, there exist positive constants \(m_{f}\) and \(M_{f}\) such that \(m_{f} \leq f ( t ) \leq M_{f}\) for all \(t \geq 0\).

Hence our nonlinear initial value problem for our population model reads

$$ N^{\prime } ( t ) = A \cdot N ( t ) + B \cdot f ( t ) \cdot N^{2} ( t ) - C \cdot N^{3} ( t ), \quad N ( 0 ) = N_{0} > 0, $$
(13)

for \(t \geq 0\). In the rest of this section, we want to prove some important properties of our model.

2.3 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 } \leq N ( t ) \leq N_{\max}\) for all \(t \geq 0\). In particular, solutions to positive initial values remain positive for all \(t \geq 0\).

2.4 Global existence

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

Theorem 2.6

Solutions to our population growth model (13) exist for all \(t \geq 0\).

2.5 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 \geq 0\).

2.6 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, \(\lVert \cdot \rVert _{\infty }\) denotes the maximum norm on a finite-dimensional Euclidean space.

Theorem 2.8

Let \(N_{a} \colon [ 0, T ] \longrightarrow [ 0, \infty )\) be a solution of (13) with initial value \(N_{a} ( 0 ) = N_{a, 0} > 0\), and let \(N_{b} \colon [ 0, T ] \longrightarrow [ 0, \infty )\) be solution of (13) with initial value \(N_{b} ( 0 ) = N_{b, 0} > 0\) with finite time \(T > 0\). We have the stability bound

$$ \bigl\lVert N_{a} ( t ) - N_{b} ( t ) \bigr\rVert _{ \infty } \leq \alpha _{1} \cdot \exp ( \alpha _{2} \cdot t ) $$
(14)

for all \(t \in [ 0, T ]\), where

$$ \alpha _{1} := \vert N_{a, 0} - N_{b, 0} \vert $$

and

$$ \alpha _{2} := A + 2 \cdot B \cdot M_{f} \cdot \bigl( \max \{ N_{a, \max}; N_{b, \max} \} \bigr) + 3 \cdot C \cdot \bigl( \max \{ N_{a, \max}; N_{b, \max} \} \bigr)^{2}. $$

Now, we want to extend this stability by additionally assuming perturbations with respect to our prescribed data.

Theorem 2.9

Let \(A_{a}, A_{b} > 0\), \(B_{a}, B_{b} \geq 0\), \(C_{a}, C_{b} > 0\), and let \(f_{a}, f_{b} \colon [ 0, T ] \longrightarrow ( 0, \infty )\) be continuous bounded functions, that is, there exist positive constants \(m_{f, a}\), \(m_{f, b}\), \(M_{f, a}\), and \(M_{f, b}\) such that \(m_{f, a} \leq f_{a} ( t ) \leq M_{f, a}\) and \(m_{f, b} \leq f_{b} ( t ) \leq M_{f, b}\) for all \(t \in [ 0, T ]\). Additionally, assume that \(N_{a} \colon [ 0, T ] \longrightarrow [ 0, \infty )\) solves

$$ N_{a}^{\prime } ( t ) = A_{a} \cdot N_{a} ( t ) + B_{a} \cdot f_{a} ( t ) \cdot N_{a}^{2} ( t ) - C_{a} \cdot N_{a}^{3} ( t ) ; \quad N_{a} ( 0 ) = N_{a, 0} > 0 $$

for all \(t \in [ 0, T ]\) and that \(N_{b} \colon [ 0, T ] \longrightarrow [ 0, \infty )\) solves

$$ N_{b}^{\prime } ( t ) = A_{b} \cdot N_{b} ( t ) + B_{b} \cdot f_{b} ( t ) \cdot N_{b}^{2} ( t ) - C_{b} \cdot N_{b}^{3} ( t ) ; \quad N_{b} ( 0 ) = N_{b, 0} > 0 $$

for all \(t \in [ 0, T ]\). Then we have the stability bound

$$ \bigl\lVert N_{a} ( t ) - N_{b} ( t ) \bigr\rVert _{ \infty } \leq \beta _{1} ( t ) \cdot \exp ( \beta _{2} \cdot t ) $$
(15)

for all \(t \in [ 0, T ]\), where

$$\begin{aligned} \beta _{1} ( t ) := & \vert N_{a, 0} - N_{b, 0} \vert + \bigl( \max \{ N_{a, \max}; N_{b, \max} \} \bigr) \cdot t \cdot \vert A_{a} - A_{b} \vert \\ &{} + B_{a} \cdot \bigl( \max \{ N_{a, \max}; N_{b, \max} \} \bigr)^{2} \cdot t \cdot \bigl\lVert f_{a} ( t ) - f_{b} ( t ) \bigr\rVert _{\infty } \\ &{} + \max \{ M_{f, a}; M_{f, b} \} \cdot \bigl( \max \{ N_{a, \max}; N_{b, \max} \} \bigr)^{2} \cdot t \cdot \vert B_{a} - B_{b} \vert \\ &{} + \bigl( \max \{ N_{a, \max}; N_{b, \max} \} \bigr)^{3} \cdot t \cdot \vert C_{a} - C_{b} \vert \end{aligned}$$

is a time-dependent coefficient, and

$$\begin{aligned} \beta _{2} := & \max \{ A_{a}; A_{b} \} + 2 \cdot B_{a} \cdot \max \{ M_{f, a}; M_{f, b} \} \cdot \max \{ N_{a, \max}; N_{b, \max} \} \\ &{} + 3 \cdot C_{a} \cdot \max \bigl( \{ N_{a, \max}; N_{b, \max} \} \bigr)^{2}. \end{aligned}$$

3 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.

3.1 Time-discrete problem formulation

Let \([ 0, T ]\) be the time interval for model simulations. Let \(\{ t_{j} \} _{j = 1}^{M}\) be a strictly increasing sequence of time points, that is, \(t_{1} = 0 < t_{2} < \cdots < t_{M - 1} < t_{M} = T\). Our time-discrete problem formulation reads

$$ \frac{N_{n + 1} - N_{n}}{\Delta _{n + 1}} = A \cdot N_{n} + B \cdot f_{n + 1} \cdot N_{n}^{2} - C \cdot N_{n + 1} \cdot N_{n}^{2}, \quad N_{1} = N_{0} > 0, $$
(16)

with \(\Delta _{n + 1} = t_{n + 1} - t_{n}\) for all \(n \in \{ 1, \ldots , 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].

3.2 Solvability

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

Theorem 3.1

The time-discrete explicit–implicit population growth model (16) is uniquely solvable.

Proof

We first observe from (16) that

$$ N_{n + 1} = N_{n} + A \cdot \Delta _{n + 1} \cdot N_{n} + B \cdot \Delta _{n + 1} \cdot f_{n + 1} \cdot N_{n}^{2} - C \cdot \Delta _{n + 1} \cdot N_{n}^{2} \cdot N_{n + 1}. $$

Rearranging yields

$$ N_{n + 1} \cdot \bigl\{ 1 + C \cdot \Delta _{n + 1} \cdot N_{n}^{2} \bigr\} = N_{n} \cdot \{ 1 + A \cdot \Delta _{n + 1} + B \cdot \Delta _{n + 1} \cdot f_{n + 1} \cdot N_{n} \} . $$

This implies

$$ N_{n + 1} = N_{n} \cdot \frac{1 + A \cdot \Delta _{n + 1} + B \cdot \Delta _{n + 1} \cdot f_{n + 1} \cdot N_{n}}{1 + C \cdot \Delta _{n + 1} \cdot N_{n}^{2}}. $$
(17)

We conclude unique solvability for all \(n \in \{ 1, \ldots , M - 1 \} \). □

3.3 Boundedness

Theorem 3.2

The unique solution of our time-discrete population growth model (16) is bounded, that is, there exist constants \(N_{\mathrm{num},\min}\) and \(N_{\mathrm{num},\max}\) such that

$$ N_{\mathrm{num},\min} \leq N_{n} \leq N_{\mathrm{num},\max} $$
(18)

for all \(n \in \{ 1, \ldots , M \} \).

Proof

1) We first show that \(N_{\mathrm{num},\min} > 0\). By (17), we observe by induction that

$$ N_{n + 1} = N_{n} \cdot \frac{1 + A \cdot \Delta _{n + 1} + B \cdot \Delta _{n + 1} \cdot f_{n + 1} \cdot N_{n}}{1 + C \cdot \Delta _{n + 1} \cdot N_{n}^{2}} > 0 $$

for all \(n \in \{ 1, \ldots , M - 1 \} \) because of \(N_{1} = N_{0} > 0\). Hence we obtain \(N_{\mathrm{num},\min} > 0\). This means that our time-discrete solution is nonnegative for all \(n \in \{ 1, \ldots , M \} \).

Additionally, we see that our solution sequence is monotonically decreasing if and only if

$$ \frac{1 + A \cdot \Delta _{n + 1} + B \cdot \Delta _{n + 1} \cdot f_{n + 1} \cdot N_{n}}{1 + C \cdot \Delta _{n + 1} \cdot N_{n}^{2}} \leq 1. $$
(19)

This yields

$$ A + B \cdot f_{n + 1} \cdot N_{n} \leq C \cdot N_{n}^{2}. $$

Square addition as in the time-continuous case results in

$$ \biggl( N_{n} - \frac{B \cdot f_{n + 1}}{2 \cdot C} \biggr)^{2} \geq \frac{4 \cdot A \cdot C + B^{2} \cdot f_{n + 1}^{2}}{4 \cdot C^{2}}. $$
(20)

From (20) and the nonnegativity, we see that our time-discrete solution is monotonically decreasing if and only if

$$ N_{n} \geq \frac{B \cdot f_{n + 1}}{2 \cdot C} + \sqrt{ \frac{4 \cdot A \cdot C + B^{2} \cdot f_{n + 1}^{2}}{4 \cdot C^{2}}} $$
(21)

as in the time-continuous case. Since the function f is bounded below by zero, this implies

$$ N_{\mathrm{num},\min} = \min \biggl\{ N_{1}, \sqrt{ \frac{A}{C}}, \frac{B \cdot m_{f}}{2 \cdot C} \biggr\} $$
(22)

for our lower bound.

2) Similarly, we can conclude

$$ N_{\mathrm{num},\max} = \max \biggl\{ N_{1}, \frac{B \cdot M_{f}}{2 \cdot C} + \sqrt{ \frac{4 \cdot A \cdot C + B^{2} \cdot M_{f}^{2}}{4 \cdot C^{2}}} \biggr\} $$
(23)

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.

3.4 Convergence

Theorem 3.3

Let us assume that the solution of our time-continuous population growth model (13) is twice continuously differentiable and that the function f is continuously differentiable. Additionally, we assume that \(\Delta _{p} \leq 1\) for all \(p \in \mathbb{N}\) and \(\lim_{p \to \infty } \Delta _{p + 1} = 0\). This implies that the solution of our time-discrete population growth model (16) converges linearly toward the solution of the time-continuous population growth model on a time interval \([ 0, T ]\) for arbitrary \(T > 0\).

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.

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 \in \{ 1, \ldots , M - 1 \} \). Here, we only consider one time step and denote the time-discrete solution at \(t_{p + 1}\) by \(\widetilde{N_{p + 1}}\). By (17), we have

$$ N_{p + 1} = N ( t_{p} ) \cdot \frac{1 + A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N ( t_{p} )}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}}. $$

Application of the triangle inequality yields

$$\begin{aligned}& \bigl\vert N ( t_{p + 1} ) - \widetilde{N_{p + 1}} \bigr\vert \\& \quad = \biggl\vert N ( t_{p + 1} ) - N ( t_{p} ) \cdot \frac{1 + A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N ( t_{p} )}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad = \biggl\vert N ( t_{p + 1} ) - N ( t_{p} ) \cdot \frac{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} + A \cdot \Delta _{p + 1}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \\& \qquad {} - N ( t_{p} ) \cdot \frac{B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N ( t_{p} ) - C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad = \biggl\vert N ( t_{p + 1} ) - N ( t_{p} ) - N ( t_{p} ) \cdot \frac{A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N ( t_{p} ) - C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad = \biggl\vert \int _{t_{p}}^{t_{p + 1}} N^{\prime } ( \tau ) \,\mathrm{d}\tau - \Delta _{p + 1} \cdot \frac{A \cdot N ( t_{p} ) + B \cdot f_{p + 1} \cdot ( N ( t_{p} ) )^{2} - C \cdot ( N ( t_{p} ) )^{3}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad = \biggl\vert \int _{t_{p}}^{t_{p + 1}} N^{\prime } ( \tau ) \,\mathrm{d}\tau - \Delta _{p + 1} \cdot \frac{A \cdot N ( t_{p} ) + B \cdot ( f_{p} + f_{p + 1} - f_{p} ) \cdot ( N ( t_{p} ) )^{2}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \\& \qquad {}- \Delta _{p + 1} \cdot \frac{- C \cdot ( N ( t_{p} ) )^{3}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad = \biggl\vert \int _{t_{p}}^{t_{p + 1}} N^{\prime } ( \tau ) \,\mathrm{d}\tau - \Delta _{p + 1} \cdot \frac{A \cdot N ( t_{p} ) + B \cdot f_{p} \cdot ( N ( t_{p} ) )^{2} - C \cdot ( N ( t_{p} ) )^{3}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \\& \qquad {}- \Delta _{p + 1} \cdot \frac{B \cdot ( f_{p + 1} - f_{p} ) \cdot ( N ( t_{p} ) )^{2}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad = \biggl\vert \int _{t_{p}}^{t_{p + 1}} N^{\prime } ( \tau ) \,\mathrm{d}\tau - \Delta _{p + 1} \cdot \frac{N^{\prime } ( t_{p} )}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \\& \qquad {}- \Delta _{p + 1} \cdot \frac{B \cdot ( f_{p + 1} - f_{p} ) \cdot ( N ( t_{p} ) )^{2}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad = \biggl\vert \int _{t_{p}}^{t_{p + 1}} N^{\prime } ( \tau ) \,\mathrm{d}\tau - \Delta _{p + 1} \cdot \frac{N^{\prime } ( t_{p} ) + B \cdot ( f_{p + 1} - f_{p} ) \cdot ( N ( t_{p} ) )^{2}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad = \biggl\vert \int _{t_{p}}^{t_{p + 1}} N^{\prime } ( \tau ) \,\mathrm{d}\tau - \Delta _{p + 1} \cdot N^{\prime } ( t_{p} ) \\& \qquad {} - \Delta _{p + 1} \cdot \frac{B ( f_{p + 1} - f_{p} ) \cdot ( N ( t_{p} ) )^{2} - C \cdot \Delta _{p + 1} ( N ( t_{p} ) )^{2} \cdot N^{\prime } ( t_{p} )}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad \leq \biggl\vert \int _{t_{p}}^{t_{p + 1}} \bigl( N^{\prime } ( \tau ) - N^{\prime } ( t_{p} ) \bigr) \,\mathrm{d}\tau \biggr\vert + \Delta _{p + 1}^{2} \cdot N_{\max}^{2} \cdot \biggl\vert B \cdot \frac{f_{p + 1} - f_{p}}{t_{p + 1} - t_{p}} - C \cdot N^{\prime } ( t_{p} ) \biggr\vert \end{aligned}$$

after some manipulations. By the mean value theorem of calculus, there are \(\tau _{1}, \xi _{1} \in ( t_{p}, t_{p + 1} )\) such that

$$ N^{\prime \prime } ( \tau _{1} ) = \frac{N^{\prime } ( \tau ) - N^{\prime } ( t_{p} )}{\tau - t_{p}} \quad \mbox{and} \quad f^{\prime } ( \xi _{1} ) = \frac{f_{p + 1} - f_{p}}{t_{p + 1} - t_{p}}. $$

These equations yield

$$\begin{aligned} \bigl\vert N ( t_{p + 1} ) - \widetilde{N_{p + 1}} \bigr\vert \leq & \biggl\vert \int _{t_{p}}^{t_{p + 1}} ( \tau - t_{p} ) \cdot N^{\prime \prime } ( \tau _{1} ) \,\mathrm{d}\tau \biggr\vert + \Delta _{p + 1}^{2} \cdot N_{\max}^{2} \cdot \bigl\vert B \cdot f^{\prime } ( \xi _{1} ) - C \cdot N^{ \prime } ( t_{p} ) \bigr\vert \\ \leq & \frac{1}{2} \cdot \Delta _{p + 1}^{2} \cdot \bigl\lVert N^{ \prime \prime } ( t ) \bigr\rVert _{\infty } + \Delta _{p + 1}^{2} \cdot N_{\max}^{2} \cdot \bigl( B \cdot \bigl\lVert f^{\prime } ( t ) \bigr\rVert _{\infty } + C \cdot \bigl\vert N^{\prime } ( t_{p} ) \bigr\vert \bigr) \\ \leq & \Delta _{p + 1}^{2} \cdot C_{\mathrm{loc}}, \end{aligned}$$

where we set

$$\begin{aligned} C_{\mathrm{loc}} := & \frac{1}{2} \cdot \bigl\lVert N^{\prime \prime } ( t ) \bigr\rVert _{\infty } + B \cdot N_{\max}^{2} \cdot \bigl\lVert f^{\prime } ( t ) \bigr\rVert _{\infty } \\ &{} + C \cdot N_{\max}^{2} \cdot \bigl\{ A \cdot N_{ \max} + B \cdot M_{f} \cdot N_{\max}^{2} + C \cdot N_{ \max}^{3} \bigr\} . \end{aligned}$$

We finally obtain

$$ \bigl\vert N ( t_{p + 1} ) - \widetilde{N_{p + 1}} \bigr\vert \leq C_{\mathrm{loc}} \cdot \Delta _{p + 1}^{2} $$
(24)

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 \(\vert N_{p + 1} - \widetilde{N_{p + 1}} \vert \). We observe that

$$\begin{aligned} N_{p + 1} = & N_{p} \cdot \frac{1 + A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N_{p}}{1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2}} \\ = & N_{p} + N_{p} \cdot \frac{A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N_{p} - C \cdot \Delta _{p + 1} \cdot N_{p}^{2}}{1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2}} \end{aligned}$$

and

$$\begin{aligned} \widetilde{N_{p + 1}} = & N ( t_{p} ) \cdot \frac{1 + A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N ( t_{p} )}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \\ = & N ( t_{p} ) + N ( t_{p} ) \cdot \frac{A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N ( t_{p} )}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \\ &{} - N ( t_{p} ) \cdot \frac{C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}}. \end{aligned}$$

By the triangle inequality, we obtain

$$\begin{aligned}& \vert N_{p + 1} - \widetilde{N_{p + 1}} \vert \\& \quad \leq \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert + \biggl\vert N_{p} \cdot \frac{A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N_{p} - C \cdot \Delta _{p + 1} \cdot N_{p}^{2}}{1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2}} \\& \qquad {} - N ( t_{p} ) \cdot \frac{A \cdot \Delta _{p + 1} + B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot N ( t_{p} ) - C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}}{1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2}} \biggr\vert \\& \quad \leq \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \\& \qquad {} + \biggl\vert \frac{A \cdot \Delta _{p + 1} \cdot \{ N_{p} \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} ) \} }{ ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} )} \\& \qquad {} - \frac{A \cdot \Delta _{p + 1} \cdot \{ N ( t_{p} ) \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \} }{ ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} )} \biggr\vert \\& \qquad {} + \biggl\vert \frac{B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot \{ N_{p}^{2} \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} ) \} }{ ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} )} \\& \qquad {} - \frac{B \cdot \Delta _{p + 1} \cdot f_{p + 1} \cdot \{ ( N ( t_{p} ) )^{2} \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \} }{ ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} )} \biggr\vert \\& \qquad {} + \biggl\vert \frac{C \cdot \Delta _{p + 1} \cdot \{ N_{p}^{3} \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} ) \} }{ ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} )} \\& \qquad {} - \frac{C \cdot \Delta _{p + 1} \cdot \{ ( N ( t_{p} ) )^{3} \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \} }{ ( 1 + C \cdot \Delta _{p + 1} \cdot N_{p}^{2} ) \cdot ( 1 + C \cdot \Delta _{p + 1} \cdot ( N ( t_{p} ) )^{2} )} \biggr\vert \\& \quad \leq \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert + A \cdot \Delta _{p + 1} \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \\& \qquad {} + A \cdot C \cdot \Delta _{p + 1}^{2} \cdot N_{\max} \cdot N_{\mathrm{num},\max} \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \\& \qquad {} + B \cdot \Delta _{p + 1} \cdot \bigl\lVert f ( t ) \bigr\rVert _{\infty } \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \cdot \bigl\vert N_{p} + N ( t_{p} ) \bigr\vert \\& \qquad {} + C \cdot \Delta _{p + 1} \cdot \bigl\vert N_{p}^{3} - \bigl( N ( t_{p} ) \bigr)^{3} \bigr\vert + C^{2} \cdot \Delta _{p + 1}^{2} \cdot N_{\max}^{2} \cdot N_{\mathrm{num},\max}^{2} \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert . \end{aligned}$$

Since

$$\begin{aligned} \bigl\vert N_{p}^{3} - \bigl( N ( t_{p} ) \bigr)^{3} \bigr\vert = & \bigl\vert N_{p}^{2} + N_{p} \cdot N ( t_{p} ) + \bigl( N ( t_{p} ) \bigr)^{2} \bigr\vert \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \\ \leq & \bigl( N_{\max}^{2} + N_{\max} \cdot N_{ \mathrm{num},\max} + N_{\mathrm{num},\max}^{2} \bigr) \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \end{aligned}$$

and \(\Delta _{p + 1} \leq 1\), we can conclude that

$$\begin{aligned} \vert N_{p + 1} - \widetilde{N_{p + 1}} \vert \leq & \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert + A \cdot \Delta _{p + 1} \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \\ &{} + A \cdot C \cdot \Delta _{p + 1} \cdot N_{\max} \cdot N_{ \mathrm{num},\max} \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \\ &{} + B \cdot \Delta _{p + 1} \cdot \bigl\lVert f ( t ) \bigr\rVert _{\infty } \cdot N_{\mathrm{num},\max} \cdot N_{\max} \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \\ &{} + C \cdot \Delta _{p + 1} \cdot \bigl( N_{\max}^{2} + N_{ \max} \cdot N_{\mathrm{num},\max} + N_{\mathrm{num},\max}^{2} \bigr) \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \\ &{} + C^{2} \cdot \Delta _{p + 1} \cdot N_{\max}^{2} \cdot N_{ \mathrm{num},\max}^{2} \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert . \end{aligned}$$

Let us define

$$\begin{aligned} C_{\mathrm{prop}} := & \bigl\{ A + A \cdot C + N_{\max} \cdot N_{\mathrm{num},\max} + B \cdot \bigl\lVert f ( t ) \bigr\rVert _{\infty } \cdot N_{\mathrm{num},\max} \cdot N_{\max} \\ & {}+ C \cdot \bigl( N_{\max}^{2} + N_{ \max} \cdot N_{\mathrm{num},\max} + N_{\mathrm{num},\max}^{2} \bigr) + C^{2} \cdot N_{\max}^{2} \cdot N_{ \mathrm{num}, \max}^{2} \bigr\} . \end{aligned}$$

This yields

$$ \vert N_{p + 1} - \widetilde{N_{p + 1}} \vert \leq ( 1 + C_{ \mathrm{prop}} \cdot \Delta _{p + 1} ) \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert $$
(25)

for error propagation in time from the pth time step to the \(( p + 1 )\)th time step.

3) Our proof is based on the inequality \(1 + x \leq \exp ( x )\) for \(x \geq 0\). Note that \(t_{1} = 0 < t_{2} \ldots < t_{M - 1} < t_{M} = T\).

3.1) First, we want to prove inductively that

$$ \begin{aligned}[b] \bigl\vert N_{p + 1} - N ( t_{p + 1} ) \bigr\vert \leq{}& \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p + 1} - t_{1} \} \bigr) \\ &{} + C_{\mathrm{loc}} \cdot \sum\limits _{k = 2}^{p + 1} \Delta _{k}^{2} \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p + 1} - t_{k} \} \bigr) \end{aligned} $$
(26)

for all \(p \in \{ 0, \ldots , 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

$$\begin{aligned}& \bigl\vert N_{2} - N ( t_{2} ) \bigr\vert \\& \quad = \bigl\vert N_{2} - \widetilde{N_{2}} + \widetilde{N_{2}} - N ( t_{2} ) \bigr\vert \\& \quad \leq \vert N_{2} - \widetilde{N_{2}} \vert + \bigl\vert \widetilde{N_{2}} - N ( t_{2} ) \bigr\vert \\& \quad \leq ( 1 + C_{\mathrm{prop}} \cdot \Delta _{2} ) \cdot \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert + C_{ \mathrm{loc}} \cdot \Delta _{2}^{2} \\& \quad \leq \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{2} - t_{1} \} \bigr) + C_{\mathrm{loc}} \cdot \sum _{k = 2}^{2} \Delta _{k}^{2} \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{2} - t_{k} \} \bigr), \end{aligned}$$

We now assume that

$$ \begin{aligned} \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert \leq {}& \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p} - t_{1} \} \bigr) \\ & {}+ C_{\mathrm{loc}} \cdot \sum\limits _{k = 2}^{p} \Delta _{k}^{2} \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p} - t_{k} \} \bigr). \end{aligned} $$

We now want to show (26). We obtain

$$\begin{aligned}& \bigl\vert N_{p + 1} - N ( t_{p + 1} ) \bigr\vert \\& \quad = \bigl\vert N_{p + 1} - \widetilde{N_{p + 1}} + \widetilde{N_{p + 1}} - N ( t_{p + 1} ) \bigr\vert \\& \quad \leq \vert N_{p + 1} - \widetilde{N_{p + 1}} \vert + \bigl\vert \widetilde{N_{p + 1}} - N ( t_{p + 1} ) \bigr\vert \\& \quad \leq ( 1 + C_{\mathrm{prop}} \cdot \Delta _{p + 1} ) \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert + C_{ \mathrm{loc}} \cdot \Delta _{p + 1}^{2} \\& \quad \leq \exp ( C_{\mathrm{prop}} \cdot \Delta _{p + 1} ) \cdot \bigl\vert N_{p} - N ( t_{p} ) \bigr\vert + C_{ \mathrm{loc}} \cdot \Delta _{p + 1}^{2} \\& \quad \leq \exp ( C_{\mathrm{prop}} \cdot \Delta _{p + 1} ) \cdot \Biggl\{ \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p} - t_{1} \} \bigr) \\& \qquad {} + C_{\mathrm{loc}} \cdot \sum_{k = 2}^{p} \Delta _{k}^{2} \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p} - t_{k} \} \bigr) \Biggr\} + C_{\mathrm{loc}} \cdot \Delta _{p + 1}^{2} \\& \quad = \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p + 1} - t_{1} \} \bigr) \\& \qquad {} + C_{\mathrm{loc}} \cdot \sum_{k = 2}^{p} \Delta _{k}^{2} \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p + 1} - t_{k} \} \bigr) + C_{\mathrm{loc}} \cdot \Delta _{p + 1}^{2} \\& \quad = \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p + 1} - t_{1} \} \bigr) + C_{\mathrm{loc}} \cdot \sum_{k = 2}^{p + 1} \Delta _{k}^{2} \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p + 1} - t_{k} \} \bigr), \end{aligned}$$

and this finishes our inductive proof.

3.2) We define \(\Delta := \max_{p \in \mathbb{N}} \Delta _{p}\). We infer that

$$\begin{aligned}& \bigl\vert N_{p + 1} - N ( t_{p + 1} ) \bigr\vert \\& \quad \leq \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p + 1} - t_{1} \} \bigr) \\& \qquad {} + C_{\mathrm{loc}} \cdot \sum_{k = 2}^{p + 1} \Delta _{k}^{2} \cdot \exp \bigl( C_{\mathrm{prop}} \cdot \{ t_{p + 1} - t_{k} \} \bigr) \\& \quad \leq \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp ( C_{\mathrm{prop}} \cdot T ) + C_{ \mathrm{loc}} \cdot \Delta \cdot \sum_{k = 2}^{p + 1} \Delta _{k} \cdot \exp ( C_{\mathrm{prop}} \cdot T ) \\& \quad = \bigl\vert N_{1} - N ( t_{1} ) \bigr\vert \cdot \exp ( C_{\mathrm{prop}} \cdot T ) + C_{ \mathrm{loc}} \cdot \Delta \cdot T \cdot \exp ( C_{ \mathrm{prop}} \cdot T ). \end{aligned}$$

If the initial conditions of our time-continuous and time-discrete models coincide and \(\Delta \to 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. □

3.5 Algorithmic summary

We summarize our algorithmic approach in Table 1.

Table 1 Algorithmic summary of our explicit–implicit time-discretion population growth model

4 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.

4.1 Example 1: comparison of time discretization algorithms

First, we compare the classical explicit Eulerian time discretization algorithm

$$ \frac{N_{n + 1} - N_{n}}{\Delta _{n + 1}} = A \cdot N_{n} + B \cdot f_{n + 1} \cdot N_{n}^{2} - C \cdot N_{n}^{3},\quad N_{1} = N_{0} > 0, $$
(27)

to our explicit–implicit numerical solution algorithm

$$ \frac{N_{n + 1} - N_{n}}{\Delta _{n + 1}} = A \cdot N_{n} + B \cdot f_{n + 1} \cdot N_{n}^{2} - C \cdot N_{n}^{2} \cdot N_{n + 1},\quad N_{1} = N_{0} > 0, $$
(28)

by (16) for all \(n \in \{ 1, 2, \ldots , M - 2, M - 1 \} \). We consider the following inputs:

  • \(A = 1\), \(B = 0.01\), and \(C = 0.00001\);

  • Time step size \(\Delta = 1.0\);

  • Initial population size \(N_{1} = 10\);

  • Simulation time \(T = 6\);

  • Function \(f ( t ) \equiv 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 explicit–implicit 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.

Figure 1
figure 1

Results of our first example from Sect. 4.1

4.2 Example 2: monotonically increasing solution

We present an example that is monotonically increasing. We consider the following inputs:

  • \(A = 0.1\), \(B = 0.02\), and \(C = 0.00002\);

  • Time step size \(\Delta = 1.0\);

  • Initial population size \(N_{1} = 100\);

  • Simulation time \(T = 1000\);

  • Function \(f ( t ) = 0.05 \cdot \frac{0.001 + \exp ( -0.1 \cdot t )}{0.0009 + \exp ( -0.1 \cdot t )}\);

  • Lower bound \(m_{f} = 0.05\) and upper bound \(M_{f} = \frac{5}{90}\) of f.

In Fig. 2, we see that we can get results that look similar to monotonically increasing solutions from Verhulst’s logistic model.

Figure 2
figure 2

Results of our second example from Sect. 4.2

4.3 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:

  • \(A = 0.1\), \(B = 0.04\), and \(C = 0.00002\);

  • Time step size \(\Delta = 1.0\);

  • Initial population size \(N_{1} = 100\);

  • Simulation time \(T = 1000\);

  • Function \(f ( t ) = 0.1 + \frac{ ( 0.01 \cdot t^{2} + 0.1 \cdot t + 1 ) \cdot \exp ( -0.02 \cdot t )}{10 + \exp ( -0.02 \cdot t )}\);

  • Lower bound \(m_{f} = 0.1\) and upper bound \(M_{f} = 1.585\) of f.

We observe from Fig. 3 that it is possible to obtain solutions that first monotonically increase and later monotonically decrease in time.

Figure 3
figure 3

Results of our third example from Sect. 4.3

4.4 Example 4: oscillating solution

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

  • \(A = 0.1\), \(B = 0.02\), and \(C = 0.00002\);

  • Time step size \(\Delta = 1.0\);

  • Initial population size \(N_{1} = 100\);

  • Simulation time \(T = 1000\);

  • Function \(f ( t ) = 0.05 + 0.02 \cdot \sin ( \frac{t}{10} )\);

  • Lower bound \(m_{f} = 0.03\) and upper bound \(M_{f} = 0.07\) for f.

The results for this example with an oscillating solution are shown in Fig. 4.

Figure 4
figure 4

Results of our fourth example from Sect. 4.4

5 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.

5.1 Parameter estimation

We consider the proposed explicit–implicit time discretization method

$$ \frac{N_{n + 1} - N_{n}}{\Delta _{n + 1}} = A \cdot N_{n} + \widetilde{f_{n + 1}} \cdot N_{n}^{2} - C \cdot N_{n}^{2} \cdot N_{n + 1} $$
(29)

for all \(n \in \{ 1, 2, \ldots , M - 2, M - 1 \} \), where \(\widetilde{f_{n + 1}} = B \cdot f_{n + 1}\). Here, we assume that time points

$$ ( t_{1}, N_{1} ), \ldots , ( t_{M}, N_{M} ) $$

are given and, for simplicity and due to real-life data, \(\Delta _{n + 1} \equiv 1\). Finding the nonnegative coefficients

$$ A, \widetilde{f_{2}}, \ldots , \widetilde{f_{M}}, C $$

is our main goal. Hence, we want to solve the optimization problem

$$ \min_{A, \widetilde{f_{2}}, \ldots , \widetilde{f_{M}}, C} \left \lVert \mathbf{Z} \cdot \begin{pmatrix} A \\ \widetilde{f_{2}} \\ \vdots \\ \widetilde{f_{M}} \\ C \end{pmatrix} - \begin{pmatrix} N_{2} - N_{1} \\ N_{3} - N_{2} \\ \vdots \\ N_{M - 1} - N_{M - 2} \\ N_{M} - N_{M - 1} \end{pmatrix} \right \rVert _{2}, $$
(30)

where we set

Z:= ( N 1 N 2 N M 2 N M 1 | N 1 2 0 0 0 N 2 2 0 0 N M 2 2 0 0 0 N M 1 2 | N 1 2 N 2 N 2 2 N 3 N M 2 2 N M 1 N M 1 2 N M )

with respect to the nonnegativity constraints on the sought parameters, where \(\lVert \cdot \rVert _{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].

5.2 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:

  • \(A_{1} = 0.0225\), \(B_{1} = 0.005\), and \(C_{1} = 0.00025\);

  • Time step size \(\Delta = 1.0\) year;

  • Initial population size \(N_{1} = 2.536431\) billion people;

  • Simulation time \(T_{1} = 500\) years;

  • Function \(f_{1} ( t ) = \frac{0.02 \cdot ( 0.01 + \exp ( -0.02 \cdot t ) )}{0.05 + \exp ( - 0.02 \cdot t )}\);

  • Lower bound \(m_{f_{1}} = 0.004\) and upper bound \(M_{f_{1}} = 0.02\) of f.

We use the following parameters in our second user-chosen setting:

  • \(A_{2} = 0.0180\), \(B_{2} = 0.0009\), and \(C_{2} = 0.0002\);

  • Time step size \(\Delta = 1.0\) year;

  • Initial population size \(N_{1} = 2.536431\) billion people;

  • Simulation time \(T_{2} = 500\) years;

  • Function \(f_{2} ( t ) = \frac{ ( 0.2 \cdot t^{3} + 15.0 \cdot t^{2} + 50.0 \cdot t ) \cdot \exp ( -0.01 \cdot t )}{20.0 \cdot t^{2} + 1.0} + 0.01\);

  • Lower bound \(m_{f_{2}} = 0.01\) and upper bound \(M_{f_{2}} = 5.975\) of f.

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 \(\Delta = 1.0\) year;

  • Initial population size \(N_{1} = 2.536431\) billion people;

  • Function \(f_{3} ( t ) = 1\).

Figure 5
figure 5

The results of our nonnegative linear least-squares optimization can be seen as the blue graph. All three different settings are depicted in different colors as given in the legend

All chosen functions \(g_{j} ( t ) = B_{j} \cdot f_{j} ( t )\) for \(j \in \{ 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.

Figure 6
figure 6

Results for all three settings of total world 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

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.

Figure 7
figure 7

Results for all three settings of total world 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

5.3 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:

  • \(A = 0.0000001\) (since our approach states \(A \approx 0\), we stay with a small positive chosen number), \(B = 0.0002050\), and \(C = 0.00000019\);

  • Time step size \(\Delta = 1.0\) year;

  • Initial population size \(N_{1} = 82.802\) million people;

  • Simulation time \(T = 200\) years;

  • Function \(f ( t ) = \exp ( -0.03 \cdot t )\).

Figure 8
figure 8

The results of our nonnegative linear least-squares optimization can be seen as the blue graph. Our user-chosen simplified time-dependent coefficient function is given by the black graph

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.

Figure 9
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

6 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 explicit–implicit 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 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].