The Verhulst-Like Equations: Integrable OΔE and ODE with Chaotic Behavior

Abstract

In this paper, we study various variants of Verhulst-like ordinary differential equations (ODE) and ordinary difference equations (O$\Delta$E). Usually Verhulst ODE serves as an example of a deterministic system and discrete logistic equation is a classic example of a simple system with very complicated (chaotic) behavior. In our paper we present examples of deterministic discretization and chaotic continualization. Continualization procedure is based on Padé approximants. To correctly characterize the dynamics of obtained ODE we measured such characteristic parameters of chaotic dynamical systems as the Lyapunov exponents and the Lyapunov dimensions. Discretization and continualization lead to a change in the symmetry of the mathematical model (i.e., group properties of the original ODE and O$\Delta$E). This aspect of the problem is the aim of further research.

1. Introduction

Verhulst is credited with formulating the ordinary differential equation (ODE)

$$\frac{dN}{dt}=rN\left(1-\frac{N}{K}\right)$$

(1)

and coining the name logistic [1,2]. Later investigators proposed variations on the Verhulst equation (e.g., difference logistic equation), sometimes continuing to refer to these as logistic models. The first application of ODE (1) was connected with population problems, and more generally, problems in ecology. If the Verhulst model is used for describing change in population size N over time t, then in Equation (1) $r$ is the Malthusian parameter (rate of maximum population growth) and $K$ is the carrying capacity (i.e., the maximum sustainable population). Equation (1) is widely used in problems of ecology, economics, chemistry, medicine, pharmacology, epidemiology, etc. [3,4,5,6,7,8]. As a rule, this model is oversimplified for quantitative estimations, but reflects the key qualitative features of processes under consideration.

Equation (1) can be reduced to the form

$$\frac{dx}{dt}=rx(1-x)$$

(2)

where $x=N/K.$

The initial condition is

$$x(0)=a.$$

(3)

The Cauchy problem (Equations (2) and (3)) has the exact solution

$$x=\frac{a}{a+(1-a)e^{-\hspace{0.17em}rt}}$$

(4)

The discrete logistic equation can be written as follows

$$x_{n+1}=R{x}_n(1-x_n)$$

(5)

where parameter $R$ characterizes the rate of reproduction (growth) of the population; $R=rh$, parameter $h$ defines the time between consecutive measurements.

The nonlinear difference Equation (5) exhibits period doubling to chaos [5,6].

We will analyze a slightly different discrete logistic equation

$$x_{n+1}-x_n=R{x}_n(1-x_n)$$

(6)

Difference Equation (6) is obtained from ODE (2) using a forward difference scheme for the derivative with the step of discretization h.

Ordinary difference equation (O$\Delta$E) (6) is close enough in behavior to solutions to the equation [4,5]

$$x_{n+1}=x_n\exp [R(1-x_n)]$$

(7)

O$\Delta$E (7) is close enough in behavior to solutions to Equation (6) for $x_n\approx 1$.

Initial condition for O$\Delta$Es (5), (6) or (7) is

$$x_0=a$$

(8)

The discrete Cauchy problems (5) and (8), (6) and (8) or (7) and (8) for sufficiently large values of the parameter $R$ describe the complex, chaotic behavior of the system [3,4,5,6,7,8]. For O$\Delta$E (6), as it is shown in [3] for $R=2.3$ the solution starts to oscillate periodically around the value $x=1$. This solution is stable as long as $R<\sqrt{6}\approx 2.449$. For $R=2.500$ the process comes to steady periodic oscillations with a period of four.

It can be mentioned that the chaotic threshold for O$\Delta$E (7) is 2.6824 [4,5].

After describing the objects of our research, we proceed to the formulation of its goals. Both continuous and discrete logistic equations have been extensively investigated. The results of these studies are described in a large number of research and review books and papers [4,5,6,7,8,9]. Our study focuses on the discretization and continualization of nonlinear ODE and O$\Delta$E, while the logistic equations serve only as convenient and simple examples.

First problem: is it possible to discretize ODE (6) in such a way that the resulting O$\Delta$E has only regular solutions? This practically important issue has been studied quite well [10,11,12,13], and we, in fact, only briefly present the known results.

On the other hand, many researchers point that discrete logistic models are more adequate to the essence of the physical, economic or biological processes precisely because they have chaotic regimes [3]. In this regard, our second problem is: can such a continualization of the original O$\Delta$E be proposed that the resulting ODE has a chaotic solution? The practical value of these equations can be considered controversial, but the construction of such continualization procedure is interesting.

It is difficult to expect that standard continualization, based on the Taylor series, will provide the desired result. It remains to be hoped for the techniques based on the use of Padé approximants [14].

The paper is organized as follows: In Section 2 we propose discretization of Verhulst ODE without chaotic behavior. Continualization of the discrete Verhulst equation using Padé approximants is described in Section 3. Then, in Section 4 we display and discuss some results of numerical integrations. Finally, in Section 5 some conclusions are presented.

2. Integrable O$\Delta$E

As it is mentioned in [9], non-invertible maps, such as the logistic map, may display chaos. Of interest is the transformation of the original discrete logistic equation into a form leading to deterministic solutions.

Following the ideas of [10], we rewrite O$\Delta$E (6), in the following form

$$x_{n+1}-x_n=R{x}_n-R{x}_{n+1}{x}_n$$

(9)

This presentation makes it possible to express $x_{n+1}$ not a polynomial, but a fractional rational function $x_n$. Equation (9) with initial condition (8) has the exact solution of the form

$$x_n=\frac{a}{a+(1-a){(1+R)}^{-n}}$$

(10)

Thus, representation (9) allows one to obtain a difference scheme without chaotic behavior. Of course, the transition from Equation (6) to Equation (9) is based on the ideas of preserving the Lie symmetry of the original ODE [10,11,12,13]. However, in our article we do not develop this topic.

3. Continualization with Padé Approximants

Let us try to construct a continuous model (i.e., ODE), describing the chaotic behavior like original O$\Delta$E. As it is mentioned in [9], for generating chaotic behavior nonlinear ODE must have dimension $D\ge 3$

To construct a logistic-like ODE with chaotic behavior, additional modifications were introduced—piecewise constant argument, delay, and fractional derivative [15,16,17,18]. We use continualization based on Maclaurin expansion and Padé approximants.

For continualization of O$\Delta$E (6) let us introduce the continuous coordinate $x$ scaled in such a way that $x_n=x(nh)$. Suppose $x(t)$ slightly changing function we use Maclaurin expansion

$$x_{n+1}-x_n=h{x}_t+\frac{h^2}{2}{x}_{tt}+\frac{h^3}{6}{x}_{ttt}+\dots$$

(11)

The third order equation thus obtained

$$h^3{x}_{ttt}+3h^2{x}_{tt}+6h{x}_t-6Rx+6R{x}_{}^2=0$$

(12)

describes completely deterministic trajectories.

Let us consider fluctuations around the second equilibrium position, $x_n=1$. Using changing of variables

$$x_n=1+y_n;\left|y_n\right|<<1$$

(13)

one obtains

$$y_{n+1}-y_n=-R{y}_n(1+y_n)$$

(14)

The initial condition for this difference equation can be written as follows

$$y_0=\mathsf{\alpha },0<\mathsf{\alpha }<<1.$$

(15)

Suppose $y(t)$ slightly changing function we use Maclaurin expansion

$$y_{n+1}-y_n=h{y}_t+\frac{h^2}{2}{y}_{tt}+\frac{h^3}{6}{y}_{ttt}+\dots$$

The fifth order ODE can be written as follows

$$h\frac{d}{dt}\left[1+\frac{h}{2}\frac{d}{dt}+\frac{h^2}{6}\frac{d^2}{d{t}^2}+\frac{h^3}{24}\frac{d^3}{d{t}^3}+\frac{h^4}{120}\frac{d^4}{d{t}^4}\right]y=-Ry(1+y)$$

(16)

Equation (16) describes deterministic trajectories.

Transform the differential operator in square brackets of Equation (16) into the diagonal Padé approximation [2,2], we obtain:

$$1+\frac{h}{2}\frac{d}{dt}+\frac{h^2}{6}\frac{d^2}{d{t}^2}+\frac{h^3}{24}\frac{d^3}{d{t}^3}+\frac{h^4}{120}\frac{d^4}{d{t}^4}=\frac{1}{3}\hspace{0.17em}\frac{h^2\frac{d^2}{d{t}^2}+6h\frac{d}{dt}+60}{h^2\frac{d^2}{d{t}^2}-8h\frac{d}{dt}+20}$$

Then

$$h^3\frac{d^3}{d{t}^3}+6h^2\frac{d^2}{d{t}^2}+60h\frac{d}{dt}=-3R\hspace{0.17em}\left(h^2\frac{d^2}{d{t}^2}-8h\frac{d}{dt}+20\right)\hspace{0.17em}y\left(1+y\right)$$

and after routine transformations one obtains:

$$\begin{array}{l}{h}^3\frac{d^{3}y}{d{t}^3}+3h^2(2+R(1+2y))\frac{d^{2}y}{d{t}^2}+12h(5-2R)\frac{dy}{dt}-\\ 48hRy\frac{dy}{dt}+6R{h}^2{\left(\frac{dy}{dt}\right)}^2+60Ry(y+1)=0\end{array}$$

(17)

Point after (17) must be omitted

Let us formulate initial conditions for the third order ODE (17). From the initial condition for original difference Equation (15) one obtains

$$y(0)=\mathsf{\alpha }$$

(18)

Additional initial conditions for Equation (17) we choose in the following form:

$$\begin{array}{l}t=0:\hspace{1em}y=\mathsf{\alpha };\hspace{0.17em}\hspace{0.17em}\frac{dy}{dt}\approx \frac{y_1-y_0}{h}=-\frac{R\mathsf{\alpha }\left(1+\mathsf{\alpha }\right)}{h}\hspace{0.17em}\hspace{0.17em}\\ \frac{d^{2}y}{d{t}^2}\approx \frac{y_2-2y_1+y_0}{h^2}=\frac{R^2\mathsf{\alpha }\left(1+\mathsf{\alpha }\right)\left(\mathsf{\alpha }+\left(1+\mathsf{\alpha }\right)\left(1-R\mathsf{\alpha }\right)\right)}{h^2}\end{array}$$

(19)

Point after (19) must be omitted

Initial value problem (Equations (17) and (19)) can be transformed by going to “dimensionless time” $T=t/h$, to the following form:

$$\frac{d^{3}y}{d{T}^3}+3(2+R(1+2y))\frac{d^{2}y}{d{T}^2}+12(5-2R)\frac{dy}{dT}-48Ry\frac{dy}{dT}+6R{\left(\frac{dy}{dT}\right)}^2+60Ry(y+1)=0$$

(20)

$$T=0:\hspace{1em}y=\mathsf{\alpha };\hspace{0.17em}\hspace{0.17em}\frac{dy}{dT}=-R\mathsf{\alpha }\left(1+\mathsf{\alpha }\right);\hspace{0.17em}\hspace{0.17em}\frac{d^{2}y}{d{T}^2}=R^2\mathsf{\alpha }\left(1+\mathsf{\alpha }\right)\left(\mathsf{\alpha }+\left(1+\mathsf{\alpha }\right)\left(1-R\mathsf{\alpha }\right)\right)$$

(21)

4. Numerical Results

Numerical integration of Cauchy problem (Equations (20) and (21)) is carried out using the Adams–Bashforth–Moulton method (a predictor-corrector method), where $f\left(t_{n+1},y_{n+1}\right)$ is found by first applying the Adams–Bashforth method (the predictor), then using the Adams–Bashforth–Moulton method (the corrector).

The calculations were performed in the Maple 18 computing environment. The graphs below are depicted in the original variable $t=T\cdot h$.

The presented numerical results can be divided into three groups: describing periodic oscillations, periodic oscillations with subharmonics, and the chaotic oscillations. For $2.5\le R\le 2.88$ one obtains periodic oscillations (Figure 1, Figure 2 and Figure 3).

For the values $2.5\le R\le 2.88$, a small change in the initial conditions does not lead to a radical change in the behavior of the system. For $2.89\le R\le 3.0$ in a periodic solution appears subharmonics (Figure 4, Figure 5 and Figure 6).

In this case, with small changes in the initial conditions, the nature of oscillations of the system now undergoes significant changes (Figure 6c).

For $R>3.0$ the behavior of the system becomes chaotic (Figure 7, Figure 8 and Figure 9).

With the chaotic behavior of the system, small changes in the initial conditions lead to significant changes in the oscillations (Figure 9c).

At the same time, phase trajectories with small changes in the initial conditions have a characteristic structure. If we exclude points that correspond to the initial mode of establishing oscillations, the structure of the phase trajectories does not depend on the initial conditions (Figure 10 and Figure 11).

To correctly characterize the dynamics and confirmations of the chaotic behavior of the system we calculated for different values of R Lyapunov exponents and Lyapunov dimensions. If the system has at least one positive Lyapunov exponent, then it is chaotic [5,8,9]. The multi-paradigm numerical computing environment MATLAB was used to calculate Lyapunov exponents. Figure 12 and Figure 13 show the dynamics of Lyapunov exponents for R = 2.5 and R = 3.0, respectively. Since all Lyapunov exponents are negative, the system is not chaotic.

At R > 3.0, the largest Lyapunov exponent becomes positive, which indicates the appearance of chaos in the system (Figure 14 and Figure 15).

The obtained values of Lyapunov exponents confirm the earlier conclusions about the absence or presence of chaos in the system.

The Lyapunov dimension $D_L$ can be calculated using the formula:

$$D_L=j+{\displaystyle \sum _{i=1}^j\frac{{\mathsf{\lambda }}_i}{\left|\hspace{0.17em}{\mathsf{\lambda }}_{j+1}\hspace{0.17em}\right|}}$$

where $j$ is defined from the conditions:

$${\displaystyle \sum _{i=1}^j{\mathsf{\lambda }}_i}>0 \mathrm{and} {\displaystyle \sum _{i=1}^{j+1}{\mathsf{\lambda }}_i}<0$$

For R = 3.05 and R = 3.1 one obtains:

$$D_L\hspace{0.17em}|\underset{\underset{R=3.05}{}}{}=2.0102; D_L\hspace{0.17em}|\underset{\underset{R=3.1}{}}{}=2.0316,$$

which are consistent with that of a third order chaotic system [19].

5. Conclusions

Differential and difference equations are the main tools for mathematical modeling of physical, economic, environmental, social processes, and phenomena [3,4,7]. To study differential equations, the entire arsenal of Calculus and Functional Analysis is used. Difference equations are the standard tools for numerical analysis. The relationship between these procedures is nontrivial, discretization or continualization often fundamentally change the nature of the solution. Therefore, the study of approaches that allow these operations to preserve the basic properties of the original systems is interesting. The question of discretization of differential equations preserving their properties has been well studied and continues to be studied. In our article, we briefly touch on this problem by the example of conservative discretization of the logistic equation. This question is closely related to group properties of the original ODE and O$\Delta$E. We do not use such a deep theory, remaining within the framework of the phenomenological approach.

Our second task was to construct a continuous approximation of a discrete logistic equation with chaotic behavior. Continualization procedure was based on Padé approximants. To correctly characterize the dynamics of obtained ODE we measured such characteristic parameters of chaotic dynamical systems as the Lyapunov exponents and the Lyapunov dimensions.