Numerical solution of non-linear Bratu-type boundary value problems via quintic B-spline collocation method

1.   Introduction

In this work, we consider the one-dimensional Bratu-type BVPs:

where $\lambda$ is a positive physical parameter and $W\left( u \right)$ is non-linear function. The Bratu-type (BT) problem appears in a large variety of application regions, for example, fuel ignition model of thermal combustion, thermal reaction, chemical reactor theory and nanotechnology, radiative heat transfer, the Chandrasekhar model of the expansion of the universe ^[1,2,3,4]. These models play a significant role in physical phenomena and arises non-linear differential equations ^{[5,6,7,8,9,10]}. The non-linear differential equations are complicated to solve either analytically or numerical. Much attention has recently been devoted to search for better and more efficient solution methods for finding the approximate solution of non-linear models. Numerical methods play a significant role for finding the approximate solution of non-linear problems.

The BVPs defined by Eq (1) with boundary conditions (2) have substantial literature on numerical and analytical solutions. Wazwaz ^[11] has presented Adomian decomposition method to find the reliable solution of BT equation. Caglar et al. ^[12] have presented three-degree B-spline collocation method (BSCM) for finding the numerical treatment of BT problems. Roul and Thula ^[13] have discussed a fourth order BSCM for finding the numerical treatment of BT and Lane-Emden problems. Feng et al. ^[14] have discussed the Homotopy perturbation method for obtaining the solution of BT problem. Variational iteration method has been presented by Batiha ^[15] to calculate the approximate solution of BT equation. Aksoy and Pakdemirli ^[16] have developed the perturbation iteration theory for linear and non-linear second-order differential equation and implemented on BT equation to obtain the solution. Venkatesh et al. ^[17] have discussed the Haar wavelet approach for solving initial and boundary value problem of BT. Abukhaled et al. ^[18] have presented the cubic BSCM, an adaptive spline collocation, and optimal collocation techniques for calculating the numerical treatment of BT equations. Raja and Ahmad ^[19] have discussed the artificial neural networks to find the numerical solution of one-dimensional BT equation with the use of three type of transforms function including log-sigmoid, radialbasis, and tan-sigmoid transfer functions. Hohsen ^[20] has done the survey of different numerical techniques for one-dimensional and two-dimensional BT problems. Ghazanfari and Sepahvandzadeh ^[21] have discussed the solution of fractional BT equations by using the Adomain decomposition method based on Taylor series. Darwish and Kashkari ^[22] have developed the optimal Homotopy asymptotic method to obtain the approximate solution of second order BT initial value problems. Abd-Elhameed et al. ^[23] have proposed operational matrix method based on shifted Legendre polynomials to find the numerical solution of second order linear and non-linear boundary value problems. Babolian et al. ^[24] have investigated reproducing kernel method for solving the fractional order BT equations. Wazwaz ^[25] has used successive differentiation method for obtaining the solution of non-linear BT problems. Raja et al. ^[26] have introduced bio-inspired computing technique for solving BT equations. They took advantage of artificial neural networks' strengths to find the solution of one-dimensional BT problem. The validity of the model is dependent on the weights of the artificial neural networks being properly adjusted. Masood et al. ^[27] have applied mexican hat wavelet based neural network for the approximate solution of non-linear BT equation. Grover and Tomer ^[28] have presented the differential transform method to compute the numerical solution of the fractional BT equations and caputo derivative is used to represent the fractional derivative. Keshavarz et al. ^[29] have discussed an effective numerical method based on the Taylor wavelets for solving initial and boundary value problems of BT equations. Sakar et al. ^[30] have developed reproducing kernel approach based on Legendre polynomials to obtain numerical treatment of fractional BT BVPs. Tomar and Pandey ^[31] have introduced an efficient analytical iterative method to obtain the solution of initial and boundary value BT problems.

The B-spline function is widely used with combination of other effective numerical methods to find numerical solution of different types of problems which are arises in various physical phenomena such as nonlinear Schrödinger equation ^[32], complex modified Korteweg-de Vries (mkdv) equation ^[33], kdv-mkdv equation ^[34], modified equal width wave equation ^[35] and modified Kawahara equation ^[36], etc. Başhan ^[37,38] has presented the Crank–Nicolson scheme and differential quadrature method (DQM) based on quintic B-spline basis function to find the numerical solution of mkdv equation and coupled kdv equations. Mirzaee and Alipour ^{[39,40,41,42,43]} have presented the Bicubic, cubic, and quintic B-spline collocation method for finding the numerical treatment of two-dimensional weakly singular stochastic integral equation, dimensional stochastic itô-Volterra integral equation, fractional-order linear stochastic integro-differential equation, and multidimensional non-linear stochastic quadratic integral equation. Mirzaee and Alipour ^[44] have used fractional-order orthogonal Bernstein polynomial for solving non-linear fractional partial Volterra integro-differential equations.

The remaining part of the article is organized as follows. Derivation of the QBSCM is discussed in section 2. Error analysis of the QBSCM is described in section 3. Numerical results are exhibited in section 4 which shows the efficiency and accuracy of the QBSCM. Conclusions of the method are summarised in section 5.

2.   Derivation of quintic B-spline Collocation method

In this section, QBSCM is described to find the numerical treatment of Eqs (1) and (2). Let $P = \left\{ {0 = {t_0} < {t_1} < \ldots < {t_{N + 1}} < {t_N} = 1} \right\}$ be the partition of the interval $\left[ {0, \, \, \, 1} \right]$ and generate the uniform mesh with equally spaced mesh points ${t_i} = ih, \, \, i = 0, 1, 2, \ldots , N,$ where $h = {1 \mathord{\left/ {\vphantom {1 N}} \right. } N}.$ We define the quintic spline space ${S_{5, \, \, P}}$ as ${S_{5, \, \, P}} = \left\{ {q\left( t \right):\, \, q\left( t \right) \in {C^4}\left[ {0, \, 1} \right]\, \, } \right\},$ where $q\left( t \right)$ is a quintic polynomial on the partition $P.$ The standard B-spline is defined ^[45,46] as

where ${\delta ^{n + 1}}$ represents the ${\left( {n + 1} \right)^{th}}$ order central difference with unit step and the function $x_ + ^n$ defined as

Putting $n = 5, \, \, x = \frac{{\, t - {t_i}}}{h} - \frac{1}{2}$ in Eq (3), we obtain the following quintic B-spline (QBS) function ${Q_i}\left( t \right), \, \, \, i = - 2, - 1, , \ldots , \, N + 1, \, \, N + 2$ for the space ${S_{5, \, \, P}}$

Each basis function ${Q_i}\left( t \right)$ is fourth times continuous differentiable on $\left[ {{t_{i - 3}}, \, \, \, {t_{i + 3}}} \right].$ From Eq (5), we can see that $\sum\limits_{i = - 2}^{N + 2} {{Q_i}\left( t \right)} = 1$ on $\left[ {0, \, \, 1} \right]$ and each function ${Q_i}\left( t \right)$ are an element of ${S_{5, \, P}}.$ To enable the B-spline basis functions, we introduce six additional mesh points outside of the interval $\left[ {0, \, \, 1} \right].$ The six extra knots of the partition $P$ as ${t_{ - 3}} < \, \, {t_{ - 2}}, < \, {t_{ - 1}}$ and ${t_{N + 1}} < \, \, {t_{N + 2}}, < \, {t_{N + 3}}.$ Let $K = \left\{ {{Q_{ - 2}}, \, \, \, {Q_{ - 1}}, \, \, \, \ldots \, , \, \, \, {Q_{N + 2}}, \, \, \, {Q_{N + 3}}} \right\}$ is the set of QBS basis function and $\phi \left( P \right) = {\text{span}}\, K.$ Since the basis function ${Q_{ - 2}}, \, \, \, {Q_{ - 1}}, \, \ldots \, , \, \, {Q_{N + 2}}, \, \, \, {Q_{N + 3}}$ are linearly independent of the interval $\left[ {0, \, \, 1} \right].$ Therefore, $\phi \left( P \right)$ is a $\left( {N + 5} \right) - {\text{dimension}}$ QBS. It is to be seen that $\phi \left( P \right) = {S_{5, \, \, P}}$ ^[47]. Let $S\left( t \right)$ is the QBS interpolation function of $u\left( t \right)$ at the knots ${t_i}$. Then the QBS is defined as

where, $S\left( t \right)$ satisfies the following interpolation formula

We represent $S_i^{\left( m \right)} = {S^{\left( m \right)}}\left( {{t_i}} \right)$ for $i = 0, 1, \ldots , N$ and ${S^{\left( m \right)}}\left( t \right) = \sum\limits_{i = - 2}^{N + 2} {{\alpha _i}Q_i^{\left( m \right)}\left( t \right)}$ be the approximation to ${u^{\left( m \right)}}\left( t \right), \, \, \, \left( {m = 0, 1, 2, 3, 4} \right).$ The values of ${Q_i}\left( t \right), \, \, {Q'_i}\left( t \right), \, {Q_i}^{\prime \prime }\left( t \right)\,$ ${Q_i}^{\left( 3 \right)}\left( t \right)$ and ${Q_i}^{\left( 4 \right)}\left( t \right)$ at the mesh points are given in the Table 1.

From Eqs (1) and (6), we get

with the boundary conditions (2), we have

and

The above Eqs (9)–(11) are providing $\left( {N + 3} \right)$ equations with $\left( {N + 5} \right)$ unknowns. To find the unique solution of these systems of equations. We required two more equations. For finding the two more equations, differentiating Eq (1) with respect to t, we have

At the point $t = {t_0}\, {\text{and}}\, \, t = {t_N},$ we obtain

Substituting the quintic B-spline solution from Eq (6) in Eqs (13) and (14), we get

Thus, Eqs (9)–(11), (15) and (16) yield a system of $\left( {N + 5} \right)$ equations with $\left( {N + 5} \right)$ unknowns. We write these equations in the matrix form:

where, $A$ is the coefficient matrix and $X = {\left[ {{\alpha _{ - 2}}, {\alpha _{ - 1}}, {\alpha _0}, \ldots , {\alpha _N}, {\alpha _{N + 1}}, {\alpha _{N + 2}}} \right]^T}$ with the right-hand side $B = \lambda \, \, {\left[ {0, \frac{d}{{dt}}W\left( {{S_0}} \right), \, W\left( {{S_0}} \right), W\left( {{S_1}} \right), \ldots , W\left( {{S_N}} \right), \frac{d}{{dt}}W\left( {{S_N}} \right), 0} \right]^T}$ and $A$ is given below:

Solving the system defined in Eq (17) by using the Newton's method for the unknowns ${\alpha _{ - 2}}, {\alpha _{ - 1}}, {\alpha _0}, \ldots {\alpha _N}, {\alpha _{N + 1}}, {\alpha _{N + 2}}$ and substituting these values in Eq (6), we get the quintic B-spline solution of Eqs (1) and (2).

3.   Error analysis

Error bound is studied for QBS $S\left( t \right)$ and its derivatives up to order four at knots ${t_i}, \, \, \, \, i = 0, 1, \ldots , N.$ From Eq (6) and Table 1, we establish the following consistency relations ^[48]:

where $\gamma$ represents the discrete operator which is defined as follows:

We express the following two lemmas for deriving the error bound of quintic spline $S.$

Lemma 1. If $S\left( t \right)$ is the quintic spline (QS) interpolation function of $u\left( t \right)$ and $u\left( t \right) \in {C^8}\left[ {0, \, \, 1} \right].$ Then we have

Proof. Utilizing the interpolation condition (7) in the RHS Eq (18) and expanding $u$ by Taylor's series, we can acquire Eq (22). In the similar manner, we can prove the reimaging relations of Lemma 1.

Lemma 2. Let the coefficient matrix of the system $\gamma {y_i} = 0, \, \, \, i = 2, 3, \cdots , N - 2$ is denoted by $P$ and ${y_0} = {y_1} = {y_{n - 1}} = {y_n} = 0$ then prove that ${\left\| {{P^{ - 1}}} \right\|_\infty } \leqslant \frac{1}{{12}},$ where

Proof: From the given matrix $P$, we can see that, it is strictly diagonally dominant and non-singular matrix. Therefore, from the Lemma 4 in Lucas ^[49], we can obtain ${\left\| {{P^{ - 1}}} \right\|_\infty } \leqslant \frac{1}{{12}}.$

We estimate the error bound for the QS and its derivative in the following theorem.

Theorem 1. Let $S\left( t \right)$ is the QS interpolation of $u\left( t \right) \in {C^8}\left[ {0, \, \, 1} \right]$ and satisfies the interpolation conditions (7) and (8). Then we obtain the following relations for $i = 0, 1, 2, \ldots , N.$

Solution: We have

Expand the R.H.S. of Eq (30) by using the Taylor's series in terms of $y$, we obtain

We consider ${y_i} = u_i^{\left( 4 \right)} - \frac{1}{{12}}{h^2}u_i^{\left( 6 \right)} + \frac{1}{{240}}{h^4}u_i^{\left( 8 \right)}, \, \, i = 2, 3, \ldots , N - 2$ and substituting it in the Eq (30). We get

Subtracting Eq (32) from Eq (25), we have

Define ${d_i} \equiv S_i^{\left( 4 \right)} - u_i^{\left( 4 \right)} + \frac{1}{{12}}{h^2}u_i^{\left( 6 \right)} - \frac{1}{{240}}{h^4}u_i^{\left( 8 \right)}, \, \, \, \, \, i = 0, 1, \ldots , N.$ From Eqs (7) and (33), we get

The inverse of the coefficient matrix $P$ of the Eq (34) has a bounded norm and using Lemma 2, we obtain

Hence, the Eq (29) is proved. Next, we have to prove Eq (28). From the consistency relation of quintic spline, we get

Using Eqs (7) and (29) in Eq (36), we obtain

Expand the above equation by Taylor's series expansion, we get

Hence Eq (28) is proved for $i = 1, 2, \ldots , N - 2.$ We still need to prove Eq (28) for $i = 0, \, \, N - 1, \, \, N.$ The following quintic B-spline consistency relations are considered:

Using Eqs (7), (29) and (38) in Eqs (39) and (40). After substitution, expand the equations by Taylor's series and we get the Eq (28) for $i = 0, \, \, N - 1, \, \, N.$ Similarly, we can prove the Eqs (26) and (27).

Now, we drive the global error bounds in the following theorem.

Theorem 2. Consider $S\left( t \right)$ is the QS solution of $u\left( t \right)$ and $u\left( t \right) \in {C^\infty }\left[ {0, \, \, 1} \right]$ then we have the global error bounds ${\left\| {{u^{\left( m \right)}}\left( t \right) - {S^{\left( m \right)}}\left( t \right)} \right\|_\infty } = O\left( {{h^{6 - m}}} \right),$ for $m = 0, 1, 2, 3, 4.$

Proof. The QS interpolation error in defined by $e\left( t \right) = u\left( t \right) - S\left( t \right)$ and first we prove

Since ${S^{\left( 4 \right)}}\left( t \right)$ is piecewise continuous linear function in the interval $\left[ {0, \, \, 1} \right]$ with respect to the partition. For $i = 1, 2, 3, \ldots , N,$ we get

Let ${\hat u^{\left( 4 \right)}}\left( t \right)$ is piecewise continuous linear interpolating function to ${u^{\left( 4 \right)}}\left( t \right).$ For $i = 1, 2, 3, \ldots , N,$ we have

From Eqs (29), (42) and (43), we obtain

Therefore,

Using the piecewise linear interpolation theory from ^[50,51,52], we get

From Eqs (45) and (46), we have

Now, we have to prove

Let $S_i^{\left( 3 \right)}\left( t \right)$ represents the restriction of ${S^{\left( 3 \right)}}\left( t \right), \, \, t \in \left[ {{t_{i - 1}}, \, \, \, {t_i}} \right], \, \, i = 1, 2, 3, \ldots , N.$ By using the Newton-Leibniz formula, Eqs (28) and (48) for $\, t \in \left[ {{t_{i - 1}}, \, \, \, {t_i}} \right],$ we get

Hence ${\left\| {{u^{\left( 3 \right)}}\left( t \right) - {S^{\left( 3 \right)}}\left( t \right)} \right\|_\infty } = O\left( {{h^3}} \right)$ holds.

Similarly, we can prove, ${\left\| {{u^{\left( 2 \right)}}\left( t \right) - {S^{\left( 2 \right)}}\left( t \right)} \right\|_\infty } = O\left( {{h^4}} \right), \, \, \, \, \, \, \, \, \, {\left\| {{u^{\left( 1 \right)}}\left( t \right) - {S^{\left( 1 \right)}}\left( t \right)} \right\|_\infty } = O\left( {{h^5}} \right)$ and ${\left\| {u\left( t \right) - S\left( t \right)} \right\|_\infty } = O\left( {{h^6}} \right).$

4.   Numerical results

In this section, QBSCM is applied on two non-linear BT BVPs to exhibit the accuracy and efficiency of the method. The efficiency of the method is assessed by calculating the maximum absolute error (MAE) and order of convergence. The MAE is defined by

where $u\left( t \right)$ represents the exact solution and $S\left( t \right)$ denotes the approximate solution. The order of convergences is defined as

Example 1. Consider the following non-linear BT BVP ^[12,13]:

with the boundary conditions (BCs):

The exact solution of the Example 1 is given by

where $\theta$ is the solution of $\theta - \sqrt {2\lambda } \cosh \left( {\frac{\theta }{4}} \right) = 0.$

Example 2. Consider the following non-linear BT BVP:

with the BCs:

The exact solution of the Example 2 is $u\left( t \right) = - \ln \left[ {1 - \cos \left( {\left( {0.5 + t} \right)\pi } \right)} \right].$

The MAE and order of convergence of Example 1 are presented in Tables 2–4 and compared with the existing methods for $\lambda = 1, \, 2, \, 3.51$ and different values of $N.$ Tables demonstrate that the proposed method gives a good approximate solution. Point wise error and MAE for $\lambda = 2$ and $N = 20, \, \, 60\, \, {\text{and}}\, \, 90$ are displaying in Tables 5–7 respectively. Obtained results are compared with existing methods ^[12,13] at the same number of mesh points. Table 8 exhibits the MAE and order of convergences of Example 2 for different values of $N.$ Numerical results have been calculated through the MATLAB (R2021b) software and execution time (in seconds) is given in Tables 2–4 and 8. Figures 1 and 2 depict the comparison between the exact and numerical solutions of Examples 1 and 2 for $N = 60$ and $N = 30$, respectively.

5.   Conclusions

This work presents QBSCM to obtain the numerical solution of non-linear BT problems. The quintic B-spline interpolation error analysis has been studied, and it gives fourth-order convergence results. The applicability and efficiency of the proposed method have been estimated through numerical examples and obtained numerical results have been compared with the exact solution and other existing methods at the same number of mesh points. The proposed method provides better numerical results than the methods ^[12,13].

Conflict of interest

The authors declare that they have no competing interests.