Restarted Adomian Decomposition Method for Solving Volterra’s Population Model

Mariam Al-Mazmumy¹, Safa Otyuan Almuhalbedi²

¹Department of Mathematics, Faculty of Science-AL Faisaliah, King Abdulaziz University, Jeddah, KSA.
²Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, KSA.

DOI: 10.4236/ajcm.2017.72016                

Abstract

In this paper, we used an efficient algorithm to obtain an analytic approximation for Volterra’s model for population growth of a species within a closed system, called the Restarted Adomian decomposition method (RADM) to solve the model. The numerical results illustrate that RADM has the good accuracy.

Keywords

Adomian Decomposition Method, Restarted Adomian Method, Integro-Differential Equations, Volterra’s Population Model

1. Introduction

Integro-differential equations arise in many areas of mathematics and sciences, such as biology, ecology, medicine, physics and technology. This class of equations arises while modeling various engineering and natural science problems, and hence it attracts much attention in numerical computation and analysis. Recently, many attempts have been made to develop analytic and approximate methods to solve the Volterra’s population model, such as, Euler method [1] , the modified Euler method [1] , the classical fourth-order Runge-Kutta method [1] , Runge-Kutta-Fehlberg method [1] , Pade approximation [2] , Adomian decomposition method [3] , Sinc-Galerkin method [3] , composite spectral functions approximations [4] , Rational Chebyshev and Hermite functions collocation [5] , Homotopy perturbation method [6] and Variation iteration method [7] . In this research, we solve Volterra’s population model by Restarted Adomian decomposition method.

2. Description of the Volterra’s Population Model

Volterra proposed this model for a population [8] [9] . The simplest case, in characterizing the population dynamics of an isolated species, is to consider an asexually reproducing organism for which age is irrelevant and behavior does not change with time or with the number of the organisms. The number, u, must be sufficiently large so that the process can be well approximated by a deterministic treatment and by real, rather than integer, numbers. It may mean number of individuals, in which case it is an integer, but it may also mean total weight, weight of certain parts, total metabolism, or some other measure of quantity of life.

Under these assumptions the change in u is given by the Malthusian equation in [8] :

$\frac{\text{d}u}{\text{d}t}=au,a>0$ (1)

where a is birth rate, by integration, a geometrical law of increase (or decrease, if $a<0$ ) is obtained:

$u=u_0{\text{e}}^{at}$ (2)

Volterra’s practice discussed deeply the restrictive assumptions under which his formulas were derived. Then, he proceeded to remove one or two of the assumptions at a time, in the best Baconian fashion. The above assumptions imply, for example, that unlimited environmental resources are available to the species. One can easily allow for a finite environmental capacity by taking for a decreasing function of u. By assuming that a decreases linearly with u, one obtains the Verhlust-Pearl equation

$\frac{\text{d}u}{\text{d}t}=\left(a-bu\right)u,\left(a,b>0\right)$ (3)

where b is crowding coefficient, the integral, often called the “logistic curve,” is widely used even outside ecology. By taking into account specific mechanisms affecting reproduction or mortality lead to a much more complex functional relationship. For example, by taking a population living in a completely closed environment, such as some microorganisms confined to a test tube, the amount of nutrients available decreases with time in proportion to the total amount of “metabolism” that takes place in the tube from the beginning of the experiment. Total metabolism also determines the concentration of toxic waste in the medium. For simplicity, it is assumed that the metabolic activity of the population is directly proportional to the number of individuals and that its total amount affects linearly the coefficient of self-increase. Hence, the system can be repre- sented by the integro-differential equation

$u^{\prime }\left(t\right)=au\left(t\right)-b{u}^2\left(t\right)-c{\displaystyle {\int }_0^{t}u\left(t\right)u\left(\tau \right)\text{d}\tau },u\left(0\right)=\alpha .$ (4)

If the integral term is missing, the well-known logistic equation with birth rate $a>0$ and crowding coefficient $b>0$ is obtained. The last term, containing the integral, indicates the “total metabolism” or total amount of toxins produced since time zero. The individual death rate is proportional to this integral, so the population death rate due to toxicity must include a factor u. The presence of the toxic term, by considering the system being always closed, causes the population level to fall to zero in the long run. The relative size of the sensitivity to toxins, c, determines the manner in which the population evolves before its fated decay. By introducing the non-dimensional variables the time and population scales are obtained as:

$T=\frac{tc}{b}$ and $P=\frac{bu}{a}$ (5)

And the non-dimensional problem takes the form:

$k{P}^{\prime }\left(T\right)=P\left(T\right)-P^2\left(T\right)-{\displaystyle {\int }_0^{t}P\left(T\right)P\left(\tau \right)\text{d}\tau },P\left(0\right)=P_0.$ (6)

where $P=P\left(T\right)$ is the scaled population of identical individuals at a time T,

and the non-dimensional parameter $K=\frac{c}{ab}$ is a prescribed parameter. The

nondimensional parameter K plays a great role in the behavior of $P\left(T\right)$ concerning the rapid rise to a certain amplitude followed by an exponential decay to extinction. It is important to point out that for K small, the population is not sensitive to toxins, whereas the population is strongly sensitive to toxins for large K.

3. General Description of the Adomian Decomposition Method

Equation (4) may be written as

$Lu\left(t\right)=au\left(t\right)-b{u}^2\left(t\right)-c{\displaystyle {\int }_0^{t}u\left(t\right)u\left(\tau \right)\text{d}\tau }$ (7)

With initial condition: $u\left(0\right)=\alpha .$

Let $L=\frac{\text{d}}{\text{d}t}$ , so $L^{-1}(.)={\displaystyle {\int }_0^t(.)\text{d}t}$ , applying $L^{-1}$ of both sides in Equation (7),

and using the initial conditions, we obtain

$u\left(t\right)=\alpha +a{L}^{-1}u\left(t\right)-b{L}^{-1}{u}^2\left(t\right)-c{L}^{-1}{\displaystyle {\int }_0^{t}u\left(t\right)u\left(\tau \right)\text{d}\tau }$ (8)

The unknown solution function u assumed to be given by a series of the form

$u\left(t\right)={\displaystyle {\sum }_{n=0}^{\infty }{u}_n\left(t\right)}$ (9)

The nonlinear term usually represented by an infinite series of the so-called Adomian polynomials $A_n\left(t\right)$ and $B_n\left(t,\tau \right)$ , respectively, i.e., we set

$\{\begin{array}{l}{u}^2\left(t\right)={\displaystyle {\sum }_{n=0}^{\infty }{A}_n\left(t\right)}\\ u\left(t\right)u\left(\tau \right)={\displaystyle {\sum }_{n=0}^{\infty }{B}_n\left(t,\tau \right)}\end{array}$ (10)

Substituting (9) and (10) in (8) given the recursive relation

$\{\begin{array}{l}{u}_0=\alpha \\ u_{k+1}=a{L}^{-1}{u}_k-b{L}^{-1}{A}_k-c{L}^{-1}{\displaystyle {\int }_0^t{B}_k\left(t,\tau \right)\text{d}s},k\ge 0\end{array}$ (11)

From this recursive relation, we can compute $u_0,u_1,u_2,\cdots .$

The solution of Equation (4) is now determined. However, in practice series ${\displaystyle {\sum }_{n=0}^{\infty }{u}_n}$ must be truncated to the series ${\phi }_n={{\displaystyle \sum }}_{i=0}^n{u}_i$ with ${\lim }_{n\to \infty }{\phi }_n=u$ .

4. Restarted Adomian Decomposition Method (RADM)

In 2003, E. Babolian, et al. [10] introduced a new algorithm called “Restarted Adomian Method”, to improve the accuracy dramatically. This new method depends on adding a term to both sides of Equation (8). Let G be the proper term, which is determined next; then

$u\left(t\right)+G=\alpha +G+a{L}^{-1}u\left(t\right)-b{L}^{-1}{u}^2\left(t\right)-c{L}^{-1}{\displaystyle {\int }_0^{t}u\left(t\right)u\left(\tau \right)\text{d}\tau }$ (12)

By applying the “Modified Adomian Decomposition Method” on Equation (12), we obtain:

$\{\begin{array}{l}{u}_0=G\\ u_1=\alpha -G+a{L}^{-1}{u}_0-b{L}^{-1}{A}_0-c{L}^{-1}{\displaystyle {\int }_0^t{B}_0\left(t,\tau \right)\text{d}s}\\ u_{k+1}=a{L}^{-1}{u}_k-b{L}^{-1}{A}_k-c{L}^{-1}{\displaystyle {\int }_0^t{B}_k\left(t,\tau \right)\text{d}s},k\ge 1\end{array}$ (13)

Hence, the following algorithm is presented.

・ The algorithm

Choose small natural numbers m, n.

Step 1: Apply the Adomian method on Equation (12) and calculate $u_0,u_1,u_2,\cdots$ . Set

$\omega =u_0+u_1+\cdots +u_n$

Step 2: For $i=2:m$ , $G={\omega }^{i-1}$

$\{\begin{array}{l}{u}_0=G\\ u_1=\alpha -G+a{L}^{-1}{u}_0-b{L}^{-1}{A}_0-c{L}^{-1}{\displaystyle {\int }_0^t{B}_0\left(t,\tau \right)\text{d}s}\\ u_{k+1}=a{L}^{-1}{u}_k-b{L}^{-1}{A}_k-c{L}^{-1}{\displaystyle {\int }_0^t{B}_k\left(t,\tau \right)\text{d}s},k\ge 1\end{array}$ (14)

Set

$\omega =u_0+u_1+\cdots +u_n$

end.

Remarks:

1) ${\omega }^m$ can be considered as the approximate solution of Equation (4).

2) The Adomian Decomposition Method’ usually gives the sum of the first few terms, and, consequently, gives an approximation of u. In the new algorithm (RADM), $u_0$ are updated, while the terms with large index are not calculated. Therefore, m and n are considered to be small, say, m = 3 and n = 2.

5. Computation Results and Analysis

Example 1

a = 1, b = 1, c = 0.1, α = 0.1

$u^{\prime }\left(t\right)=u\left(t\right)-u^2\left(t\right)-0.1{\displaystyle {\int }_0^{t}u\left(t\right)u\left(\tau \right)\text{d}\tau },u\left(0\right)=0.1$

Solution:

Applying $L^{-1}(.)={\displaystyle {\int }_0^t(.)\text{d}t}$ in both sides given,

$u\left(t\right)=0.1+L^{-1}u\left(t\right)-L^{-1}{u}^2\left(t\right)-0.1L^{-1}{\displaystyle {\int }_0^{t}u\left(t\right)u\left(\tau \right)\text{d}\tau }$

By (ADM) the recursive relation is

$\{\begin{array}{l}{u}_0=0.1\\ u_{k+1}=L^{-1}{u}_k\left(t\right)-L^{-1}{A}_k\left(t\right)-0.1L^{-1}{\displaystyle {\int }_0^t{B}_k\left(t,\tau \right)\text{d}\tau },k\ge 0\end{array}$

Or

$\{\begin{array}{l}{u}_1=0.09x-0.0005x^2\\ u_2=0.036x^2-0.00058333x^3+0.00000167x^4\\ ⋮\end{array}$

The series solutions are

$u\left(t\right)=u_0+u_1+u_2+\cdots$

By applying the new algorithm (RADM) with $n=2$ and $m=3$ , we obtain:

Step 1

$\{\begin{array}{l}{u}_0=0.1\\ u_{k+1}=L^{-1}{u}_k\left(t\right)-L^{-1}{A}_k\left(t\right)-0.1L^{-1}{\displaystyle {\int }_0^t{B}_k\left(t,\tau \right)\text{d}\tau },k\ge 0,\end{array}$

Or

$\{\begin{array}{l}{u}_1=0.09x-0.0005x^2\\ u_2=0.036x^2-0.00058333x^3+0.00000167x^4\end{array}$

${\omega }^1=u_0+u_1+u_2=0.1+0.09x+0.0355x^2-0.00058333x^3+0.00000167x^4$

Step 2

$\{\begin{array}{l}{u}_0=0.1+0.09x+0.0355x^2-0.00058333x^3+0.00000167x^4\\ u_1=0.1-\left(\text{0}.1+0.09x+0.0355x^2-0.00058333x^3+0.00000167x^4\right)\\ +L^{-1}{u}_0\left(t\right)-L^{-1}{A}_0\left(t\right)-0.1L^{-1}{\displaystyle {\int }_0^t{B}_0\left(t,\tau \right)\text{d}\tau }\\ u_2=L^{-1}{u}_1\left(t\right)-L^{-1}{A}_1\left(t\right)-0.1L^{-1}{\displaystyle {\int }_0^t{B}_1\left(t,\tau \right)\text{d}\tau }\end{array}$

Or

$\{\begin{array}{l}{u}_1=0.0068999999x^3-0.001935x^4-0.00028258x^5+\cdots \\ u_2=0.00137999x^4-0.00057532x^5-0.00006516x^6+\cdots \end{array}$

${\omega }^2=u_0+u_1+u_2=0.1+0.09x+0.0355x^2+0.00631667x^3+\cdots$

Step 3

$\{\begin{array}{l}{u}_0=0.1+0.09x+0.0355x^2+0.00631667x^3+\cdots \\ u_1=0.1-\left(.1+0.09x+0.0355x^2+0.00631667x^3+\cdots \right)\\ +L^{-1}{u}_0\left(t\right)-L^{-1}{A}_0\left(t\right)-0.1L^{-1}{\displaystyle {\int }_0^t{B}_0\left(t,\tau \right)\text{d}\tau }\\ u_2=L^{-1}{u}_1\left(t\right)-L^{-1}{A}_1\left(t\right)-0.1L^{-1}{\displaystyle {\int }_0^t{B}_1\left(t,\tau \right)\text{d}\tau }\end{array}$

Or

$\{\begin{array}{l}{u}_1=0.00034499x^4+0.00046423x^5-0.00007783x^6+\cdots \\ u_2=-0.00034499x^4-0.00011923x^5+0.00005931x^6+\cdots \end{array}$

$\begin{array}{l}{\omega }^3=u_0+u_1+u_2=0.1+0.09x+0.0355x^2+0.00631667x^3-0.00055375x^4\\ -0.00051289x^5-0.00083165x^6+\cdots \end{array}$

The results produced by the present method with only few components (m = 5) are in a very good agreement with the best of the results of the methods listed in Table 1.

Example 2

a = 1, b = 1, c = 1, α = 0.1.

$u^{\prime }\left(t\right)=u\left(t\right)-u^2\left(t\right)-{\displaystyle {\int }_0^{t}u\left(t\right)u\left(\tau \right)\text{d}\tau },u\left(0\right)=0.1$

Solution:

Applying $L^{-1}(.)={\displaystyle {\int }_0^t(.)\text{d}t}$ in both sides given,

$u\left(t\right)=0.1+L^{-1}u\left(t\right)-L^{-1}{u}^2\left(t\right)-L^{-1}{\displaystyle {\int }_0^{t}u\left(t\right)u\left(\tau \right)\text{d}\tau }$

By (ADM) the recursive relation is

$\{\begin{array}{l}{u}_0=0.1\\ u_{k+1}=L^{-1}{u}_k\left(t\right)-L^{-1}{A}_k\left(t\right)-L^{-1}{\displaystyle {\int }_0^t{B}_k\left(t,\tau \right)\text{d}\tau },k\ge 0\end{array}$

Or

$\{\begin{array}{l}{u}_1=0.09x-0.005x^2\\ u_2=0.036x^2-0.00583333x^3+0.00016667x^4\\ ⋮\end{array}$

The series solutions are

$u\left(t\right)=u_0+u_1+u_2+\cdots$

Table 1. Error values obtained by the Adomian decomposition method and Restarted Adomian decomposition method for u(t).

By applying the new algorithm (RADM) with $n=2$ and $m=3$ , we obtain:

Step 1

$\{\begin{array}{l}{u}_0=0.1\\ u_{k+1}=L^{-1}{u}_k\left(t\right)-L^{-1}{A}_k\left(t\right)-L^{-1}{\displaystyle {\int }_0^t{B}_k\left(t,\tau \right)\text{d}\tau },k\ge 0,\end{array}$

Or

$\{\begin{array}{l}{u}_1=0.09x-0.005x^2\\ u_2=0.036x^2-0.00583333x^3+0.00016667x^4\\ ⋮\end{array}$

${\omega }^1=u_0+u_1+u_2=0.1+0.09x+0.031x^2-0.00583333x^3+0.00016667x^4$

Step 2

$\{\begin{array}{l}{u}_0=0.1+0.09x+0.031x^2-0.00583333x^3+0.00016667x^4\\ u_1=0.1-\left(0.1+0.09x+0.031x^2-0.00583333x^3+0.00016667x^4\right)\\ +L^{-1}{u}_0\left(t\right)-L^{-1}{A}_0\left(t\right)-L^{-1}{\displaystyle {\int }_0^t{B}_0\left(t,\tau \right)\text{d}\tau }\\ u_2=L^{-1}{u}_1\left(t\right)-L^{-1}{A}_1\left(t\right)-L^{-1}{\displaystyle {\int }_0^t{B}_1\left(t,\tau \right)\text{d}\tau }\end{array}$

${\omega }^2=u_0+u_1+u_2=0.1+0.09x+0.031x^2+0.00106667x^3+\cdots$

Step 3

$\{\begin{array}{l}{u}_0=0.1+0.09x+0.031x^2+0.00106667x^3+\cdots \\ u_1=0.1-\left(0.1+0.09x+0.031x^2+0.00106667x^3+\cdots \right)\\ +L^{-1}{u}_0\left(t\right)-L^{-1}{A}_0\left(t\right)-L^{-1}{\displaystyle {\int }_0^t{B}_0\left(t,\tau \right)\text{d}\tau }\\ u_2=L^{-1}{u}_1\left(t\right)-L^{-1}{A}_1\left(t\right)-L^{-1}{\displaystyle {\int }_0^t{B}_1\left(t,\tau \right)\text{d}\tau }\end{array}$

${\omega }^3=u_0+u_1+u_2=0.1+0.09x+0.031x^2+0.00106667x^3+\cdots$

The results produced by the present method with only few components (m = 5) are in a very good agreement with the best of the results of the methods listed in Table 2.

Table 2. Error values obtained by the Adomian decomposition method and Restarted Adomian decomposition method for u(t).

6. Conclusion

In this paper, we have applied Restarted Adomian method in solving nonlinear integro-differential equations. The numerical results show that RADM is more accurate than Adomian decomposition method of the solution Volterra’s population model in Table 1 and Table 2.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Te Beest, K. (1997) Numerical and Analytical Solutions of Volterras Population Model. SIAM Review, 39,484-493. https://doi.org/10.1137/S0036144595294850
[2]	Wazwaz, A.M. (1999) Analytical Approximations and Pad Approximants for Volterras Population Model. Applied Mathematics and Computation, 100, 3-25. https://doi.org/10.1016/S0096-3003(98)00018-6
[3]	Al-Khaled, K. (2005) Numerical Approximations for Population Growth Models. Applied Mathematics and Computation, 160, 865-873. https://doi.org/10.1016/j.amc.2003.12.005
[4]	Ramezani, M., Razzaghi, M. and Dehghan, M. (2007) Composite Spectral Functions for Solving Volterras Population Model. Chaos, Solitons & Fractals, 34, 588-593. https://doi.org/10.1016/j.chaos.2006.03.067
[5]	Parand, K., Rezaei, A. and Taghavi, A. (2010) Numerical Approximations for Population Growth Model by Rational Chebyshev and Hermite Functions Collocation Approach: A Comparison. Mathematical Methods in the Applied Sciences, 33, 2076-2086. https://doi.org/10.1002/mma.1318
[6]	Mohyud-Din, S.T., Yildirim, A. and Gulkanat, Y. (2010) Analytical Solution of Volterra’s Population Model. Journal of King Saud University (Science), 22, 247-250. https://doi.org/10.1016/j.jksus.2010.05.005
[7]	Al-Wesabi, Y., Dahawi, A.A., Daniel, Y. and Murid, A.H.M. (2014) Analytical Solution of Volterra’s Population Model Using Variation Iteration Method (VIM). 1st International Conference of Recent Trends in Information and Communication Technologies, 413-423.
[8]	Scudo, F. (1971) Vito Volterra and Theoretical Ecology. Theoretical Population Biology, 2, 1-23. https://doi.org/10.1016/0040-5809(71)90002-5
[9]	Small, R.D. (1989) Population Growth in a Closed System and Mathematical Modeling. Classroom Notes in Applied Mathematics, SIAM, Philadelphia, 317-320.
[10]	Babolian, E. and Javadi, Sh. (2003) Restarted Adomian Method for Algebraic Equations. Applied Mathematics and Computation, 146, 533-541. https://doi.org/10.1016/S0096-3003(02)00603-3