An Effective Approximation Algorithm for Second-Order Singular Functional Differential Equations

Abstract

In this research study, a novel computational algorithm for solving a second-order singular functional differential equation as a generalization of the well-known Lane–Emden and differential-difference equations is presented by using the Bessel bases. This technique depends on transforming the problem into a system of algebraic equations and by solving this system the unknown Bessel coefficients are determined and the solution will be known. The method is tested on several test examples and proves to provide accurate results as compared to other existing methods from the literature. The simplicity and robustness of the proposed technique drive us to investigate more of their applications to several similar problems in the future.

1. Introduction

The primary concern of this research work is to develop a computationally effective technique, which relies on novel Bessel polynomials and a set of collocation points to find the solutions of the following second-order singular functional differential equations (SFDEs)

$$\frac{d^2}{d{x}^2}w(ax+b)+\frac{\mu }{x}\frac{d}{dx}w(cx+d)+r\left(x\right)w(ex+f)=p\left(x\right),0<x\le \ell ,$$

(1)

where the constants $a,b,c,d,e,f$, and $\mu$ are in $\mathbb{R}$. Moreover, ℓ is a positive real number and $r\left(x\right),p\left(x\right)$ are given real-valued functions. The above SFDEs are accompanied with the initial conditions

$$w\left(0\right)=w_0,\frac{d}{dx}w\left(0\right)=w_1.$$

(2)

The study of the singular functional differential equations (SFDEs) is one of the most important areas of study with a variety of applications in engineering and science topics. The researcher community is continuously investigating the possible applications of these types of equations and numerous fields are reported including electrodynamics [1], models of infectious diseases [2], population growth models [3], the simulation of tumor growth [4], the processing of chemical systems [5], understanding the gene system [6], and viral infectious models [7]. These models have attracted the attention of many scientists with their singularity at the origin or other points. One of the most important models that has this type of equation is the well-known Lane–Emden type equations (LEE) that have been used since they was named after the famous astrophysicists Jonathan Homer Lane and Robert Emden back in 1870. The Lane–Emden equations have some applications in the field of models of thermal explosion [8], models of isothermal gas spheres [9], stellar structure [10], and the study of thermionic current. These potential and important applications have motivated researchers to investigate the solution of these models more.

Finding the approximate solution of the LEE is one of the most interesting subjects for scientists and researchers. For obtaining the results for the solutions of these types of problems, some numerical and analytical techniques have been adapted, acquiring good results. For example, Mirazaee et al. in [11] investigated a Fibonacci polynomials based method for solving this problem. Moreover, Kadalbajoo et al. [12] adapted a method with Taylor series expansion for solving a similar type of equation. The asymptotic solutions are analyzed in [13] for a class of nonlinear singular perturbed equations. Sabir et al. [14] performed a neuro-swarm intelligent computing algorithm for solving a second-order type equation. Other reported methods can be found in [15,16,17,18,19] and references therein with multiple other applications of this problem. All of the above-mentioned methods, either analytical or numerical, have some advantages and disadvantages in terms of errors or computational costs. Thus, we searched for a suitable collocation based method to adapt and acquire more accurate results which urges us to use the Bessel collocation method.

The appearance of new Bessel functions related to the Bessel functions of the first type was systematically shown in a seminal paper [20]. Since then, many research papers have been devoted to discovering many characteristics of these polynomials from the algebraic point of view; for more details, see [20,21,22,23,24,25,26]. In recent years, the study of different models with the aid of these polynomials has witnessed a large increase of research due to their simplicity and the ability to provide good results. For the recent applications of the Bessel polynomials, we draw your attention to some recent works we have done in [27,28,29]. The main goal of this research work is to propose a spectral approach based on a combination of novel Bessel bases as well as some appropriate collocation points for an approximate treatment of the SFDEs in (1). Supposedly, the underlying model problem has a solution in terms of the Bessel series expansion on $[0,\ell ]$, the proper representations of the unknowns and their derivatives yield to find the unknown series coefficients, through solving an algebraic system of equations. Indeed, the explicit and original representation of the Bessel functions is

$${\mathcal{B}}_m\left(\xi \right)=\sum _{\sigma =0}^m\frac{2^{-\sigma }(m+\sigma )!}{\sigma !(m-\sigma )!}{\xi }^{\sigma },m=0,1,\dots .$$

(3)

It can be obviously observed that the coefficients of ${\mathcal{B}}_m\left(\xi \right)$ are all positive. Contrary to the Bessel functions of the first kind [30], the novel Bessel functions ${\mathcal{B}}_m\left(\xi \right)$ are the unique solutions of the following differential equation [27,28,29]

$${\xi }^2{\mathcal{B}}_m^{″}\left(\xi \right)+(2\xi +2){\mathcal{B}}_m^{\prime }\left(\xi \right)=m(m+1){\mathcal{B}}_m\left(\xi \right).$$

The main capability of of the proposed Bessel spectral algorithm is that it converts the SFDEs (1) to a system of algebraic equations while reducing computational complexity. Usages of the proposed technique but with different bases such as Legendre, Chebyshev, Chelyshkov, alternative Bessel, and Jacobi functions can be found in [31,32,33,34,35,36,37,38,39,40].

The outline of this study is structured as follows. The next section provides the methodology of the Bessel matrix procedure for the SFDEs in an elegant manner. In Section 3, three test examples are solved in order to evaluate the reliability and accuracy of the presented matrix technique. In Section 4, we present a summary and conclusion.

2. The Bessel Matrix Technique

In order to utilize the Bessel functions on $[0,\ell ]$, we first make a change of variable in (3). We take $\xi =x/\ell$ where $\ell >0$ in (3). Henceforth, the shifted Bessel functions will be denoted by ${\mathcal{B}}_{m,\ell }\left(x\right)$, which are orthogonal with respect to $g_{\ell }\left(x\right):=exp(-2\ell /x)$; see [20,27]. Suppose that the unknown solution $w\left(x\right)$ of (1) can be expanded in the Bessel polynomials form

$$w_N\left(x\right)=\sum _{m=0}^N{c}_m{\mathcal{B}}_{m,\ell }\left(x\right),0\le x\le \ell .$$

(4)

Now, the ultimate goal is to seek the coefficients $c_m$ for $m=0,1,\dots ,N$. To proceed, we introduce the unknown vector

$${\mathit{C}}_N={[c_0{c}_1\dots c_N]}^T,$$

and the vector containing the shifted Bessel functions of order $m=0$ to $m=N$ as

$${\mathit{V}}_N\left(x\right)=\left[{\mathcal{B}}_{0,\ell }\left(x\right){\mathcal{B}}_{1,\ell }\left(x\right)\dots {\mathcal{B}}_{N,\ell }\left(x\right)\right].$$

Thus, on the other hand, we have the following expression for $w_N\left(x\right)$ in (4) in the form

$$w_N\left(x\right)={\mathit{V}}_N\left(x\right){\mathit{C}}_N.$$

(5)

Next, we introduce the following Bessel matrix L, which has a low-triangular structure and is of size $(N+1)\times (N+1)$

$$\mathit{L}={\left(\begin{array}{cccccc}1& 0& 0& \dots & 0& 0\\ 1& 1& 0& \dots & 0& 0\\ 1& 3& 3& \dots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & \ddots & ⋮\\ 1& \frac{2^{-1}N!}{{\ell }^1(N-2)!1!}& \frac{2^{-2}(N+1)!}{{\ell }^2(N-3)!2!}& \dots & \frac{2^{1-N}(2N-2)!}{{\ell }^{N-1}0!(N-1)!}& 0\\ 1& \frac{2^{-1}(N+1)!}{{\ell }^1(N-1)!1!}& \frac{2^{-2}(N+2)!}{{\ell }^2(N-2)!2!}& \dots & \frac{2^{1-N}(2N-1)!}{{\ell }^{N-1}1!(N-1)!}& {\displaystyle \frac{2^{-N}\left(2N\right)!}{{\ell }^{N}0!N!}}\end{array}\right)}_{(N+1)\times (N+1)}.$$

This allows us to write further the vector ${\mathit{V}}_N\left(x\right)$ in (5) in the product representation form

$${\mathit{V}}_N\left(x\right)={\mathbf{\Xi }}_N\left(x\right){\mathit{L}}^T,$$

(6)

where ${\mathbf{\Xi }}_N\left(x\right)=\left[1x{x}^2\dots x^N\right]$ stands for the vector of monomials.

Finally, we place relation (6) into (5) to arrive at the following form of approximate solution $w_N\left(x\right)$ in (4) as

$$w_N\left(x\right)={\mathit{V}}_N\left(x\right){\mathit{C}}_N={\mathbf{\Xi }}_N\left(x\right){\mathit{L}}^T{\mathit{C}}_N.$$

(7)

We now mention a pertinent result about the convergence of Bessel functions. The following theorem asserts that the approximate solution $w_N\left(x\right)$ is exponentially convergent (in the weighed $L_2$ norm) to the exact solution $w\left(x\right)$ if we let N tend to infinity.

Theorem 1.

Let us denote by $w_N\left(x\right)={\mathit{V}}_N\left(x\right){\mathit{C}}_N$ the best square approximation to $w\left(x\right)$ and also we have $w\left(x\right)\in C^{N+1}(0,\ell ]$. Then, the following error estimate holds

$$\begin{array}{c}\hfill {\parallel w\left(x\right)-w_N\left(x\right)\parallel }_{g_{\ell }}\le \frac{{\ell }^{N+\frac{3}{2}}}{\sqrt{2N+3}}\frac{e^{-1}{K}^{\infty }}{(N+1)!},\end{array}$$

where $K^{\infty }:={max}_{x\in (0,\ell ]}|w^{(N+1)}\left(x\right)|$.

Proof. 

For a similar proof, we refer interested readers to [27,30].    □

We now associate a sequence of collocation points ${\left\{x_q\right\}}_{q=0}^N$ on $[0,\ell ]$ to our approximation algorithm. In this respect, the following grid points are employed

$$x_q=\xi +\frac{(\ell -\xi )q}{N},q=0,1,\dots ,N,0<\xi <\ell .$$

(8)

Once the preceding collocation points are placed into the relation (7), we get the matrix expression for the unknown solution itself as

$$\mathit{W}=\chi {\mathit{L}}^T{\mathit{C}}_N,\chi =\left(\begin{array}{c}{\mathbf{\Xi }}_N\left(x_0\right)\\ {\mathbf{\Xi }}_N\left(x_1\right)\\ ⋮\\ {\mathbf{\Xi }}_M\left(x_N\right)\end{array}\right),\mathit{W}=\left(\begin{array}{c}{w}_N\left(x_0\right)\\ w_N\left(x_1\right)\\ ⋮\\ w_N\left(x_N\right)\end{array}\right).$$

(9)

In order to express $w(\alpha x+\beta )$ and its derivatives in the matrix forms, we state and prove the next theorem.

Theorem 2.

For any constants α and β, the matrix representations of $w(\alpha x+\beta )$, $\frac{d}{dx}w(\alpha x+\beta )$, and $\frac{d^2}{d{x}^2}w(\alpha x+\beta )$ at the collocation points (8) can be represented as

$$\begin{array}{ccc}\hfill {\mathit{W}}_{\alpha ,\beta }& =\hfill & \chi {\mathit{H}}_{\alpha ,\beta }^T{\mathit{L}}^T{\mathit{C}}_N,\end{array}$$

(10)

$$\begin{array}{ccc}{\dot{\mathit{W}}}_{\alpha ,\beta }& =\hfill & \chi {\mathit{H}}_{\alpha ,\beta }^T{\mathit{B}}^T{\mathit{L}}^T{\mathit{C}}_N,\end{array}$$

(11)

$$\begin{array}{ccc}{\ddot{\mathit{W}}}_{\alpha ,\beta }& =\hfill & \chi {\mathit{H}}_{\alpha ,\beta }^T{\left({\mathit{B}}^T\right)}^2{\mathit{L}}^T{\mathit{C}}_N.\end{array}$$

(12)

Here, the matrix $\mathit{B}$ is defined in (15) and the matrix ${\mathit{H}}_{\alpha ,\beta }$ is given at (14). Moreover, we have

$${\mathit{W}}_{\alpha ,\beta }=\left(\begin{array}{c}{w}_N(\alpha x_0+\beta )\\ w_N(\alpha x_1+\beta )\\ ⋮\\ w_N(\alpha x_N+\beta )\end{array}\right),{\dot{\mathit{W}}}_{\alpha ,\beta }=\left(\begin{array}{c}{w}_N^{\prime }(\alpha x_0+\beta )\\ w_N^{\prime }(\alpha x_1+\beta )\\ ⋮\\ w_N^{\prime }(\alpha x_N+\beta )\end{array}\right),{\ddot{\mathit{W}}}_{\alpha ,\beta }=\left(\begin{array}{c}{w}_N^{\prime \prime }(\alpha x_0+\beta )\\ w_N^{\prime \prime }(\alpha x_1+\beta )\\ ⋮\\ w_N^{\prime \prime }(\alpha x_N+\beta )\end{array}\right).$$

Proof. 

According to (7), we may write

$$w_N(\alpha x+\beta )={\mathbf{\Xi }}_N(\alpha x+\beta ){\mathit{L}}^T{\mathit{C}}_N.$$

(13)

Our aim is to express ${\mathbf{\Xi }}_N(\alpha x+\beta )$ in terms of ${\mathbf{\Xi }}_N\left(x\right)$, which is defined in (6). With the help of the binomial expansion,

$${(\alpha x+\beta )}^n=\sum _{j=0}^n\left(\genfrac{}{}{0pt}{}{n}{j}\right){\alpha }^j{\beta }^{n-j}{x}^j,$$

we can represent the vector ${\mathbf{\Xi }}_N(\alpha x+\beta )$ as follows

$${\mathbf{\Xi }}_N(\alpha x+\beta )={\mathbf{\Xi }}_N\left(x\right){\mathit{H}}_{\alpha ,\beta }^T,$$

(14)

where the matrix ${\mathit{H}}_{\alpha ,\beta }^T$ is dependent on two parameters $\alpha$ and $\beta$ and is defined as

$${\mathit{H}}_{\alpha ,\beta }^T=\left(\begin{array}{ccccc}\left(\genfrac{}{}{0pt}{}{0}{0}\right){\alpha }^0{\beta }^0& \left(\genfrac{}{}{0pt}{}{1}{0}\right){\alpha }^0{\beta }^1& \left(\genfrac{}{}{0pt}{}{2}{0}\right){\alpha }^0{\beta }^2& \dots & \left(\genfrac{}{}{0pt}{}{N}{0}\right){\alpha }^0{\beta }^N\\ 0& \left(\genfrac{}{}{0pt}{}{1}{1}\right){\alpha }^1{\beta }^0& \left(\genfrac{}{}{0pt}{}{2}{1}\right){\alpha }^1{\beta }^1& \dots & \left(\genfrac{}{}{0pt}{}{N}{1}\right){\alpha }^1{\beta }^{N-1}\\ 0& 0& \left(\genfrac{}{}{0pt}{}{2}{2}\right){\alpha }^2{\beta }^0& \dots & \left(\genfrac{}{}{0pt}{}{N}{2}\right){\alpha }^2{\beta }^{N-2}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & \left(\genfrac{}{}{0pt}{}{N}{N}\right){\alpha }^N{\beta }^0\end{array}\right).$$

Now, we combine the relations (13) and (14) to obtain

$$w_N(\alpha x+\beta )={\mathbf{\Xi }}_N\left(x\right){\mathit{H}}_{\alpha ,\beta }^T{\mathit{L}}^T{\mathit{C}}_N.$$

We are now ready to put the collocation points (8) into the preceding equation followed by utilizing the relation (9); the proof of (10) is done.

We then find a relationship between ${\mathbf{\Xi }}_N\left(x\right)$ and $\frac{d^s}{d{x}^s}{\mathbf{\Xi }}_N\left(x\right)$ for $s=1,2$. For this purpose, we define the differentiation matrix $\mathit{B}$ through defining

$${\mathit{B}}^T={\left(\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 2& \dots & 0\\ ⋮& ⋮& 0& ⋮& ⋮\\ 0& 0& 0& \ddots & N\\ 0& 0& 0& \dots & 0\end{array}\right)}_{(N+1)\times (N+1)}.$$

It can be easily seen that [30]

$$\begin{array}{c}\hfill \frac{d}{dx}{\mathbf{\Xi }}_N\left(x\right)={\mathbf{\Xi }}_N\left(x\right){\mathit{B}}^T,\end{array}$$

(15)

Differentiating the relation (15) once more, we get

$$\begin{array}{c}\hfill \frac{d^2}{d{x}^2}{\mathbf{\Xi }}_N\left(x\right)={\mathbf{\Xi }}_N\left(x\right){\left({\mathit{B}}^T\right)}^2.\end{array}$$

(16)

Our next aim is to differentiate (7) with respect to the variable x and utilize (15). Therefore, we will get the following approximation for $\frac{d}{dx}{w}_N(\alpha x+\beta )$ as

$$\begin{array}{c}\hfill \frac{d}{dx}{w}_N(\alpha x+\beta )\approx \frac{d}{dx}{w}_N\left(x\right){|}_{x\to (\alpha x+\beta )}={\mathbf{\Xi }}_N\left(x\right){|}_{x\to (\alpha x+\beta )}{\mathit{B}}^T{\mathit{L}}^T{\mathit{C}}_N.\end{array}$$

(17)

In the same manner, after using (16) we get an approximation for $\frac{d^2}{d{x}^2}{w}_N(\alpha x+\beta )$ as follows

$$\begin{array}{c}\hfill \frac{d^2}{d{x}^2}{w}_N(\alpha x+\beta )\approx \frac{d^2}{d{x}^2}{w}_N\left(x\right){|}_{x\to (\alpha x+\beta )}={\mathbf{\Xi }}_N\left(x\right){|}_{x\to (\alpha x+\beta )}{\left({\mathit{B}}^T\right)}^2{\mathit{L}}^T{\mathit{C}}_N.\end{array}$$

(18)

We now replace the vector ${\mathbf{\Xi }}_N(\alpha x+\beta )$ in (17) and (18) via (14). Thus, we get the following matrix representation forms

$$\begin{array}{cc}\hfill \frac{d}{dx}{w}_N(\alpha x+\beta )& ={\mathbf{\Xi }}_N\left(x\right){\mathit{H}}_{\alpha ,\beta }^T{\mathit{B}}^T{\mathit{L}}^T{\mathit{C}}_N,\hfill \\ \hfill \frac{d^2}{d{x}^2}{w}_N(\alpha x+\beta )& ={\mathbf{\Xi }}_N\left(x\right){\mathit{H}}_{\alpha ,\beta }^T{\left({\mathit{B}}^T\right)}^2{\mathit{L}}^T{\mathit{C}}_N.\hfill \end{array}$$

(19)

The proofs of (11) and (12) are straightforward after inserting the collocation points (8) into the foregoing (19).    □

Remark 1.

It is worth mentioning the cost of calculating ${\mathbf{\Xi }}_N(\alpha x+\beta )$ as well as $\frac{d^s}{d{x}^s}{\mathbf{\Xi }}_N(\alpha x+\beta )$ for $s=1,2$ in (19) is $\mathcal{O}\left({(N+1)}^2\right)$. For practical applications, however, it would be very useful to consider an algorithm with linear complexity $\mathcal{O}(N+1)$. This task can be accomplished by taking direct differentiation of the monomials ${\mathbf{\Xi }}_N\left(x\right)$. Upon calling Algorithm 1, the s-order derivatives $(s\ge 1,s\in \mathbb{N})$ of ${\mathbf{\Xi }}_N\left(x\right)$ can be achieved directly. Let us assume that Algorithm 1 takes $(N=4,s)$ as inputs. The outputs with $s=1$ and $s=2$ are as follows, respectively:

$$\begin{array}{c}\hfill \frac{d}{dx}{\mathbf{\Xi }}_4\left(x\right)=\left[012x3x^{2}4x^3\right],\frac{d^2}{d{x}^2}{\mathbf{\Xi }}_4\left(x\right)=\left[0026x12x^2\right].\end{array}$$

In this work, we are particularly interested in computing the first and second derivatives of ${\mathbf{\Xi }}_N\left(x\right)$ by Algorithm 1. This enables us to compute $\frac{d^s}{d{x}^s}{\mathbf{\Xi }}_N(\alpha x+\beta )$ for $s=0,1,2$ after invoking Algorithm 1 followed by substituting $x\to (\alpha x+\beta )$ in them.

Algorithm 1: The computation of s-derivative of the vector ${\mathbf{\Xi }}_N\left(x\right)$.

procedure$\left[{\mathbf{\Xi }}_N^{\left(s\right)}\right]$= compute_DerX$(N,s)$ ${\mathbf{\Xi }}_N^{\left(s\right)}\left[1\right]:=0$; for $j:=1,\dots ,N$do  if ($j-s<0$) then   ${\mathbf{\Xi }}_N^{\left(s\right)}[j+1]:=0$;  else   $\phantom{(}$$\phantom{(}$${\displaystyle {\mathbf{\Xi }}_N^{\left(s\right)}[j+1]:=\frac{j!}{(j-s)!}{x}^{j-s};}$  end if end for end;

The considered SFDE problem (1) will be collocated at the set of collocation points (8) to arrive at

$$\frac{d^2}{d{x}^2}w(a{x}_q+b)+\frac{\mu }{x_q}\frac{d}{dx}w(c{x}_q+d)+r\left(x_q\right)w(e{x}_q+f)=p\left(x_q\right),q=0,1,\dots ,N.$$

(20)

To express the former $N+1$ equations in a matrix representation, we exploit the results of Theorem 2 to have

$${\ddot{\mathit{W}}}_{a,b}+\mathit{Q}{\dot{\mathit{W}}}_{c,d}+\mathit{R}{\mathit{W}}_{e,f}=\mathit{P}.$$

(21)

Two matrices $\mathit{Q},\mathit{R}$ and the vector $\mathit{P}$ are

$$\mathit{Q}=\left(\begin{array}{cccc}\frac{\mu }{x_0}& 0& \dots & 0\\ 0& \frac{\mu }{x_1}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \frac{\mu }{x_N}\end{array}\right),\mathit{R}=\left(\begin{array}{cccc}r\left(x_0\right)& 0& \dots & 0\\ 0& r\left(x_1\right)& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & r\left(x_N\right)\end{array}\right),\mathit{P}=\left(\begin{array}{c}p\left(x_0\right)\\ p\left(x_1\right)\\ ⋮\\ p\left(x_N\right),\end{array}\right),$$

(22)

We finally obtain the so-called fundamental matrix equation for the underlying model (1).

Lemma 1.

Suppose that the solution of model problem (1) can be expanded in Bessel function form (4) or (7). Then we have

$$\underset{\mathit{Z}}{\underbrace{\left(\chi {\mathit{H}}_{a,b}^T{\left({\mathit{B}}^T\right)}^2+\mathit{Q}\chi {\mathit{H}}_{c,d}^T{\mathit{B}}^T+\mathit{R}\chi {\mathit{H}}_{e,f}^T\right){\mathit{L}}^T}}{\mathit{C}}_N=\mathit{P},or[\mathit{Z};\mathit{P}].$$

(23)

Proof. 

To conclude the fundamental matrix Equation (23), it is sufficient to place relations (10), (11), and (12) into (21). □

Let us emphasize that the algebraic matrix Equation (23) is linear and the unknown coefficient ${\mathit{C}}_N$ can be calculated through solving it. However, the initial conditions (2) must be utilized in implementation of the fundamental matrix Equation (23). This aim will be considered below.

Initial Conditions in the Matrix Form

Analogously, we are able to approximate the initial conditions (2) in the matrix form, which allows us to find the solution of (1) via solving the fundamental matrix Equation (23). We first convert $w\left(0\right)=w_0$ into a matrix representation. To this end, let us tend $t\to 0$ in (7) to arrive at

$${\widehat{\mathit{Z}}}_0:={\mathbf{\Xi }}_N\left(0\right){\mathit{L}}^T{\mathit{C}}_N=w_0,\mathrm{or}[{\widehat{\mathit{Z}}}_0;w_0].$$

For the second condition $\frac{d}{dx}w\left(0\right)=w_1$, we first differentiate (7). Hence, we combine the resultant equation with (15) to get

$$\frac{d}{dx}{w}_N\left(x\right)=\frac{d}{dx}{\mathbf{\Xi }}_N\left(x\right){\mathit{L}}^T{\mathit{C}}_N={\mathbf{\Xi }}_N\left(x\right){\mathit{B}}^T{\mathit{L}}^T{\mathit{C}}_N.$$

Now, let $t\to 0$ in the foregoing equation to obtain

$${\widehat{\mathit{Z}}}_1:={\mathbf{\Xi }}_N\left(0\right){\mathit{B}}^T{\mathit{L}}^T{\mathit{C}}_N=w_1,\mathrm{or}[{\widehat{\mathit{Z}}}_1;w_1].$$

Now, to the fundamental matrix Equation (23), we add the initial conditions (2). For this purpose, the replacements of the first and the second row of the matrix $[\mathit{Z};\mathit{P}]$ are done by the row matrices $[{\widehat{\mathit{Z}}}_0;w_0]$ and $[{\widehat{\mathit{Z}}}_1;w_1]$. Let us denote by $[\widehat{\mathit{Z}};\widehat{\mathit{P}}]$ the modified version of the fundamental matrix. After solving this modified form, the unknown coefficients $c_m$, $m=0,1,\dots ,N$ will be calculated and thus the desired approximation $w_N\left(x\right)$ of SFDEs (1) will be determined.

3. Computational Simulations

Let us illustrate the practicability of our Bessel matrix approach through numerical simulations. In this respect, computational results of three test examples are performed to show the reliability and validity of the proposed numerical model. Thus, we show that our suggested approximation algorithm can produce improved results compared to existing available computational procedures. For numerical simulations, we use MATlAB software version 2017a for programming and visualization. Moreover, the value $\xi =0.01$ is taken in the set of collocation points (8) in order to ensure that the zero point is excluded from this set. In order to evaluate the accuracy as well as the convergence of the proposed Bessel matrix technique, we define

$$\begin{array}{cc}\hfill & {\mathbb{E}}_N\left(x\right):=|w_N\left(x\right)-w\left(x\right)|,x\in [0,\ell ],\hfill \\ \hfill & {\mathbb{R}}_N\left(x\right):={\mathbb{E}}_N\left(x\right)/\left|w\right(x\left)\right|,x\in [0,\ell ].\hfill \end{array}$$

Example 1.

We firstly consider a non-homogeneous singular differential difference model problem of the form

$$\frac{d^2}{d{x}^2}w(3x-1)+\frac{2}{x}\frac{d}{dx}w\left(2x\right)+xw(x+1)=-6+44x+3x^2+3x^3+x^4,0\le x\le \ell ,$$

with initial conditions $w\left(0\right)=1$, $\frac{d}{dx}w\left(0\right)=0$. A straightforward calculation shows that true exact solution is given by $w\left(x\right)=1+x^3$.

For this test problem, we set $\ell =10$, and utilize $N=3$, which is sufficient to obtain an accurate solution. To do so, we express the solution of (1) in terms of Bessel bases as considered in (4). Afterwards, through solving the fundamental matrix Equation (23), we obtain

$$\begin{array}{c}\hfill C_3={\left[-997/3600-1000/3200/3\right]}^T.\end{array}$$

By exploiting the first four Bessel basis functions and multiply them by $C_3$, we get the approximate solution $w_3\left(x\right)$ on $0\le x\le 10$ as follows

$$\begin{array}{c}\hfill w_3\left(x\right)=\left[1{\mathcal{B}}_{1,\ell }\left(x\right){\mathcal{B}}_{2,\ell }\left(x\right){\mathcal{B}}_{3,\ell }\left(x\right)\right]C_3=1+x^3.\end{array}$$

Clearly this is the exact solution of (1). In the next experiments, the results of absolute errors utilizing diverse values of $\ell =1,5,20$ are computed.

Table 1. The comparison of (absolute) errors in the Bessel matrix approach in Example 1 using $N=3$, $\ell =1,5,20$, and diverse $x\in [0,\ell ]$.

	Bessel ($\mathit{N}=3$)			ANNs ($\mathit{\ell }=1$) [17]
x	$\mathit{\ell }=1$	$\mathit{\ell }=5$	$\mathit{\ell }=20$	Min	Mean	S.D
$0.1\ell$	${9.2006}_{-19}$	${5.4584}_{-16}$	${1.1589}_{-14}$	${6.20}_{-09}$	${1.42}_{-05}$	${2.03}_{-05}$
$0.2\ell$	${3.3134}_{-18}$	${1.9357}_{-15}$	${3.3996}_{-14}$	${1.60}_{-07}$	${7.36}_{-05}$	${4.17}_{-04}$
$0.3\ell$	${6.6298}_{-18}$	${3.7979}_{-15}$	${4.8678}_{-14}$	${5.21}_{-07}$	${1.44}_{-04}$	${8.90}_{-04}$
$0.4\ell$	${1.0319}_{-17}$	${5.7610}_{-15}$	${3.7092}_{-14}$	${2.86}_{-07}$	${2.16}_{-04}$	${1.37}_{-03}$
$0.5\ell$	${1.3831}_{-17}$	${7.4534}_{-15}$	${1.9303}_{-14}$	${1.34}_{-07}$	${2.80}_{-04}$	${1.80}_{-03}$
$0.6\ell$	${1.6615}_{-17}$	${8.5035}_{-15}$	${1.3905}_{-13}$	${1.72}_{-07}$	${3.24}_{-04}$	${2.11}_{-03}$
$0.7\ell$	${1.8121}_{-17}$	${8.5398}_{-15}$	${3.4069}_{-13}$	${1.96}_{-07}$	${3.48}_{-04}$	${2.30}_{-03}$
$0.8\ell$	${1.7799}_{-17}$	${7.1907}_{-15}$	${6.4277}_{-13}$	${6.20}_{-07}$	${3.58}_{-04}$	${2.40}_{-03}$
$0.9\ell$	${1.5099}_{-17}$	${4.0846}_{-15}$	${1.0638}_{-12}$	${9.39}_{-07}$	${3.67}_{-04}$	${2.47}_{-03}$
$1.0\ell$	${9.4700}_{-18}$	${1.1500}_{-15}$	${1.6224}_{-12}$	${5.82}_{-07}$	${3.79}_{-04}$	${2.56}_{-03}$

tabulates the results of ${\mathbb{E}}_N\left(x\right)$ which are calculated at some points $x_j=j\ell /10$ and j varies from 1 to 10. A comparison with the outcomes of artificial neural networks (ANNs) reported in [17] on $[0,1]$ and with 10 neurons is carried out in Table 1 to testify to the validity of our numerical results. An obvious conclusion can be made from Table 1 that our numerical model results are not only highly accurate on the unit interval but also have sufficient accuracy for the larger values of ℓ.

Utilizing $\ell =1,5,20$, the obtained approximate solutions by the Bessel matrix approach are given by

$$\begin{array}{cc}\hfill w_3\left(x\right)& =1.0x^3-1.011761533\times {10}^{-16}{x}^2+1.0,x\in [0,1],\hfill \\ \hfill w_3\left(x\right)& =1.0x^3-2.43108357\times {10}^{-15}{x}^2-5.450188595\times {10}^{-107}x+1.0,x\in [0,5],\hfill \\ \hfill w_3\left(x\right)& =1.0x^3+3.669958207\times {10}^{-15}{x}^2+1.0,x\in [0,20].\hfill \end{array}$$

Finally, we consider $N=4$ and $\ell =50$ for this example. The obtained approximate solution $w_4\left(x\right)$ for $x\in [0,\ell ]$ takes the form

$$\begin{array}{cc}\hfill w_4\left(x\right)& =-7.132116273\times {10}^{-18}{x}^4+1.0x^3-7.257031883\times {10}^{-15}{x}^2\hfill \\ \hfill & -3.966717388\times {10}^{-105}x+1.0.\hfill \end{array}$$

The above approximation together with the exact solutions are visualized in Figure 1. The resulting absolute errors are also presented in this figure, but on the right part.

Example 2.

As the second test problem, let us consider a Lane–Emden differential difference model problem of the form

$$\frac{d^2}{d{x}^2}w(2x-1)+\frac{2}{x}\frac{d}{dx}w\left(3x\right)+xw(x+1)=e^{2x-1}+\frac{2}{x}{e}^{3x}+x{e}^{x+1},0\le x\le \ell ,$$

with initial conditions $w\left(0\right)=1$, $\frac{d}{dx}w\left(0\right)=1$. It is not a difficult task to show that the exact solution of (2) is $w\left(x\right)=e^x$.

Firstly, $N=6$ is used for this example. Considering $\ell =1$, the resulting approximate solution via the Bessel matrix technique takes the following form on $x\in [0,\ell ]$

$$\begin{array}{cc}\hfill & w_6\left(x\right)=0.00464175985x^6-0.001670260859x^5+0.04740032168x^4+0.1803886514x^3\hfill \\ \hfill & +0.4870266818x^2+1.0x+1.0.\hfill \end{array}$$

Let us check that our proposed method yields a reasonable result in comparison with the exact solution. In this respect, the first seven terms of its series expansion are written as

$$e^x\approx 1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\frac{x^5}{5!}+\frac{x^6}{6!}.$$

However, more accurate results based on the Bessel matrix approach will be achieved if one increases N. For instance, the approximated solution $w_{10}\left(x\right)$ can be analogous on $[0,1]$ as

$$\begin{array}{cc}\hfill & w_{10}\left(x\right)=0.000001068044011x^{10}-0.000002908206541x^9+0.00004176809694x^8\hfill \\ \hfill & +0.0001888482682x^7+0.001313363417x^6+0.008508813772x^5\hfill \\ \hfill & +0.04165172897x^4+0.1663795908x^3+0.5001750161x^2+1.0x+1.0.\hfill \end{array}$$

Additionally, the graphical representation of the numerical results using different values of $N=6,10$ and $N=15$ are shown in Figure 2. In this plot, the absolute errors ${\mathbb{E}}_N\left(x\right)$ for $x\in [0,1]$ and for these values of N are also depicted.

In

Table 2. The comparison of absolute errors in the Bessel matrix approach in Example 2 utilizing $N=10,20,30$, $\ell =1,2,5$, and diverse $x\in [0,\ell ].$

	Bessel			ANNs ($\mathit{\ell }=1$) [17]
x	$\mathit{\ell }=1,\mathit{N}=10$	$\mathit{\ell }=2,\mathit{N}=20$	$\mathit{\ell }=5,\mathit{N}=30$	Min	Mean	S.D
$0.1\ell$	${1.4633}_{-6}$	${1.9321}_{-9}$	${1.7737}_{-7}$	${3.30}_{-8}$	${4.51}_{-5}$	${8.65}_{-5}$
$0.2\ell$	${4.7314}_{-6}$	${9.8570}_{-9}$	${2.8830}_{-7}$	${1.09}_{-6}$	${1.81}_{-4}$	${3.67}_{-4}$
$0.3\ell$	${8.2497}_{-6}$	${2.3617}_{-8}$	${9.0792}_{-8}$	${3.79}_{-7}$	${4.22}_{-4}$	${8.68}_{-4}$
$0.4\ell$	${1.0729}_{-5}$	${4.0088}_{-8}$	${1.1621}_{-7}$	${2.04}_{-6}$	${7.52}_{-4}$	${1.59}_{-3}$
$0.5\ell$	${1.1221}_{-5}$	${5.5145}_{-8}$	${3.1451}_{-8}$	${1.05}_{-5}$	${1.16}_{-3}$	${2.50}_{-3}$
$0.6\ell$	${9.1480}_{-6}$	${6.5255}_{-8}$	${3.1422}_{-7}$	${1.51}_{-5}$	${1.58}_{-3}$	${3.50}_{-3}$
$0.7\ell$	${4.2955}_{-6}$	${6.8403}_{-8}$	${6.4486}_{-7}$	${2.80}_{-7}$	${1.96}_{-3}$	${4.47}_{-3}$
$0.8\ell$	${3.2249}_{-6}$	${6.4385}_{-8}$	${6.9379}_{-7}$	${2.02}_{-5}$	${2.26}_{-3}$	${5.30}_{-3}$
$0.9\ell$	${1.3031}_{-5}$	${5.4593}_{-8}$	${3.7913}_{-7}$	${3.87}_{-6}$	${2.43}_{-3}$	${5.85}_{-3}$
$1.0\ell$	${2.4539}_{-5}$	${4.1499}_{-8}$	${1.7989}_{-7}$	${5.39}_{-6}$	${2.45}_{-3}$	${6.05}_{-3}$

, the computed values of the values of the absolute errors ${\mathbb{E}}_N\left(x\right)$ at some points $x\in [0,\ell ]$ with $\ell =1$, $\ell =2$ as well as $\ell =5$ for Example 2 are shown. The corresponding number of basis functions are $N=10,20$ and $N=30$, respectively. Furthermore, the statistical results for this test example utilizing the ANN approach for $\ell =1$ and with 10 neurons that was proposed in [17] are reported in Table 2 for comparison.

In the next and last experiment performed for the test problem (2), we consider a relatively large domain of computation by choosing $\ell =10$. To get a reasonable accuracy, we take $N=60$. In Figure 3, the curve of approximate solution $w_{60}\left(x\right)$ as well as the related relative error ${\mathbb{R}}_N\left(x\right)$ are shown.

Example 3.

The last and third test case is devoted to the following singular differential difference with trigonometric functions defined on $0\le x\le \ell$

$$\frac{d^2}{d{x}^2}w(3x-1)+\frac{2}{x}\frac{d}{dx}w\left(3x\right)+xw(x+1)=-sin(3x-1)+\frac{2}{x}cos\left(3x\right)+xsin(x+1),$$

with initial conditions $w\left(0\right)=0$, $\frac{d}{dx}w\left(0\right)=1$. An easy calculation shows that $w\left(x\right)=sinx$ is the exact solution of this model problem.

Utilizing $N=5,10$ in the Bessel matrix procedure, we get the following polynomial forms for the approximate solutions on $x\in [0,1]$ as follows

$$\begin{array}{cc}\hfill w_5\left(x\right)& =0.004432261297x^5+0.008818459263x^4-0.1704928679x^3\hfill \\ \hfill & -0.002798292238x^2+1.0x-6.38693976\times {10}^{-109},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill w_{10}\left(x\right)& =-1.8339\times {10}^{-7}{x}^{10}+3.1712\times {10}^{-6}{x}^9+1.2818\times {10}^{-7}{x}^8-1.9967\times {10}^{-4}{x}^7\hfill \\ \hfill & -6.0962\times {10}^{-7}{x}^6+0.008337076448x^5+2.0318\times {10}^{-7}{x}^4-0.1666714768x^3\hfill \\ \hfill & +1.7925\times {10}^{-6}{x}^2+1.0x-9.337157052\times {10}^{-109}.\hfill \end{array}$$

Let us consider the series form of the exact solution, i.e., $sinx\approx x-\frac{x^3}{3!}+\frac{x^5}{5!}-\dots +\frac{x^9}{9!}$. A comparison between the achieved approximations and the exact one indicates the good alignment between them. The absolute errors ${\mathbb{E}}_N\left(x\right)$ utilizing diverse values of $N=5,10,15$ in the approximate solutions are presented in Figure 4. The curves of $w_N\left(x\right)$ for $N=5,10$ are also visualized in Figure 4.

We next examine the benefits of the presented Bessel matrix scheme and validate our results for $\ell =1,\pi ,2\pi$. We also utilize various $N=10,20$, and $N=30$ for each computational domain $[0,\ell ]$. In this respect and in terms of accuracy, some comparisons are performed in Example 3 in

Table 3. The comparison of absolute errors in the Bessel matrix method in Example 3 utilizing $N=10,20,30$, $\ell =1,\pi ,2\pi$, and diverse $x\in [0,\ell ].$

	Bessel			ANNs ($\mathit{\ell }=1$) [17]
x	$\mathit{\ell }=1,\mathit{N}=10$	$\mathit{\ell }=\mathit{\pi },\mathit{N}=20$	$\mathit{\ell }=2\mathit{\pi },\mathit{N}=30$	Min	Mean	S.D
$0.1\ell$	${1.3172}_{-8}$	${4.9934}_{-12}$	${5.9707}_{-11}$	${3.51}_{-8}$	${1.03}_{-6}$	${1.29}_{-6}$
$0.2\ell$	${3.4686}_{-8}$	${1.3423}_{-11}$	${7.7946}_{-11}$	${7.83}_{-9}$	${2.37}_{-6}$	${2.34}_{-6}$
$0.3\ell$	${4.1486}_{-8}$	${1.7895}_{-11}$	${8.3353}_{-11}$	${8.63}_{-8}$	${3.90}_{-6}$	${3.94}_{-6}$
$0.4\ell$	${1.8096}_{-8}$	${2.0670}_{-11}$	${7.7780}_{-11}$	${6.81}_{-9}$	${5.52}_{-6}$	${5.42}_{-6}$
$0.5\ell$	${4.1676}_{-8}$	${2.2442}_{-11}$	${6.2475}_{-11}$	${8.48}_{-8}$	${6.83}_{-6}$	${6.74}_{-6}$
$0.6\ell$	${1.3460}_{-7}$	${2.3120}_{-11}$	${3.7118}_{-11}$	${6.03}_{-8}$	${7.52}_{-6}$	${7.82}_{-6}$
$0.7\ell$	${2.4942}_{-7}$	${2.3020}_{-11}$	${1.3218}_{-12}$	${1.37}_{-7}$	${7.74}_{-6}$	${8.37}_{-6}$
$0.8\ell$	${3.6947}_{-7}$	${2.2240}_{-11}$	${4.5159}_{-11}$	${1.28}_{-7}$	${7.61}_{-6}$	${8.51}_{-6}$
$0.9\ell$	${4.7614}_{-7}$	${2.0646}_{-11}$	${1.0211}_{-10}$	${8.94}_{-9}$	${7.53}_{-6}$	${8.38}_{-6}$
$1.0\ell$	${5.5254}_{-7}$	${1.8200}_{-11}$	${1.6874}_{-10}$	${2.68}_{-8}$	${7.70}_{-6}$	${8.27}_{-6}$

. Besides the numerical results reported at some x in $[0,\ell ]$, the outcomes of the achieved errors obtained via ANNs with 10 neurons are further presented in Table 3 as has been done in the previous test examples. Comparing our numerical achievements in Table 1, Table 2 and Table 3 with the outcomes of ANNs reveals that our approach is more accurate and can be easily applied on long computational domains.

Finally, we take $\ell =5\pi$. Moreover, the value of $N=60$ is utilized as the number of the basis of functions. The numerical solution $w_{60}\left(x\right)$ is depicted in Figure 5. The graphical representation of the absolute error ${\mathbb{E}}_N\left(x\right)$ is also shown in Figure 5. It can be obviously observed that, to keep the accuracy of the proposed method up on a long domain of computation, one has to adjust the number of basis functions accordingly.

4. Conclusions

In this manuscript, an accurate and reliable numerical technique is designed and implemented to obtain an approximate solution for a class of singular second-order functional differential models using the Bessel polynomials. The proposed approach is viably developed for solving various test examples of the second-order singular functional differential model problems. A precise and accurate performance is witnessed for this technique consistently achieving high accuracy from the existing exact results for the given test cases based upon the second-order singular functional differential model. The presented approach looks proficient and promising for solving similar applicable problems in the future.