#  Vieta-Lucas Wavelet based schemes for the numerical solution of the
singular models

Rakesh Kumar  School of Mathematics, Computers and Information Sciences,
Central University of Himachal Pradesh, Shahpur Campus, Shahpur 176206, India
Shivani Aeri e-mail<EMAIL_ADDRESS>School of Mathematics, Computers
and Information Sciences, Central University of Himachal Pradesh, Shahpur
Campus, Shahpur 176206, India  Dumitru Baleanu  Department of Mathematics,
Cankaya University, 06810 Ankara, Turkey  Institute of Space Sciences,
Magurele-Bucharest, Romania  Kottakkaran Soopy Nisar  Department of
Mathematics, College of Arts Sciences, Prince Sattam Bin Abdulaziz University,
Saudi Arabia

###### Abstract

In this paper, numerical methods based on Vieta-Lucas wavelets are proposed
for solving a class of singular differential equations. The operational matrix
of the derivative for Vieta-Lucas wavelets is derived. It is employed to
reduce the differential equations into the system of algebraic equations by
applying the ideas of the collocation scheme, Tau scheme, and Galerkin scheme
respectively. Furthermore, the convergence analysis and error estimates for
Vieta-Lucas wavelets are performed. In the numerical section, the comparative
analysis is presented among the different versions of the proposed Vieta-Lucas
wavelet methods, and the accuracy of the approaches is evaluated by computing
the errors and comparing them to the existing findings.

Keywords: Vieta-Lucas wavelets, generating function, Rodrigues’ formula,
collocation method, Galerkin method, Tau method.

##  1 Introduction

Singular differential equations (SDEs) have been attracting applied
mathematicians for many years because of their applicability in various
branches of science, engineering, and technology [1, 2]. We consider the
following SDEs [3, 4]:

$\rm{Y^{\prime\prime}(\zeta)+\frac{\mu}{\zeta}Y^{\prime}(\zeta)+f(\zeta,Y(\zeta))=g(\zeta),\leavevmode\nobreak\
\leavevmode\nobreak\ \zeta\in(0,L),}$ (1)

with initial conditions

$\rm{Y(\zeta)\mid_{\zeta=0}=\alpha_{0},\>Y^{\prime}(\zeta)\mid_{\zeta=0}=\alpha_{1}},$
(2)

or boundary conditions

$\rm{Y(\zeta)\mid_{\zeta=0}=\beta_{0},\>Y(\zeta)\mid_{\zeta=L}=\beta_{1}},$
(3)

where $\rm{g(\zeta)}$ is an analytical function, $\rm{f(\zeta,Y(\zeta))}$ is a
continuous real valued function, and $\alpha_{0}$, $\alpha_{1}$, $\beta_{0}$,
$\beta_{1}$ are arbitrary constants. Several researchers have been interested
in singular models since the available singularity in these models makes
computational efforts more tedious. The few well-known singular models which
have wider applications in thermodynamics, astrophysics, and atomic physics,
are named as Emden-Fowler model, Lane-Emden model, and Thomas-Fermi model [5].
The Runge-Kutta technique, Euler’s method, Adams method, Milne predictor-
corrector approach, and other traditional methods fail to analyze the singular
models efficiently [6]. SDEs generally appear with variable coefficients, and
the solution to these equations can be obtained by the power series method,
continued fractions, and integral transforms. The solution to these equations
leads to the generation of several special functions [7]. The existence and
uniqueness of the solutions for the SDEs have been investigated by Ford and
Pennline [8]. The singularities occurring in these equations have encouraged
numerous researchers to explore the solutions by proposing novel numerical
schemes. To address singular nonlinear problems, Bender et al. [9] presented a
perturbation approach. Russell and Shampine [10] used patch bases, collocation
method, and finite difference method for the treatment of SBVPs. Wazwaz [11,
12] respectively employed Adomian decomposition to solve SIVPs and variational
iteration method to solve the nonlinear SBVPs. By incorporating a quasi-
linearization approach, El-Gebeily and O’Regan[13] solved second-order
nonlinear SDEs. Yildirim and Ozis [14] employed the homotopy perturbation
method to solve SIVPs. Pandey [15] used the finite difference method to solve
a class of SBVPs. Sabir et al. used artificial neural networking based
algorithms to solve various SDEs [16, 17, 18]. Researchers have also utilized
Chebyshev polynomials [19], Legendre polynomials [20], Laguerre polynomials
[21, 22], Jacobi polynomials [23] and Hermite polynomials [24] to derive the
solutions of SDEs. However, less emphasis has been placed on the Vieta-Lucas
polynomials (VLPs). Recently, VLPs based schemes are used by Agarwal [25, 26,
27, 28] and El-Sayed to solve fractional advection-dispersion equation and
Heydari et. al [29] to solve variable order fractional Ginzburg-landau
equations.
In recent decades, wavelets have piqued the interest of many academicians due
to their ability in solving differential equations with high precision and
minimal processing effort [30, 31]. Recently, Chebyshev wavelets [32, 33],
Legendre wavelets[34], Hermite wavelets [35], Jacobi wavelets [36], Gegenbauer
wavelets [37] and Lucas wavelets[38, 39] have been successfully utilized in
literature for the simulation of various phenomena of scientific and technical
importance. Vieta-Lucas wavelets have not been used in past to solve singular
differential equations yet. Recently, some work is reported in the literature
on Vieta-Lucas polynomials, but not much attention has been given to Vieta-
Lucas wavelets. The wavelets are localized in nature and have compact support
so to fill the literature gap, we choose to construct the Vieta-Lucas
wavelets.
The main concern of the paper is to propose a novel class of wavelets called
as Vieta-Lucas wavelets that can handle various differential equations
efficiently. The novelty of the work includes the derivation of generating
function and Rodrigues’ formula for VLPs. The shifted form of VLPs is used to
prepare the Vieta-Lucas wavelets. The operational matrix (OM) of derivatives
is proposed to formulate the numerical scheme. The Emden-Fowler and Lane-Emden
type SDEs are also solved by the proposed approaches and perform well near the
singularity. The accuracy of the methods is analyzed by computing the errors
in $L^{2}$ and $L^{\infty}$ norms.
The remaining part is organized as: A brief overview of VLPs is given in
section 2 by including generating function, Rodrigues’ formula, and other
significant properties. In order to solve the problems over $[0,L]$, shifted
version of VLPs is given in section 3. In section 4, Vieta-Lucas wavelets and
their functional approximations are presented, and their operational matrix of
the derivative is derived in section 5. In section 6 Vieta-Lucas wavelets-
based numerical schemes are provided. Section 7 deals with the convergence and
error estimation for Vieta-Lucas wavelets. In section 8, numerous
illustrations are solved by the proposed approaches, and results are provided
in section 9 that validate the efficiency and reliability of the proposed
methods. Finally, the conclusion of the proposed work is demonstrated in
section 10.

##  2 Vieta-Lucas polynomials and properties

In this section, we give brief introduction about Vieta-Lucas polynomials,
recurrence relation, orthogonality property, generating function, Vieta-Lucas
differential equation and Rodrigues’ formula.

######  Definition 2.1.

The Vieta-Lucas polynomials $\rm{VL_{m}(\zeta)}$ of degree $\rm{m}$
($\rm{m\in\mathbb{N}\cup\\{0\\}}$) can be defined as [40, 41]:

$\rm{VL_{m}(\zeta)=2\cos(m\delta)},$ (4)

where $\delta=\arccos{(\frac{\zeta}{2})}$ and $|\zeta|\in[-2,2]$,
$\delta\in[0,\pi].$

######  Proposition 2.1.

The recurrence relation for Vieta-Lucas polynomials $\rm{VL_{m}}(\zeta)$ is
given by [40]:

$\rm{VL_{m}(\zeta)=\zeta VL_{m-1}(\zeta)-VL_{m-2}(\zeta)},\leavevmode\nobreak\
\leavevmode\nobreak\ m\geq 2,$ (5)

with $\rm{VL_{0}(\zeta)=2}$ and $\rm{VL_{1}(\zeta)=\zeta}$.

######  Proposition 2.2.

The Vieta-Lucas polynomials can be represented in terms of power series
expansion as [40]:

$\rm{VL_{m}(\zeta)=\sum_{i=0}^{\lfloor
m/2\rfloor}(-1)^{i}\frac{m(m-i-1)!\>}{i!(m-2i)!}\zeta^{m-2i}},\leavevmode\nobreak\
\leavevmode\nobreak\ m\geq 1.$ (6)

The first few Vieta-Lucas polynomials are given as:

$\begin{split}\rm{VL_{0}(\zeta)}&=2,\\\ \rm{VL_{1}(\zeta)}&=\zeta,\\\
\rm{VL_{2}(\zeta)}&=\zeta^{2}-2,\\\ \rm{VL_{3}(\zeta)}&=\zeta^{3}-3\zeta,\\\
\rm{VL_{4}(\zeta)}&=\zeta^{4}-4\zeta^{2}+2,\\\
\rm{VL_{5}(\zeta)}&=\zeta^{5}-5\zeta^{3}+5\zeta,\\\
\rm{VL_{6}(\zeta)}&=\zeta^{6}-6\zeta^{4}+9\zeta^{2}-2.\end{split}$

######  Proposition 2.3.

The Vieta-Lucas polynomials $\rm{VL_{n}(\zeta)}$ and $\rm{VL_{m}(\zeta)}$
defined over $[-2,2]$ are orthogonal in weighted sense with the weight
function $\rm{w(\zeta)=\frac{1}{\sqrt{4-\zeta^{2}}}}$ and satisfy the
following condition [25]:

$\rm{\langle
VL_{n}(\zeta),VL_{m}(\zeta)\rangle_{w(\zeta)}=\int_{-2}^{2}VL_{n}(\zeta)VL_{m}(\zeta)w(\zeta)\,d\zeta=\begin{cases}4\pi,&n=m=0,\\\
2\pi,&n=m\neq 0,\\\ 0,&n\neq m.\end{cases}}$ (7)

######  Proposition 2.4.

The generating function for Vieta-Lucas polynomials is defined as:

$\rm{\sum_{m=0}^{\infty}VL_{m}(\zeta)t^{m}=\frac{2-\zeta t}{1-\zeta
t+t^{2}}}.$ (8)

###### Proof.

Since,

$\rm{\sum_{m=0}^{\infty}\rm{VL}_{m}(\zeta)t^{m}=\rm{VL}_{0}(\zeta)t^{0}+\rm{VL}_{1}(\zeta)t^{1}+\sum_{m=2}^{\infty}\rm{VL}_{m}(\zeta)t^{m}}.$

Therefore,

$\begin{split}\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m}}&=\rm{2+\zeta
t}+\sum_{m=2}^{\infty}\rm{[\zeta\rm{VL}_{m-1}(\zeta)-\rm{VL}_{m-2}(\zeta)]t^{m}},\\\
\implies\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m}}&=\rm{2+\zeta
t+\zeta}\sum_{m=1}^{\infty}\rm{VL_{m}(\zeta)t^{m+1}}-\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m+2}},\\\
\implies\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m}}&=\rm{2+\zeta t+\zeta
t}[\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m}-\rm{VL}_{0}(\zeta)]-t^{2}}\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m}},\\\
\implies\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m}}&=\rm{2-\zeta t+(\zeta
t-t^{2})}\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m}},\\\
\implies\sum_{m=0}^{\infty}\rm{VL_{m}(\zeta)t^{m}}&=\rm{\frac{2-\zeta
t}{1-\zeta t+t^{2}}}.\end{split}$

∎

######  Proposition 2.5.

The Vieta-Lucas differential equation is defined as [42]:

$\rm{(4-\zeta^{2})\diff[2]Y\zeta-\zeta\diff
Y\zeta+m^{2}Y=0},\leavevmode\nobreak\ \leavevmode\nobreak\ m\in\mathbb{N}.$
(9)

######  Theorem 2.1.

The Rodrigues’ formula for Vieta-Lucas polynomials $\rm{VL}_{m}(\zeta)$ can be
obtained as:

$\rm{VL_{m}(\zeta)=(-1)^{m}\>2\,\frac{m!}{(2m)!}(4-{\zeta}^{2})^{\frac{1}{2}}\frac{d^{m}}{d\zeta^{m}}\Big{\\{}(4-{\zeta}^{2})^{m-\frac{1}{2}}\Big{\\}}}.$
(10)

###### Proof.

Consider

$\rm{f(\zeta)=(4-{\zeta}^{2})^{m-\frac{1}{2}}}.$

On differentiating, we obtain

$\rm{f^{\prime}(\zeta)=-(2m-1)\zeta(4-{\zeta}^{2})^{m-\frac{3}{2}}},$
$\Rightarrow\rm{(4-{\zeta}^{2})f^{\prime}(\zeta)+(2m-1)\zeta f(\zeta)=0}.$

On differentiating it (m+1) times, we have

$\rm{(4-{\zeta}^{2})\>D^{m+2}f(\zeta)-3\,\zeta\>D^{m+1}f(\zeta)+(m^{2}-1)\>D^{m}f(\zeta)=0,}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ where\leavevmode\nobreak\ D\equiv\frac{d}{d\zeta}.$ (11)

Let $\rm{g(\zeta)=(4-{\zeta}^{2})^{\frac{1}{2}}\>D^{m}f(\zeta)\>}$, then (11)
reduces to

$\rm{(4-\zeta^{2})\,g^{\prime\prime}(\zeta)-\zeta\,g^{\prime}(\zeta)+m^{2}g(\zeta)=0},$

which is equivalent to Vieta-Lucas differential equation. Thus, we can choose
$\rm{VL}_{m}(\zeta)=C_{m}\,g(\zeta)$, where $\rm C_{m}$ is a constant to be
determined. To find $\rm C_{m}$, the coefficients of $\rm{\zeta^{m}}$ in
$\rm{VL}_{m}(\zeta)$ and $\rm{g(\zeta)}$ are required to be compared. The
coefficient of $\rm{\zeta^{m}}$ in $\rm{VL}_{m}(\zeta)$ is 1.
Since,

$\begin{split}\rm{g(\zeta)}&=\rm{(4-{\zeta}^{2})^{\frac{1}{2}}\>D^{m}\,(4-{\zeta}^{2})^{m-\frac{1}{2}}},\\\
&=\rm{(4-{\zeta}^{2})^{\frac{1}{2}}\>\sum_{j=0}^{m}\binom{m}{j}\>D^{m-j}\,(2+\zeta)^{m-\frac{1}{2}}\>D^{j}\,(2-\zeta)^{m-\frac{1}{2}}},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{ (by Leibniz's rule)}\\\
&=\rm{(-1)^{m}\>m!\>\sum_{j=0}^{m}\binom{m-\frac{1}{2}}{j}\binom{m-\frac{1}{2}}{m-j}(\zeta-2)^{m-j}\>(\zeta+2)^{j}}.\end{split}$

Therefore, the coefficient of $\rm{\zeta^{m}}$ in $\rm{g(\zeta)}$ is

$\rm{(-1)^{m}\>m!\>\sum_{j=0}^{m}\binom{m-\frac{1}{2}}{j}\binom{m-\frac{1}{2}}{m-j}}.$

The Vandermonde’s convolution modifies the above coefficient of
$\rm{\zeta^{m}}$ to

$\rm{\frac{(-1)^{m}}{2}\frac{(2m)!}{m!}}.$

Thus, $\rm{C_{m}=(-1)^{m}\>2\;\frac{m!}{(2m)!}}$. Hence the Rodrigues’
formula. ∎

######  Remark 2.1.

The zeroes and extreme points of $\rm{VL}_{m}(\zeta)$ in $[-2,2]$ are
respectively presented as
$\zeta_{j}=2\cos{\left((j-\frac{1}{2})\frac{\pi}{m}\right)}$ ,
$\zeta_{j}=2\cos{\left(j\frac{\pi}{m}\right)}$, where j=1,2,3,…, m [40].

######  Proposition 2.6.

The function $\zeta^{m}$ can be expressed in terms of $\rm{VL}_{m}(\zeta)$ as:

$\rm{\zeta^{m}}=\sum_{j=0}^{\lfloor m/2\rfloor^{\star}}{m\choose
j}\rm{VL}_{m-2j}(\zeta),$ (12)

where ’${\star}$’ denotes that the last term in the sum is to be halved when m
is even.

###### Proof.

Since,

$\zeta^{m}=(2\cos{\delta})^{m}=(e^{i\delta}+e^{-i\delta})^{m}.$ (13)

Therefore, the use of binomial expansion gives

$\begin{split}\zeta^{m}&=e^{im\delta}+{m\choose
1}e^{i(m-2)\delta}+\dots+{m\choose m-1}e^{-i(m-2)\delta}+e^{-im\delta},\\\
&=(e^{im\delta}+e^{-im\delta})+{m\choose
1}(e^{i(m-2)\delta}+e^{-i(m-2)\delta})+\dots,\\\ &=2\cos{m\delta}+{m\choose
1}2\cos{(m-2)\delta}+{m\choose 2}2\cos{(m-4)\delta}+\dots,\\\
&=\rm{VL}_{m}(\zeta)+{m\choose 1}\rm{VL}_{m-2}(\zeta)+{m\choose
2}\rm{VL}_{m-4}(\zeta)+\dots.\end{split}$

Hence the result. ∎

##  3 Shifted Vieta-Lucas polynomials

######  Definition 3.1.

The shifted Vieta-Lucas polynomials defined over $[0,2]$ are denoted by
$\rm{VL^{*}_{m}(\zeta)}$ of degree $\rm{m\in\mathbb{N}\cup\\{0\\}}$ as [29]:

$\rm{VL^{*}_{m}(\zeta)=VL_{m}(2\zeta-2)}.$ (14)

The recurrence relation for shifted Vieta-Lucas polynomials
$\rm{VL^{*}_{m}(\zeta)}$ is [29]:

$\rm{VL^{*}_{m}(\zeta)=(2\zeta-2)VL^{*}_{m-1}(\zeta)-VL^{*}_{m-2}(\zeta)},$
(15)

provided $\rm{VL^{*}_{0}(\zeta)=2}$ and $\rm{VL^{*}_{1}(\zeta)=2\zeta-2}.$
The shifted Vieta-Lucas polynomials satisfy the following orthogonality
property [29]:

$\rm{\langle
VL^{*}_{n}(\zeta),VL^{*}_{m}(\zeta)\rangle_{w^{*}(\zeta)}=\int_{0}^{2}VL^{*}_{n}(\zeta)VL^{*}_{m}(\zeta)w^{*}(\zeta)\,d\zeta=\begin{cases}4\pi,&n=m=0,\\\
2\pi,&n=m\neq 0,\\\ 0,&n\neq m,\end{cases}}$ (16)

where $\rm{w^{*}(\zeta)=\frac{1}{\sqrt{2\zeta-\zeta^{2}}}}$ is the weight
function of shifted Vieta-Lucas polynomials.

##  4 Vieta-Lucas wavelets and function approximation

Wavelet is a type of function that is derived through the dilation and
translation of the mother wavelet. The continuous wavelet family with dilation
parameter $h$ and translation parameter $r$ is defined as [43]:

$\gamma_{h,r}(\zeta)=|h|^{-1/2}\gamma(\frac{\zeta-r}{h}),\hskip
14.22636pth,r\in\mathbb{R},h\neq 0.$ (17)

If $h=h_{0}^{-k},r=sr_{0}h_{0}^{-k},h_{0}>1,r_{0}>0$, then the discrete
wavelet family consists of the following members:

$\gamma_{k,s}(\zeta)=|h_{0}|^{k/2}\gamma(h_{0}\zeta-sr_{0}),\hskip
14.22636ptk,s\in\mathbb{Z}^{+},$ (18)

where $\gamma_{k,s}(x)$ constitutes the wavelet basis for $L_{2}(\mathbb{R})$.
For a specific choice of $h_{0}=2$ and $r_{0}=1$, $\gamma_{k,s}(x)$
constitutes an orthonormal basis.

######  Definition 4.1.

The Vieta-Lucas wavelets $\rm{\Upsilon_{s,m}(\zeta)=\Upsilon(k,s,m,\zeta)}$ is
defined on the interval [0,2) as:

$\rm{\Upsilon_{s,m}}(\zeta)=\begin{cases}\rm{2^{\frac{k}{2}}\>\widehat{\rm{VL}}_{m}(\>2^{k}\zeta-\hat{s}\>),}&\rm{\frac{\hat{s}-2}{2^{k}}\leq\zeta<\frac{\hat{s}+2}{2^{k}}},\\\
0,&\text{Otherwise},\end{cases}$ (19)

where

$\widehat{\rm{VL}}_{m}(\zeta)=\begin{cases}\frac{1}{\sqrt{\pi}},&m=0,\\\
\frac{1}{\sqrt{2\pi}}\rm{VL_{m}}(\zeta),&m\geq 1,\end{cases}$ (20)

$m=0,1,2,\dots,M-1;$ M be the maximum order of Vieta-Lucas polynomials,
$s=1,2,\dots,2^{k-1}$; $k=1,2,3,\dots;$ $\hat{s}=2(2s-1)$.

######  Remark 4.1.

Vieta-Lucas wavelets form an orthogonal set with respect to the weight
functions $\rm{w_{s}(\zeta)=w(2^{k}\zeta-\hat{s})}$.

######  Definition 4.2.

A function $\rm{Y(\zeta)}$ defined over ${L^{2}}_{w_{s}}[0,2]$ can be written
in terms of Vieta-Lucas wavelets series as:

$\rm{Y(\zeta)=\sum_{s=1}^{\infty}\sum_{m=0}^{\infty}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)},$
(21)

with

$\rm{\Lambda_{s,m}=\langle
Y(\zeta),\Upsilon_{s,m}(\zeta)\rangle_{w_{s}(\zeta)}=\int_{0}^{2}Y(\zeta)\Upsilon_{s,m}(\zeta)w_{s}(\zeta)\,d\zeta},$
(22)

where $\langle*,*\rangle$ denotes the inner product.

The truncated form of Vieta Lucas wavelet series can be written as:

$\rm{\bar{Y}(\zeta)}\cong\rm{\sum_{s=1}^{2^{k-1}}\sum_{m=0}^{M-1}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)=\rm{\Lambda^{T}\Upsilon(\zeta)}},$
(23)

where $\Lambda$ and $\Upsilon(\zeta)$ are $\eta\times 1$ matrices in the
following form for $\eta=2^{k-1}M$, and

$\displaystyle\rm{\Lambda}$
$\displaystyle=\rm{[\Lambda_{1,0},\Lambda_{1,1},\dots,\Lambda_{1,M-1},\Lambda_{2,0},\dots,\Lambda_{2,M-1},\dots,\Lambda_{2^{k-1},0},\dots,\Lambda_{2^{k-1},M-1}]^{T}},$
(24) $\displaystyle\rm{\Upsilon(\zeta)}$
$\displaystyle=\rm{[\Upsilon_{1,0}(\zeta),\dots,\Upsilon_{1,M-1}(\zeta),\Upsilon_{2,0}(\zeta),\dots,\Upsilon_{2,M-1}(\zeta),\dots,\Upsilon_{2^{k-1},0}(\zeta),\dots,\Upsilon_{2^{k-1},M-1}(\zeta)]^{T}}.$
(25)

##  5 Vieta-Lucas wavelet based operational matrix

######  Theorem 5.1.

Let $\Upsilon(\zeta)$ be the Vieta-Lucas wavelets vector defined in (25). Then
the derivative of the vector $\Upsilon(\zeta)$ can be expressed as:

$\frac{\rm{d\Upsilon(\zeta)}}{\rm{d\zeta}}=\rm{D\Upsilon(\zeta)},$ (26)

where $D$ is the $\eta\times\eta$ matrix given by

$\rm{D}=\begin{pmatrix}\rm{F}&0&\dots&0\\\ 0&\rm{F}&\dots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&0&\rm{F}\end{pmatrix}$

in which $F$ is $M\times M$ square matrix whose $(u,v)^{th}$ element is
defined as:

$\rm{F_{u,v}}=\begin{cases}\rm{\frac{2^{k}(u-1)}{\sqrt{\alpha_{u-1}}\sqrt{\alpha_{v-1}}}},&\rm{u=2,\dots,M;}\hskip
2.84544pt\rm{v=1,2,\dots,u-1;}\hskip 2.84544pt\rm{(u+v)}=\text{odd},\\\
0,&\text{Otherwise}.\end{cases}$ (27)

###### Proof.

Using shifted Vieta-Lucas polynomials vector, $\Upsilon(\zeta)$ can be
rewritten as

$\rm{\Upsilon_{u}(\zeta)=\Upsilon_{s,m}(\zeta)=2^{\frac{k}{2}}\sqrt{\frac{1}{2\pi\alpha_{m}}}VL^{*}_{m}(2^{k-1}\zeta-2s+2)\chi_{[\frac{s-1}{2^{k-2}},\frac{s}{2^{k-2}}]}},$
(28)

where $\rm{u=(s-1)M+(m+1)}$; $\rm{s=1,2,3,....2^{k-1}}$,
$\rm{m=0,1,2,....M-1}$ and

$\rm{\chi_{[\frac{s-1}{2^{k-2}},\frac{s}{2^{k-2}}]}}=\begin{cases}1,&\rm{\zeta\in[\frac{s-1}{2^{k-2}},\frac{s}{2^{k-2}}]},\\\
0,&\text{Otherwise}.\end{cases}$ (29)

On differentiating (28) with respect to $\rm{\zeta}$, we obtain

$\rm{\frac{d\Upsilon_{u}(\zeta)}{d\zeta}=2^{\frac{3k}{2}-1}\sqrt{\frac{1}{2\pi\alpha_{m}}}\frac{d}{d\zeta}VL^{*}_{m}(2^{k-1}\zeta-2s+2)\chi_{[\frac{s-1}{2^{k-2}},\frac{s}{2^{k-2}}]}},$
(30)

since (30) vanishes outside the interval
$\rm{\zeta\in[\frac{s-1}{2^{k-2}},\frac{s}{2^{k-2}}]}$. Thus, the nonzero
components of Vieta-Lucas wavelets expansion exist only in the interval
$\rm{\zeta\in[\frac{s-1}{2^{k-2}},\frac{s}{2^{k-2}}]}$ i.e,
$\rm{\Upsilon_{i}(\zeta)}$ for $i=(s-1)M+1,(s-1)M+2,\dots,sM.$ Vieta-Lucas
wavelets expansion can now be written as

$\rm{\frac{d\Upsilon_{u}(\zeta)}{d\zeta}=\sum_{i=(s-1)M+1}^{sM}a_{i}\Upsilon_{i}(\zeta)}.$

Here,

$\rm{\frac{d}{d\zeta}\Upsilon_{u}(\zeta)=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{for}\hskip
5.69046ptu=1,M+1,2M+1,\dots,(2^{k-1}-1)M+1},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{because}\leavevmode\nobreak\
\rm{\frac{d}{d\zeta}VL_{0}(\zeta)=0}\leavevmode\nobreak\ \text{everywhere.}$

So, the first row of matrix $F$ is zero.
Now, the first derivative of shifted Vieta-Lucas polynomials can be expressed
as

$\rm{\frac{d}{d\zeta}VL^{*}_{m}(\zeta)=2\sum_{\begin{subarray}{c}j=0\\\
j+m=odd\end{subarray}}^{M-1}\frac{m}{\alpha_{j}}VL^{*}_{j}(\zeta)}.$ (31)

Using (31) in (30), we obtain

$\rm{\frac{d\Upsilon_{u}(\zeta)}{d\zeta}}=\rm{2^{\frac{3k}{2}}\sqrt{\frac{1}{2\pi\alpha_{m}}}\sum_{\begin{subarray}{c}j=0\\\
j+m=odd\end{subarray}}^{M-1}\rm{\frac{m}{\alpha_{j}}VL^{*}_{m}(2^{k-1}\zeta-2s+2)\chi_{[\frac{s-1}{2^{k-2}},\frac{s}{2^{k-2}}]}},}$
(32)

which can be rewritten as

$\rm{\frac{d\Upsilon_{u}(\zeta)}{d\zeta}=2^{k}\sum_{\begin{subarray}{c}v=1\\\
u+v=odd\end{subarray}}^{u-1}\frac{u-1}{\sqrt{\alpha_{u-1}}\sqrt{\alpha_{v-1}}}\Upsilon_{(s-1)M+v}(\zeta).}$

Therefore

$\rm{\frac{d\Upsilon_{u}(\zeta)}{d\zeta}=F_{u,v}\Upsilon_{(s-1)M+v}(\zeta),}$

with

$\rm{F_{u,v}}=\begin{cases}\rm{2^{k}\frac{u-1}{\sqrt{\alpha_{u-1}}\sqrt{\alpha_{v-1}}}},&\rm{u=2,\dots,M;\hskip
2.84544ptv=1,2,\dots,u-1;\hskip 2.84544pt(u+v)=\text{odd}},\\\
0,&\text{Otherwise}.\end{cases}$ (33)

which leads to the desired expression. ∎

For example, if we select k=2 and M=3, then the discrete members of shifted
Vieta-Lucas wavelets can be written as:

$\displaystyle\rm{\Upsilon_{1}(\zeta)=\Upsilon_{1,0}(\zeta)}=$
$\displaystyle\begin{cases}\rm{\frac{2}{\sqrt{\pi}}},&0\leq\zeta<1,\\\
0,&\text{Otherwise}.\end{cases}$
$\displaystyle\rm{\Upsilon_{2}(\zeta)=\Upsilon_{1,1}(\zeta)}=$
$\displaystyle\begin{cases}\rm{\frac{2\sqrt{2}}{\sqrt{\pi}}(2\zeta-1)},&0\leq\zeta<1,\\\
0,&\text{Otherwise}.\end{cases}$
$\displaystyle\rm{\Upsilon_{3}(\zeta)=\Upsilon_{1,2}(\zeta)}=$
$\displaystyle\begin{cases}\rm{\frac{2\sqrt{2}}{\sqrt{\pi}}(8\zeta^{2}-8\zeta+1)},&0\leq\zeta<1,\\\
0,&\text{Otherwise}.\end{cases}$
$\displaystyle\rm{\Upsilon_{4}(\zeta)=\Upsilon_{2,0}(\zeta)}=$
$\displaystyle\begin{cases}\rm{\frac{2}{\sqrt{\pi}}},&1\leq\zeta<2,\\\
0,&\text{Otherwise}.\end{cases}$
$\displaystyle\rm{\Upsilon_{5}(\zeta)=\Upsilon_{2,1}(\zeta)}=$
$\displaystyle\begin{cases}\rm{\frac{2\sqrt{2}}{\sqrt{\pi}}(2\zeta-3),}&1\leq\zeta<2,\\\
0,&\text{Otherwise}.\end{cases}$
$\displaystyle\rm{\Upsilon_{6}(\zeta)=\Upsilon_{2,2}(\zeta)}=$
$\displaystyle\begin{cases}\rm{\frac{2\sqrt{2}}{\sqrt{\pi}}(8\zeta^{2}-24\zeta+17)},&1\leq\zeta<2,\\\
0,&\text{Otherwise}.\end{cases}$

As a result, the first order derivatives of the shifted Vieta-Lucas wavelets
over the domain $[0,2]$ are:

$\displaystyle\rm{\frac{d\Upsilon_{1}}{d\zeta}}$ $\displaystyle=0,$
$\displaystyle\rm{\frac{d\Upsilon_{2}}{d\zeta}}$
$\displaystyle=\frac{4\sqrt{2}}{\sqrt{\pi}}=2\sqrt{2}\leavevmode\nobreak\
\Upsilon_{1},$ $\displaystyle\rm{\frac{d\Upsilon_{3}}{d\zeta}}$
$\displaystyle=\frac{16\sqrt{2}}{\sqrt{\pi}}(2\zeta-1)=8\leavevmode\nobreak\
\Upsilon_{2},$ $\displaystyle\rm{\frac{d\Upsilon_{4}}{d\zeta}}$
$\displaystyle=0,$ $\displaystyle\rm{\frac{d\Upsilon_{5}}{d\zeta}}$
$\displaystyle=\frac{4\sqrt{2}}{\sqrt{\pi}}=2\sqrt{2}\leavevmode\nobreak\
\Upsilon_{4},$ $\displaystyle\rm{\frac{d\Upsilon_{6}}{d\zeta}}$
$\displaystyle=\frac{16\sqrt{2}}{\sqrt{\pi}}(2\zeta-3)=8\leavevmode\nobreak\
\Upsilon_{5}.$

So, the matrix D is as follows:

$\rm{D}=\begin{pmatrix}0&0&0&0&0&0\\\ 2\sqrt{2}&0&0&0&0&0\\\ 0&8&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&0&0&2\sqrt{2}&0&0\\\ 0&0&0&0&8&0\\\ \end{pmatrix}$.

######  Corollary 5.1.

The $m^{th}$ order differential OM can be achieved as:

$\rm{\frac{d^{(m)}\Upsilon(\zeta)}{d\zeta^{(m)}}=D^{(m)}\Upsilon(\zeta)},$
(34)

and $D^{(m)}$ denotes the $m^{th}$ order derivative of $D$.

##  6 Numerical Scheme

Consider the most general appearance of second order differential equations as

$\rm{Y^{\prime\prime}(\zeta)=F(\zeta,Y(\zeta),Y^{\prime}(\zeta)),\hskip
14.22636pt0\leq\zeta\leq L}.$ (35)

with

$\rm{Y(\zeta)|_{\zeta=0}=\alpha_{0},\>Y^{\prime}(\zeta)|_{\zeta=0}=\alpha_{1}},$
(36)

or

$\rm{Y(\zeta)|_{\zeta=0}=\beta_{0},\>Y(\zeta)|_{\zeta=L}=\beta_{1}}.$ (37)

Let $\rm{\bar{Y}(\zeta)}$ be the Vieta-Lucas wavelet series approximation to
the solution of equations (35)-(37)

$\rm{\bar{Y}(\zeta)=\sum_{s=1}^{2^{k-1}}\sum_{m=0}^{M-1}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)=\Lambda^{T}\Upsilon(\zeta)}.$
(38)

Now by using (34), we obtain

$\rm{\bar{Y}^{\prime}(\zeta)=\Lambda^{\rm{T}}D\Upsilon(\zeta)}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\rm{\bar{Y}^{\prime\prime}(\zeta)=\Lambda^{\rm{T}}D^{(2)}\Upsilon(\zeta)}.$
(39)

The residual function $\rm{R(\zeta)}$ can be obtained by using equations (38)
and (39) in equation (35) as

$\rm{R(\zeta)=\Lambda^{\rm{T}}D^{(2)}\Upsilon(\zeta)-F(\zeta,\Lambda^{\rm{T}}\Upsilon(\zeta),\Lambda^{\rm{T}}D\Upsilon(\zeta))}.$
(40)

The corresponding updation in equations (35) and (36) are as follows

$\rm{\Lambda^{T}\Upsilon(0)=\alpha_{0},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\Lambda^{T}D\Upsilon(0)=\alpha_{1}},$ (41)

or

$\rm{\Lambda^{T}\Upsilon(0)=\beta_{0},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \Lambda^{T}\Upsilon(l)=\beta_{1}}.$
(42)

When $\rm{R(\zeta)}$ is zero, the exact solution is obtained, although it is
quite difficult to get $\rm{R(\zeta)}$ identically zero. So our primary focus
is to make the residual value as small as possible. Thus, we utilize the
weighted residual methods in order to minimize the residual function. Now, the
weighted residual equation is written as

$\rm{\langle
W,R(\zeta)\rangle}=\int_{0}^{2}WR(\zeta)\,d\zeta=0,\leavevmode\nobreak\
\leavevmode\nobreak\ $ (43)

where $\rm{W}$ is the weighted function in integral sense.
Thus, three weighted residual approaches are presented which are based on the
different choices of $\rm{W}$.

Approach - I : Collocation Approach
In this approach, the weighted function is chosen as Dirac delta($\delta$)
function which vanishes everywhere except at the collocation points, this
yields

$\rm{\langle\delta(\zeta-\zeta_{i}),R(\zeta)\rangle=\int_{0}^{2}\delta(\zeta-\zeta_{i})R(\zeta)\,d\zeta,\leavevmode\nobreak\
\leavevmode\nobreak\ i=1,2,3,\dots,2^{k-1}M-2,}$

which gives

$\rm{R(\zeta_{i})}=0,\hskip 14.22636pti=1,2,3,\dots,2^{k-1}M-2.$ (44)

Here, the extrema of Vieta-Lucas polynomials are chosen as the collocation
points. Substituting the collocation points in equation (44) leads to the
following system of algebraic equations as

$\rm{\Lambda^{\rm{T}}D^{(2)}\Upsilon(\zeta_{i})-F(\zeta_{i},\Lambda^{\rm{T}}\Upsilon(\zeta_{i}),\Lambda^{\rm{T}}D\Upsilon(\zeta_{i}))=0},\hskip
14.22636pti=1,2,3,\dots,2^{k-1}M-2.$ (45)

Now, the $\rm{2^{k-1}M}$ system of equations (equation (45) with conditions
(41) or (42)) yield the values of unknown coefficients, and thus the required
solution. The solution achieved in this sense will be called as
$\rm{Y_{VLWC}(\zeta)}$ solution.

Approach - II : Tau Approach
In this approach, the weighted function is chosen to be the same as the Vieta-
Lucas wavelets $\Upsilon_{i}$, which yields $\rm{(2^{k-1}M-2)}$ nonlinear
equations as

$\rm{\langle\Upsilon_{i},R(\zeta)\rangle=\int_{0}^{2}\Upsilon_{i}\leavevmode\nobreak\
R(\zeta)\,d\zeta=0,\hskip 19.91684pti=1,2,3,\dots,2^{k-1}M-2.}$ (46)

Thus, solving $\rm{2^{k-1}M}$ equations (equation (46) with (41) or (42)), we
get the unknown coefficients which leads to the appropriate solution, called
as $\rm{Y_{VLWT}(\zeta)}$.

Approach - III : Galerkin Approach
The main idea behind this approach is to expand the solutions not only in
terms of usual Vieta-Lucas wavelets expansion, but with some combinations of
Vieta-Lucas wavelets which satisfy the boundary requirements. In this
approach, the weighted functions are taken as the trial series solution, which
gives

$\rm{\langle\Psi_{i},R(\zeta)\rangle=\int_{0}^{2}\Psi_{i}\leavevmode\nobreak\
{R}(\zeta)\,d\zeta=0,\hskip 19.91684pti=1,2,3,\dots,2^{k-1}M,}$ (47)

where the trial series solution for Galerkin approach is written as

$\rm{{Y_{VLWG}}(\zeta)=\Lambda^{T}\Psi(\zeta)},$ (48)

and $\rm{\Psi(\zeta)=\nu(\zeta)\Upsilon(\zeta)}$, where $\rm{\nu(\zeta)}$ be
the trial function. Now, solving $\rm{2^{k-1}M}$ system of equations from (47)
provide the values of unknown coefficients, and thus we get the required
solution $\rm{{Y_{VLWG}}(\zeta)}$.

##  7 Convergence and error bound estimation

######  Theorem 7.1.

Let $\rm{Y(\zeta)}\in L^{2}_{w_{s}}[0,2]$ then the Vieta-Lucas series defined
in (21) converges to $\rm{Y(\zeta)}$ by the use of Vieta-Lucas wavelets i.e.

$\rm{Y(\zeta)=\sum_{s=1}^{\infty}\sum_{m=0}^{\infty}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)}.$

###### Proof.

Suppose $\rm{Y(\zeta)}\in L^{2}_{w_{s}}[0,2]$, where $L^{2}_{w_{s}}[0,2]$ be
the Hilbert space and $\rm{\Upsilon_{s,m}}$ defined in (19) forms an
orthonormal basis with respect to weight function
$w_{s}(\zeta)=w(2^{k}\zeta-\hat{s}).$
Consider, $\rm{Y(\zeta)=\sum_{m=0}^{M-1}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)}$
where $\rm{\Lambda_{s,m}=\langle
Y(\zeta),\Upsilon_{s,m}(\zeta)\rangle_{w_{s}(\zeta)}}$ for a fixed $\rm{s}$
and $\rm{\Upsilon_{s,m}(\zeta)=\Upsilon_{j}(\zeta)}$,
$\rm{\varsigma_{j}=\langle Y(\zeta),\Upsilon_{j}(\zeta)\rangle}$. The sequence
of partial sum $\\{\rm{\rho_{s}}\\}$ of
$\\{\rm{\varsigma_{j}\Upsilon_{j}(\zeta)}\\}_{j\geq 1}$ is defined as
$\rm{\rho_{s}=\sum_{j=1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta)}.$
Now,

$\rm{\langle Y(\zeta),\rho_{s}\rangle=\langle
Y(\zeta),\sum_{j=1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta)\rangle=\sum_{j=1}^{s}|\varsigma_{j}|^{2}}.$

For s $>$ m, we assert that

$\rm{\left\lVert\rho_{s}-\rho_{m}\right\rVert^{2}=\sum_{j=m+1}^{s}|\varsigma_{j}|^{2}}.$

Now,

$\rm{\left\lVert\sum_{j=m+1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta)\right\rVert^{2}=\langle\sum_{j=m+1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta),\sum_{j=m+1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta)\rangle=\sum_{j=m+1}^{s}|\varsigma_{j}|^{2},\hskip
5.69046pts>m.}$

Therefore,

$\rm{\left\lVert\sum_{j=m+1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta)\right\rVert^{2}=\sum_{j=m+1}^{s}|\varsigma_{j}|^{2},\hskip
5.69046pts>m.}$

From Bessel’s inequality, we know
$\rm{\sum_{j=1}^{\infty}|\varsigma_{j}|^{2}}$ is convergent.
Thus we have

$\rm{\left\lVert\sum_{j=m+1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta)\right\rVert^{2}}\rightarrow
0\hskip 14.22636pt\text{as}\hskip 5.69046pts\rightarrow\infty.$

Which implies

$\rm{\left\lVert\sum_{j=m+1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta)\right\rVert}\rightarrow
0,$

and $\\{\rm{\rho_{s}}\\}$ is a Cauchy sequence that converges to
$\rm{\rho}$(say).
Thus,

$\begin{split}\rm{\langle\rho-Y(\zeta),\Upsilon(\zeta)\rangle}&=\rm{\langle\rho,\Upsilon(\zeta)\rangle-\langle
Y(\zeta),\Upsilon(\zeta)\rangle,}\\\
&=\rm{\langle\lim_{s\rightarrow\infty}\rho_{s},\Upsilon(\zeta)\rangle-\varsigma_{j},}\\\
&=\rm{\varsigma_{j}-\varsigma_{j}}.\\\ \end{split}$

Which implies

$\rm{\langle\rho-Y(\zeta),\Upsilon(\zeta)\rangle=0}.$

Therefore, $\rm{Y(\zeta)=\rho}$ and
$\rm{\sum_{j=1}^{s}\varsigma_{j}\Upsilon_{j}(\zeta)}$ converges to
$\rm{Y(\zeta)}$ for $s\rightarrow\infty$. ∎

######  Theorem 7.2.

Let $Y(\zeta)$ be a second order square integrable function defined over
$[0,2]$ with bounded second order derivative say
$|Y^{\prime\prime}(\zeta)|\leq H$ for some constant H. Then $Y(\zeta)$ can be
expanded as an infinite sum of Vieta Lucas wavelets and the series converges
to $Y(\zeta)$ uniformly, that is

$\rm{Y(\zeta)=\sum_{s=1}^{\infty}\sum_{m=0}^{\infty}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)},$

where $\rm{\Lambda_{s,m}=\langle
Y(\zeta),\Upsilon_{s,m}(\zeta)\rangle_{{L^{2}}_{w_{s}}[0,2]}}.$

###### Proof.

From (22), we have

$\displaystyle\rm{\Lambda_{s,m}=\langle
Y(\zeta),\Upsilon_{s,m}(\zeta)\rangle_{{L^{2}}_{w_{s}}[0,2]}=\int_{0}^{2}Y(\zeta)\Upsilon_{s,m}(\zeta)w_{s}(\zeta)\,d\zeta},$
$\displaystyle\rm{=\int_{\frac{\hat{s}-2}{2^{k}}}^{\frac{\hat{s}+2}{2^{k}}}Y(\zeta)2^{\frac{k}{2}}\sqrt{\frac{1}{2\pi}}{\rm{VL_{m}}}(2^{k}\zeta-\hat{s})\,w_{s}(2^{k}\zeta-\hat{s})\,d\zeta},$

where $\hat{s}=2(2s-1)$. Substituting $2^{k}\zeta-\hat{s}=2\cos{\delta}$ and
from the definition of Vieta-Lucas polynomial, we obtain

$\displaystyle\rm{\Lambda_{s,m}}=$
$\displaystyle\rm{2^{\frac{-k}{2}}\sqrt{\frac{1}{2\pi}}\int_{0}^{\pi}Y\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\rm{VL_{m}}}(\delta)\,w_{s}(\delta)\,d\delta,$
$\displaystyle=$
$\displaystyle\rm{2^{\frac{-k}{2}}\sqrt{\frac{1}{2\pi}}\int_{0}^{\pi}Y\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)2\cos{m\delta}\,\frac{1}{\sqrt{4-(2\cos{\delta})^{2}}}\,2\sin{\delta}\,d\delta},$
$\displaystyle=$
$\displaystyle\rm{2^{\frac{-k}{2}}\sqrt{\frac{1}{2\pi}}\int_{0}^{\pi}Y\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)2\cos{m\delta}\,\frac{1}{2\sqrt{1-\cos^{2}{\delta}}}\,2\sin{\delta}\,d\delta},$
$\displaystyle=$
$\displaystyle\rm{2^{\frac{-k+1}{2}}\sqrt{\frac{1}{\pi}}\int_{0}^{\pi}Y\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\cos{m\delta}\,d\delta},$

Using the integration by parts, we get

$\displaystyle\rm{\Lambda_{s,m}}=$
$\displaystyle\rm{\frac{2^{\frac{-3k+1}{2}}}{m}\sqrt{\frac{1}{\pi}}\int_{0}^{\pi}Y^{\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)2\sin{m\delta}\sin{\delta}\,d\delta},$
$\displaystyle=$
$\displaystyle\rm{\frac{2^{\frac{-3k+1}{2}}}{m}\sqrt{\frac{1}{\pi}}\int_{0}^{\pi}Y^{\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)(\cos{(m-1)\delta}-\cos{(m+1)\delta})\,d\delta},$
$\displaystyle=$
$\displaystyle\rm{\frac{2^{\frac{-3k+1}{2}}}{m}\sqrt{\frac{1}{\pi}}\left(\int_{0}^{\pi}Y^{\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\cos{(m-1)\delta}\,d\delta-\int_{0}^{\pi}Y^{\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\cos{(m+1)\delta}\,d\delta\right),}$
$\implies\rm{|\Lambda_{s,m}|}\leq\leavevmode\nobreak\
\rm{\frac{2^{\frac{-3k+1}{2}}}{m}\sqrt{\frac{1}{\pi}}(|I_{1}|+|I_{2}|)},$ (49)

where
$\rm{I_{1}}=\int_{0}^{\pi}Y^{\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\cos{(m-1)\delta}\,d\delta$
and
$\rm{I_{2}}=\int_{0}^{\pi}Y^{\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\cos{(m+1)\delta}\,d\delta$.
Next we estimate $\rm{I_{1}}$ and $\rm{I_{2}}$ respectively. On Integrating
$\rm{I_{1}}$ by parts, we have

$\displaystyle\rm{I_{1}}=$
$\displaystyle\rm{\frac{2^{-k+1}}{m-1}\int_{0}^{\pi}Y^{\prime\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\sin{(m-1)\delta}\sin{\delta}\,d\delta},$
$\displaystyle\rm{|I_{1}|}\leq$
$\displaystyle\rm{\frac{2^{-k+1}}{m-1}\int_{0}^{\pi}\left|Y^{\prime\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\right|\,d\delta},$

which gives

$\rm{|I_{1}|}\leq\rm{\frac{H\pi 2^{-k+1}}{m-1}}.$ (50)

Similarly on Integrating $\rm{I_{2}}$ , we obtain

$\displaystyle\rm{I_{2}}=$
$\displaystyle\rm{\frac{2^{-k+1}}{m+1}\int_{0}^{\pi}Y^{\prime\prime}\left(\frac{2\cos{\delta}+\hat{s}}{2^{k}}\right)\sin{(m+1)\delta}\sin{\delta}\,d\delta},$

which gives

$\rm{|I_{2}|}\leq\rm{\frac{H\pi 2^{-k+1}}{m+1}}.$ (51)

By using (50) and (51) in (49), we obtain

$\displaystyle\rm{|\Lambda_{s,m}|\leq\frac{H\sqrt{\pi}\,2^{\frac{-5k+5}{2}}}{m^{2}-1}},\hskip
14.22636ptm>1.$

Since $\rm{s\leq 2^{k-1}}$ and $\rm{s\geq 1}$. Therefore we get

$\rm{|\Lambda_{s,m}|\leq\frac{H\sqrt{\pi}\,}{s^{\frac{5}{2}}(m^{2}-1)}}.$ (52)

For m = 1, we have

$\displaystyle\rm{|\Lambda_{s,1}|\leq\frac{\sqrt{\pi}\,}{2\leavevmode\nobreak\
s^{\frac{3}{2}}(m^{2}-1)}\max_{0\leq\zeta\leq 2}|\rm{Y^{\prime}(\zeta)}|}.$

Also,

$\displaystyle\left|\rm{\sum_{s=1}^{\infty}\sum_{m=0}^{\infty}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)}\right|$
$\displaystyle\leq\rm{\left|\sum_{s=1}^{\infty}\Lambda_{s,0}\Upsilon_{s,0}(\zeta)\right|+\sum_{s=1}^{\infty}\sum_{m=1}^{\infty}|\Lambda_{s,m}||\Upsilon_{s,m}(\zeta)|}$
$\displaystyle\leq\rm{\left|\sum_{s=1}^{\infty}\Lambda_{s,0}\Upsilon_{s,0}(\zeta)\right|+\sum_{s=1}^{\infty}\sum_{m=1}^{\infty}|\Lambda_{s,m}|}<\infty.$

which implies that the series
$\rm{\sum_{s=1}^{2^{k-1}}\sum_{m=0}^{M-1}\Lambda_{s,m}}$ is absolutely
convergent. Thus the series
$\rm{\sum_{s=1}^{2^{k-1}}\sum_{m=0}^{M-1}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)}$
converges to $\rm{Y(\zeta)}$ uniformly. ∎

######  Lemma 7.3.

Let $f(\zeta)$ be a continuous, positive, decreasing function for $\zeta\geq
m$ if $f(k)=\Lambda_{k}$, provided that $\sum\Lambda_{s}$ is convergent and
$R_{s}=\sum_{k=s+1}^{\infty}\Lambda_{k}$, then
$R_{s}\leq\int_{s}^{\infty}f(\zeta)d\zeta$ [44].

######  Theorem 7.4.

For the Vieta-Lucas wavelets expansion, if $Y(\zeta)$ satisfies the theorem
(7.2), then the following error estimate holds for $M>2$,

$\rm{\left\lVert
Y(\zeta)-\bar{Y}(\zeta)\right\rVert_{w_{s}}}<H\sqrt{\rm{\pi\frac{(M^{2}-2M)ln(M)-M^{2}ln(M-2)+(2ln(M-2)-2)M+2}{(2^{5(2^{k-1}-1)}5\ln{2})\leavevmode\nobreak\
4M(M-2)}}}.$

###### Proof.

Considering the Vieta-Lucas wavelets expansion as

$\rm{\bar{Y}(\zeta)=\sum_{s=1}^{2^{k-1}}\sum_{m=0}^{M-1}\Lambda_{s,m}\Upsilon_{s,m}(\zeta)}.$

From (19), by using the orthogonality property of Vieta-Lucas wavelets with
respect to weight function $w_{s}(\zeta)$, we obtain

$\rm{\left\lVert
Y(\zeta)-\bar{Y}(\zeta)\right\rVert^{2}_{w_{s}}=\sum_{s=2^{k-1}+1}^{\infty}\sum_{m=M}^{\infty}\Lambda^{2}_{s,m}}.$

By using theorem (7.2), it can be expressed as

$\rm{\left\lVert
Y(\zeta)-\bar{Y}(\zeta)\right\rVert^{2}_{w_{s}}<H^{2}\pi\sum_{s=2^{k-1}+1}^{\infty}\sum_{m=M}^{\infty}\frac{1}{2^{5k-5}(m^{2}-1)^{2}}}.$

From lemma (7.3), we obtain

$\begin{split}\rm{\left\lVert
Y(\zeta)-\bar{Y}(\zeta)\right\rVert^{2}_{w_{s}}}&<\rm{H^{2}\pi\int_{s=2^{k-1}}^{\infty}\int_{m=M-1}^{\infty}\frac{1}{2^{5\zeta-5}(z^{2}-1)^{2}}d\zeta
dz,}\\\
&=\rm{H^{2}\pi\left(\int_{s=2^{k-1}}^{\infty}\frac{1}{2^{5\zeta-5}}d\zeta\right)\left(\int_{m=M-1}^{\infty}\frac{1}{(z^{2}-1)^{2}}dz\right),}\\\
&=\rm{H^{2}\pi\left(\frac{1}{2^{5(2^{k-1}-1)}5\ln{2}}\right)\left(\frac{(M^{2}-2M)ln(M)-M^{2}ln(M-2)+(2ln(M-2)-2)M+2}{4M(M-2)}\right)}.\end{split}$

which completes the proof. ∎

##  8 Numerical Examples

This section contains numerical examples that demonstrate the efficiency and
reliability of the presented schemes.

######  Example 8.1.

Consider the non-homogeneous SDE [46]:

$\rm{Y^{\prime\prime}(\zeta)+\frac{1}{\zeta}Y^{\prime}(\zeta)={\left(\frac{8}{8-\zeta^{2}}\right)}^{2},\hskip
14.22636pt0\leq\zeta\leq 1,}$ (53)

with

$\rm{Y(\zeta)|_{\zeta=0}=0,\hskip
11.38092ptY^{\prime}(\zeta)|_{\zeta=0}=0.}\\\ $ (54)

The exact solution is

$\rm{Y_{Exact}(\zeta)=2\log{\left(\frac{8}{8-\zeta^{2}}\right)}}.$

The following solutions are obtained by using the proposed numerical schemes
at $\eta$ = 12:

$\begin{split}\rm{Y_{VLWC}(\zeta)}&=(1.60\times 10^{-16})-(1.18\times
10^{-16})\zeta+\dots+(1.82\times 10^{-5})\zeta^{11},\\\
\rm{Y_{VLWT}(\zeta)}&=(1.83\times 10^{-16})+(4.51\times
10^{-17})\zeta+\dots+(8.79\times 10^{-4})\zeta^{11},\\\
\rm{Y_{VLWG}(\zeta)}&=0.25\zeta^{2}-0.08\zeta^{3}+0.40\zeta^{4}+\dots-(9.75\times
10^{-4}\zeta^{13}).\end{split}$

######  Example 8.2.

Take the nonlinear SDE as [47]:

$\rm{Y^{\prime\prime}(\zeta)+\pi^{3}\frac{(Y^{2}(\zeta))}{\sin{(\pi\zeta)}}=0,\hskip
14.22636pt0<\zeta<1,}$ (55)

with boundary restrictions

$\rm{Y(\zeta)|_{\zeta=0}=0,\hskip 11.38092ptY(\zeta)|_{\zeta=1}=0.}$ (56)

The exact solution is

$\rm{Y_{Exact}(\zeta)=\frac{\sin{(\pi\zeta)}}{\pi}}.$

The following solutions are obtained for $\eta$ = 6:

$\begin{split}\rm{Y_{VLWC}(\zeta)}&=(1.38\times
10^{-16})+0.99\zeta+0.10\zeta^{2}-2.20\zeta^{3}+1.10\zeta^{4}-(7.04\times
10^{-13})\zeta^{5},\\\ \rm{Y_{VLWT}(\zeta)}&=(-2.77\times
10^{-17})+0.93\zeta+0.64\zeta^{2}-3.62\zeta^{3}+2.56\zeta^{4}-0.51\zeta^{5},\\\
\rm{Y_{VLWG}(\zeta)}&=1.00\zeta-0.03\zeta^{2}-1.39\zeta^{3}-0.75\zeta^{4}+1.96\zeta^{5}-0.90\zeta^{6}+0.12\zeta^{7}.\end{split}$

######  Example 8.3.

Consider the nonhomogeneous Emden-Fowler type SDE [11]:

$\rm{Y^{\prime\prime}(\zeta)+\frac{8}{\zeta}Y^{\prime}(\zeta)+\zeta
Y(\zeta)=\zeta^{5}-\zeta^{4}+44\zeta^{2}-30\zeta,\hskip 14.22636pt\zeta\geq
0,}$ (57)

with

$\rm{Y(\zeta)|_{\zeta=0}=0,\hskip 11.38092ptY^{\prime}(\zeta)|_{\zeta=0}=0.}$
(58)

The following solutions are obtained for $\eta$ = 5:

$\begin{split}\rm{Y_{VLWC}(\zeta)}&=\rm{(5.55\times 10^{-16})+(9.32\times
10^{-15})\zeta^{2}-\zeta^{3}+\zeta^{4}},\\\
\rm{Y_{VLWT}(\zeta)}&=\rm{\zeta^{4}-\zeta^{3}},\\\
\rm{Y_{VLWG}(\zeta)}&=\rm{\zeta^{4}-\zeta^{3}}.\end{split}$

The above solutions are computed in the interval $0\leq\zeta\leq 2$ and it is
observed that both $\rm{Y_{VLWT}}$ and $\rm{Y_{VLWG}}$ yield the exact
solution.

######  Example 8.4.

Consider the nonlinear Emden-Fowler type SDE [48]:

$\rm{Y^{\prime\prime}(\zeta)+\frac{8}{\zeta}Y^{\prime}(\zeta)+18Y(\zeta)=-4Y(\zeta)\ln{(Y(\zeta})),\hskip
14.22636pt0<\zeta\leq 1,}$ (59)

with initial conditions

$\rm{Y(\zeta)|_{\zeta=0}=1,\hskip 11.38092ptY^{\prime}(\zeta)|_{\zeta=0}=0.}$
(60)

The approximated solutions obtained by using proposed schemes for $\eta$ = 3 :

$\displaystyle\rm{Y_{VLWC}(\zeta)}$ $\displaystyle=\rm{e^{(-3.33067\times
10^{-16})-\zeta^{2}}},$ $\displaystyle\rm{Y_{VLWT}(\zeta)}$
$\displaystyle=\rm{e^{-\zeta^{2}}},$ $\displaystyle\rm{Y_{VLWG}(\zeta)}$
$\displaystyle=\rm{e^{-\zeta^{2}}}.$

Here, $\rm{Y_{VLWT}}$ and $\rm{Y_{VLWG}}$ leads to the exact solution.

######  Example 8.5.

Let us consider the nonlinear Lane-Emden type SDE [11]:

$\rm{Y^{\prime\prime}(\zeta)+\frac{2}{\zeta}Y^{\prime}(\zeta)-6Y(\zeta)=4Y(\zeta)\ln{(Y(\zeta))},\hskip
14.22636pt0<\zeta\leq 1,}$ (61)

with initial conditions

$\rm{Y(\zeta)|_{\zeta=0}=1,\hskip 11.38092ptY^{\prime}(\zeta)|_{\zeta=0}=0.}$
(62)

The approximated solutions obtained by using proposed schemes at $\eta$ = 3:

$\displaystyle\rm{Y_{VLWC}(\zeta)}$ $\displaystyle=\rm{e^{(3.8857\times
10^{-16})+\zeta^{2}}},$ $\displaystyle\rm{Y_{VLWT}(\zeta)}$
$\displaystyle=\rm{e^{\zeta^{2}}},$ $\displaystyle\rm{Y_{VLWG}(\zeta)}$
$\displaystyle=\rm{e^{\zeta^{2}}}.$

The solutions obtained by $\rm{Y_{VLWT}}$ and $\rm{Y_{VLWG}}$ yield the exact
solution.

##  9 Results and Discussions

(a) $\eta=6$

(b) $\eta=12$

(c) $\eta=6$

(d) $\eta=8$

Figure 1:  Solution curves for (a)-(b) Example 8.1 and (c)-(d) Example 8.2.

(a) $\eta=5$

(b) $\eta=3$

(c) $\eta=3$

Figure 2:  Solution curves for (a) Example 8.3 (b) Example 8.4 and (c) Example
8.5.

(a) $\rm{Y_{VLWC}(\zeta)}$

(b) $\rm{Y_{VLWT}(\zeta)}$

(c) $\rm{Y_{VLWG}(\zeta)}$

Figure 3:  Comparison of logarithmic values of absolute errors at different
resolutions for Example 8.1 by proposed approaches.

(a) $\eta=8$

(b) $\eta=10$

(c) $\eta=12$

Figure 4:  Comparison of logarithmic values of absolute errors of proposed
schemes for Example 8.1 at $\eta=8,10$ and 12.

Figure 1 demonstrates the solution curves of the exact solution and
approximate solutions by the proposed numerical approaches for Example 8.1 and
8.2 at different resolutions. It is observed from the zoomed profile of figure
1 that as the resolution increases the approximate solutions become more
accurate and overlap the exact solution curve. The solution plots for Example
8.3, 8.4 and 8.5 are shown in figure 2. It can be seen from figure 2 that the
proposed approaches provide accurate results at the smaller resolutions which
proves that all three proposed numerical approaches are reliable. Figure 3
depicts a comparison of logarithmic values of absolute errors at various
resolutions for Example 8.1, as obtained by proposed numerical approaches. It
demonstrates that the errors are highest for $\eta=6$ and least for $\eta=12$,
indicating that as the resolution increases, the error reduces for all three
proposed numerical schemes. Figure 4 shows a comparative analysis of the
presented numerical schemes at various resolutions, and it is observed that
the errors are bounded for all of the proposed schemes, demonstrating the
efficiency and accuracy of the schemes.
The absolute error obtained by the proposed methods for Example 8.1 and
Example 8.2 are compared in Table 1. In Example 8.3, 8.4 and 8.5,
$\rm{Y_{VLWT}(\zeta)}$ and $\rm{Y_{VLWG}(\zeta)}$ produce the exact solution.
Table 2 compares the absolute error obtained by $\rm{Y_{VLWC}(\zeta)}$ with
the existing findings. The solutions obtained from the proposed approaches are
in good agreement with the existing results, demonstrating the reliability of
the proposed schemes.

Table 1:  Absolute error comparisons for Example 8.1 and Example 8.2. Example
8.1
---
$\zeta$ |  Exact Solution |  $\rm{\mid Y(\zeta)-Y_{VLWC}(\zeta)\mid}$ |  $\rm{\mid Y(\zeta)-Y_{VLWT}(\zeta)\mid}$ |  $\rm{\mid Y(\zeta)-Y_{VLWG}(\zeta)\mid}$
0.1 |  0.0025015 |  1.3701E-12 |  9.4083E-07 |  4.1019E-05
0.2 |  0.0100250 |  9.0015E-13 |  3.9899E-07 |  7.1956E-05
0.3 |  0.0226275 |  1.1882E-12 |  5.3283E-08 |  8.3708E-05
0.4 |  0.0404054 |  1.4496E-12 |  4.7147E-07 |  9.2989E-05
0.5 |  0.0634973 |  1.5533E-12 |  1.0288E-06 |  1.0235E-04
0.6 |  0.0920878 |  2.0826E-12 |  8.6287E-07 |  1.0947E-04
0.7 |  0.1264121 |  8.8321E-13 |  2.4217E-07 |  1.1425E-04
0.8 |  0.1667632 |  7.3194E-13 |  9.5880E-08 |  1.1836E-04
0.9 |  0.2135007 |  1.8588E-12 |  1.9422E-07 |  1.2276E-04
1.0 |  0.2670627 |  9.3192E-13 |  7.7329E-07 |  1.2706E-04
Example 8.2
0.1 |  0.0983631 |  1.7106E-04 |  1.7790E-03 |  2.0126E-04
0.2 |  0.1870978 |  3.9765E-04 |  6.3775E-04 |  3.4647E-04
0.3 |  0.2575181 |  1.2307E-04 |  2.6134E-03 |  5.0466E-04
0.4 |  0.3027306 |  3.6135E-04 |  2.6338E-03 |  5.4599E-04
0.5 |  0.3183098 |  5.9103E-04 |  1.0177E-03 |  4.0975E-04
0.6 |  0.3027306 |  3.6135E-04 |  1.1302E-03 |  1.5710E-04
0.7 |  0.2575181 |  1.2307E-04 |  2.6834E-03 |  8.3777E-05
0.8 |  0.1870978 |  3.9765E-04 |  2.9582E-03 |  2.0029E-04
0.9 |  0.0983631 |  1.7106E-04 |  1.8955E-03 |  1.5510E-04
Table 2:  Absolute error comparisons for Example 8.3, Example 8.4 and Example
8.5.  Example 8.3
---
$\zeta$ |  Exact Solution |  HFC Error[24] |  $\rm{\mid Y(\zeta)-Y_{VLWC}(\zeta)\mid}$
0.01 |  -0.0000009 |  5.7E-08 |  5.5603E-16
0.10 |  -0.0009000 |  8.4E-08 |  6.3881E-16
0.50 |  -0.0625000 |  2.2E-06 |  1.8318E-15
1.00 |  0.0000000 |  8.2E-07 |  2.7755E-15
2.00 |  8.0000000 |  1.7E-07 |  1.7763E-15
Example 8.4
$\zeta$ |  Exact Solution |  LWM Error[45] |  $\rm{\mid Y(\zeta)-Y_{VLWC}(\zeta)\mid}$
0.1 |  0.9900498 |  4.8E-06 |  3.3306E-16
0.2 |  0.9607894 |  6.8E-06 |  3.3306E-16
0.3 |  0.9139311 |  8.0E-07 |  3.3306E-16
0.4 |  0.8521437 |  8.3E-06 |  3.3306E-16
0.5 |  0.7788007 |  1.2E-05 |  2.2204E-16
0.6 |  0.6976763 |  5.3E-05 |  1.1102E-16
0.7 |  0.6126263 |  2.0E-04 |  1.1102E-16
0.8 |  0.5272924 |  5.9E-04 |  2.2204E-16
0.9 |  0.4448580 |  1.4E-03 |  5.5511E-17
1.0 |  0.3678794 |  3.0E-03 |  5.6511E-17
Example 8.5
$\zeta$ |  Exact Solution |  HFC Error[24] |  $\rm{\mid Y(\zeta)-Y_{VLWC}(\zeta)\mid}$
0.01 |  1.0001000 |  2.2E-08 |  4.4408E-16
0.02 |  1.0004000 |  1.5E-08 |  2.2204E-16
0.05 |  1.0025031 |  2.1E-08 |  2.2204E-16
0.1 |  1.0100501 |  1.7E-08 |  4.4408E-16
0.2 |  1.0408107 |  2.1E-08 |  4.4408E-16
0.5 |  1.2840254 |  3.0E-08 |  4.4408E-16
0.7 |  1.6323162 |  4.2E-08 |  6.6613E-16
0.8 |  1.8964808 |  5.1E-08 |  2.2204E-16
0.9 |  2.2479079 |  9.2E-08 |  4.4408E-16
1.0 |  2.7182818 |  8.8E-08 |  8.8817E-16

##  10 Conclusion

In this work, the second-order SDEs are analyzed by using the novel numerical
schemes which are based on Vieta-Lucas wavelets. The operational matrix
approach combined with three weighted residual approaches is used to solve
singular IVPs and singular BVPs. Simultaneously two convergence theorems are
proved, and the upper bound for the tolerable error is explored. The numerical
examples are taken to demonstrate the efficacy of the proposed approaches. The
approximate solutions derived through the proposed methods are found to be
extremely accurate, which proves the accuracy and applicability.

##  11 Declaration

Conflict of interest: Authors have no conflict of interest.
Availability of data and material: Not applicable.
Compliance with ethical standards.

## References

  * [1] Stakgold I (2000) Boundary Value Problems of Mathematical Physics: Volume 1. Society for Industrial and Applied Mathematics
  * [2] Gatica JA, Oliker V, Waltman P (1989) Singular nonlinear boundary value problems for second-order ordinary differential equations. Journal of Differential Equations 79(1):62-78
  * [3] Sabir Z, Baleanu D, Shoaib M, Raja MA (2021) Design of stochastic numerical solver for the solution of singular three-point second-order boundary value problems. Neural Computing and Applications 33:2427-43
  * [4] Sabir Z, Günerhan H, Guirao JL (2020) On a new model based on third-order nonlinear multi singular functional differential equations. Mathematical Problems in Engineering:1-9
  * [5] Sabir Z, Manzar MA, Raja MA, Sheraz M, Wazwaz AM (2018) Neuro-heuristics for nonlinear singular Thomas-Fermi systems. Applied Soft Computing 65:152-69
  * [6] Sabir Z, Wahab HA, Javeed S, Baskonus HM (2021) An efficient stochastic numerical computing framework for the nonlinear higher order singular models. Fractal and Fractional 5(4):176
  * [7] da Costa Campos LM (2019) Singular Differential Equations and Special Functions. CRC Press
  * [8] Ford WF, Pennline JA (2009) Singular non-linear two-point boundary value problems: Existence and uniqueness. Nonlinear Analysis: Theory, Methods & Applications 71(3-4):1059-72
  * [9] Bender CM, Milton KA, Pinsky SS, Simmons Jr LM (1989) A new perturbative approach to nonlinear problems. Journal of Mathematical Physics 30(7):1447-55
  * [10] Russell RD, Shampine LF (1975) Numerical methods for singular boundary value problems. SIAM Journal on Numerical Analysis 12(1):13-36
  * [11] Wazwaz AM (2002) A new method for solving singular initial value problems in the second-order ordinary differential equations. Applied Mathematics and Computation 128(1):45-57
  * [12] Wazwaz AM (2011) The variational iteration method for solving nonlinear singular boundary value problems arising in various physical models. Communications in Nonlinear Science and Numerical Simulation 16(10):3881-6
  * [13] El-Gebeily M, O’Regan D (2007) A quasilinearization method for a class of second order singular nonlinear differential equations with nonlinear boundary conditions. Nonlinear Analysis: Real World Applications 8(1):174-86
  * [14] Yıldırım A, Öziş T (2007) Solutions of singular IVPs of Lane–Emden type by homotopy perturbation method. Physics Letters A 369(1-2):70-6
  * [15] Pandey RK (1997) A finite difference method for a class of singular two point boundary value problems arising in physiology. International journal of computer mathematics 65(1-2):131-40
  * [16] Guirao JL, Sabir Z, Saeed T (2020) Design and numerical solutions of a novel third-order nonlinear Emden–Fowler delay differential model. Mathematical Problems in Engineering:1-9
  * [17] Sabir Z, Amin F, Pohl D, Guirao JL (2020) Intelligence computing approach for solving second order system of Emden–Fowler model. Journal of Intelligent & Fuzzy Systems 38(6):7391-406
  * [18] Sabir Z, Khalique CM, Raja MA, Baleanu D (2021) Evolutionary computing for nonlinear singular boundary value problems using neural network, genetic algorithm and active-set algorithm. The European Physical Journal Plus 136:1-9
  * [19] Kadalbajoo MK, Aggarwal VK (2005) Numerical solution of singular boundary value problems via Chebyshev polynomial and B-spline. Applied mathematics and computation 160(3):851-63
  * [20] El‐Sayed AA, Agarwal P (2019) Numerical solution of multiterm variable‐order fractional differential equations via shifted Legendre polynomials. Mathematical Methods in the Applied Sciences 42(11):3978-91
  * [21] Zhou F, Xu X (2016) Numerical solutions for the linear and nonlinear singular boundary value problems using Laguerre wavelets. Advances in Difference Equations 2016(1):1-5
  * [22] Baishya C, Veeresha P (2021) Laguerre polynomial-based operational matrix of integration for solving fractional differential equations with non-singular kernel. Proceedings of the Royal Society A 477(2253):20210438
  * [23] El-Sayed AA, Baleanu D, Agarwal P (2020) A novel Jacobi operational matrix for numerical solution of multi-term variable-order fractional differential equations. Journal of Taibah University for Science 14(1):963-74
  * [24] Parand K, Dehghan M, Rezaei AR, Ghaderi SM (2010) An approximation algorithm for the solution of the nonlinear Lane–Emden type equations arising in astrophysics using Hermite functions collocation method. Computer Physics Communications 181(6):1096-108
  * [25] Agarwal P, El-Sayed AA (2020) Vieta–Lucas polynomials for solving a fractional-order mathematical physics model. Advances in Difference Equations 2020(1):1-8
  * [26] Agarwal P, Qi F, Chand M, Singh G (2019) Some fractional differential equations involving generalized hypergeometric functions. Journal of Applied Analysis 25(1):37-44
  * [27] Zhou H, Yang L, Agarwal P (2017) Solvability for fractional p-Laplacian differential equations with multipoint boundary conditions at resonance on infinite interval. Journal of Applied Mathematics and Computing 53:51-76
  * [28] Agarwal P, Al-Mdallal Q, Cho YJ, Jain S (2018) Fractional differential equations for the generalized Mittag-Leffler function. Advances in difference equations:1-8
  * [29] Heydari MH, Avazzadeh Z, Razzaghi M (2021) Vieta-Lucas polynomials for the coupled nonlinear variable-order fractional Ginzburg-Landau equations. Applied Numerical Mathematics 165:442-58
  * [30] Torrence C, Compo GP (1998) A practical guide to wavelet analysis. Bulletin of the American Meteorological Society 79(1):61-78
  * [31] Mehandiratta V, Mehra M, Leugering G (2021) An approach based on Haar wavelet for the approximation of fractional calculus with application to initial and boundary value problems. Mathematical Methods in the Applied Sciences 44(4):3195-213
  * [32] Yuanlu LI (2010) Solving a nonlinear fractional differential equation using Chebyshev wavelets. Communications in Nonlinear Science and Numerical Simulation 15(9):2284-92
  * [33] Rostami M, Hashemizadeh E, Heidari M (2012) A comparative study of numerical integration based on block-pulse and sinc functions and Chebyshev wavelet. Mathematical Sciences 6:1-7
  * [34] Razzaghi M, Yousefi S (2001) The Legendre wavelets operational matrix of integration. International Journal of Systems Science 32(4):495-502
  * [35] Kumar S, Ghosh S, Kumar R, Jleli M (2021) A fractional model for population dynamics of two interacting species by using spectral and Hermite wavelets methods. Numerical Methods for Partial Differential Equations 37(2):1652-72
  * [36] Eslahchi MR, Kavoosi M (2018) The use of Jacobi wavelets for constrained approximation of rational Bézier curves. Computational and Applied Mathematics 37:3951-66
  * [37] Usman M, Hamid M, Ul Haq R, Wang W (2018) An efficient algorithm based on Gegenbauer wavelets for the solutions of variable-order fractional differential equations. The European Physical Journal Plus 133:1-6
  * [38] Kumar R, Koundal R, Srivastava K, Baleanu D (2020) Normalized Lucas wavelets: an application to Lane–Emden and pantograph differential equations. The European Physical Journal Plus 135:1-24
  * [39] Koundal R, Kumar R, Srivastava K, Baleanu D (2022) Lucas wavelet scheme for fractional Bagley–Torvik equations: Gauss–Jacobi approach. International Journal of Applied and Computational Mathematics 8(1):3
  * [40] Horadam AF (2002) Vieta polynomials. Fibonacci Quarterly 40(3):223-32
  * [41] Robbins N (1991) Vieta’s triangular array and a related family of polynomials. International Journal of Mathematics and Mathematical Sciences 14(2):239-44
  * [42] Izadi M, Yüzbaşı Ş, Ansari KJ (2021) Application of Vieta–Lucas series to solve a class of multi-pantograph delay differential equations with singularity. Symmetry 13(12):2370
  * [43] Ray SS, Gupta AK (2018) Wavelet methods for solving partial differential equations and fractional differential equations. CRC Press
  * [44] Stewart J (2012) Single variable essential calculus: early transcendentals. Cengage Learning
  * [45] Aminikhah H, Moradian S (2013) Numerical solution of singular Lane-Emden equation. International Scholarly Research Notices
  * [46] Nasab AK, Kılıçman A, Babolian E, Atabakan ZP (2013) Wavelet analysis method for solving linear and nonlinear singular boundary value problems. Applied Mathematical Modelling 37(8):5876-86
  * [47] Wazwaz AM (2001) A reliable algorithm for obtaining positive solutions for nonlinear boundary value problems. Computers & Mathematics with Applications 41(10-11):1237-44
  * [48] Wazwaz AM (2005) Adomian decomposition method for a reliable treatment of the Emden–Fowler equation. Applied Mathematics and Computation 161(2):543-60