A Discretization Approach for the Nonlinear Fractional Logistic Equation

The present study aimed to develop and investigate the local discontinuous Galerkin method for the numerical solution of the fractional logistic differential equation, occurring in many biological and social science phenomena. The fractional derivative is described in the sense of Liouville-Caputo. Using the upwind numerical fluxes, the numerical stability of the method is proved in the L∞ norm. With the aid of the shifted Legendre polynomials, the weak form is reduced into a system of the algebraic equations to be solved in each subinterval. Furthermore, to handle the nonlinear term, the technique of product approximation is utilized. The utility of the present discretization technique and some well-known standard schemes is checked through numerical calculations on a range of linear and nonlinear problems with analytical solutions.


Introduction
In studies of elementary population dynamics the simplest model for the growth of a population is known as rate equation and structured by Malthus in (1798) [1] $ & % dMptq dt " r Mptq, t ą 0, where Mptq denotes the population at time t, the non-zero parameter r equals to r " β´α, where β and α are the per capita birth and death rates respectively. Here, M 0 is the population at time t " 0. The exact analytical solution of Malthus population model (1) is explained the constant population growth rate Mptq " M 0 e rt . The Maithusian grow model is unrealistic over long times due to the fact that the solution of the rate equation is not included two main factors such as spread of diseases and the limitation on food supply. To model the effects of these factors in a population model, the logistic equation was considered by P. R. Verhulst in 1838 [2] dNptq dt " rNptqˆ1´N ptq K˙, " X 0 . ( where X 0 " Mp0q{M max . The exact solution of this equation can be easily obtained as Xptq " X 0 X 0`p 1´X 0 qe´σ t . In the last decades, many efforts have been devoted to extend the integer-order models to the corresponding fractional-order models, which are more descriptive and can provide a powerful and valuable instrument for the explanation of hereditary and memory properties of several materials and process [3,4]. Replacing the classical derivative operator in (2) by a fractional one, the fractional logistic equation will be obtained. This model of population growth has been found applications in numerous disciplines of science and engineering. For instance, the growth of tumors in medicine [5] can be modelled as the fractional logistic equation (FLE). In addition, the milstone of various important mathematical models is based on the fractional logistic equation such as two models in Radar signals [6] and electroanalytical chemistry [7]. Several variations of the population growth model have been studied in the literature [8]. In the present study, we are going to investigate the following logistic population model of fractional order in the form $ & % LC a D ν t Xptq " σ Xptq´1´Xptq¯": σ Xptq gpXptqq, t ą 0, Xp0q " X 0 , where the symbol LC a D ν t denotes the fractional derivative operator of Liouville-Caputo type and ν P p0, 1s. It should be emphasize that in (3) we have used the function gpsq " 1´s, which corresponds to the nonlinear logistic equation. However, to address the linear counterpart of this equation we also consider gpsq " 1. The issue of existence and the uniqueness of the solution of (3) is discussed in details in Reference [9].
It is known that for most fractional differential equations there is no possibility to find the exact solutions analytically. Consequently, exploring an approximate or numerical technique is of primary interest for such fractional equations. Many efforts have been made toward the exact analytical solution of the problem (3). The first one is proposed by West [10], which is based on the Carleman embedding technique. Later, it is shown that in Reference [11] the this analytical function is only very close to the numerical solutions of the FLE. The other analytical methods for the FLE include the fractional Taylor expansion method [12], a method based on Euler's numbers [13], and the varational iterative method [14]. Besides the analytical investigations, numerous computational approaches have been proposed for the nonlinear FLE. Let us mention the predictor-corrector approaches [9,15], the finite difference schemes [14,16], the spectral methods [17,18], the Bessel collocation method [19], the Chebyshev wavelet method [20], the Laguerre collocation method [21], and the fractional spline collocation method [22].
Many other numerical and approximation methods as well as computational approaches have been developed and applied for the FDEs which are based upon various closely-related models of real-world problems. For example, Baleanu et al. [23] made use of a Chebyshev spectral method based on operational matrices, a remarkable survey of numerical methods can be found in [24], a study of the fractional-order Bessel, Chelyshkov, and Legendre collocation schemes for the fractional Riccati equation was presented in [25], an operational matrix of fractional-order derivatives of Fibonacci polynomials was developed in [26], an introductory overview and recent developments involving FDEs was presented in [27], efficiency of the spectral collocation method in the dynamic simulation of the fractional-order epidemiological model of the Ebola virus was investigated in [28], the Jacobi collocation method and a spectral tau method based on shifted second-kind Chebyshev polynomilas for the approximate solution of some families of the fractional-order Riccati differential equations were discussed in [29,30], computational approaches to FDEs for the biological population model were discussed in [31], the generalized Chebyshev and Bessel colllocation approaches for fractional BVPs and multi-order FDEs were considered in [32,33], and a general wavelet quasi-linearization method for solving fractional-order population growth model was developed and applied in [34].
In this work, we take a further step towards proposing a numerical method for solving the FLE. We utilize a discontinuous finite element approach, i.e, the local discontinuous Galerkin (LDG) discretization approach for the FLE (3). To apply the LDG scheme, we must rewrite a given FDEs as a system of first-order ordinary differential equations (ODEs) with together a fractional integral. Hence, the discontinuous Galerkin (DG) method is employed to discretize the resulting system as well as the fractional integral. The first DG method was introduced by Reed and Hill [35] in 1973 for numerically solving neutron transport, that is, a time-independent linear hyperbolic equation. Since then the DG schemes have been well implemented for the classical ODEs was started by the work [36]. DG schemes as a subclass of finite element methods (FEMs) allow us to exploit discontinuous discrete basis functions. These local basis functions are usually selected as piecewise polynomials. Exploiting completely discontinuous basis functions offers great opportunities compared to traditional FEMs when used to discretize differential equations. In summary, the main gains of the DG methods are in terms of flexibility, accuracy as well as parallelizability, see cf. Reference [37].
To the best of our knowledge, the LDG approaches for the ODEs of fractiona-order including one-term and multi-terms were first discussed in Reference [38] and then have been applied to many model problems [39][40][41]. It is worth mentioning that the success of LDG methods is based on the designing of appropriate numerical fluxes at the interface elements. In this work, we utilize the upwind numerical flux as natural choice for the FLE. By choosing the upwind fluxes we arel able to prove the numerical stability of the LDG scheme.
The rest of this paper is organized as follows. In the next Section, we review some fractional calculus preliminaries and state some of their properties that will be used later on. The formulation of the LDG scheme for the logistic equation is established in Section 3. Hence, the algebraic form of the LDG scheme is obtained with the aid of shifted Legendre basis functions. The technique of product approximation is also applied to deal with the nonlinear term in the weak formulation. In Section 4 we establish the numerical stability of the scheme in the linear case and a discussion about the error estimation is made. In Section 5, the applicability and utility of the present numerical schemes are verified by performing several simulations on two linear and nonlinear population growth and logistic model problems. Finally, a conclusion is drawn in Section 6.

Fractional Calculus
Now, we present some fundamental and mathematical preliminaries of the fractional calculus theory to be utilized in our subsequent sections, see References [3,4,27].

Definition 1.
Let ν ě 0 is given. The (left) Riemann-Liouville fractional integral operator of order ν is given by The integral operator I ν has many properties. Among others, we make use of the followings (1) I ν I β f ptq " I ν`β f ptq, (2) I ν pc 1 f ptq`c 2 gptqq " c 1 I ν f ptq`c 2 I ν gptq, c 1 , c 2 P R, The corresponding definition of the right Riemann-Liouville fractional integral on the interval rt, bs instead of ra, ts is given by 1 Γpνq ν " m, m P N.
Here, we have used the ceiling and floor functions rνs, tνu respectively. It should be noted that, two operators I ν and D ν are related through the following expression

Discretized LDG Formulation
In order to formulate the LDG method for the logistic equation in (3), some basic notations will first be introduced.
Let us consider (3) on L " p0, Tq for some given T ą 0. To rewrite (3) as a first-order system, we introduce two new variables z 0 ptq " Xptq and z 1 ptq " dXptq dt and use the relation (6) to get being ν P p0, 1s and t P L. By ∆ we denote a partitioning of the interval L into into J subintervals L l " pt l´1 , t l q for l " 1, . . . , J. The grid points of ∆ will be denoted as By h l we mean the length of each L l , that is, h l " t l´tl´1 for l " 1, 2, . . . , N. The maximum length of these element is taken as h :" max J l"1 h l . We associate the mesh ∆ with the broken Sobolev spaces H 1 ∆ " tw : L Ñ Rˇˇw| L l P H 1 pL l q, l " 1, 2, . . . , Ju. and S ∆ " tw : L Ñ Rˇˇw| L l P L 2 pL l q, l " 1, 2, . . . , Ju, By using these function spaces, let assume that the solutions of system (7) belong to corresponding spaces´z It should be noted that the elements of space H 1 ∆ may be discontinuous in t at discrete time level t 1 . In this respect, at the mesh grid points, defining the left-sided as well as the right-sided limits of a function w is necessary, where w : L Ñ R is a piecewise continuous function. By wń and wǹ , we let the left-and right-sided limits of w at t l wl " w`pt l q " wptl q :" lim sÑ0`w pt n`s q, wĺ " w´pt l q " wptĺ q :" lim tÑ0´w pt n`s q.
For any positive integer number r, we denote by P r pL l q the space of polynomials of degree less or equal than r on the element L l P ∆. We then let the approximate solutions z 0 ptq, z 1 ptq belong to a subspace V prq Ă H 1 ∆ , which is a finite dimensional space. This subspace is defined as the space of discontinuous and piecewise polynomial functions V prq " tw : L Ñ Rˇˇw| L l P P r pL l q, l " 1, 2, . . . , Ju.
We further define Z 0 ptq and Z 1 ptq as the DG approximations to the exact solutions z 0 ptq and z 1 ptq of the system (7) respectively on the element L l . Below, we make use of the following notations pw, vq l :" xw, wy l .
For obtaining the weak DG formulation, we first multiply the first equation in (7) by a test function w 0 P V prq and integrate over L l . By applying the integrating by parts we get Hence, the second integral equation in (7) will multiplied by a test function w 1 P V prq and integrate over L l . To advance the solution in time, we replace Z 0 ptl´1q by the upwind flux Z 0 ptĺ´1q in (8). Thus, the discrete formulation for finding Z 0 , Z 1 P V prq takes the following form for all w 0 , w 1 P V prq , and l " 1, 2, . . . , J It should be noted that, to start computations on the first element L 1 " pt 0 , t 1 q we use the given initial condition Z 0 pt0 q " X 0 . Hence, by utilizing the upwind flux as natural choice, we are able to solve the resultant equations element by element on each subinterval L l for l " 1, 2, . . . , J. In each element, we just need to invert a local matrix of size pr`1qˆpr`1q in place of a global matrix of size Jpr`1qˆJpr`1q.

Algebraic Formulation
Since the functions in V prq may be discontinuous across interfaces of the element, various local bases can be selected for finite element approximation in (9). Let us choose a basis in the space P r pL l q Entropy 2020, 22, 1328 6 of 17 formed by functions φ l 0 , φ l 1 , . . . , φ l r . Thus the numerical approximations Z 0 of z 0 and Z 1 of z 1 in every element L l can be expressed as Here, the coefficients α l i , β l i , i " 0, . . . , r denote the degrees of freedom to be sought in each L l ,l " 1, . . . , J. To proceed, we take the test functions in each element L l in the form w j " φ l j ptq for j " 0, 1, . . . , r and l " 0, 1, . . . , J. Now, by specifying the basis functions as we done below, the discrete LDG formulation (9) is reduced to a algebraic system of equations.
For practical implementation of the LDG scheme (9) for the FLE (3), we use the set of orthogonal Legendre polynomials for the space V prq . Let us recall that, the i'th degree Legendre polynomials P i psq can be generated by the well-known Rodriguez formula The Legendre polynomials satisfy the following relations [17] ż 1 where δ ij denotes the Kronecker delta. The first property shows that these set of orthogonal polynomials are orthogonal with respect to weighting function wptq " 1 on p´1, 1q. The Legendre polynomial P i psq of degree i can be explicitly expressed as follows where M i " i{2 or pi´1q{2, whichever is an integer. Due to the fact that these polynomials are orthogonal on r´1, 1s, we map them onto the element L l by using the following change of variable Let the resultant shifted Legendre polynomials denoted by L i ptq. Thus, the explicit form of L i ptq of degree i takes the form By means of the binomial formula, one can further simplify the last expression as follows where the coefficients C ikm are defined as Now, we choose φ l i ptq " L i ptq in (10) for l " 1, 2, . . . , J, where L i is the shifted Legendre polynomial of degree i in t defined in (12). With this transformation, the unknown values α l i , β l i in (10) can be interpreted as the Legendre coefficients of the expansion of Z 0 , Z 1 . Hence, by the virtue of the Legendre properties (11) and inserting (10) into the discrete formulation (9) we have for l " 1, . . . , J as for j " 0, . . . , r. To proceed, we need to deal with two main difficulties involving the integral and nonlinear terms that appear in (13). To tackle the integral term, the properties (1)-(3) of fractional integration in 2 is used to obtain Next, the explicit form (12) is utilized for L j ptq and then 0 I p1´νq t L i ptq will be inserted into the inner product. Now, by integration over L l we obtain with the coefficients C 2 ikmjk 1 m 1 :" C 1 ikm C jk 1 m 1 {pm`m 1`2´ν q.
Th nonlinear term in (13) can be computed using the Legendre polynomials. For instance, if r " 1 we may write it as a product of two vectors for j " 0, 1. Therefore, it is not a difficult task to calculate n DC j by direct computation (D.C.) using the shifted Legendre polynomials on each L l for different j. Of course one may exploit the symbolic toolbox in Matlab to facilitate the process of integration of these polynomials. Alternatively, to handle the nonlinear term in (13), the product approximation (P.A.) technique [42] is used in the following manner This technique enables us to write the nonlinear part as n PA i,j :"´Z 2 0 ptq, L j ptq¯l " Now, it suffices to calculate the two first terms in (13). To this end, we compute the elements of the mass matrix as m i,j :"´L i ptq, L j ptq¯l " Finally, the entries of the stiffness matrix need to be calculated. In the new coordinate, we recursively employ the Legendre property (11b) to derive By applying the orthogonality relation (11a) to the preceding equation and then simplifying the involved integral in s i,j , we finally get 2, if i ą j and pi`jq is even, 0, otherwise.
Note in (18)  Clearly, the value of Z 0 ptĺ´1q is already known from the preceding time interval L l´1 . Obviously this value at the first time interval is computed as X 0 , which known as the initial condition. Also, the obtained system (18) is a nonlinear algebraic system of equations have to be solved in each L l for l " 1, . . . , J. This system can be solved for example, via Newton type schemes. It is known that this method converges quadratically whenever the approximation is close to the actual solution of the given nonlinear system. Using the D.C. approach, we also arrive at a nonlinear system of equation in the general form F F Fpα α α l , β β β l q " 0 0 0 to be solved in each interval L l . As we show in the numerical experiments, this approach is more accurate than the corresponding P.A. approach.

Numerical Stability and Error Estimates
Now, we are going to establish the stability of proposed LDG scheme when applied to the logistic equation in the linear case by considering gptq " 1 in (3). In this case we have Without loss of generality, let us assume that σ ă 0. The numerical scheme of (19) is to find Z 0 , Z 1 P V prq such that for all w 0 , w 1 P V prq , and l " 1, 2, . . . , J. Let us state the next lemma, which based on the semigroup properties of fractional integral operators and will be used below, a proof of which can be found in Reference [38].

Lemma 1.
Suppose that ν P p0, 1q, then we have Let us assume that r Z 0 , r Z 1 P V prq be the approximate solutions of Z 0 , Z 1 respectively. Now, the numerical errors are defined as E X i :" r Z i´Zi for i " 0, 1. It can be seen that r Z 0 and r Z 1 both satisfy (20). If we subtract Equation (20) from the same equations with r Z 0 and r Z 1 , the following error equations will be obtained which holds for all w 0 , w 1 P V prq . Taking w 0 " E X 0 and w 1 " E X 1 in (21) followed by collecting these two equations, we conclude that To deal with the third term, we utilize the identity´u, du dt¯l " pu 2 ptĺ q´u 2 ptl´1qq{2 with u " E X 0 . Hence, we multiply the preceding equation by two. Adding and subtracting E 2 X 0 ptĺ´1q to the modified equation and rearranging the terms to obtaiń By summing over elements for l " 1, . . . , J, we get By using Lemma 1, we have established the following stability of the LDG in the L 8 norm for (20) (see also References [38,40]: Lemma 2. We have the following L 8 stability of the LDG scheme (20) and for the numerical errors hold We close this section by pointing out some facts about the order of convergence of the proposed LDG scheme. In Reference [38] it is shown that the solution can be calculated with optimal order of convergence pr`1q in the L 2 norm. In this work the mechanism of superconvergence is also discussed. The authors observed the superconvergence of order pr`1q`mintr, νu at downwind point of each element.

Numerical Results and Discussions
In this section, we present some results of computations using the proposed LDG scheme described in the preceding sections to test their accuracy and efficiency when applied to the logistic equation. To assess the accuracy of the present numerical algorithms, we calculate the difference between the true exact and numerical solutions whenever the exact solution is available. For this purpose, we also consider a linear fractional population model and then we solve the fractional logistic equation numerically.
In order to asses the numerical scheme more qualitatively, by EOC we denote the estimated order of convergence calculated through defining where E a phq is the absolute error corresponding to the step-size h. Moreover, to test the validity and accuracy of proposed LDG method and to make a comparison between our numerical model results with the results of other existing methods, we employ the predictor-corrector PECE method of Adams-Bashforth-Moulton type considered in Reference [43] as well as the implicit product integration of trapezoidal type described in Reference [24]. All experimental computations have been done by using MATLAB R2017a.

Linear Model
In this section, we consider a linear test problem to show the effectiveness of the proposed LDG approach. For this purpose, we consider the fractional population growth where 0 ă ν ď 1 and σ ą 0. This model problem is previously studied in Reference [22] and can be considered as a generalization of the Malthusian model (1) to the fractional-order derivative. By the aid of the Laplace transform, the exact analytical solution of the initial-value problem can be obtained in terms of well-knwon Mittag-Leffler function [10] Xptq " X 0 E ν pσ ν t ν q, E ν pzq " Note that by taking ν " 1 the exact solution becomes Xptq " X 0 e σ t .
To start computation, we take σ " 1 for simplicity and set X 0 " 3{4. By considering ν " 1 and J " 1, the approximate solutions for r " 3, 6, and r " 9 on the interval 0 ď t ď 2 are obtained as follows These approximations together with the corresponding absolute errors are depicted in Figure 1. Clearly, as r increased, more accurate results will be obtained. Note, in all cases, the step size is taken as h " 2. Moreover, we emphasize that numerical solutions for this model problem based on the fractional spline collocation scheme have been proposed in Reference [22] with achieved absolute errors larger than 1ˆ10´4, see Figure 2 in this paper. The parameters used in this approach related to ν " 1 were M 1 " 2 6 , 2 7 , 2 8 , N 1 " 37, 69, 133, which obviously are much more greater than our used parameters. Absolute-error r = 3 r = 6 r = 9 Figure 1. The approximated LDG with exact solutions (left) and the corresponding absolute errors (right) for J " 1, ν " 1, σ " 1, X 0 " 0.75, and different r " 3, 6, 9. Additionally, to justify our numerical model results, a comparison in Table 1 has been performed between the previous work on PECE [15,43] in terms of the number of (sub)intervals J is used in the computation. In this comparison, we compute the numerical solutions corresponding to Xp2q as well as absolute errors |Xp2q´Z 0 p2q| in these methods via different values of J " 2 i for i " 0, 1, . . . 7. For our LDG method we take r " 2 and ν " 1. The last column in each method reports the corresponding EOC. The exact value of Xp2q up to 30 digits is Xp2q " 5.54179207419798736111715697916. Table 1. Comparison of absolute errors in LDG with r " 2 and PECE for different number of interval J and ν " 1. Numbers in bold show that the correct digits are obtained by the LDG. The observed EOC seen for PECE in Table 1 is approximately 2 as was proved in Reference [43]. However, the spuperconvergence EOC about 5 («2r + 1) is clearly achieved for our results. This comparison indicates the thoroughness of the proposed method. The numerical solutions for various values of ν " 0.65, 0.75, 0.85, 0.95 using r " 5 and J " 1 are depicted in Figure 2, left plot. In all plots, the exact solutions are indicated by a solid line while the numerical counterpart are visualized by (coloured) dotted, dashed, and dash-dotted curves. Note that the computational domain is r0, 1s, which implies that the time step is h " 1. It can be seen from Figure 2 that the numerical solution obtained by the present LDG scheme has a good accuracy even using a relatively large time step and a low degree of the approximating polynomials. Furthermore, an appropriate choice of these computational parameters can improve the approximation accuracy. Finally, for the linear model problem (23), we investigate the standard L1 approximation method [44] and its variant known as the fast L1 method [45]. To implement these approaches, we use a uniform mesh with the step size h " 1{1000 on the interval r0, 1s. In the LDG scheme, we utilize J " 1 or h " 1 and r " 5 as the results shown in Figure 2. The numerical model results are presented in Table 2 for ν " 0.75 and ν " 0.5. For each ν, the corresponding exact solutions are also reported in the last column.

Nonlinear Model
We now consider the FLE (3) on r0, 1s with the initial condition given by X 0 " 1{2 and the parameter σ " 1{2. Using ν " 1, the analytical exact solution of the logistic equation is given by The simulation results for this example can be found in Figures 3 and 4 for the number of elements equals to J " 5 and the polynomial degree r " 2. In Figure 3, we take ν " 1 to compare the numerical results to the exact solution. Furthermore, we also use different approaches to treat the nonlinear term Entropy 2020, 22, 1328 13 of 17 in the weak formulation, that is, the D.C. and P.A., which are utilized to compute n DC j and n PA i,j in (15). As one can see that from Figure 3 that a slightly more accurate result is obtained by means of direct computation rather than product approximation, however, as mentioned it is more time-consuming. In order to observe the behaviour of numerical solutions more closely, a magnification of these solutions at t " 0.4 is done in Figure 3. The exact solution is depicted by a solid line.  In the next experiment, we plot the absolute errors when utilizing two approaches D.C. and P.A., as one observes in Figure 4. The computational parameters are the same as those applied for Figure 3. In Figure 4, the left plot corresponds to the D.C. and the right plot is when we use P.A. technique. Note that in all plots we have divided further each interval L l into ten subinterval uniformly to see the behaviour of the corresponding curves more precisely.  Let us interpret the numerical errors depicted in Figure 4. On the left picture in which the P.A. technique is used, the smallest errors are obtained at upwind points. Almost the same magnitude of errors is achieved at downwind points. On the contrary, on the right picture without using the P.A. this process is reversed. This implies that the minimum values of absolute errors are achieved at downwind points and there exist considerable difference between them and the errors obtained at upwind points in each L l . In the next experiments, we compare the numerical errors achieved at the final point T " 1.0, which is clearly a downwind point in the first approach.
In Tables 3 and 4, we summarize the numerical results related to Xp1q and its numerical approximation Z 0 p1q are obtained by the LDG procedure (9). Here, we use r " 1, 2 and a different choice of the number of grid points J " 1, 2, 4, 8 and 16 are utilized. In these tables, we further compare the performance of two different D.C. and P.A. approaches. All calculations are shown with 10 decimal places of accuracy. In the last column of each table, the estimated order of convergence (EOC) is given. The exact value is Xp1q " 0.622459331201855. Table 3. Comparison of absolute errors in LDG with r " 1 using P.A. and D.C. for different number of interval J and ν " 1. Numbers in bold show that the correct digits are obtained by the LDG.

P.A.
D.C. It can be seen from Tables 3 and 4 that using r " 1 and r " 2 in the D.C. approach, the results are accurate respectively to 6 and 10 decimal places for only J " 4 intervals. In other words, achieving an order of accuracy equal to 3 and 5 is possible if one uses the LDG scheme with r " 1, 2 degree of polynomials and for a small number of elements. These EOC are also confirmed the superconvergence order at downwind points previously reported in Reference [38]. Note that by utilizing the P.A. technique, the obtained EOC is equal to 2. We emphasize also that using the scheme PECE for the nonlinear logistic equation the EOC at most 2 will be achieved and of course a larger number of intervals J is required. In the next plot, we examine the behaviour of the absolute errors in the log scale for various polynomial degrees as well as with respect to the number of elements J, see Figure 5. In the next experiment we show the impact of the fractional derivative on the approximated obtained solutions. In Figure 6 we present the approximated solutions at J " 4, r " 3 with different values of the fractional derivatives ν " 0.65, 0.75, 0.85, 0.95 as well as ν " 1.0. In these plots, we also compare the performance of two P.A. and D.C. approaches for these values of ν. In each case, for ν " 1.0 the exact solution is also shown by a solid line. To justify our computed results, the implicit product-integration of trapezoidal (IPIT) rule with the step size h " 1{256 is used [24].
From both depictions in Figure 6, one can observe that the numerical solutions for ν P p0, 1q are approaching to the solutions correspond to ν " 1 for which the exact solution is known. Of course, more reliable results is obtained through the D.C. as previously tested for ν " 1 in Tables 3 and 4.

Conclusions
In this work, an approximation algorithm based on the LDG scheme is developed for the fractional-order logistic equation occurring in many biological and chemical phenomena. To be more precise, our numerical scheme based on discontinuous Galerkin finite element concept with Legendre basis functions yields to a set of nonlinear equations to be solved in each subinterval. The numerical stability in the linear case is proved and the order of convergence is also discussed. Beside the direct computation of the nonlinear term, the technique of product approximation is also utilized and then their performance are compared for various J, r and ν. We have tested the performance of the LDG scheme on the linear as well as nonlinear growth and logistic differential equations of fractional order. Comparing our numerical results with the PECE indicates that the present approaches produce an accurate approximation for the underlying model problems.