Numerical experiments with the Bratu equation in one, two and three dimensions

The non-linear Bratu equation is solved numerically in one, two and three dimensions. The numerical results are obtained with help of three different numerical methods namely pseudospectral method, finite difference method and radial basis functions method. All methods give approximately the same results. One-dimensional case is used for testing numerical methods as the analytical solution is known in this case. In two dimensions, the existing results are generalized to larger area of eigenvalues. Three-dimensional case is solved for the first time and this is the most important result of the paper.


Introduction
The aim of the present paper is to solve numerically the Bratu equation (see below) in three space dimensions. This equation coupled with appropriate boundary conditions is used to model the temperature distribution in combustion models (Bebernes and Eberly 1989). The three-dimensional model was considered in investigations on the suns core temperatures, see Chandrasekhar (1967); Buckmire (2003). The Bratu model stimulates also a thermal reaction process in a rigid material where the process depends on the balance between chemically generated heat and heat transfer by conduction (Wazwaz 2005;Noor and Mohyud-Din 2007a,b). We have chosen three different methods to achieve our goal: pseudospectral method, finite difference method and radial basis functions (RBF) method. Let us note that the RBF method we have used is original as it generalizes the classical approach on the boundary (see Sect. 3.4 for details). These methods are first intensively tested in the one-dimensional case in which the analytical solution is known. This case was also investigated numerically with help of many different methods: spline-based treatment (Abukhaled et al. 2012), (modified) variational iteration method (Jin 2010;Mohyud-Din 2008), Adomian's decomposition method (Wazwaz 2005;Syam and Hamdan 2006), the Lie-group shooting method (Abbasbandy et al. 2011), one-point pseudospectral collocation (Boyd 2011) and others. Our methods work fine and are in good agreement with analytical formulae and numerical results. The two-dimensional case is more complicated and analytical solutions are not known. It was studied numerically (Boyd 1986;Misirli and Gurefe 2011;Kapania 1990) and we verify and generalize the existing results. The three-dimensional case is examined, to our knowledge, for the first time and our results are novel. It is much more difficult than one-and twodimensional cases. It requires parallel computations to save time, and the starting point of the linearization method must be very carefully chosen to achieve convergence of the whole procedure-this will be discussed later. Let us also stress that the qualitative behaviour of the three-dimensional bifurcation curve is very different from that in one and two dimensions.
The paper is organized as follows. In Sect. 2 we describe the Bratu problem in detail. Section 3 is devoted to numerical methods. We review the Newton's method used to linearize our problem and the details of the numerical procedures mentioned above. Numerical results concerning all three examined cases are presented in Sect. 4. Section 5 involves short conclusions.
In the present paper, we solve numerically the above equation for n = 1, 2, 3. That is we examine how many solutions exist for a given eigenvalue λ > 0 and compute these solutions with the help of different numerical techniques. We have used the simple finite difference method, the pseudospectral method and the RBF method.
In the simplest, one-dimensional case, the Bratu equation reduces to in the interval |x| ≤ 1, with boundary conditions This case lets us test our computational methods, as the analytical solution is known up to some numerical constant z(λ) which can be obtained from the transcendental equation (Boyd 2011) This equation has two real branches meeting at the limit point λ limit = 0.878457679781290 3015, z limit = 1.8101705806989772753, above which there are no solutions. The exact solution of Eq.
(1) reads Thus ϕ(x; λ) has also two branches meeting at the bifurcation point. The analytical solution of the two-dimensional Bratu-equation is not known and the numerical study of it (Boyd 1986) shows also the existence of two branches which were (partly) computed only approximately. The present paper is, to our knowledge, the first (numerical) study of the three-dimensional case.

Newton's linearization method
Let us first discuss a more general problem of the form where L is some linear operator (Laplace operator in our case) and F a non-linear function (−λ · exp for the Bratu equation). Assume that we know ϕ n : the n-th approximation to the solution of our problem. We seek the next (n + 1)-th approximation in the form where δϕ n is a small correction. Linearizing our problem we get the following equation for the δϕ n Because we finally obtain This is a linear problem with the linear operator L − F (ϕ n ) and the source term on the right hand side of the above formula. The procedure requires the starting pair (ϕ 0, ϕ 1 )-ϕ 0 should be carefully chosen and satisfy the desired boundary conditions, ϕ 1 in turn can be calculated from Eq. (2) for n = 0. In the case of the Bratu problem, the Eq. (3) takes the form It must be solved on the area [−1, 1] dim with zero boundary conditions and the whole procedure is stopped when the norm |δϕ n | is less than some small positive constant eps.

Finite difference method
In this method each interval [−1, 1] is divided into n subintervals of equal length h = 2/n. The grid in each dimension consists of n + 1 equidistant points x 0 , x 1 , . . . , x n and the second derivative of the function u is approximated by where i, j, k are the node numbers in x, y, z direction, respectively. Inserting this formula into Eq. (4) we obtain the symmetric, sparse system of linear equations for δϕ n (i, j, k). There are several excellent numerical packages solving such systems-in our practice we have used IBM Watson Sparse Matrix Package (WSMP) https://researcher.ibm.com/researcher/ view_project.php?id=1426 and the Parallel Direct Sparse Solver (PARDISO) included in INTEL MKL library http://software.intel.com/en-us/articles/intel-mkl/. Let us stress that these packages work fine in multiprocessor (multicore) environment.

Pseudospectral method
In this method, we choose the Chebyshev points of the second kind to form the grid in the interval [−1, 1]. These points are given by The second derivative in the given node is approximated by the second derivative of the Lagrange polynomial interpolant in that point. The appropriate formula is linear in the node values of the function u(i) and is usually written with help of the differentiation matrix (Berrut et al. 2005;Berrut and Trefethen 2004) where D (2) is the second differentiation matrix with elements and D (1) is the first differentiation matrix The so-called barycentric weights w j are given by Applying the pseudospectral method to the Eq. (4), we obtain a sparse, unsymmetric system of linear equations, which can also be solved by WSMP or PARDISO package as these packages include unsymmetric solvers as well.

RBF method
A RBF is a real-valued function whose value depends only on the distance from some point c, called a centre, so that In this method, we choose the centres c i , i = 1, . . . , N lying in the computation area and the shape of the function f . The RBF method becomes very popular nowadays and is extensively studied both from theoretical and numerical point of view [see Boyd (2010a,b);   Table 1. [see for instance Larsson and Fornberg (2003)] Let us note that the RBFs shown in the Table 1 can depend on the shape parameter ε and we denote r = ||x − c||. The simplest RBF interpolant is a linear combination of RBFs centred at the scattered points c i where the second sum is over all centres lying on the boundary. The reason is that we want our Laplace equations to be valid also at the boundary centres where the (Dirichlet) boundary conditions must be fulfilled as well. Thus, we need two RBFs for each boundary centre. This Let us now apply the RBF interpolant for δϕ n to the Bratu problem. We must solve Eq. (4) subject to zero boundary conditions. We require this equation to be strictly valid at all centres and boundary conditions must be fulfilled at boundary centres. We have As a result we must solve the dense linear system for N + N B variables. The previously mentioned numerical packages WSMP and PARDISO are not suitable for this case as they have been derived for sparse linear systems and we have used Parallel Linear Algebra Software for Multicore Archtectures (PLASMA) package http://icl.cs.utk.edu/plasma/ which turned out to be excellent for solving dense linear systems on multicore processors.

The one-dimensional case
As was mentioned earlier the analytical solution is known in this case, so we compare our numerical results with this solution. Let us first show two real branches of the solutions meeting at the bifurcation point λ bif = 0.8784576797812903015.
They are presented in Fig. 1. It is interesting to show the errors of the numerical solutions depending on the eigenvalue λ. The error is defined as and we present the logarithm of the errors for our three numerical methods in Fig. 2. The number of grid points is 2,000 for the finite difference method, 25 for the pseudospectral method and 200 for the RBF method. It is evident that the best results we have obtained with the pseudospectral method (with the exception of the area of small λ for the upper (right) branch).
To be complete we present also the dependence of the numerical error on the number on grid points. Figure 3 presents this dependence for the finite difference method.

The two-dimensional case
This case was examined long ago by Boyd (1986). We present in Fig. 4 the Boyd's bifurcation curve obtained both by numerical and approximated analytical methods and compare it with our purely numerical bifurcation curve obtained by three different numerical methods described above. Let us stress that our curve obeys larger area although the upper branch slightly differs for different numerical methods especially for small values of the eigenvalue λ (Fig. 4).
We have obtained our numerical results for the following numbers of grid points: 25 × 25 for pseudospectral method, 100 × 100 for finite difference method and 1,850 for rbf method. An example of the two-dimensional numerical ϕ-surface (both for lower (left) and upper (right) branch) is presented in Fig. 5 4.3 The three-dimensional case The most important aim of the numerical study of this case was to get the bifurcation curve of the solutions in the largest possible area. As in the previous cases we applied three numerical 3.4 for details). All three methods worked fine for the lower branch but for the upper branch they failed for λ approximately equal to 1.3. Therefore, the neighbourhood of this point required very careful investigations. After doing them it turned out that the qualitative behaviour of the bifurcation curve in three dimensions is quite different than in one and two dimensions. We have found two new bifurcation points: λ 2 ≈ 1.321 (local minimum) and λ 3 ≈ 1.556 (local maximum). The first bifurcation point λ 1 ≈ 2.475 is the global maximum of our bifurcation curve which is presented in Fig. 6. Thus, starting from λ ≈ 0.0 we have two solutions for a given λ (up to λ = λ 2 ). For λ = λ 2 and λ = λ 3 , three solutions are present and between these values four solutions exist. Above λ 3 , we have two solutions up to λ = λ 1 and for λ greater than this value there are no solutions at all. Let us note that obtaining upper branches of the three-dimensional bifurcation curve was a difficult numerical task as this process was very sensitive to the choice of the starting point of the Newton's procedure (see Sect. 3.1). We present also two examples of the solutions for λ = 2.25, one from the lower branch and the second from the first upper branch. Figure 7 shows their sections by the hyperplane z = 0.

Conclusions
In the present paper, we have presented numerical solutions of the Bratu equation in one, two and three dimensions. Our results in the one dimension show that the numerical procedures we have applied work fine and our numerics agrees with analytical formulae with high precision. Therefore, it was tempting to apply them in higher dimensions. In the two-dimensional case, the numerical results we have obtained confirm qualitatively the shape of the bifurcation curve presented by Boyd, but our bifurcation curve covers a wider range (it is valid also for small values of the eigenvalue λ in the upper branch).
The three-dimensional case is numerically much more complicated. The lower branch of the bifurcation curve was rather easy to obtain numerically, but for the upper branch the numerical procedures failed for λ ≈ 1.3. The reason of such behaviour became clear for us after very careful investigation of the neighbourhood of this point. It turned out that two new bifurcation points exist (the local minimum and the local maximum) and the qualitative behaviour of the three-dimensional bifurcation curve is quite different than in one and two dimensions (see Fig. 6). Let us add that finding the additional branches of the bifurcation curve was a very hard task. First, we had to get one point of a given branch. Then, we obtained other points by changing infinitesimally the parameter λ → λ + δλ and taking the solution for λ as the starting point of the Newton's procedure (see Sect. 3.1) that yielded the solution for λ + δλ. We repeated this procedure for all values of λ on the given branch. It is worth mentioning that the pseudospectral method seems to be the quickest in three dimensions. Our results are summarized in Fig. 8, where the all three bifurcation curves are presented.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.