- Email: [email protected]
Rational Chebyshev tau method for solving Volterra’s population model
Applied Mathematics and Computation 149 (2004) 893–900 www.elsevier.com/locate/amc
Rational Chebyshev tau method for solving VolterraÕs population model K. Parand a, M. Razzaghi a
a,b,*
Department of Applied Mathematics, Amirkabir University of Technology, Tehran 15914, Iran b Department of Mathematics and Statistics, Mississippi State University, Mississippi State, MS 39762, USA
Abstract An approximate method for solving VolterraÕs population model for population growth of a species in a closed system is proposed. VolterraÕs model is a nonlinear integro-differential equation where the integral term represents the effect of toxin. The approach is based on a rational Chebyshev tau method. The VolterraÕs population model is first converted to a nonlinear ordinary differential equation. The operational matrices of derivative and product of rational Chebyshev functions are presented. These matrices together with the tau method are then utilized to reduce the solution of the VolterraÕs model to the solution of a system of algebraic equations. Illustrative examples are included to demonstrate the validity and applicability of the technique and a comparison is made with existing results. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Rational Chebyshev; Tau method; VolterraÕs population model
1. Introduction The Volterra model for population growth of a species within a closed system is given in [1,2] as Z t du j ¼ u u2 u uðxÞ dx; uð0Þ ¼ u0 ; ð1Þ dt 0
*
Corresponding author. E-mail address: [email protected] (M. Razzaghi).
0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.09.006
894
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
where uðtÞ is the scaled population of identical individuals at a time t, and j ¼ c=ab is a prescribed nondimensional parameter. Moreover, a > 0 is the birth rate coefficient, b > 0 is the crowding coefficient, and c > 0 is the toxicity coefficient. The coefficient c indicates the essential behavior of the population evolution before its level falls to zero in the long term. One may show that the only equilibrium solution of (1) is the trivial solution uðtÞ ¼ 0. Furthermore, the analytical solution [3] Z t Z s 1 uðtÞ ¼ u0 exp 1 uðsÞ uðxÞ dx ds ð2Þ j 0 0 shows that uðtÞ > 0 for all t if u0 > 0. The solution of Eq. (1) has been of considerable concern. Although a closed form solution has been achieved in [3,4], but it was formally shown that the closed form solution cannot lead to any insight into the behavior of the population evolution [1]. Some approximate and numerical solutions for VolterraÕs population model have been reported. In Ref. [2], the successive approximations method was offered for the solution of Eq. (1), but was not implemented. In [3], the singular perturbation method for VolterraÕs population model is considered. It is shown in [3] that for the case j 1, where populations are weakly sensitive to toxins, a rapid rise occurs along the logistic curve that will reach a peak and be followed by a slow exponential decay. And, for j large, where populations are strongly sensitive to toxins, the solution is proportional to sech2 ðtÞ. In this case the solution uðtÞ has a smaller amplitude compared to the amplitude of uðtÞ for the case j 1. In [4], three numerical algorithms namely the Euler method, the modified Euler method and the fourth order Runge–Kutta method for the solution of Eq. (1) are obtained. Moreover, a phase-plane analysis are implemented. In [4], the numerical results are correlated to give an insight of the problem and its solution without using perturbation techniques. However, the performance of the traditional numerical techniques is well known in that it provides grid points only, and in addition, it requires a large amounts of calculations. In [4], 20 iteration steps is carried out to obtain the grid points. In [1], the series solution method and the decomposition method are implemented independently to Eq. (1) and to a related ordinary differential equation (ODE). Furthermore, the Pade approximations are used in the analysis to capture the essential behavior of the population uðtÞ of identical individuals. In the present paper, we introduce a new computational method for solving Eq. (1). The VolterraÕs population model is first converted to an equivalent nonlinear ODE, the solution of which is then approximated by a rational Chebyshev functions with unknown coefficients. The operational matrices of derivative and product of rational Chebyshev functions are given. These matrices together with the tau method are then utilized to evaluate the un-
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
895
known coefficients and find approximate solutions for uðtÞ. The tau method was invented by Lanczos in 1938 [5]. The method is based on expanding the required approximate solution as the elements of a complete set of orthogonal functions. In tau method, unlike the Galerkin approximation, the expansion functions are not required to satisfy the boundary constraint individually [6,7]. The sections of this paper are organized as follows: in Section 2 we describe the formulation of rational Chebyshev functions required for our subsequent development. Section 3 converts VolterraÕs population to a nonlinear ODE. Section 4 summarizes the application of the rational Chebyshev tau method to the solution of nonlinear ODE for the Volterra population model. As a result a set of nonlinear algebraic equations is formed, and a solution of the considered ODE is introduced. In Section 5 the proposed method is applied to the numerical example, and a comparison is made with existing methods in the literature.
2. Properties of rational Chebyshev functions 2.1. Rational Chebyshev functions The well-known Chebyshev polynomials are orthogonal inffi the interval pffiffiffiffiffiffiffiffiffiffiffiffi ½1; 1 with respect to the weight function wðyÞ ¼ 1= 1 y 2 and can be determined with the aid of the following recurrence formula: T0 ðyÞ ¼ 1;
T1 ðyÞ ¼ y;
Tnþ1 ðyÞ ¼ 2yTn ðyÞ Tn1 ðyÞ;
n P 1:
The rational Chebyshev (RC) functions are defined in [8] by x1 Rn ðxÞ ¼ Tn : xþ1 Thus RC functions satisfy: R0 ðxÞ ¼ 1;
R1 ðxÞ ¼
Rnþ1 ðxÞ ¼ 2
x1 ; xþ1
x1 Rn ðxÞ Rn1 ðxÞ; xþ1
ð3Þ n P 1:
ð4Þ
Rational Chebyshev pffiffiffi functions are orthogonal with respect to the weight function wðxÞ ¼ 1= xðx þ 1Þ on the interval ½0; þ1Þ, with the orthogonality property: Z 1 cm p dnm Rn ðxÞRm ðxÞwðxÞ dx ¼ 2 0
896
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
with cm ¼
m ¼ 0; m 6¼ 1;
2; 1;
where dnm is the Kronecker function. 2.2. Function approximation A function f ðxÞ defined over the interval ½0; þ1Þ may be expanded as f ðxÞ ¼
1 X
ai Ri ðxÞ;
ð5Þ
i¼0
where ai ¼
2 ci p
Z
1
Ri ðxÞf ðxÞwðxÞ dx: 0
If f ðxÞ in Eq. (5) is truncated up to the N th terms, then Eq. (5) can be written as f ðxÞ ¼
N 1 X
ai Ri ðxÞ ¼ AT RðxÞ
i¼0
with T
A ¼ ½a0 ; a1 ; . . . ; aN 1 ;
ð6Þ T
RðxÞ ¼ ½R0 ðxÞ; R1 ðxÞ; . . . ; RN 1 ðxÞ :
ð7Þ
2.3. Operational matrix of derivative The derivative of the vector RðxÞ defined in Eq. (7) can be expressed by R0 ðxÞ ¼
dRðxÞ ’ DRðxÞ; dx
ð8Þ
where D is N N operational matrix for derivative. Differentiating Eq. (4) we get: R0nþ1 ðxÞ
x1 ¼ Rn ðxÞ þ 2 R0n ðxÞ R0n1 ðxÞ; xþ1 ðx þ 1Þ2 4
n P 1:
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
897
By using 2
3 1 ¼ R0 ðxÞ R1 ðxÞ þ R2 ðxÞ ¼ R01 ðxÞ; 4 4 ðx þ 1Þ 2
the elements dij , of the matrix D can be obtained from 9 0 R0nþ1 ðxÞ ¼ 2ðR1 ðxÞ Rn ðxÞÞ R0n1 ðxÞ; n P 1; > = R01 ðxÞ ¼ 34 R0 ðxÞ R1 ðxÞ þ 14 R2 ðxÞ; > ; R00 ðxÞ ¼ 0:
ð9Þ
The general form of the matrix D is a lower-Heisenberg matrix. The matrix D can be expressed as D ¼ D1 þ D2 , where D1 is a tridiagonal matrix which is obtained from 7 i i; i; D1 ¼ diag: ; i ¼ 0; . . . ; N 1 4 4 and the dij elements of matrix D2 are obtained from d10 ¼ 1 and 2; j P i 1; dij ¼ kicj ; j < i 1;
ð10Þ
where k ¼ ð1Þiþjþ1 , c0 ¼ 1 and cj ¼ 2 for j > 1. For N ¼ 6 we have 1 0 0 0 0 0 0 0 B 3=4 1 1=4 0 0 0 C C B B 2 7=2 2 1=2 0 0 C C: B B 3 6 21=4 3 3=4 0 C C B @ 4 8 8 7 4 1 A 5 10 10 10 35=4 5 2.4. The product operational matrix The following property of the product of two rational Chebyshev function vectors will also be used. Since Rm ðxÞRn ðxÞ ¼ 12½Rmþn ðxÞ þ Rjnmj ðxÞ; if we truncate the terms higher than RN 1 , then we have ~ T ðxÞ; RðxÞRT ðxÞA ¼ AR
ð11Þ
~ is a N N product operational matrix for the vector A and is given by where A
898
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
0 B B B B B ~¼B A B B B B @
a0
a1
1 a 2 1 1 a 2 2
a0 þ 12 a2 1 ða1 þ a3 Þ 2
1 ða1 2
1 a 2 3
1 ða2 2
1 ða1 2
.. . 1 a 2 N 1
þ a4 Þ .. .
a2
a3
þ a3 Þ a0 þ 12 a4
1 a 2 N 2
þ a5 Þ .. .
1 a 2 N 3
1 ða2 2 1 ða1 2
þ a4 Þ þ a5 Þ
a0 þ 12 a6 .. .
1 a 2 N 4
aN 1
1
C C C C C C: 1 a 2 N 4 C C .. C . C A a0 1 a 2 N 2 1 a 2 N 3
3. Converting Volterra’s population model to a nonlinear ODE equation In this section we convert VolterraÕs population model in Eq. (1) to an equivalent nonlinear ordinary differential equation. Let Z x yðxÞ ¼ uðtÞ dt: ð12Þ 0
This leads to y 0 ðxÞ ¼ uðxÞ; 00
ð13Þ
0
y ðxÞ ¼ u ðxÞ:
ð14Þ
Inserting Eqs. (12)–(14) into Eq. (1) yields the nonlinear differential equation jy 00 ðxÞ ¼ y 0 ðxÞ ðy 0 ðxÞÞ2 yðxÞy 0 ðxÞ
ð15Þ
with the initial conditions yð0Þ ¼ 0;
ð16Þ
0
y ð0Þ ¼ u0 ;
ð17Þ
obtained by using Eqs. (12) and (13) respectively.
4. Solution of nonlinear ODE by rational Chebyshev functions To solve Eq. (15) with initial conditions in Eqs. (16) and (17), we express yðxÞ, and y ðjÞ ðxÞ, j ¼ 1; 2, by rational Chebyshev functions as yðxÞ ¼
N 1 X
ai Ri ðxÞ ¼ AT RðxÞ;
ð18Þ
i¼0
y ðjÞ ðxÞ ¼
N 1 X i¼0
ðjÞ
ai Ri ðxÞ ¼ AT Dj RðxÞ;
j ¼ 1; 2;
ð19Þ
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
899
where Dj is the jth power of the matrix D given in Eq. (8). Using Eqs. (18) and (19) we get ~ y 0 ðxÞyðxÞ ¼ AT DRðxÞRT ðxÞA ¼ AT DARðxÞ;
ð20Þ
2 y 0 ðxÞ ¼ AT DRðxÞRT ðxÞDT A ¼ AT DRðxÞRT ðxÞF ¼ AT D Fe RðxÞ;
ð21Þ
~ and Fe are calculated similarly to Eq. (11). where F ¼ DT A and the matrices A Using Eqs. (18)–(21) the residual ResðxÞ for Eq. (15) can be written as i h ~ RðxÞ: ResðxÞ ¼ jAT D2 AT D þ AT D Fe þ AT DA As in a typical tau method [6] we generate ðN 2Þ algebraic equations by applying Z þ1 hResðxÞ; Rk ðxÞi ¼ ResðxÞRk ðxÞwðxÞ dx ¼ 0; k ¼ 0; 1; . . . ; N 3: 0
Using Eqs. (16)–(19) we get yð0Þ ¼ AT Rð0Þ ¼ 0;
ð22Þ
y 0 ð0Þ ¼ AT DRð0Þ ¼ u0 :
ð23Þ
Eqs. (22) and (23) generates two set of algebraic equations. Consequently, the unknown coefficients of the vector A in Eq. (18) can be calculated.
5. Illustrative example We applied the method presented in this paper to examine the mathematical structure of uðtÞ. In particular, we seek to study the rapid growth along the logistic curve that will reach a peak, then followed by the slow exponential decay where uðtÞ ! 0 as t ! 1. The mathematical behavior so defined was introduced by Scudo [2] and justified by Small [3] by using singular perturbation methods for the inner and outer solutions. Further, these properties were also confirmed by TeBeest [4] upon using a phase-plane analysis and [1] by using Pade approximations. We applied the method presented in this paper and solved Eq. (15) for uð0Þ ¼ 0:1 and j ¼ 0:02, 0.04, 0.1, 0.2, and 0.5 with N ¼ 10 and then evaluated umax and tcritical , which also are evaluated in [1]. In Table 1, the resulting values using the present method with N ¼ 10 together with the results obtained in [1] and the exact values obtained by using j umax ¼ 1 þ jln ; 1 þ j u0 evaluated in [4] are presented.
900
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
Table 1 A comparison between method in [1] and the present method (N ¼ 10) with the exact values for umax j
Critical t
Method in [1]
Present method
Exact umax
0.02 0.04 0.1 0.2 0.5
0.111845 0.210246 0.464476 0.816858 1.626711
0.903838 0.861240 0.765113 0.657912 0.485282
0.923463 0.873708 0.769734 0.659045 0.485188
0.923447 0.873719 0.769741 0.659050 0.485190
6. Conclusion The operational matrices of derivative and product of rational Chebyshev functions together with the tau method are used for solving the Volterra model for population growth of a species within a closed system. The stability and convergence of the tau approximations (see [7]) make this approach very attractive and contributed to the good agreement between approximate and exact values for umax in the numerical example.
References [1] A.M. Wazwaz, Analytical approximation and Pade approximation for VolterraÕs population model, Appl. Math. Comput. 100 (1999) 13–25. [2] F.M. Scudo, Volterra and theoretical ecology, Theoret. Popul. Biol. 2 (1971) 1–23. [3] R.D. Small, Population Growth in a Closed System and Mathematical Modelling, in: Classroom Notes in Applied Mathematics, SIAM, Philadelphia, PA, 1989, pp. 317–320. [4] K.G. TeBeest, Numerical and analytical solutions of VolterraÕs population model, SIAM Rev. 39 (1997) 484–493. [5] C. Lanczos, Applied Analysis, Englewood Cliffs, Prentice-Hall, NJ, 1956. [6] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamic, Springer, New York, 1988. [7] D. Gottlieb, M. Hussaini, S. Orszg, Theory and Applications of Spectral Methods in Spectral Methods for Partial Differential Equations, SIAM, Philadelphia, 1984. [8] B.Y. Guo, J. Shen, Z.Q. Wang, Chebyshev rational spectral and pseudospectral methods on a semi-infinite interval, Int. J. Numer. Math. Engng. 53 (2002) 65–84.
Rational Chebyshev tau method for solving VolterraÕs population model K. Parand a, M. Razzaghi a
a,b,*
Department of Applied Mathematics, Amirkabir University of Technology, Tehran 15914, Iran b Department of Mathematics and Statistics, Mississippi State University, Mississippi State, MS 39762, USA
Abstract An approximate method for solving VolterraÕs population model for population growth of a species in a closed system is proposed. VolterraÕs model is a nonlinear integro-differential equation where the integral term represents the effect of toxin. The approach is based on a rational Chebyshev tau method. The VolterraÕs population model is first converted to a nonlinear ordinary differential equation. The operational matrices of derivative and product of rational Chebyshev functions are presented. These matrices together with the tau method are then utilized to reduce the solution of the VolterraÕs model to the solution of a system of algebraic equations. Illustrative examples are included to demonstrate the validity and applicability of the technique and a comparison is made with existing results. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Rational Chebyshev; Tau method; VolterraÕs population model
1. Introduction The Volterra model for population growth of a species within a closed system is given in [1,2] as Z t du j ¼ u u2 u uðxÞ dx; uð0Þ ¼ u0 ; ð1Þ dt 0
*
Corresponding author. E-mail address: [email protected] (M. Razzaghi).
0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.09.006
894
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
where uðtÞ is the scaled population of identical individuals at a time t, and j ¼ c=ab is a prescribed nondimensional parameter. Moreover, a > 0 is the birth rate coefficient, b > 0 is the crowding coefficient, and c > 0 is the toxicity coefficient. The coefficient c indicates the essential behavior of the population evolution before its level falls to zero in the long term. One may show that the only equilibrium solution of (1) is the trivial solution uðtÞ ¼ 0. Furthermore, the analytical solution [3] Z t Z s 1 uðtÞ ¼ u0 exp 1 uðsÞ uðxÞ dx ds ð2Þ j 0 0 shows that uðtÞ > 0 for all t if u0 > 0. The solution of Eq. (1) has been of considerable concern. Although a closed form solution has been achieved in [3,4], but it was formally shown that the closed form solution cannot lead to any insight into the behavior of the population evolution [1]. Some approximate and numerical solutions for VolterraÕs population model have been reported. In Ref. [2], the successive approximations method was offered for the solution of Eq. (1), but was not implemented. In [3], the singular perturbation method for VolterraÕs population model is considered. It is shown in [3] that for the case j 1, where populations are weakly sensitive to toxins, a rapid rise occurs along the logistic curve that will reach a peak and be followed by a slow exponential decay. And, for j large, where populations are strongly sensitive to toxins, the solution is proportional to sech2 ðtÞ. In this case the solution uðtÞ has a smaller amplitude compared to the amplitude of uðtÞ for the case j 1. In [4], three numerical algorithms namely the Euler method, the modified Euler method and the fourth order Runge–Kutta method for the solution of Eq. (1) are obtained. Moreover, a phase-plane analysis are implemented. In [4], the numerical results are correlated to give an insight of the problem and its solution without using perturbation techniques. However, the performance of the traditional numerical techniques is well known in that it provides grid points only, and in addition, it requires a large amounts of calculations. In [4], 20 iteration steps is carried out to obtain the grid points. In [1], the series solution method and the decomposition method are implemented independently to Eq. (1) and to a related ordinary differential equation (ODE). Furthermore, the Pade approximations are used in the analysis to capture the essential behavior of the population uðtÞ of identical individuals. In the present paper, we introduce a new computational method for solving Eq. (1). The VolterraÕs population model is first converted to an equivalent nonlinear ODE, the solution of which is then approximated by a rational Chebyshev functions with unknown coefficients. The operational matrices of derivative and product of rational Chebyshev functions are given. These matrices together with the tau method are then utilized to evaluate the un-
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
895
known coefficients and find approximate solutions for uðtÞ. The tau method was invented by Lanczos in 1938 [5]. The method is based on expanding the required approximate solution as the elements of a complete set of orthogonal functions. In tau method, unlike the Galerkin approximation, the expansion functions are not required to satisfy the boundary constraint individually [6,7]. The sections of this paper are organized as follows: in Section 2 we describe the formulation of rational Chebyshev functions required for our subsequent development. Section 3 converts VolterraÕs population to a nonlinear ODE. Section 4 summarizes the application of the rational Chebyshev tau method to the solution of nonlinear ODE for the Volterra population model. As a result a set of nonlinear algebraic equations is formed, and a solution of the considered ODE is introduced. In Section 5 the proposed method is applied to the numerical example, and a comparison is made with existing methods in the literature.
2. Properties of rational Chebyshev functions 2.1. Rational Chebyshev functions The well-known Chebyshev polynomials are orthogonal inffi the interval pffiffiffiffiffiffiffiffiffiffiffiffi ½1; 1 with respect to the weight function wðyÞ ¼ 1= 1 y 2 and can be determined with the aid of the following recurrence formula: T0 ðyÞ ¼ 1;
T1 ðyÞ ¼ y;
Tnþ1 ðyÞ ¼ 2yTn ðyÞ Tn1 ðyÞ;
n P 1:
The rational Chebyshev (RC) functions are defined in [8] by x1 Rn ðxÞ ¼ Tn : xþ1 Thus RC functions satisfy: R0 ðxÞ ¼ 1;
R1 ðxÞ ¼
Rnþ1 ðxÞ ¼ 2
x1 ; xþ1
x1 Rn ðxÞ Rn1 ðxÞ; xþ1
ð3Þ n P 1:
ð4Þ
Rational Chebyshev pffiffiffi functions are orthogonal with respect to the weight function wðxÞ ¼ 1= xðx þ 1Þ on the interval ½0; þ1Þ, with the orthogonality property: Z 1 cm p dnm Rn ðxÞRm ðxÞwðxÞ dx ¼ 2 0
896
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
with cm ¼
m ¼ 0; m 6¼ 1;
2; 1;
where dnm is the Kronecker function. 2.2. Function approximation A function f ðxÞ defined over the interval ½0; þ1Þ may be expanded as f ðxÞ ¼
1 X
ai Ri ðxÞ;
ð5Þ
i¼0
where ai ¼
2 ci p
Z
1
Ri ðxÞf ðxÞwðxÞ dx: 0
If f ðxÞ in Eq. (5) is truncated up to the N th terms, then Eq. (5) can be written as f ðxÞ ¼
N 1 X
ai Ri ðxÞ ¼ AT RðxÞ
i¼0
with T
A ¼ ½a0 ; a1 ; . . . ; aN 1 ;
ð6Þ T
RðxÞ ¼ ½R0 ðxÞ; R1 ðxÞ; . . . ; RN 1 ðxÞ :
ð7Þ
2.3. Operational matrix of derivative The derivative of the vector RðxÞ defined in Eq. (7) can be expressed by R0 ðxÞ ¼
dRðxÞ ’ DRðxÞ; dx
ð8Þ
where D is N N operational matrix for derivative. Differentiating Eq. (4) we get: R0nþ1 ðxÞ
x1 ¼ Rn ðxÞ þ 2 R0n ðxÞ R0n1 ðxÞ; xþ1 ðx þ 1Þ2 4
n P 1:
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
897
By using 2
3 1 ¼ R0 ðxÞ R1 ðxÞ þ R2 ðxÞ ¼ R01 ðxÞ; 4 4 ðx þ 1Þ 2
the elements dij , of the matrix D can be obtained from 9 0 R0nþ1 ðxÞ ¼ 2ðR1 ðxÞ Rn ðxÞÞ R0n1 ðxÞ; n P 1; > = R01 ðxÞ ¼ 34 R0 ðxÞ R1 ðxÞ þ 14 R2 ðxÞ; > ; R00 ðxÞ ¼ 0:
ð9Þ
The general form of the matrix D is a lower-Heisenberg matrix. The matrix D can be expressed as D ¼ D1 þ D2 , where D1 is a tridiagonal matrix which is obtained from 7 i i; i; D1 ¼ diag: ; i ¼ 0; . . . ; N 1 4 4 and the dij elements of matrix D2 are obtained from d10 ¼ 1 and 2; j P i 1; dij ¼ kicj ; j < i 1;
ð10Þ
where k ¼ ð1Þiþjþ1 , c0 ¼ 1 and cj ¼ 2 for j > 1. For N ¼ 6 we have 1 0 0 0 0 0 0 0 B 3=4 1 1=4 0 0 0 C C B B 2 7=2 2 1=2 0 0 C C: B B 3 6 21=4 3 3=4 0 C C B @ 4 8 8 7 4 1 A 5 10 10 10 35=4 5 2.4. The product operational matrix The following property of the product of two rational Chebyshev function vectors will also be used. Since Rm ðxÞRn ðxÞ ¼ 12½Rmþn ðxÞ þ Rjnmj ðxÞ; if we truncate the terms higher than RN 1 , then we have ~ T ðxÞ; RðxÞRT ðxÞA ¼ AR
ð11Þ
~ is a N N product operational matrix for the vector A and is given by where A
898
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
0 B B B B B ~¼B A B B B B @
a0
a1
1 a 2 1 1 a 2 2
a0 þ 12 a2 1 ða1 þ a3 Þ 2
1 ða1 2
1 a 2 3
1 ða2 2
1 ða1 2
.. . 1 a 2 N 1
þ a4 Þ .. .
a2
a3
þ a3 Þ a0 þ 12 a4
1 a 2 N 2
þ a5 Þ .. .
1 a 2 N 3
1 ða2 2 1 ða1 2
þ a4 Þ þ a5 Þ
a0 þ 12 a6 .. .
1 a 2 N 4
aN 1
1
C C C C C C: 1 a 2 N 4 C C .. C . C A a0 1 a 2 N 2 1 a 2 N 3
3. Converting Volterra’s population model to a nonlinear ODE equation In this section we convert VolterraÕs population model in Eq. (1) to an equivalent nonlinear ordinary differential equation. Let Z x yðxÞ ¼ uðtÞ dt: ð12Þ 0
This leads to y 0 ðxÞ ¼ uðxÞ; 00
ð13Þ
0
y ðxÞ ¼ u ðxÞ:
ð14Þ
Inserting Eqs. (12)–(14) into Eq. (1) yields the nonlinear differential equation jy 00 ðxÞ ¼ y 0 ðxÞ ðy 0 ðxÞÞ2 yðxÞy 0 ðxÞ
ð15Þ
with the initial conditions yð0Þ ¼ 0;
ð16Þ
0
y ð0Þ ¼ u0 ;
ð17Þ
obtained by using Eqs. (12) and (13) respectively.
4. Solution of nonlinear ODE by rational Chebyshev functions To solve Eq. (15) with initial conditions in Eqs. (16) and (17), we express yðxÞ, and y ðjÞ ðxÞ, j ¼ 1; 2, by rational Chebyshev functions as yðxÞ ¼
N 1 X
ai Ri ðxÞ ¼ AT RðxÞ;
ð18Þ
i¼0
y ðjÞ ðxÞ ¼
N 1 X i¼0
ðjÞ
ai Ri ðxÞ ¼ AT Dj RðxÞ;
j ¼ 1; 2;
ð19Þ
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
899
where Dj is the jth power of the matrix D given in Eq. (8). Using Eqs. (18) and (19) we get ~ y 0 ðxÞyðxÞ ¼ AT DRðxÞRT ðxÞA ¼ AT DARðxÞ;
ð20Þ
2 y 0 ðxÞ ¼ AT DRðxÞRT ðxÞDT A ¼ AT DRðxÞRT ðxÞF ¼ AT D Fe RðxÞ;
ð21Þ
~ and Fe are calculated similarly to Eq. (11). where F ¼ DT A and the matrices A Using Eqs. (18)–(21) the residual ResðxÞ for Eq. (15) can be written as i h ~ RðxÞ: ResðxÞ ¼ jAT D2 AT D þ AT D Fe þ AT DA As in a typical tau method [6] we generate ðN 2Þ algebraic equations by applying Z þ1 hResðxÞ; Rk ðxÞi ¼ ResðxÞRk ðxÞwðxÞ dx ¼ 0; k ¼ 0; 1; . . . ; N 3: 0
Using Eqs. (16)–(19) we get yð0Þ ¼ AT Rð0Þ ¼ 0;
ð22Þ
y 0 ð0Þ ¼ AT DRð0Þ ¼ u0 :
ð23Þ
Eqs. (22) and (23) generates two set of algebraic equations. Consequently, the unknown coefficients of the vector A in Eq. (18) can be calculated.
5. Illustrative example We applied the method presented in this paper to examine the mathematical structure of uðtÞ. In particular, we seek to study the rapid growth along the logistic curve that will reach a peak, then followed by the slow exponential decay where uðtÞ ! 0 as t ! 1. The mathematical behavior so defined was introduced by Scudo [2] and justified by Small [3] by using singular perturbation methods for the inner and outer solutions. Further, these properties were also confirmed by TeBeest [4] upon using a phase-plane analysis and [1] by using Pade approximations. We applied the method presented in this paper and solved Eq. (15) for uð0Þ ¼ 0:1 and j ¼ 0:02, 0.04, 0.1, 0.2, and 0.5 with N ¼ 10 and then evaluated umax and tcritical , which also are evaluated in [1]. In Table 1, the resulting values using the present method with N ¼ 10 together with the results obtained in [1] and the exact values obtained by using j umax ¼ 1 þ jln ; 1 þ j u0 evaluated in [4] are presented.
900
K. Parand, M. Razzaghi / Appl. Math. Comput. 149 (2004) 893–900
Table 1 A comparison between method in [1] and the present method (N ¼ 10) with the exact values for umax j
Critical t
Method in [1]
Present method
Exact umax
0.02 0.04 0.1 0.2 0.5
0.111845 0.210246 0.464476 0.816858 1.626711
0.903838 0.861240 0.765113 0.657912 0.485282
0.923463 0.873708 0.769734 0.659045 0.485188
0.923447 0.873719 0.769741 0.659050 0.485190
6. Conclusion The operational matrices of derivative and product of rational Chebyshev functions together with the tau method are used for solving the Volterra model for population growth of a species within a closed system. The stability and convergence of the tau approximations (see [7]) make this approach very attractive and contributed to the good agreement between approximate and exact values for umax in the numerical example.
References [1] A.M. Wazwaz, Analytical approximation and Pade approximation for VolterraÕs population model, Appl. Math. Comput. 100 (1999) 13–25. [2] F.M. Scudo, Volterra and theoretical ecology, Theoret. Popul. Biol. 2 (1971) 1–23. [3] R.D. Small, Population Growth in a Closed System and Mathematical Modelling, in: Classroom Notes in Applied Mathematics, SIAM, Philadelphia, PA, 1989, pp. 317–320. [4] K.G. TeBeest, Numerical and analytical solutions of VolterraÕs population model, SIAM Rev. 39 (1997) 484–493. [5] C. Lanczos, Applied Analysis, Englewood Cliffs, Prentice-Hall, NJ, 1956. [6] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamic, Springer, New York, 1988. [7] D. Gottlieb, M. Hussaini, S. Orszg, Theory and Applications of Spectral Methods in Spectral Methods for Partial Differential Equations, SIAM, Philadelphia, 1984. [8] B.Y. Guo, J. Shen, Z.Q. Wang, Chebyshev rational spectral and pseudospectral methods on a semi-infinite interval, Int. J. Numer. Math. Engng. 53 (2002) 65–84.