An Efﬁcient Scheme for Time-Dependent Emden-Fowler Type Equations Based on Two-Dimensional Bernstein Polynomials

: In this study, we introduce an efﬁcient computational method to obtain an approximate solution of the time-dependent Emden-Fowler type equations. The method is based on the 2D-Bernstein polynomials (2D-BPs) and their operational matrices. In the cases of time-dependent Lane–Emden type problems and wave-type equations which are the special cases of the problem, the method converts the problem to a linear system of algebraic equations. If the problem has a nonlinear part, the ﬁnal system is nonlinear. We analyzed the error and give a theorem for the convergence. To estimate the error for the numerical solutions and then obtain more accurate approximate solutions, we give the residual correction procedure for the method. To show the effectiveness of the method, we apply the method to some test examples. The method gives more accurate results whenever increasing n , m for linear problems. For the nonlinear problems, the method also works well. For linear and nonlinear cases, the residual correction procedure estimates the error and yields the corrected approximations that give good approximation results. We compare the results with the results of the methods, the homotopy analysis method, homotopy perturbation method, Adomian decomposition method, and variational iteration method, on the nodes. Numerical results reveal that the method using 2D-BPs is more effective and simple for obtaining approximate solutions of the time-dependent Emden-Fowler type equations and the method presents a good accuracy.


Introduction
The heat equation which expresses the diffusion of heat can be given by the following form [1][2][3]: (x, t) + a f (x, t) g(y(x, t)) + h(x, t) = ∂y ∂t (x, t), 0 < x ≤ L, 0 < t < T, subject to y(0, t) = α, y x (0, t) = 0, where y(x, t) represents the temperature and t the time, r > 0 and a is an integer, α is a constant and the term f (x, t)g(y) + h(x, t) is the nonlinear heat source, where f , g and h are some smooth functions. For the steady case with h(x, t) = 0 and r = 2, Equation (1) is the Emden-Fowler equation [4][5][6], that is, In addition, if f (x) = 1, then Equation (3) becomes the Lane-Emden equation of mathematical physics [7,8].
In this work, we study the numerical solution, which is a simply applied and effective method, of the time-dependent Emden-Fowler type equations. First, the heat-type Equation (1) is considered. To solve Equation (1) by the proposed numerical technique, an operational matrix based on two-dimensional Bernstein polynomials is constructed. Second, with the same notations, we obtain the approximate solutions of the following wave-type equation: subject to y(0, t) = α, y x (0, t) = 0.
Since the operational matrices method can achieve the singularity behavior at the point x = 0, it has been applied to some singular problems. Yousefi and Behroozifar [24] applied the Bernstein operational matrix to solve the Emden-Fowler equation for the steady-state case. The same problem was considered by Gupta and Sharma [25] and solved by the Taylor series method. Another work to solve the steady-state problem was given by Wazwaz et al. [26]. They applied VIM to the problem. When only the Lane-Emden equation is considered, there are many numerical methods that depend on ADM [27,28], VIM [29][30][31][32], HPM [8,33], and operational matrices method [34] or Bernstein collocation method [35] were used to get analytic or numerical solutions.
The methods based on operational matrix of differentiation were given in various forms in the literature. Some of them were constituted by Chebyshev polynomials [36], Legendre polynomials [37], and Bernstein polynomials [38][39][40][41][42][43][44]. One of the common properties of these polynomials is the basis of polynomial spaces. On the other hand, because of Bernstein polynomials are dense in L 2 (Ω) [45] and thus yields good approximation results, they have been used often recently to solve the problems in two-dimensions. The Bernstein operational matrices method is also known as the Bernstein matrix method [44] and Bernstein series solution [35].
In this work, we shall develop a method based Bernstein operational matrices to solve the nonlinear problem 1. The method is applied easily and presents a good accuracy. We first give Bernstein polynomials in 1D and 2D in Section 2. The matrix representations of the approximate functions and their operational matrices for both dimensions are given next. Section 3 presents the method of solution for heat-type and wave-type equations. In Section 4, we constitute the residual correction procedure both to estimate the absolute error and to get the corrected approximate solutions. The most important property of the procedure can be applied even if the exact solution is unknown. Section 5 contains some examples including linear and nonlinear models to demonstrate the applicability of the method. First, we apply the method to time-dependent Lane-Emden type problems and singular wave-type equations. We perform the method for different number of nodes. The numerical results show that increasing number of nodes gives a sequence of approximations, which converges to the exact solution, for linear problems. The method also gives good results for nonlinear problems. Residual correction procedure estimates the error in general. On the other hand, the corrected approximate solutions can be obtained and its error is less than the approximate solution. We compare the approximate solutions obtained by the method with the approximate solutions obtained by ADM, HPM, HAM, and VIM by calculating the error on some nodes in [0, L] × [0, T]. In the last section, we summarize the results.

2D Bernstein Polynomials (2D-BPs) and Their Operational Matrices
In this section, because of defining two-dimensional Bernstein polynomials easily, we will give one-dimensional, Bernstein polynomials first. The Bernstein polynomials of degree m (1D-BPs) on [0, 1] are defined as [46] where (6) can be rewritten as Let y : I → R be any real-valued continuous function. The truncated Bernstein polynomials approximation of y is where To ease computational complexity, the matrix representations of are to be determined. Let (A i+1 ) 1,j and T m (x) be 1 × (m + 1) and (m + 1) × 1 matrices respectively for j = 0, 1, . . . , m as follows: Then, the polynomial B i,m (x) is given by which gives the identity: where Note that and where By using (10), we write y m (x) as where D = ALA −1 . Similarly, for higher-order derivatives of y m (x), one has The Bernstein polynomials of degree mn (2D-BPs) on [0, 1] × [0, 1] are defined as [47,48] which can be rewritten as [45]  Now, assume that y is a real-valued function defined on Ω ⊂ R 2 . We want to approximate y(x, t) by the truncated 2D-BPs series, i.e., If the identity in [47] is applied to (17), we get Then, y m,n in (17) can be found using where The partial derivatives ∂y/∂x and ∂y/∂t can be approximated as Upon using (9), (13), and (15), one gets the identities where Hence, we get For the higher-order derivatives, we have

Solving Heat-Type and Wave-Type Equations by 2D-BPs
To approximate the solution of (1) subject to the initial condition (2), one substitutes the identities in (25), (26), and (27) into (1). The residual Re(x, t) are given by Taking the nodes in (28), we obtain (m − 1)(n + 1) equations The conditions (2) then yield Thus, we have (m + 1)(n + 1) equations. By solving these equations by the command 'fsolve' which uses Newton's method for the unknown coefficients C i,j , we can find the matrix C so that y m,n (x, t) is obtained.

Error Analysis
In this section, we first give the Bernstein series that converges for the function in C 1 . The residual correction procedure is given to estimate the absolute error and to obtain corrected approximate solutions. To determine the convergence of the method, we define the constant M n+1 as follows.
where all partial derivatives of f of order n + 1 are bounded in magnitude M n+1 , which is Now, the main theorem of function approximation using 2D-BPs are stated as follows.
where f m,n is the best approximation out of Y m,n and supposing that then lim Proof. Letf Therefore, by considering Definition 1 and by taking η = max{m, n}, we get By considering Taylor's expansion formula in Definition 1, one has Hence, we get lim This concludes the proof. Now, we constitute residual correction procedure for the problem. Let y(x, t) be the exact and y m,n (x, t) the 2D-BPs approximate solutions of Equation (1), respectively. Writing e m,n (x, t) = y(x, t) − y m,n (x, t) as the error, then Substituting the approximate solution along with the constructed operational matrix in (42) yields the residual equation Re err (x, t) for the error as To solve Equation (43) numerically, we compute it at the following collocation points: Applying the proposed method to (43) with the given initial conditions gives us the coefficient matrix C err so that we get an approximate solutionẽ s,p for the error. Hence, y m,n +ẽ s,p is another approximate solution, which is called the corrected 2D-BPs approximate solution for the problem. In practice, selecting s ≥ m, p ≥ n yields better approximation results.
The solution y s,p m,n := y m,n +ẽ s,p improves the approximation y m,n provided that e m,n −ẽ s,p = y − (y m,n +ẽ s,p ) ≤ y s,p m,n .
Similar definitions and results can be given for the wave-type Equation (4).

Application to Several Test Problems
The applicability of the proposed 2D-BPS method is demonstrated via several test problems. To simulate the results, three types of equations are selected. Our results are compared with the methods of ADM, HPM, VIM, and HAM with h = −1. We use the Newton-Cotes points (32) and Maple 15 to simulate the numerical outlets with 50-digit precision.
The exact solution is [1][2][3][4]: To solve the above problem via the proposed method, we get from Equations (30) and (31) and Equations (46) and (47) Re Substituting (32) into Equations (49) and (50) for m = n = 3 yields a system of linear equations for the unknown coefficients C i,j , which once computed numerically gives We perform the method to the problem for various n = m, s = p = 10 and s = p = 20. The results are given in Tables 1 and 2. As seen from the tables, increasing n gives the decreasing error sequence while the computational time increases. On the other hand, residual correction procedure estimates the error well for n = p. The corrected solutions are more accurate than the solutions for n < p. We give the computational costs of the method in Table 3. Increasing n = m yields more computational costs which means that the method needs more computational times, numbers of multiplications, summations, assignments, and storages for big n, m. As a result, we can say that the numbers n, m are selected as not too big or not too small. It seems to be suitable for selecting n = m around 10 for Example 1.
The graphs of |e 12,12 |, |ẽ 15,15 | and |y − (y 12,12 +ẽ 15,15 )| are plotted in Figure 2 which also depicts the absolute errors of the 5th-order approximate solutions of ADM, HPM, HAM and VIM [1][2][3][4]. We also give a comparison for the method with HAM and VIM in Table 4 for different nodes number. It can be seen from Figure 2 and Table 4 that the corrected 2D-BPs approximate solution yields better approximate results than the results obtained by ADM, HPM, HAM, and VIM. Moreover, the residual correction procedure estimates the error more accurately.    Table 3. The computational costs of the obtaining Re(x i , t j ) of the method for 0 ≤ i, j ≤ n in the Maple code for Example 1. 5  7827  217  3901  217  10  40,955  727  5503  727  15  309,207  1537  184,281  1537  20 851,768 2647 549,717 2647

Singular Wave-Type Equations
Example 2. We consider now where the exact solution is y(x, t) = x 3 + e x 2 −t .
The solution in series form by ADM [1], HAM [2], and HPM [3] is We applied the method to the problem for various m = n. We can say from Table 5 that the norm of the absolute error decreases when m, n increases. The procedure estimates the error well and the corrected solutions are better than the solutions in the norm. The absolute error and corrected absolute error for s = p = 18 are shown in Figure 3. We also give the solutions of the problem by ADM, HPM, and HAM [1][2][3] to make a comparison in the same figure. We can conclude from this figure that the corrected absolute error is better than the absolute error, also better than the absolute error of ADM, HPM, and HAM.

Nonlinear Models
Example 3. Consider y xx + 5 x y x + (24t + 16t 2 x 2 )e y − 2x 2 e y/2 = y t , The exact solution of the problem [1][2][3][4] is y(x, t) = −2 ln(1 + tx 2 ). We apply the method to the linear problem for various m = n and apply the procedure for s = p = 9. The results are given in Tables 6 and 7. As seen from the tables, increasing m = n yields more accurate results, whereas the computational times increase. On the other hand, the proposed method yields a much better approximate solution than the solutions obtained by ADM, HAM, and HPM [1][2][3], of order 5. We plot the error, the estimation of the error, and the corrected error for m = n = 8 and ms = p = 9 in Figure 4.
We can say from the figures that the estimation by using the procedure fits the error well and the corrected error is smaller than the error.

Conclusions
In this study, we have proposed a numerical method to solve time-dependent Emden-Fowler type equations. In the numerical scheme, the 2D-BPs are utilized. The residual correction procedure is applied to estimate the error of the approximate solution. By using the suggested procedure, we obtain a new approximate solution which is more accurate than the 2D-BPs approximate solution. The error analysis as well as convergence of the proposed method have been investigated. We applied the method to several test problems including linear and nonlinear problems and also compared with some other methods to show the efficiency of the new method. As seen from the results of text examples, increasing the node numbers n, m yields a sequence which converges to the exact solution for all examples in the norm. The computational times and necessary operations in the code increase if n, m increases. Hence, the optimum n, m values are observed around 10. The procedure estimates the error with a good accuracy for n, m = s, p. In addition, by using the procedure, more accurate results are obtained. We can conclude that the obtained results are consistent with the results of ADM, HPM, HAM, and VIM. Especially for the nonlinear problem, the method gives better approximate solutions with respect to the semi-analytic methods.