A nonstandard compact finite difference method for a truncated Bratu–Picard model

1.   Introduction

This paper considers the truncated Bratu$-$Picard (tBP) model in a one-dimensional case as follows:

In Model (1.1), if we consider $M = \infty$ and $\lambda > 0$, then, the classical Bratu model will be obtained where its exact solution is known. This case has many applications in science and engineering. The Bratu type equation is also used in a large variety of applied fields, such as modeling thermal reaction processes in combustible non-deformable materials, including the solid fuel ignition model, the electrospinning process for the production of ultra-fine polymer fibers, modeling some chemical reaction-diffusion, questions in geometry and relativity about the Chandrasekhar model, radiative heat transfer, and nanotechnology ^{[1,2,3,4,5,6,7,8,9]}. Several numerical methods have been developed to approximate the solution of the Bratu equation ^{[10,11,12,13,14]}. Most existing methods yield one of the two solutions to the equation (the lower solution), for example, a Laplace transform decomposition algorithm ^[15], the direct shooting and Lie-group shooting methods ^[16], the perturbation iterations, parameter perturbations, splines methods ^[17], finite difference methods, and multigrid methods ^[18]. The lower and upper solutions in the case of multiple solutions are obtained using Boyd's approach ^[19]. In ^[1], Mickens' nonstandard finite difference method (NSFD) has been used to solve the Bratu$-$Gelfand equation, and a comparison with the standard finite difference method has shown that the results of the NSFD method are more accurate. Also, the NSFD scheme converges to both lower and upper solutions. Buckmire ^[1] applied Mickens' nonstandard finite difference method (NSFD) and compared the performances of the Adomian decomposition method, Boyd's pseudospectral method, the nonlinear shooting method, standard finite difference (SFD), and NSFD methods. Buckmire reported that the NSFD method may converge to both solutions (the lower and the upper one) and is more accurate than SFD. A smart NSFD scheme for the second-order nonlinear boundary value problem has been discussed in Erdogan ^[20]. A more general compact exponentially fitted method is used in ^[21] and SFD and NSFD approaches are considered as special cases. Recently, a numerical method has been presented to solve the Bratu-type equation based on the compact finite difference method (CFD) ^[22]; this method converges to lower and upper solutions and is more accurate than the finite difference approach.

In this paper, we intend to present the nonstandard compact finite difference method (NSCFD) to study the tBP model. Most previous articles considered only positive solutions but we obtained all smooth solutions using our proposed method, where some of them are periodic and others are semi-periodic. We also show theoretically and numerically that there exists a unique solution for $\lambda \leq 0$. We observe that NSFD has a similar simplicity as an SFD approximation but it is slightly more accurate, in most cases. In addition, the NSFD method preserves some qualitative features of the continuous-time model such as boundedness and positivity. The most important weakness of the used method is that there is no specific method to find the best denominator function in the nonstandard finite difference method.

We organize the rest of the article as follows: In Section 2, we describe the solutions in different cases of the tBP model. In Section 3, the NSCFD method is presented for solving the tBP model. In Section 4, the convergence analysis of the NSCFD method is investigated. In Section 5, the numerical results obtained by the methods of this study are presented. Finally, the conclusion is drawn in Section 6.

2.   Description of the solutions of the tBP model for different values of $M$

As we know, in Eq (1.1), $\quad M \in \mathbb{N}\cup \{0\}$. We treat the cases of $M = 0$, $M = 1$, and $M = \infty$ separately.

For the other values of $M$, we consider the following two subsets of $\mathbb{N}$.

In the following section, we describe the behaviors of five different states separately.

2.1.   Case of $M=0$

In this case, we have the following model ^[23]:

The exact solution of this case can be obtained directly,

Also, the maximum value of the solution is:

2.2.   Case of $M=1$

For this case, Eq (1.1) is converted to ^[24]:

Based on the value of $\lambda$, the solution of Model (1.1) is as follows:

(1) If $\lambda = 0$, the equation has only trivial solution $u_{\lambda}^{x} = 0$.

(2) If $\lambda < 0$, we have the following solution:

As $\lambda\rightarrow -\infty$, the solution has a horizontal asymptote, $\|u_{\lambda}^{1}\|_{\infty}\uparrow 1$.\\ (3) For $\lambda > 0$, $\lambda\neq (m\pi)^{2}$, the following solution is obtained:

But for $\lambda = (m\pi)^{2}$, $(m\in \{2, 4, 6, ...\})$, two cases are distinguished. For $\lambda = (m\pi)^{2}$, $(m\in \{2, 4, 6, ...\})$, Eq (1.1) has many solutions in the following form:

and for $\lambda = (m\pi)^{2}$, $(m\in \{3, 5, 7, ...\})$, there is no continuous solution for Eq (1.1).

2.3.   Case of $M=\infty$

In this case, Eq (1.1) is converted to the Gelfand$-$Bratu equation in standard form ^[23]:

We have the following two cases based on the value of $\lambda$:

(1) The exact solution for $\lambda > 0$ has the following form:

Equation (2.7) has two solutions for $0 < \lambda < \lambda_{c}$, one solution for $\lambda = \lambda_{c}$, and no solutions for $\lambda > \lambda_{c}$. From the relation $1 = \frac{1}{4}\sqrt{2\lambda_{c}}\sinh\left(\frac{\theta}{4}\right)$, the critical value of $\lambda_{c}$ is obtained as $\lambda_{c}\approx3.513830719$.

Also, the exact solution of (2.8) is symmetric around $x = \frac{1}{2}$.

(2) The exact solution for $\lambda < 0$ is unique and has the following form:

2.4.   Case of $M\in \mathcal{N}_{2}$

For example, we consider the case of $M = 2$ and we have the following equation ^[23]:

Model (2.10) has two solutions for $0 < \lambda < \lambda_{c}$, one solution for $\lambda = \lambda_{c}$, and no solution for $\lambda > \lambda_{c}$, which $\lambda_{c}\approx 3.96$. For $\lambda < 0$, there is a unique solution. Also, a similar result can be obtained for the other values in $\mathcal{N}_{2}$.

2.5.   Case of $M\in \mathcal{N}_{3}$

For example, we consider the case of $M = 3$ and $M = 5$, and we obtain the following equations ^[23]:

Equations (2.11) and (2.12) have infinitely many solutions for $\lambda > 0$ and a unique solution for $\lambda < 0$.

3.   The solution method

In this section, we use the CFD and NSCFD methods to approximate the solutions of the truncated Bratu$-$Picard Problem (1.1). To compute the numerical solution, we first subdivide the range of $[0, 1]$ into $N$ subintervals of width $h = \frac{1}{N}$, thus node points $x_{j} = jh$, $j = 0, 1, ..., N$, are obtained. Consider the following notations:

For the second derivative, we have the following compact finite difference scheme ^[24]:

The truncation error for System (3.1) is $O(h^{4})$ and the matrix form of (3.1) is as follows:

where

3.1.   Implement of the CFD scheme for the tBP model

From Eq (1.1), we have $u''_{i} = -\lambda \sum_{n = 0}^{M}\frac{(u_{i})^{n}}{n!}$, $i = 1, ..., N-1$. By inserting this relation into System (3.1), the following nonlinear system is obtained:

3.2.   Implement of the NSCFD scheme for the tBP model

The nonstandard finite-difference scheme (NSFD) is well-developed by Mickens ^[25,26]. It has many advantages that have been shown by many researchers ^[20,21,27]. One of the critical points in NSFD schemes is that the second derivative can be approximated by the following general form:

instead of the following standard form approximation:

where the function $\phi(h)$ satisfies:

In fact, as $h\rightarrow 0$, the standard finite difference and nonstandard finite difference schemes are identical.

For obtaining the nonstandard compact finite difference method, we replace the function $h^{2}$ in System (3.1) with the nonlinear function $\phi(h) = 2\ln \left(\cosh(h)\right)$ such that it satisfies Property (3.5).

Therefore, we have the following relations:

The matrix form of (3.6) is as follows:

where matrices $A$ and $B$ have been defined in Relation (3.2).

By inserting the Relation $u''_{i} = -\lambda \sum_{n = 0}^{M}\frac{(u_{i})^{n}}{n!}$, $i = 1, ..., N-1$, in System (3.6), the following nonlinear system is obtained:

By solving the nonlinear System (3.7), the numerical solution of Model (1.1) is obtained.

Algorithm 1 shows the steps of solving Eq (1.1) using the NSCFD method.

4.   Convergence analysis of the nonstandard compact finite difference method

In this section, we discuss the issue of convergence. For this purpose, let $\bar{U} = \left[u(x_{1}), ..., u(x_{N-1})\right]^{T}$ be the vector of the exact solution. Also, consider, $U = \left[u_{1}, ..., u_{N-1}\right]^{T}$ as the vector of the numerical solution. Moreover, consider $\|.\|$ as $\|.\|_{\infty}$ and $E = \bar{U}-U$.

Lemma 4.1. Let $T = \left[t_{1}, ..., t_{N-1}\right]^{T}$ be the vector of the local truncation error to (3.6). Then, we have

Proof. For $i = 2, ..., N-2$, the local truncation error is obtained as:

By using the Taylor expansion, it can be written as

By replacing Relation (4.3) in Eq (4.2), we have

Similarly, for $i = 1$ and $i = N-1$, we have

Therefore, we have

□

Theorem 4.2. Let $\bar{U} = \left[u(x_{1}), ..., u(x_{N-1})\right]^{T}$ be the vectors to the exact solution of the boundary-value Problem (1.1), $U = \left[u_{1}, ..., u_{N-1}\right]^{T}$ be the obtained numerical solution by solving the nonlinear System (3.7), and $E = \bar{U}-U$. Then, provided $M = \infty$, and\\

we have

Proof. By replacing Taylor expansion $2\ln\left(\cosh(h)\right) = h^{2}-\frac{1}{6}h^{4}+\frac{2}{45}h^{6}+O(h^{8})$ with the matrix form of Eq (3.6), we can write

By replacing Relation $U'' = -\lambda e^{U}$ in Eq (4.7), we have

For the exact solution, we have

where $T = O(h^4)$ is the local truncation error of (3.6). By using Relations (4.8) and (4.9), we have

where $E = \bar{U}-U$, $e^{\bar{U}} -e^{U} = JE+O(h^{2})$, $J = diag\{\frac{\partial e^{u(x)}}{\partial x}:x = x_{i}, i = 1, ..., N-1\}$ is a diagonal matrix of order $N-1$, and $T' = O(h^{4})$. Because the matrix $B_{2}$ is invertible, Relation (4.10) can be written as follows:

Now if $\lambda h^{2}\left(1-\frac{1}{6}h^{2}+\frac{2}{45}h^{4}\right)\left\|B_{2}^{-1}\right\|\|A_{2}\|\|J\|\leq 1$, then

is invertible and we have

Thus,

By using the geometric series theorem, it follows that

Therefore, we have

□

5.   Numerical solution

In this section, we apply the presented numerical method for the numerical solution of Model (1.1) and obtain the solutions for different values of $\lambda$ and $M$. To solve the nonlinear System (3.7), we use a simple approach similar to that used by Boyd ^[19]. We consider $u_{0} (x) = A\sin(k \pi x)$ as an initial guess because it satisfies the boundary conditions.

The values of parameters $A$ and $k$ depend on $M$ and $\lambda$. Specially, in the case of $\lambda > 0$, we need the condition $A < u_{max}$ to obtain the lower solutions, and also, we need $A > u_{max}$ to obtain the upper solutions, where $u_{max}$ is an approximation for the maximum value of the solution.

For $M = 3$, we consider $A = 9$ and $k = 2$ for the first periodic solution. Also, for the case of $M = 5$, we consider $A = 13$ and $k = 3$ for the first semi-periodic solution.

We use Maple 17 for obtaining numerical results and the fsolve command to solve the nonlinear system of equations. We test three cases, obtain the numerical results in each case, and show the problem behavior in all cases with a diagram. Also, the convergence order of the proposed method is calculated using the following formula:

where $E_{new}$ and $E_{old}$ are the maximum absolute errors corresponding to the new mesh size ($h_{new}$) and old mesh size ($h_{old}$), respectively. Also, we report the central processing unit (CPU) time for our numerical results.

5.1.   Case of $M=\infty$

In this section, the upper and lower numerical solutions of the nonlinear System (3.7) for $\lambda = -1, 0.0001,$ $0.001, 0.01, 0.1, 1, 2, 3, 3.51$ and $N = 5, 11, 21, 41, 81, 161, 321$ are compared with the exact solution of the problem and the maximum error is obtained. The results are shown in Tables 1–9. The bifurcated nature of the computed solution for different values of $\lambda$ has been plotted in Figure 1. In Table 10, we compare the upper solutions of the CFD method with the NSCFD method for $\lambda = 1$ and $N = 11, 21, 41, 81, 161$. Also, in Table 11, we compare the lower solutions of the CFD method with the NSCFD method for $\lambda = 2$ and $N = 11, 21, 41, 81, 161$.

5.2.   Case of $M=1$

In this part, we consider Example (2.3) in Section 2 and obtain numerical results using System (3.7). Figure 2 shows the numerical solutions for $\lambda = -1$, $\lambda = 1$, and $\lambda = 4{\pi}^{2}$ in this case.

5.3.   Case of $M=2$

In this part, we consider Example (2.10) in Section 2 and obtain numerical results using System (3.7). Figure 3 shows the upper and lower solutions for $\lambda = 1$ in this case. The bifurcation behavior of the upper and lower solutions for the positive $\lambda$ is shown in Figure 4.

5.4.   Cases of $M=3$ and $M=5$

In this section, we consider Examples (2.11) and (2.12) in Section 2 and obtain numerical results using System (3.7). In this case, there are an infinite number of solutions. Figure 5 shows five semi-periodic solutions for $M = 3$ and $\lambda = 1$. Figure 6 shows numerical solutions for $M = 3$ and $\lambda = -1$. Figure 7 shows the convergence behavior of the solutions using the NSCFD method for $M = 3$ and $M = 5$. Figure 8 shows seven periodic solutions for $M = 5$ and $\lambda = 1$. Finally, Figure 9 shows numerical solutions for $M = 5$ and $\lambda = -1$.

6.   Conclusions

In this paper, we obtained numerical solutions for the truncated Bratu$-$Picard model using the nonstandard compact finite difference method. The solutions are presented for different values of $\lambda$ and $M$, and the graph of each case is plotted. Numerical results showed the existence of two, one, and zero solutions for $\lambda > 0$ and $M \in \mathcal{N}_{2}$, which is similar to the $M = \infty$ case. In Figures 4 and 6, we showed that there are infinite numbers of solutions for the $M = 3$ and $M = 5$ cases. These solutions are either periodic or semi-periodic. Finally, we presented the bifurcating nature of Model (1.1) for each case. Previous articles considered only positive solutions but we obtained all smooth solutions using our proposed method, where some of them are periodic and others are semi-periodic. We also show theoretically and numerically that there exists a unique solution for $\lambda \leq 0$. We observe that NSFD has a similar simplicity as an SFD approximation but it is slightly more accurate, in most cases. In addition, the NSFD method preserves some qualitative features of the continuous-time model such as boundedness and positivity.

The most important weakness of the used method is that there is no specific method to find the best denominator function in the nonstandard finite difference method. Also, the most important disadvantage of the compact finite difference formulation is that compact schemes are implicit and require solving a matrix system for the evaluation of solutions or derivatives at the grid points. But, we accept this limitation of the compact finite difference approach due to its excellent stability properties. We intend to use the proposed method for solving the fractional order version of Model (1.1) in the future.

Author contributions

All authors of this article have been contributed equally. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

The authors extend their appreciation to Taif University, Saudi Arabia, for supporting this work through project number (TU-DSPP-2024-46).

Conflict of interest

The authors declare no conflicts of interest.