Exact solution to fractional logistic equation

16 downloads 214708 Views 473KB Size Report
Feb 21, 2015 - The linear representation of the FLE is exact and its solution involves the ... functions that reduces to the well known analytic solution in the limit ...
Physica A 429 (2015) 103–108

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Exact solution to fractional logistic equation Bruce J. West Information Sciences Directorate, Army Research Office, Research Triangle Park, NC 27709, United States

highlights • • • •

The logistic equation is generalized to include memory using the fractional calculus. The fractional logistic equation (FLE) is solved using an infinite-order linear representation. The linear representation of the FLE is exact and its solution involves the inversion of an infinite-order matrix. The technique is applicable to a large class of polynomial nonlinear fractional rate equations.

article

info

Article history: Received 18 January 2015 Available online 21 February 2015 Keywords: Fractional calculus Logistic Nonlinear Exact solution

abstract The logistic equation is one of the most familiar nonlinear differential equations in the biological and social sciences. Herein we provide an exact solution to an extension of this equation to incorporate memory through the use of fractional derivatives in time. The solution to the fractional logistic equation (FLE) is obtained using the Carleman embedding technique that allows the nonlinear equation to be replaced by an infinite-order set of linear equations, which we then solve exactly. The formal series expansion for the initial value solution of the FLE is shown to be expressed in terms of a series of weighted Mittag-Leffler functions that reduces to the well known analytic solution in the limit where the fractional index for the derivative approaches unity. The numerical integration to the FLE provides an excellent fit to the analytic solution. We propose this approach as a general technique for solving a class of nonlinear fractional differential equations. Published by Elsevier B.V.

1. Introduction The simplest ordinary differential equation that has found explanatory value in phenomena ranging from the decay of radioactive particles to the growth of populations is the rate equation, where the rate of change in the number of entities of interest N (t ) is determined by dN (t ) dt

= kN (t ).

(1)

The solution to the rate equation is the simple exponential N (t ) = N (0)ekt , with N (0) being the initial population. The number of remaining radioactive particles shrinks exponentially in time since the rate constant is negative (k < 0), whereas in animal population studies the rate constant is positive (k > 0) so that the population grows exponentially, as conjectured for human population by the eighteenth century cleric Malthus [1]. A fundamental assumption underlying Eq. (1) is that for the process being discussed time is a homogeneous passive quantity that indexes the changing events through the regular ticking of a clock. However this is not necessarily the case for

E-mail address: [email protected]. http://dx.doi.org/10.1016/j.physa.2015.02.073 0378-4371/Published by Elsevier B.V.

104

B.J. West / Physica A 429 (2015) 103–108

complex phenomena wherein the time derivative is replaced by a fractional derivative to yield the fractional rate equation (FRE) [2] Dαt [u (t )] = kα u(t ).

(2)

where the rate of growth has been raised to a power of the order of the fractional derivative to retain dimensional consistency. We restrict the index to the domain 1 ≥ α > 0 in order to retain the spirit of a growth process and we interpret the fractional time derivative in the Caputo sense. The Laplace transform of the Caputo fractional derivative is given by LT Dαt [f (t )] ; s ≡









e−st Dαt [f (t )] = sα f (s) − sα−1 f (0);

1≥α>0

(3)

0

and  f (s) is the Laplace transform of f (t ). The solution to the fractional rate equation in terms of the Laplace variable is

 u (s) =

sα−1

u (0) sα + k whose inverse Laplace transform yields [3,2] u (t ) = Eα (kt α ) u (0)

(4)

(5)

where the Mittag-Leffler function (MLF) Eα (·) is defined by the series Eα ( z ) =

∞ 

zn

n =0

Γ (nα + 1)

.

(6)

In the limit α = 1 Eq. (5) simplifies to the exponential solution to the linear rate equation. The solutions to such FREs have provided excellent fits to empirical data in a wide range of applications from biological reactions in the making of tequila [4] and other bioengineering applications [5], to allometry relations characterizing the size-dependence of the time experienced by elephants, humming birds and all the animals in between Ref. [6]. A nonlinear growth equation was introduced into population dynamics by Verhulst [7] to quench the unbounded growth in human population proposed by Malthus [1] and mitigate the utterly dismal predictions of the fate of mankind entailed by Malthus’ prediction. Verhulst [7] introduced the nonlinear term into the rate equation and obtained what subsequently became known as the logistic equation du(t )

= ku(t ) [1 − u(t )] (7) dt where the variable u(t ) is the population N (t ) normalized to its maximum attainable value Nmax ; u (t ) = N (t ) /Nmax . The logistic equation is one of the few nonlinear equations that has a known exact closed form solution: u0 (8) u(t ) = u0 + (1 − u0 )e−kt where the initial (t = 0) state is u0 = N (0)/Nmax . In the asymptotic limit the normalized population approaches the value unity. This kind of saturation behavior has been used to explain all manner of complex phenomena where the initial exponential growth of a population is eventually curtailed by limited resources becoming exhausted. The sigmoidal behavior of the solution to the logistic equation in addition to applications in ecology has been observed in the product concentration of autocatalytic reactions, which has been argued by some to play a major role in the process of life itself [8,9]. There has been increasing application of the logistic curve in medicine where it has been used to model the growth of tumors [10]. It has also been used to model the social dynamics of replacement technologies by Fisher and Pry [11] and the adaptability of society to innovation. Examples of the latter in the transportation of goods appear as, canals being replaced by railroads, which were ceded to intercoastal highway systems, that in turn were replaced by international airfreight, each in its turn being represented by its own sigmoidal growth [12]. The dual aspects of complexity, fractional growth and nonlinearly induced saturation, dovetail into the fractional logistic equation (FLE): Dαt [u (t )] = kα u(t ) [1 − u(t )] .

(9)

Of course one cannot use the Laplace transform method to directly solve such a nonlinear fractional dynamic equation and a number of creative approximation techniques have been devised for its solution [13–15]. Kowalski [16] pointed out that in 1932 Carleman [17], following the ideas of Poincaré and Fredholm, demonstrated that a nonlinear system of ordinary differential equations with polynomial nonlinearities can be reduced to an infinite-order system of linear differential equations, without approximation. One can then adopt their favorite method for solving the set of linear equations and extract the solution to the original nonlinear equation. Herein we construct an exact solution to the FLE using the Carleman embedding technique. In Section 2 we demonstrate how to apply the technique to construct an infinite-order system of linear fractional differential equations equivalent to the FLE. The infinite-order linear system of equations is solved exactly using Laplace transforms and matrix methods to obtain a solution in terms of a weighted sum over MLF’s. In Section 3 the analytic solution to the FLE is shown to yield the logistic form of Eq. (8) when α = 1 and to agree with the direct numerical integration of the FLE for a variety of values of α . We draw some conclusions in Section 4.

B.J. West / Physica A 429 (2015) 103–108

105

2. Carleman embedding technique In this section we apply the Carleman technique to solve the FLE. To accomplish this we introduce a function of powers of the normalized population u:

φ n = un

(10)

where n ∈ Z+ is contained in the set of non-negative integers. Consequently, the logistic equation is replaced with the linear differential-difference equation dφn dτ

=



Cnm φm .

(11)

m≤n

and we have introduced the scaled time τ = kt. Eq. (11) represents an infinite-order linear system of rate equations in which the dynamics of Eq. (7) are embedded at n = 1. Note that the replacement of the solution to the nonlinear dynamic equation u(t ) with that for φn (t ) can be generalized to multiple dimensions without changing the conclusions drawn herein. Montroll [18] demonstrated how to solve a rate equation with a quadratic nonlinearity using the Carleman embedding technique to obtain the solution from an infinite set of linear equations. We are guided by his analysis and consider the set of equations generated by Eq. (10) to obtain dφn dτ

dun

du = nun−1 dτ dτ = nφn − nφn+1

=

(12)

reducing the coefficients in Eq. (11) to just two terms for each n. In this way just as the logistic equation can be replaced with the FLE the equivalent linear system of equations for the powers of the population can be replaced by the fractional differential-difference equations: Dατ [φn ] = nφn − nφn+1 .

(13)

This latter system is an infinite set of linear equations that has embedded in it the fractional logistic equation (Eq. (9)) at n = 1 and which could not have been obtained by taking the fractional derivative of Eq. (10) [19]. Given the linear form of Eq. (13) Laplace transforms can be used to obtain a solution. The Laplace transform of the solution  φn (s) and the initial conditions φ1 (0) = u0 ≡ C ; φn (0) = C n can be used in conjunction with the Laplace transform of the Caputo derivative

  LT Dατ [φn ] ; s = sα φn (s) − sα−1 C n ,

(14)

to construct the system of linear equations for the Laplace variable:

φn (s) + n φn+1 (s) = sα−1 C n . (sα − n) 

(15)

In matrix form we now have A φ (s) = sα−1 C,

(16)

where A is the infinite dimensional matrix of diagonal form that couples sequential powers of the population and whose explicit form is sketched in the Appendix. In the Appendix the algebra for inverting the infinite-order coefficient matrix in Eq. (16) is outlined. The first row of the inverse matrix allows us to obtain an infinite series for the Laplace variable  φ1 (s), which can be summed and inverse Laplace transformed to yield the analytic solution to the initial value problem in terms of MLFs u(t ) =

n ∞   u0 − 1 n =0

u0

Eα (−nkα t α ) .

(17)

Eq. (17) has the form of an expansion over eigenfunctions, where the exponentials that appear in the linear case are replaced by MLFs in the nonlinear case. The eigenvalues are given by −nkα . 3. Discussion Of course, we want to test the validity of this new solution. However we do not have a second solution to the FLE obtained using another technique with which to compare Eq. (17). Lacking such a comparison we have two ways to test the validity of the solution: (1) determine that it goes to the appropriate limits asymptotically in time and (2) compare it with the direct numerical integration of the FLE. The series solution Eq. (17) with α = 1 is shown to be equivalent to Eq. (8) using the property of MLF’s: lim Eα (−nkα t α ) = e−nkt .

α→1

(18)

106

B.J. West / Physica A 429 (2015) 103–108

Fig. 1. The analytic solution to the FLE given by Eq. (17) (dashed) is compared to the numerical solutions to Eq. (9) for values of the fractional order

α = 0.9, 0.75, 0.5, 0.25 for u(0) = 0.75, k = 01 and 1t = 0.001.

Inserting this exponential into Eq. (17) results in the explicitly summable series u (t ) =

n ∞   u0 − 1 u0

n =0

e−nkt =

u0 u0 + (1 − u0 )e−kt

,

(19)

which is the exact solution to the logistic equation. The asymptotic properties u0 = u(0) and u(∞) = 1 are shared by the solution to the FLE given by Eq. (17), although the rate of approach to these asymptotic values are quite different in the two cases; being exponential in time for the logistic equation, but inverse power law in time for the MLF. The FLE can also be integrated numerically. The integration of the FLE is compared with the analytic solution for four values of the fractional-order parameter α in Fig. 1. It is clear that there is excellent agreement between the analytic solution and the numerical integration of the FLE at early times for all values of the fractional order parameter. At late times there is a modest deviation from the analytic solution for α ≤ 0.5, which we attribute to the fact that we have used the same integration parameters in all the calculations. This is compelling evidence that we have obtained the exact solution to the initial value problem for the FLE. 4. Conclusion The Carleman embedding approach has successfully determined the solution to many nonlinear differential equations as presented and discussed by Kowalski and Steeb [20]. Herein we adopted that technique to solve the FLE, noting that this approach has not previously been applied to a nonlinear FRE. The Carleman embedding technique is straight forward to implement for a polynomial function on the right hand side of the rate equation for a single variable. However, the more orders in the polynomial the more off-diagonal terms that appear in the coupling matrix. Consequently, it would be advantageous to find an alternative way to implement the linearization ideas of Carleman. We are presently working on a linearization technique, using an operator method to replace the matrix approach used herein, and the preliminary results are encouraging. Acknowledgment The author thanks the Army Research Office for supporting this research and M. Turalska for supplying the figure. Appendix The coefficient matrix in Eq. (16) can be written explicitly as



(sα − 1)

0  0 0 ·

1 (sα − 2) 0 0

·

0 2 (sα − 3) 0 0

0 0 3 (sα − 4) 0

C 0 ϕ1 C 2  0 φ2      0 φ3  = sα−1 C 3  .    ·  4 ·

 

·

·

 

·

(20)

B.J. West / Physica A 429 (2015) 103–108

107

The inverse of the coefficient matrix is not difficult to obtain but does require some effort



A −1

1

(sα

(sα

− 1)

  0     0 =    0   0

−1 − 1) (sα − 2)

2

(sα

2) (sα

− −2 (sα − 2) (sα − 3)

− 3)

−2 · 3 − − 2) (sα − 3) (sα − 4) 2·3 (sα − 2) (sα − 3) (sα − 4) −3 α (s − 3) (sα − 4) (sα

1) (sα

·

·



0

(sα

0

0

(sα − 4)

·

0

0

0

·

  ·     · . (21)    ·   ·

0

0

0

0

·

1

(sα − 2)

0



1) (sα

1

− 3)

1

· ·

Multiplying Eq. (20) on the left by the inverse matrix provides a formal solution for all the powers of the Laplace variable solution. Embedded in this set of solutions is the one to our original equation, which we take from the first row sα−1 C 2 2sα−1 C 3 3 · 2sα−1 C 4 sα−1 C  +··· φ 1 ( s) = α − α + α − α α α α α s −1 (s − 1) (s − 2) (s − 1) (s − 2) (s − 3) (s − 1) (s − 2) (sα − 3) (sα − 4) ∞  sα−1 Γ (n) Γ (sα − n) = (−1)n−1 C n Γ (sα ) n=1

=

∞  (−1)n−1 C n sα−1 B (n, sα − n)

(22)

n=1

where we have used for the infinite dimensional linear system.

Γ (sα ) = (sα − 1) Γ (sα − 1) = (sα − 1) (sα − 2) Γ (sα − 2) = · · ·

= (sα − 1) (sα − 2) · · (sα − n) Γ (sα − n)

(23)

and introduced the beta function B(n, m) =

1



z n−1 (1 − z )m−1 dz = 0

Γ (n) Γ (m) . Γ (n + m)

(24)

Hence the Laplace variable solution can be re-expressed using this integral to obtain

 φ 1 ( s) =

∞ 

sα−1 (−1)n−1 C n



1

α z n−1 (1 − z )s −n−1 dz

0

n =1

in which we can explicitly evaluate the sum to obtain

 φ1 (s) = sα−1 C

1

 0

α

(1 − z )s −1 dz . 1 + [C − 1] z

(25)

Introducing the final transformation of variables (1 − z ) = e−τ and noting the change in integration limits allows us to write

 φ 1 ( s) = C

 0



α sα−1 e−s τ dτ

(26)

C + (1 − C ) e−τ

and evaluating the integral, by expanding the function in an infinite series under the condition 1 > C > 0.5 yields

 φ 1 ( s) =

n α−1 ∞   C −1 s . C n + sα n =0

(27)

The inverse Laplace transform of Eq. (27) for the time-dependent solution to the FLE

φ1 (t ) =

n ∞   C −1 n =0

C

LT −1



sα−1 n+s

 ; kt , α

(28)

which from the inverse Laplace transform [3,21] can be expressed in terms of a series over Mittag-Leffler functions given by Eq. (17). References [1] T.R. Malthus, Malthus—Population: The First Essay (1798), University of Michigan Press, Ann Arbor, MI, 1959. [2] B.J. West, P. Grigolini, Complex Webs: Anticipating the Improbable, Cambridge University Press, Cambridge, UK, 2011.

108 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

B.J. West / Physica A 429 (2015) 103–108 I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, CA, 1999. R. Toledo-Hernandez, V. Rico-Ramrez, G.A. Iglesias-Silva, U.M. Diwekar, Chem. Eng. Sci. 117 (2014) 217. R.L. Magin, Fractinal Calculus in Bioengineering, Begell House, Connecticut, 2006. D. West, B.J. West, Phys. Life Rev. 10 (2013) 210. P.F. Verhulst, Notice sur la loi que la population sint dons son accroissement, Math. Phys. 10 (1838) 113. S. Kauffman, At Home in the Universe: The Search for the Laws of Self-Organization and Complexity, Oxford University Press, UK, 1995. R. Ulanowicz, Ecology, The Ascendent Perspective, Columbia University Press, 1997. U. Forys, A. Marciniak-Czochra, Logistic equations in tumor growth modelling, Int. J. Appl. Math. Comput. Sci. 13 (2003) 317. T.C. Fisher, R.H. Fry, A simple substitution model of technological change, Technol. Forecast. Soc. Change 3 (1971) 75. B.J. West, Complex Worlds: Uncertain, Unequal and Unfair, Black Rose Writing, TX, 2012. S. Das, P.K. Gupta, K. Vishal, Approximate approach to the Das model of fractional logistic population growth, Appl. Appl. Math. 5 (2010) 605. A.M.A. El-Sayed, A.E.M. El-Mesiry, H.A.A. El-Saka, On the fractional-order logistic equation, Appl. Math. Lett. 20 (2007) 817. M.M. Khader, M.M. Babatin, On approximate solutions for fractional logistic differential equation, Math. Probl. Eng. 2013 (2013). Article ID 391901. K. Kowalski, Nonlinear dynamical systems and classical orthogonal polynomials, J. Math. Phys. 38 (1997) 2403. T. Carleman, Acta Math. 59 (1932) 63. E.W. Montroll, On the solution of nonlinear rate equations by matrix inversion, AIP Conf. Proc. 46 (1978) 337. A. Svenkeson, M.T. Beig, M. Turalska, B.J. West, P. Grigolini, Fractional trajectories: decorrelation versus friction, Physica A 393 (2013) 5663–5672. K. Kowalski, W.-H. Steeb, Nonlinear Dynamical Systems and Carleman Linearization, World Scientific, Singapore, 1991. B.J. West, M. Bologna, P. Grigolini, Physics of Fractal Operators, Cambridge Pub., New York, 2003.