A Pair of Optimized Nyström Methods with Symmetric Hybrid Points for the Numerical Solution of Second-Order Singular Boundary Value Problems

Abstract

This paper introduces an efficient approach for solving Lane–Emden–Fowler problems. Our method utilizes two Nyström schemes to perform the integration. To overcome the singularity at the left end of the interval, we combine an optimized scheme of Nyström type with a set of Nyström formulas that are used at the fist subinterval. The optimized technique is obtained after imposing the vanishing of some of the local truncation errors, which results in a set of symmetric hybrid points. By solving an algebraic system of equations, our proposed approach generates simultaneous approximations at all grid points, resulting in a highly effective technique that outperforms several existing numerical methods in the literature. To assess the efficiency and accuracy of our approach, we perform some numerical tests on diverse real-world problems, including singular boundary value problems (SBVPs) from chemical kinetics.

1. Introduction

The problem of interest is described by the differential equation

$$q″\left(x\right)+\frac{\lambda }{x}{q}^{\prime }\left(x\right)=k(x,q\left(x\right)),0<x\le x_N=1$$

(1)

with one of the following three different types of boundary conditions (BCs)

$$q\left(0\right)=q_a,q\left(1\right)=q_b,$$

(2)

$$q\left(0\right)=q_a,q\prime \left(1\right)=q_b^{\prime },$$

(3)

$$q\prime \left(0\right)=q_a^{\prime },q\left(1\right)=q_b,$$

(4)

where $\lambda ,q_a,q_b,q_a^{\prime },$ and $q_b^{\prime }$ are given, and the function $k(x,q)$ is assumed to be continuous. It is worth mentioning that the existence and uniqueness of the solution to (1) together with appropriate boundary conditions have been rigorously determined by Zhang [1].

Many researchers, such as Thula and Roul [2], Rufai and Ramos [3,4], and Tunc [5], discussed the wide applicability and the theoretical analysis of second-order singular initial/boundary value ordinary differential equations in different fields. By modeling using the system represented by (1), researchers can gain insights into the behavior and properties of many complex physical systems.

Researchers across various fields of applied sciences and engineering have shown significant interest in solving equations of the Emden–Fowler type, which are represented by (1)–(4). However, these equations pose challenges due to their singularity at $x=0$ and nonlinear characteristics, making them difficult to handle theoretically. In order to overcome these challenges and obtain meaningful solutions, numerical methods have emerged as crucial tools. By discretizing the equations and performing computations on a computational grid, these methods enable researchers to address the complexities arising from the singularity and nonlinearity, providing reasonable and practical approximate solutions.

Considerable efforts have been dedicated to obtain approximate solutions of the aforementioned problems. Various strategies have been used to tackle the challenge posed by the singularity at $x=0$ in Equation (1). Numerical analysts have proposed various numerical methods to solve the type of problem described in Equation (1). Some of these existing numerical methods include the finite difference methods proposed in [6,7], the spline methods discussed in [8,9], the Nyström methods reported in [10], the approximation methods introduced in [11,12], the hybrid block technique in [13], the high-order compact finite differences method in [14], the semi-numerical approach in [15], the pseudospectral method in [16], and the collocation method presented in [17]. The choice of these methods depends on the specific problem characteristics, such as the presence of singularities, smoothness of the solution, computational resources available, and desired accuracy.

Despite the extensive research conducted by numerous scholars to address the challenge of solving SBVPs described by Equations (1)–(4), the accuracies of many existing methods still require improvements. In order to address this issue, this work introduces a novel approach called Pair of Optimized Nyström Methods (PONM). This method is specifically designed to enhance the accuracy of numerical solutions by directly integrating second-order SBVPs. By utilizing PONM, we provide more reliable and accurate results compared to existing methods, contributing to the advancement of solving these challenging boundary value problems.

2. Pair of Optimized Nyström Methods

The method under consideration, namely, a pair of optimized Nyström methods, has been previously introduced in [3] for solving second-order singular initial value problems. The interested reader can refer to this source to better understand the method’s mathematical derivation and its characteristics. This method is represented by the following equations, which are fully explained in [3].

$$\begin{array}{ccc}\hfill & & q_{n+1}=q_n+h{q}_n^{\prime }+\frac{1}{360}{h}^2\left(7\left(\sqrt{21}+7\right)f_{n+u}+64f_{n+v}-7\left(\sqrt{21}-7\right)f_{n+w}+18f_n\right),\hfill \\ & & q_{n+1}^{\prime }=q_n^{\prime }+h\left(\frac{1}{20}{f}_n+\frac{49}{180}{f}_{n+u}+\frac{16}{45}{f}_{n+v}+\frac{49}{180}{f}_{n+w}+\frac{1}{20}{f}_{n+1}\right).\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill q_{n+u}& =& q_n-\frac{1}{14}\left(\sqrt{21}-7\right)h{q}_n^{\prime }\hfill \\ & & +\frac{h^2}{17640}\left(140f_{n+u}+32\left(31-7\sqrt{21}\right)f_{n+v}+7\left(227-49\sqrt{21}\right)f_{n+w}\right)\hfill \\ & & +\frac{h^2}{17640}\left(\left(435-63\sqrt{21}\right)f_n-6f_{n+1}\right),\hfill \\ \hfill q_{n+v}& =& q_n+\frac{h}{2}{q}_n^{\prime }+\frac{h^2}{5760}\left(7\left(8\sqrt{21}+35\right)f_{n+u}+80f_{n+v}+7\left(35-8\sqrt{21}\right)f_{n+w}\right)\hfill \\ & & +\frac{h^2}{5760}\left(147f_n+3f_{n+1}\right),\hfill \\ \hfill q_{n+w}& =& q_n+\frac{\left(\sqrt{21}+7\right)h}{14}{q}_n^{\prime }+\frac{h^2}{17640}7\left(49\sqrt{21}+227\right)f_{n+u}\hfill \\ & & +\frac{h^2}{17640}\left(32\left(7\sqrt{21}+31\right)f_{n+v}+140f_{n+w}+\left(63\sqrt{21}+435\right)f_n-6f_{n+1}\right),\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill q_{n+u}^{\prime }& =& q_n^{\prime }+\frac{h}{17640}\left(7\left(343-9\sqrt{21}\right)f_{n+u}+64\left(49-12\sqrt{21}\right)f_{n+v}+7\left(343-69\sqrt{21}\right)f_{n+w}\right.\hfill \\ & & \left.+9\left(3\sqrt{21}+119\right)f_n+27\left(\sqrt{21}-7\right)f_{n+1}\right),\hfill \\ \hfill q_{n+v}^{\prime }& =& q_n^{\prime }+\frac{h}{2880}\left(7\left(15\sqrt{21}+56\right)f_{n+u}+512f_{n+v}+7\left(56-15\sqrt{21}\right)f_{n+w}+117f_n+27f_{n+1}\right),\hfill \\ \hfill q_{n+w}^{\prime }& =& q_n^{\prime }+\frac{h}{17640}\left(7\left(69\sqrt{21}+343\right)f_{n+u}+64\left(12\sqrt{21}+49\right)f_{n+v}+7\left(9\sqrt{21}+343\right)f_{n+w}\right.\hfill \\ & & \left.+9\left(119-3\sqrt{21}\right)f_n-27\left(\sqrt{21}+7\right)f_{n+1}\right),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill q_1& =& q_0+h{q}_0^{\prime }\hfill \\ \hfill & & +h^2\left(0.2009319137389590f_{\overline{u}}+0.2292411063595862f_{\overline{v}}+0.0698269799014541f_{\overline{w}}\right),\hfill \\ \hfill q_1^{\prime }& =& q_0^{\prime }+h\left(0.2204622111767684f_{\overline{u}}+0.3881934688431719f_{\overline{v}}\right.\hfill \\ \hfill & & \left.+0.3288443199800597f_{\overline{w}}+0.0625f_1\right),\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill q_{\overline{u}}& =& q_0+0.0885879595127039h{q}_0^{\prime }\hfill \\ & & +h^2\left(0.0053826755294719f_{\overline{u}}+0.0024215917832576f_{\overline{v}}\right.\hfill \\ & & \left.+0.001564645634154f_{\overline{w}}-0.0006018160950569f_1\right),\hfill \\ \hfill q_{\overline{v}}& =& q_0+0.4094668644407347h{q}_0^{\prime }\hfill \\ & & +h^2\left(0.0695583040205626f_{\overline{u}}+0.0161202500910538f_{\overline{v}}\right.\hfill \\ & & \left.-0.002478766567991f_{\overline{w}}+0.000631768993838f_1\right),\hfill \\ \hfill q_{\overline{w}}& =& q_0+0.7876594617608470h{q}_0^{\prime }\hfill \\ & & +h^2\left(0.1545378137303791f_{\overline{u}}+0.1448580872610296f_{\overline{v}}\right.\hfill \\ & & \left.+0.0111501356039639f_{\overline{w}}-0.00034232274467f_1\right),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill q_{\overline{u}}^{\prime }& =& q_0^{\prime }+h\left(0.112999479323156f_{\overline{u}}-0.0403092207235222f_{\overline{v}}\right.\hfill \\ & & \left.+0.0258023774203363f_{\overline{w}}-0.0099046765072664f_1\right),\hfill \\ \hfill q_{\overline{v}}^{\prime }& =& q_0^{\prime }+h\left(0.2343839957474002f_{\overline{u}}+0.206892573935358f_{\overline{v}}\right.\hfill \\ & & \left.-0.0478571280485407f_{\overline{w}}+0.0160474228065162f_1\right),\hfill \\ \hfill q_{\overline{w}}^{\prime }& =& q_0^{\prime }+h\left(0.2166817846232503f_{\overline{u}}+0.4061232638673733f_{\overline{v}}\right.\hfill \\ & & \left.+0.189036518170056f_{\overline{w}}-0.024182104899832f_1\right),\hfill \end{array}$$

(10)

where

$$\begin{array}{c}\hfill u=\frac{\left(7-\sqrt{21}\right)}{14},v=\frac{1}{2},w=\frac{\left(7+\sqrt{21}\right)}{14},\end{array}$$

$$\begin{array}{ccc}\hfill \overline{u}& =& \frac{1}{7}\left(3-\sqrt{6}sin\left(\frac{1}{3}{tan}^{-1}\left(7\right)\right)-\sqrt{2}cos\left(\frac{1}{3}{tan}^{-1}\left(7\right)\right)\right)\simeq 0.0885879595127039,\hfill \\ \hfill \overline{v}& =& \frac{1}{7}\left(3+\sqrt{6}sin\left(\frac{1}{3}{tan}^{-1}\left(7\right)\right)-\sqrt{2}cos\left(\frac{1}{3}{tan}^{-1}\left(7\right)\right)\right)\simeq 0.4094668644407347,\hfill \\ \hfill \overline{w}& =& \frac{3}{7}+\frac{2}{7}\sqrt{2}cos\left(\frac{1}{3}{tan}^{-1}\left(7\right)\right)\simeq 0.7876594617608470.\hfill \end{array}$$

To solve the SBVP in (1) subject to any of the boundary conditions in (2)–(4) using PONM, we first select a mesh grid with step-size h. Then, we utilized the formulas presented in (5)–(7) for $n=1\left(1\right)N-1$, in addition to those given in (8)–(10), corresponding to the first integration step. This systematic approach enabled us to obtain a global method that could accurately approximate the solution and complete the integration process over the interval $[0,x_N]$. Through these steps, we could effectively use PONM to give numerical solutions to the challenging nature of the given SBVPs and obtain reliable numerical results.

Convergence Analysis

Firstly, we will define the concept of convergence. We will proceed to show that the PONM is convergent and provide a suitable matrix–vector representation for the formulas in (5)–(7) and (8)–(10). By doing so, we can analyze the convergence of the PONM and prove its effectiveness in approximating the solution to the considered problem. We assume that the exact solution is sufficiently smooth, as needed.

Definition 1. 

Let $q\left(x\right)$ be the exact solution of the considered problem and ${\left\{q_j\right\}}_{j=0}^N$ the approximations provided by the PONM technique. The numerical scheme is said to be convergent of order p if, when $h\to 0$, there is a constant K that is independent of h, satisfying

$$\underset{0\le j\le N}{max}|q\left(x_j\right)-q_j|\le K{h}^p.$$

Remark 1. 

Note that K denotes a generic positive real constant, which, in general, cannot be estimated. This does not affect the convergence result, as, when $h\to 0$, we have

$${\displaystyle \underset{0\le j\le N}{max}|q\left(x_j\right)-q_j|\to 0,}$$

no matter how big K is. This K depends, of course, on the specific problem at hand.

In the following lines, we will analyze the convergence of the proposed method for solving the SBVP given in (1) and (2), but we could have considered any other boundary conditions.

Theorem 1. 

Let $q\left(x\right)$ denote the exact solution of problem (1) together with the BCs given in (2), and ${\left\{q_j\right\}}_{j=0}^N$ represents the approximate solution obtained using the proposed scheme. Then, the proposed method exhibits convergence of, at least, order five.

Proof. 

The matrix M with dimensions $8N\times 8N$ can be defined as follows:

$$M=\left[\begin{array}{cccc}{M}_{1,1}& M_{1,2}& \cdots & M_{1,2N}\\ ⋮& ⋮& & ⋮\\ M_{2N,1}& M_{2N,2}& \cdots & M_{2N,2N}\end{array}\right],$$

where the matrix M comprises various sub-matrices, denoted as $M_{i,j}$, where most of them have dimensions of $4\times 4$. However, there are two exceptions: $M_{i,N+1}$ with dimensions $4\times 3$, and $M_{i,2N}$ with dimensions $4\times 5$. These sub-matrices are given as follows:

$$M_{i,i}=I,i=N+2,\cdots ,2N,I\mathrm{is}\mathrm{a}\mathrm{four}\text{-}\mathrm{dimensional}\mathrm{identity}\mathrm{matrix},$$

$$M_{N,N}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right];M_{i,i-1}=\left[\begin{array}{cccc}0& 0& 0& -1\\ 0& 0& 0& -1\\ 0& 0& 0& -1\\ 0& 0& 0& -1\end{array}\right],i=N+2,\cdots ,2N;$$

$$M_{N+1,N+1}=\left[\begin{array}{ccccc}-1& 1& 0& 0& 0\\ -1& 0& 1& 0& 0\\ -1& 0& 0& 1& 0\\ -1& 0& 0& 0& 1\end{array}\right];$$

$$M_{1,N+1}=h\left[\begin{array}{ccccc}-0.0885879595127039& 0& 0& 0& 0\\ -0.4094668644407347& 0& 0& 0& 0\\ -0.7876594617608471& 0& 0& 0& 0\\ -1.0000000000000000& 0& 0& 0& 0\end{array}\right];$$

$$M_{i,N+i}=h\left[\begin{array}{cccc}0& 0& 0& -\frac{1}{14}\left(7-\sqrt{21}\right)\\ 0& 0& 0& -{\displaystyle \frac{1}{2}}\\ 0& 0& 0& -\frac{1}{14}\left(7+\sqrt{21}\right)\\ 0& 0& 0& -1\end{array}\right],i=1\cdots ,N-1;$$

$$M_{N,2N}=h\left[\begin{array}{ccccc}0& 0& 0& 0& -\frac{1}{14}\left(7-\sqrt{21}\right)\\ 0& 0& 0& 0& -{\displaystyle \frac{1}{2}}\\ 0& 0& 0& 0& -\frac{1}{14}\left(7+\sqrt{21}\right)\\ 0& 0& 0& 0& -1\end{array}\right].$$

The remaining sub-matrices, which were not mentioned previously, were null matrices represented by $M_{i,j}=\mathbb{O}$. These null sub-matrices contributed to the structure of the overall matrix. Furthermore, we introduced another matrix called U, which had dimensions of $8N\times (4N+1)$. The matrix U played a role in the problem and complemented the matrix M in the numerical computations.

$$U=\left[\begin{array}{cccc}{U}_{1,1}& U_{1,2}& \cdots & U_{1,N}\\ ⋮& ⋮& & ⋮\\ U_{2N,1}& U_{2N,2}& \cdots & U_{2N,N}\end{array}\right],$$

where the submatrices, denoted by $U_{i,j}$, have dimensions of $4\times 4$, except for $U_{i,1}$, where i ranges from 1 to $2N$, which has a size of $4\times 5$. The explicit expressions for these submatrices are provided below:

$$U_{1,1}=h\left[\begin{array}{ccccc}0& -0.005382675529& -0.002421591783& -0.001564645635& 0.000601816095\\ 0& -0.069558304020& -0.016120250091& 0.002478766567& -0.000631768993\\ 0& -0.154537813730& -0.144858087261& -0.011150135603& 0.0003423227446\\ 0& -0.200931913738& -0.229241106359& -0.069826979901& 0.0000000000000\end{array}\right];$$

$$U_{i,i}=h\left[\begin{array}{cccc}-\frac{1}{126}& \frac{4\left(-31+7\sqrt{21}\right)}{2205}& \frac{-227+49\sqrt{21}}{2520}& \frac{1}{2940}\\ -\frac{7\left(35+8\sqrt{21}\right)}{5760}& -\frac{1}{72}& \frac{7\left(-35+8\sqrt{21}\right)}{5760}& -\frac{1}{1920}\\ \frac{-227-49\sqrt{21}}{2520}& -\frac{4\left(31+7\sqrt{21}\right)}{2205}& -\frac{1}{126}& \frac{1}{2940}\\ -\frac{7}{360}\left(7+\sqrt{21}\right)& -\frac{8}{45}& \frac{7}{360}\left(-7+\sqrt{21}\right)& 0\end{array}\right],i=2\cdots ,N;$$

$$U_{i,i-1}=h\left[\begin{array}{cccc}0& 0& 0& \frac{-145+21\sqrt{21}}{5880}\\ 0& 0& 0& -\frac{49}{1920}\\ 0& 0& 0& \frac{-145-21\sqrt{21}}{5880}\\ 0& 0& 0& -\frac{1}{20}\end{array}\right],i=3,\cdots ,N;U_{2,1}=h\left[\begin{array}{ccccc}0& 0& 0& 0& \frac{-145+21\sqrt{21}}{5880}\\ 0& 0& 0& 0& -\frac{49}{1920}\\ 0& 0& 0& 0& \frac{-145-21\sqrt{21}}{5880}\\ 0& 0& 0& 0& -\frac{1}{20}\end{array}\right];$$

$$U_{N+1,1}=\left[\begin{array}{ccccc}0& -0.112999479323& 0.040309220723& -0.025802377420& 0.009904676507\\ 0& -0.23438399574& -0.206892573935& 0.047857128048& -0.016047422806\\ 0& -0.216681784623& -0.406123263867& -0.189036518170& 0.024182104899\\ 0& -0.220462211176& -0.388193468843& -0.328844319980& -0.062500000000\end{array}\right];$$

$$U_{N+j,j}=\left[\begin{array}{cccc}\frac{-343+9\sqrt{21}}{2520}& -\frac{8}{45}+\frac{32}{35\sqrt{21}}& \frac{-343+69\sqrt{21}}{2520}& -\frac{3\left(-7+\sqrt{21}\right)}{1960}\\ -\frac{7\left(56+15\sqrt{21}\right)}{2880}& -\frac{8}{45}& \frac{7\left(-56+15\sqrt{21}\right)}{2880}& -\frac{3}{320}\\ \frac{-343-69\sqrt{21}}{2520}& -\frac{8\left(49+12\sqrt{21}\right)}{2205}& \frac{-343-9\sqrt{21}}{2520}& \frac{3\left(7+\sqrt{21}\right)}{1960}\\ -\frac{49}{180}& -\frac{16}{45}& -\frac{49}{180}& -\frac{1}{20}\end{array}\right],j=2,\cdots ,N;$$

$$U_{N+j,j-1}=\left[\begin{array}{cccc}0& 0& 0& \frac{-119-3\sqrt{21}}{1960}\\ 0& 0& 0& -\frac{13}{320}\\ 0& 0& 0& \frac{-119+3\sqrt{21}}{1960}\\ 0& 0& 0& -\frac{1}{20}\end{array}\right],j=3,\cdots ,N;U_{N+2,1}=\left[\begin{array}{ccccc}0& 0& 0& 0& \frac{-119-3\sqrt{21}}{1960}\\ 0& 0& 0& 0& -\frac{13}{320}\\ 0& 0& 0& 0& \frac{-119+3\sqrt{21}}{1960}\\ 0& 0& 0& 0& -\frac{1}{20}\end{array}\right].$$

All submatrices, $U_{i,j}$, which were not mentioned earlier, are null matrices, denoted by $U_{i,j}=\mathbb{O}$.

We can express the vectors of exact values in the following manner:

$$\begin{array}{ccc}\hfill Q_{8N}& =& {\left(q\left(x_{\overline{u}}\right),q\left(x_{\overline{v}}\right),q\left(x_{\overline{w}}\right),q\left(x_1\right),\cdots ,q\left(x_{N-1+w}\right),q\prime \left(x_0\right),q\prime \left(x_{\overline{u}}\right),\cdots ,q\prime \left(x_N\right)\right)}^T,\hfill \\ \hfill F_{4N+1}& =& \left(f(x_0,q\left(x_0\right),q\prime \left(x_0\right)),f(x_{\overline{u}},q\left(x_{\overline{u}}\right),q\prime \left(x_{\overline{u}}\right)),\cdots ,f(x_N,q\left(x_N\right),q\prime \left(x_N\right)))\right..\hfill \end{array}$$

It is important to highlight that the vector Y contains a total of $8N$ terms, but the vector F contains $(4N+1)$ components. This is due to the fact that, as indicated in (2), the boundary conditions specify that the values of $q\left(x\right)$ at $x_0$ and $x_N$ are known, with $q\left(x_0\right)=q_a$, and $q\left(x_N\right)=q_b$. The SBVP can be discretized and approximated using the following formulas:

$$M_{8N\times 8N}{Q}_{8N}+h{U}_{8N\times (4N+1)}{F}_{4N+1}+C_{8N}=LT{\left(h\right)}_{8N}.$$

(11)

In the above equation, the vector $C_{8N}$ is constructed by assembling the predetermined values, given by

$$C_{8N}={\left(-q_a,-q_a,-q_a,-q_a,0,\cdots ,0,q_b,0,\cdots ,0\right)}^T,$$

where T denotes the transpose operation. On the other hand, the vector $LT{\left(h\right)}_{8N}$ represents the Local Truncation Errors (LTEs) associated with the formulas used in the discretization process. These LTEs can be expressed as follows:

$$LT{\left(h\right)}_{8N}\simeq \left[\begin{array}{c}2.805556\times {10}^{-6}{h}^6{q}^{\left(6\right)}\left(x_0\right)+\mathcal{O}\left(h^7\right)\\ 9.161677\times {10}^{-7}{h}^6{q}^{\left(6\right)}\left(x_0\right)+\mathcal{O}\left(h^7\right)\\ -2.962403\times {10}^{-6}{h}^6{q}^{\left(6\right)}\left(x_0\right)+\mathcal{O}\left(h^7\right)\\ 1.417233\times {10}^{-7}{h}^8{q}^{\left(8\right)}\left(x_0\right)+\mathcal{O}\left(h^9\right)\\ \frac{h^7{q}^{\left(7\right)}\left(x_1\right)}{1152480\sqrt{21}}+\mathcal{O}\left(h^8\right)\\ -\frac{h^8{q}^{\left(8\right)}\left(x_1\right)}{30965760}+\mathcal{O}\left(h^9\right)\\ -\frac{h^7{q}^{\left(7\right)}\left(x_1\right)}{1152480\sqrt{21}}+\mathcal{O}\left(h^8\right)\\ \frac{h^9{q}^{\left(9\right)}\left(x_1\right)}{177811200}+\mathcal{O}\left(h^{10}\right)\\ ⋮\\ \frac{h^9{q}^{\left(9\right)}\left(x_{N-1}\right)}{177811200}+\mathcal{O}\left(h^{10}\right)\\ 4.58527\times {10}^{-5}{h}^6{q}^{\left(6\right)}\left(x_0\right)+\mathcal{O}\left(h^7\right)\\ -4.81348\times {10}^{-5}{h}^6{q}^{\left(6\right)}\left(x_0\right)+\mathcal{O}\left(h^7\right)\\ 2.60817\times {10}^{-5}{h}^6{q}^{\left(6\right)}\left(x_0\right)+\mathcal{O}\left(h^7\right)\\ -2.02462\times {10}^{-8}{h}^9{q}^{\left(9\right)}\left(x_0\right)+\mathcal{O}\left(h^{10}\right)\\ \frac{h^7{q}^{\left(7\right)}\left(x_1\right)}{493920}+\mathcal{O}\left(h^8\right)\\ -\frac{h^7{q}^{\left(7\right)}\left(x_1\right)}{322560}+\mathcal{O}\left(h^8\right)\\ \frac{h^7{q}^{\left(7\right)}\left(x_1\right)}{493920}+\mathcal{O}\left(h^8\right)\\ -\frac{h^{10}{y}^{\left(10\right)}\left(x_1\right)}{1422489600}+\mathcal{O}\left(h^{11}\right)\\ ⋮\\ -\frac{h^{10}{q}^{\left(10\right)}\left(x_{N-1}\right)}{1422489600}+\mathcal{O}\left(h^{11}\right)\end{array}\right].$$

Regarding the approximate values, they are obtained from the following system of equations:

$$M_{8N\times 8N}{\overline{Q}}_{8N}+h{U}_{8N\times (4N+1)}{\overline{F}}_{4N+1}+C_{8N}=0,$$

(12)

${\overline{Q}}_{8N}$ being an approximation of the vector $Q_{8N}$, such that:

$${\overline{Q}}_{8N}={\left(q_{\overline{u}},q_{\overline{v}},q_{\overline{w}},q_1,\cdots ,q_{N-1+w},q_0^{\prime },q_{\overline{u}}^{\prime },\cdots ,q_N^{\prime }\right)}^T,$$

and

$${\overline{F}}_{4N+1}={\left(f_0,f_{\overline{u}},f_{\overline{v}},f_{\overline{w}},f_1,f_{1+u},\cdots ,f_N\right)}^T.$$

By subtracting (12) from (11), we obtain the following equation:

$$M_{8N\times 8N}{E}_{8N}+h{U}_{8N\times (4N+1)}\left(F_{4N+1}-{\overline{F}}_{4N+1}\right)=LT{\left(h\right)}_{8N},$$

(13)

where the vector $E_{8N}$ denotes the errors at the discrete points, which are obtained by subtracting ${\overline{Q}}_{8N}$ from the vector $Q_{8N}$. It is represented as a column vector containing the errors $e_{\overline{u}},e_{\overline{v}},\cdots ,e_{N-1+w},e_0^{\prime },e_{\overline{u}}^{\prime },\cdots ,e_N^{\prime }$ at the considered points.

This equation characterizes the truncation errors resulting from the numerical approximation process. By utilizing the Mean Value Theorem, as presented in [18], we have that, for any suitable subindex i, there exists a value ${\xi }_i$ such that:

$$\begin{array}{ccc}& & f(x_i,q\left(x_i\right),q\prime \left(x_i\right))-f(x_i,q_i,q_i^{\prime })=\left(q\left(x_i\right)-q_i\right)\frac{\partial f}{\partial q}\left({\xi }_i\right)+\left(q\prime \left(x_i\right)-q_i^{\prime }\right)\frac{\partial f}{\partial q\prime }\left({\xi }_i\right),\hfill \end{array}$$

where ${\xi }_i$ represents an intermediate point lying on the line segment between $(x_i,q\left(x_i\right),q\prime \left(x_i\right))$ and $(x_i,q_i,q_i^{\prime })$. This implies that:

$$\begin{array}{ccc}& & {(F-\overline{F})}_{4N+1}=\hfill \\ & & \left(\begin{array}{cccccccc}\frac{\partial f}{\partial q}\left({\xi }_0\right)& 0& \cdots & 0& \frac{\partial f}{\partial q\prime }\left({\xi }_0\right)& 0& \cdots & 0\\ 0& \frac{\partial f}{\partial q}\left({\xi }_{\overline{u}}\right)& \cdots & 0& 0& \frac{\partial f}{\partial q\prime }\left({\xi }_{\overline{u}}\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \frac{\partial f}{\partial q}\left({\xi }_N\right)& 0& 0& \cdots & \frac{\partial f}{\partial q\prime }\left({\xi }_N\right)\end{array}\right)\left(\begin{array}{c}{e}_0\\ e_{\overline{u}}\\ ⋮\\ e_N\\ e_0^{\prime }\\ e_{\overline{u}}^{\prime }\\ ⋮\\ e_N^{\prime }\end{array}\right)=\hfill \\ & & \left(\begin{array}{cccccc}0& \cdots & 0& \frac{\partial f}{\partial q\prime }\left({\xi }_0\right)& \cdots & 0\\ \frac{\partial f}{\partial q}\left({\xi }_{\overline{u}}\right)& \cdots & 0& 0& \cdots & 0\\ ⋮& \ddots & 0& ⋮& \ddots & ⋮\\ 0& \cdots & \frac{\partial f}{\partial q}\left({\xi }_{N-1+w}\right)& 0& \cdots & 0\\ 0& \cdots & 0& 0& \cdots & \frac{\partial f}{\partial q\prime }\left({\xi }_N\right)\end{array}\right)\left(\begin{array}{c}{e}_{\overline{u}}\\ e_{\overline{v}}\\ ⋮\\ e_{N-1+w}\\ e_0^{\prime }\\ e_{\overline{u}}^{\prime }\\ ⋮\\ e_N^{\prime }\end{array}\right)=\hfill \\ & & J_{(4N+1)\times 8N}{E}_{8N}.\hfill \end{array}$$

It is important to note that the second equation utilizes the fact that $e_0=q\left(x_0\right)-q_a=0$, and $e_N=q\left(x_N\right)-q_b=0$. Taking this into account, we can reorganize the equation given in (13) in the following manner:

$$\left(M_{8N\times 8N}+h{U}_{8N\times (4N+1)}{J}_{(4N+1)\times 8N}\right)E_{8N}=LT{\left(h\right)}_{8N},$$

(14)

and, taking $D=M+hUJ$, we can write

$$D_{8N\times 8N}{E}_{8N}=LT{\left(h\right)}_{8N}.$$

(15)

For sufficiently small values of $h>0$, we can express Equation (15) as follows:

$$E=D^{-1}\left(LT\left(h\right)\right).$$

(16)

After expanding the components of $D^{-1}$ in powers of h, it is obtained that $\parallel D^{-1}\parallel =\mathcal{O}\left(h^{-1}\right)$. Considering that $q\left(x\right)$ has derivatives that are bounded within the interval $[0,x_N]$ up to the required order, we can derive the following inequality from (16) and the vector $LT\left(h\right)$:

$$\begin{array}{cc}\hfill \parallel E\parallel \le & \parallel D^{-1}\parallel \parallel LT\left(h\right)\parallel \hfill \\ \hfill =& \mathcal{O}\left(h^{-1}\right)\mathcal{O}\left(h^6\right)\hfill \\ \hfill \le & K{h}^5.\hfill \end{array}$$

Therefore, the method suggested exhibits a minimum convergence order of five. This means that the error between the exact solution and the numerical approximation decreases at a rate of at least fifth order as the step size is reduced. □

Remark 2. 

The above result shows a fifth order of convergence as a global method, i.e., considering all the points, including the intermediate ones. However, taking into account the expression of the vector of local truncation errors $LT{\left(h\right)}_{8N}$, we see that, at the grid points (with integer index), the method exhibits a superconvergence order (see [19]):

• $|e_1|=|q\left(x_1\right)-q_1|\le |O\left(h^{-1}\right)\left|\right|O\left(h^8\right)|\le K{h}^7,$
• $|e_i|=|q\left(x_i\right)-q_i|\le |O\left(h^{-1}\right)\left|\right|O\left(h^9\right)|\le K{h}^8,i=2,3,\cdots ,N.$

This behavior can be observed in the numerical experiments (see

Table 1. Order of convergence for Problem (17).

h	Proposed Method	MAE	ROC
$\frac{1}{4}$	PONM	$2.4386\times {10}^{-9}$
$\frac{1}{8}$	PONM	$1.0498\times {10}^{-11}$	$7.8597$
$\frac{1}{16}$	PONM	$4.0634\times {10}^{-14}$	$8.0132$
$\frac{1}{32}$	PONM	$2.2204\times {10}^{-16}$	$7.5157$

, where the calculated rates of convergence are close to eight).

3. Implementation

To implement the PONM, we used a block global approach. We rearranged the system in (12) as $\mathit{F}(\tilde{\mathit{Q}})=\mathbf{0}$, the unknowns being as follows:

$$\begin{array}{ccc}\hfill \tilde{\mathbf{Q}}& =& \left\{q_{\overline{u}},q_{\overline{v}},q_{\overline{w}},q_{\overline{u}}^{\prime },q_{\overline{v}}^{\prime },q_{\overline{w}}^{\prime }\right\}\bigcup {\left\{q_j\right\}}_{j=1\left(1\right)N-1}\bigcup {\left\{q_j^{\prime }\right\}}_{j=0\left(1\right)N}\hfill \\ & & \bigcup {\left\{q_{j+u},q_{j+v},q_{j+w},q_{j+u}^{\prime },q_{j+v}^{\prime },q_{j+w}^{\prime }\right\}}_{j=1\left(1\right)N-1}.\hfill \end{array}$$

To solve the resulting nonlinear equations, we employed the Modified Newton’s Method (MNM) due to the implicit nature of the PONM method. The MNM can be expressed as:

$${\tilde{\mathbf{Q}}}^{\mathbf{i}+\mathbf{1}}={\tilde{\mathbf{Q}}}^{\mathbf{i}}-{\left({\mathbf{J}}^{\mathbf{i}}\right)}^{-\mathbf{1}}{\mathbf{F}}^{\mathbf{i}},$$

where $\mathbf{J}$ is the Jacobian matrix of $\mathbf{F}$. The MNM iterations start with initial values obtained from linear interpolation based on the boundary conditions.

We provide a comprehensive summary of how the PONM is applied to obtain numerical solutions to the SBVPs:

• Select a positive integer N and define the step size h as $h=\frac{x_N-x_0}{N}$ to create the partition $P_N$. This partition consists of the points

  $$P_N=\left\{x_{\overline{u}},x_{\overline{v}},x_{\overline{w}}\right\}\bigcup {\left\{x_j\right\}}_{j=0\left(1\right)N}\bigcup {\left\{x_{j+k}\right\}}_{k=u,v,w;j=1\left(1\right)N-1}.$$
• Utilize the formulas given in (8)–(10) and the ones in (5)–(7), taking $n=1\left(1\right)N-1$ to construct a system whose unknowns are those in $\tilde{\mathbf{Q}}$.
• Combine all the equations from Step 2 within $P_N$ to form a single block matrix equation
• Then, solve the system obtained in the previous step, in order to get an efficient and accurate approximation of the SBVP solution on the grid and intermediate points over $\left[x_0,x_N\right]$.

4. Computational Examples

Here, we present the computational experiments and discussions of the proposed PONM method applied to solve singular model problems described by Equations (1) to (4). The accuracy of the PONM is evaluated using the usual formulas:

$$AE=\parallel q\left(x_j\right)-q_j\parallel ,MAE=\underset{j=0,1,\cdots ,N}{max}\parallel q\left(x_j\right)-q_j\parallel .$$

In the above formulas, absolute error (AE) measures the absolute error at each node as the difference between the theoretical solution $q\left(x_j\right)$ and the approximate solution $q_j$ obtained from the PONM. Maximum absolute error (MAE) is used to indicate the largest deviation between the theoretical and approximate solutions obtained using the PONM over the considered interval. These formulas are useful in determining the accuracy of the PONM in approximating the solutions of the singular model problems.

4.1. Numerical Example 1

Consider the following physical model SBVP problem of the isothermal gas sphere equilibrium, as described in [20]:

$$\begin{array}{c}\hfill q″\left(x\right)+\frac{2}{x}q\prime \left(x\right)+q{\left(x\right)}^5=0,q\left(0\right)=0,q\prime \left(1\right)=\sqrt{\frac{3}{4}},\end{array}$$

(17)

with analytical solution $q\left(x\right)=\sqrt{\frac{3}{3+x^2}}$.

The rate of convergence (ROC) of the proposed PONM is presented in Table 1 using the formula $ROC=-{log}_2\left(\frac{MA{E}_h}{MA{E}_{2h}}\right)$.

Table 2. Comparison of AE for Example (17).

x	PONM	Method in [21]	Method in [20]	Method in [22]
0	$1.50857\times {10}^{-12}$	$1.65000\times {10}^{-6}$	$6.32000\times {10}^{-4}$	$6.52000\times {10}^{-6}$
$0.1$	$1.75970\times {10}^{-12}$	$6.63000\times {10}^{-6}$	$6.27000\times {10}^{-4}$	$6.46000\times {10}^{-6}$
$0.2$	$9.21152\times {10}^{-13}$	$1.59000\times {10}^{-6}$	$6.11000\times {10}^{-4}$	$6.30000\times {10}^{-6}$
$0.3$	$6.06404\times {10}^{-13}$	$1.53000\times {10}^{-6}$	$5.86000\times {10}^{-4}$	$6.05000\times {10}^{-6}$
$0.4$	$4.24327\times {10}^{-13}$	$1.44000\times {10}^{-6}$	$5.28000\times {10}^{-4}$	$5.70000\times {10}^{-6}$
$0.5$	$2.99538\times {10}^{-13}$	$1.34000\times {10}^{-6}$	$5.09000\times {10}^{-4}$	$5.30000\times {10}^{-6}$
$0.6$	$2.06835\times {10}^{-13}$	$1.22000\times {10}^{-6}$	$4.53000\times {10}^{-4}$	$4.84000\times {10}^{-6}$
$0.7$	$1.35669\times {10}^{-13}$	$1.10000\times {10}^{-6}$	$3.82000\times {10}^{-4}$	$4.33000\times {10}^{-6}$
$0.8$	$7.99361\times {10}^{-14}$	$9.58000\times {10}^{-7}$	$2.88000\times {10}^{-4}$	$3.86000\times {10}^{-6}$
$0.9$	$3.57492\times {10}^{-14}$	$7.30000\times {10}^{-7}$	$1.64000\times {10}^{-4}$	$3.24000\times {10}^{-6}$
$1.0$	$0.00000$	$1.89000\times {10}^{-14}$	$0.00000$	$1.45000\times {10}^{-13}$

compares the absolute error (AE) of the methods presented in [20,21,22] with the proposed method, using $h=\frac{1}{10}$. The results in Table 2 demonstrate that the proposed PONM is more accurate than the methods proposed in [20,21,22]. Furthermore, Figure 1 is produced using $h=\frac{1}{16}$, showing good agreement between the numerical and exact solutions.

4.2. Numerical Example 2

Consider the following model SBVP [23]:

$$\begin{array}{c}\hfill q″\left(x\right)+\frac{2}{x}q\prime \left(x\right)-{\varphi }^{2}q{\left(x\right)}^n=0,q\left(1\right)=1,q\prime \left(0\right)=0,x\in [0,1].\end{array}$$

(18)

Although the general exact solution for (18) is not known, the solution for $n=1$ can be expressed as $q\left(x\right)=\frac{sinh\left(x\varphi \right)}{xsinh\left(\varphi \right)}$, where $\varphi$ represents the Thiele modulus. The value of $\varphi$ is determined by the ratio of the reaction rate at the catalyst surface to the diffusion rate through the catalyst pores.

In

Table 3. Comparison of AE for Problem (18).

x	AE with NM	AE with PONM
$0.1$	$1.87474\times {10}^{-9}$	$5.37389\times {10}^{-11}$
$0.2$	$9.16894\times {10}^{-10}$	$1.74590\times {10}^{-11}$
$0.3$	$1.91515\times {10}^{-9}$	$8.02497\times {10}^{-12}$
$0.4$	$2.61712\times {10}^{-9}$	$4.47804\times {10}^{-12}$
$0.5$	$3.29006\times {10}^{-9}$	$2.85583\times {10}^{-12}$
$0.6$	$3.95964\times {10}^{-9}$	$1.94503\times {10}^{-12}$
$0.7$	$4.50469\times {10}^{-9}$	$1.28325\times {10}^{-12}$
$0.8$	$4.59898\times {10}^{-9}$	$6.94778\times {10}^{-13}$
$0.9$	$3.55123\times {10}^{-9}$	$1.87517\times {10}^{-13}$
1	$0.00000$	$0.00000$

, the Absolute Error (AE) of the Nyström method (NM) presented in [23] is compared with the proposed method using $h=\frac{1}{10}$, $n=1$, and $\varphi =5$. The results in Table 3 indicate that the proposed PONM is more accurate than the NM method proposed in [23]. Additionally, Figure 2 is generated using $h=\frac{1}{20}$, $n=1$, and $\varphi =2$, demonstrating good agreement between the numerical solution and the exact solution.

4.3. Numerical Example 3

Consider the SBVP for the thermal explosion in a cylindrical vessel presented in [24]:

$$\begin{array}{c}\hfill q″\left(x\right)+\frac{1}{x}q\prime \left(x\right)-e^{q\left(x\right)}=0,q\left(0\right)=0,q\prime \left(1\right)=0.\end{array}$$

(19)

The exact solution of this problem is $q\left(x\right)=2log\left(\frac{2\sqrt{6}+1-5}{\left(2\sqrt{6}-5\right)x^2+1}\right)$.

We solve Problem (19) numerically using the new PONM for various stepsize $\left(h\right)$ values. The summarized numerical results and the comparisons of MAE between PONM and the method in [24] are presented in

Table 4. Comparison of MAE for Example (19).

h	Methods	MAE
$\frac{1}{8}$	PONM	$1.51282\times {10}^{-11}$
$\frac{1}{8}$	Method in [24]	$8.53810\times {10}^{-10}$
$\frac{1}{16}$	PONM	$2.33730\times {10}^{-13}$
$\frac{1}{16}$	Method in [24]	$2.19100\times {10}^{-11}$
$\frac{1}{32}$	PONM	$2.58127\times {10}^{-15}$
$\frac{1}{32}$	Method in [24]	$3.92400\times {10}^{-13}$

. It is worth noting that the PONM performs better than the technique in [24]. Additionally, Figure 3 is obtained using $h=\frac{1}{8}$, and it demonstrates a good agreement between the numerical and exact solutions.

4.4. Numerical Example 4

Consider the SBVP for the thermal explosion arising in chemistry and chemical kinetics presented in [25] and the following physical model describing it:

$$\begin{array}{c}\hfill q″\left(x\right)+\frac{\lambda }{x}q\prime \left(x\right)={\varphi }^{2}q\left(x\right)exp\left(\frac{sr(1-q(x\left)\right)}{1+c(1-q(x\left)\right)}\right),q\left(0\right)=0,q\prime \left(1\right)=0.\end{array}$$

(20)

The true solution of this SBVP is not available.

Problem (20) is solved using the PONM method for $\varphi =r=s=c=1$. The numerical results obtained by PONM, the spline method (SM), and the Adomian decomposition method (ADM) in [25,26] are presented in

Table 5. Comparison of numerical solutions for Example (20).

x	PONM, $\mathit{h}=\frac{1}{10}$	SM in [26], $\mathit{h}=\frac{1}{10}$	ADM in [25], $\mathit{h}=\frac{1}{50}$
$0.1$	$0.8383648968813362$	$0.838364878696000$	$0.83836491959750$
$0.2$	$0.8431842515033859$	$0.843184233589000$	$0.84318428772800$
$0.3$	$0.8512302074809177$	$0.851230190133000$	$0.85123026453741$
$0.4$	$0.8625224114670785$	$0.862522394979000$	$0.86252249405263$
$0.5$	$0.8770865272385724$	$0.877086511985000$	$0.87708663616863$
$0.6$	$0.8949523006665648$	$0.894952287104000$	$0.89495243141739$
$0.7$	$0.9161509382762109$	$0.916150926969000$	$0.91615107927542$
$0.8$	$0.9407117001190575$	$0.940711691749000$	$0.94071183074544$
$0.9$	$0.9686575885214240$	$0.968657583887000$	$0.96865767679753$
$1.0$	$1.0000000000000000$	$1.0000000000000000$	$1.000000000000000$

. It is evident from Table 4 that the PONM method outperforms the other methods significantly.

4.5. Numerical Example 5

Consider the following SBVP:

$$\begin{array}{c}\hfill q″\left(x\right)+\frac{0.5}{x}q\prime \left(x\right)+e^{2q\left(x\right)}-0.5e^{q\left(x\right)}=0,q\left(0\right)=log\left(2\right),q\left(1\right)=0.\end{array}$$

(21)

The exact solution of this problem is $q\left(x\right)=log\left(\frac{2}{x^2+1}\right)$.

The new PONM was used to solve Problem (21) numerically for different values of the stepsize $\left(h\right)$. Table 4 presents the numerical results and compares PONM with [14] in terms of MAE. It is evident that PONM provides better accuracy than the method in [14].

5. Conclusions

The paper proposed a Pair of Optimized Nyström Methods (PONM) to solve the SBVP of the Lane–Emden type, which is expressed in (1)–(4). The primary aim of this approach is to provide more accurate approximate solutions to some physical model SBVP problems. The numerical solutions obtained from the proposed PONM method are presented in Table 2, Table 3, Table 4, Table 5 and

Table 6. Comparison of MAE for Example (21).

h	Methods	MAE
$\frac{1}{8}$	PONM	$3.38905\times {10}^{-10}$
$\frac{1}{8}$	Method in [14]	$5.62700\times {10}^{-5}$
$\frac{1}{16}$	PONM	$1.55342\times {10}^{-12}$
$\frac{1}{16}$	Method in [14]	$1.45823\times {10}^{-7}$
$\frac{1}{32}$	PONM	$6.21724\times {10}^{-15}$
$\frac{1}{32}$	Method in [14]	$1.37249\times {10}^{-9}$

and Figure 2 and Figure 3. These results demonstrate that the proposed PONM method performs better than various existing numerical techniques used for comparison. Hence, it can be concluded that the PONM method proposed in this study is an efficient numerical method to solve SBVPs of the Lane–Emden type and other similar problems in diverse fields of science and engineering.