Department of Mathematics, Faculty of Science, Minia University, Minia, 61519, EG
Department of Mathematics, Faculty of Science, New-Valley University, El-Khargah, 72511, EG
Abstract
The present study is related to the numerical solutions of new mathematical models based on the variable order Emden-Fowler delay differential equations. The shifted fractional Gegenbauer,
CS,j(α,μ)(t),
operational matrices (OMs) of VO differentiation, in conjunction with the spectral collocation method are used to solve aforementioned models numerically. The VO operator of differentiation will be used in the Caputo sense. The proposed technique simplifies these models by reducing them to systems of algebraic equations that are easy to solve. To determine the effectiveness and accuracy of the sugested technique, the absolute errors and maximum absolute errors for four realistic models are studied and illustrated by several tables and graphs at different values of the VO and the SFG parameters;
α
and
μ.
Also numerical comparisons between the suggested technique with other numerical methods in the existing literature are held. The numerical results confirm that the suggested technique is accurate, computationally efficient and easy to implement.
Introduction
In recent years, many researchers have been interested in the study of initial value problems of some special second-order nonlinear singular ordinary differential equations (ODEs). One of these most important equations is the Emden-Fowler (EF) equation, which has a variety of applications in various fields. Robert Emden [
1
] and Fowler [
2
] have studied the EF equation in detail. The general mathematical model of EF is expressed as
y″(t)+rty′(t)+γg(t)f(y(t))=0,r≥1,
with the initial conditions (ICs)
y(0)=θ1,y′(0)=θ2.
where
g
(
t
) and
f
(
y
(
t
)) are functions of
t
and
y
(
t
), respectively,
r,γ,θ1
, and
θ2
are constants. Solving differential equations with singularity behavior is a defy. In particular, the current EF equation problem, which involves a singularity at
t=0
, is crucial in practical applications. Analytically, solving this equation is difficult, so to solve this problem various numerical techniques were used [
3
,
4
]. For
g(t)=1,
and
γ=1,
Eq. (
1.1
) moderates to the Lane-Emden (LE) equation in its usual form. This equation has appeared in the work of Homer Lane and Robert Emden work. The inner construction of polytrophic stars, the gas cloud model, cluster galaxies, and radiative cooling are all represented by these models. No one can deny the usefulness and importance of these models, which have a wide range of applications in mathematical physics, astrophysics, and celestial mechanics such as isotropic continuous media, isothermal gas sphere models [
5
], dusty fluid models [
6
], oscillating magnetic systems [
7
], catalytic diffusion reactions [
8
], stellar structure models [
9
], electromagnetic theory [
10
], gaseous star density [
11
], and classical quantum mechanics [
12
].
The nonlinear singular ODE (
1.1
) was first developed by Jonathan Homer Lane [
13
], a U.S. astrophysicist who was interested in estimating both the temperature and the density of mass on the surface and after 37 years was investigated in greater depth by Emden [
1
]. This ODE in astrophysics is a dimensionless version of Poisson’s equation for a simple stellar model’s gravitational potential [
14
]. Go back to [
15
,
16
] to get a sufficient explanation of deriving of the classic LE equation for modeling the thermal dynamics of a spherical cloud, which takes the following form:
y″(t)+2ty′(t)+ym(t)=0,t>0,0, \end{aligned}$$\end{document}]]>
with the ICs
y(0)=1,y′(0)=0.
In quantum mechanics and astrophysics, the values of m physically lie in the interval [0,5]. In different literatures, the exact solutions of problem (
1.3
–
1.4
) have only been found at
m=0,1,5
[
16
], and for other values of m, several methods have been successfully used to approximate the solution of this problem [
17
,
18
].
One of the most essential significant classes of the differential equations is the delay differential (DD) equations. Due to their contributions significantly to solve the real-life problems. These equations have numerous applications in many biological models, as well as scientific phenomena such as communication system models, dynamical population models, economic systems, engineering systems, and transport models, so they have attracted a lot of attention from the research community [
19
,
20
], and various numerical/analytical techniques are introduced to solve them [
21
,
22
].
Fractional calculus is a new topic that appears as a generalization of integer calculus, in which the order will be non-integer (arbitrary). The dynamic behavior of various phenomena in engineering sciences, physics, and other disciplines of research can be explored more precisely using non-integer order differential equations. So, there are a growing number of applications of fractional-order calculus in different fields as e.g. Long electrical lines, electrochemical processes, dielectric polarization, colored noise, viscoelastic materials, chaos, control theory, and in many other areas. Recently, the VO differential equation in which the order is a function of space and/or time have been demonstrated their ability to accurately describe real world phenomena [
23
]. In the VO model, the next state is influenced by all of its former states, not only upon its current state as in the case of integer-order derivatives. So, the VO derivative is more realistic and generalized than the standard and fractional ones. VO differential equations have gotten a lot of attention lately as a result of their suitability modeling a wide range of phenomena across many fields of science and engineering, such as anomalous diffusion [
24
], viscoelastic mechanics [
25
], control systems [
26
], petroleum engineering [
27
], and many other branches of physics and engineering [
28
].
With this respect, the importance of the VO-EF-DDEs are quite clear, due to their massive applications in mathematical physics and engineering problems. For this significance, obtaining the numerical solutions of this type of equations are required since the exact solutions are usually not easy to claim, because of the kernel of the VO operators have a changing exponent. So, the main objective of the present work is to develop an efficient and accurate numerical technique for solving the following new classes of VO-EF-DDEs:
A class of variable order Emden-Fowler delay differential equations (VO-EF-DDEs) of the following form:
Dtν(t)y(at+b)+λty′(ct+d)+g(t)y(et+f)=F(t),
with the ICs
y(0)=γ1,y′(0)=γ2,
where
λ≥1,ν(t),g(t),
and
F
(
t
) are given functions, the parameters
a,b,c,d,e,f,γ1,
and
γ2
are known constants also
Dtν(t)
denotes the Caputo fractional derivative of variable order,
1<ν(t)≤2
and
y
(
t
) is an unknown function. By taking
g(t)=1
the EF model reduces to the Lane-Emden model. The description of problem (
1.5
) in the VO operator
ν(t)
has many advantages one of them; it describes the behavior of the system and its properties in a more accurate manner. Also, it’s important to mention that the model (
1.5
–
1.6
) covers the classical and the constant fractional models when
ν(t)
is constant.
A system of variable order Emden-Fowler delay of differential equations (VO-EF-DDEs) of the following form
Dtν(t)y1(at+b)+λ1ty1′(ct+d)+g1(t)y1(et+f)=F1(t),Dtν(t)y2(at+b)+λ2ty2′(ct+d)+g2(t)y2(et+f)=F2(t),
with the following ICs
y1(0)=β1,y1′(0)=β2,y2(0)=β3,y2′(0)=β4,
where
λ1,λ2≥1,a,b,c,d,e,f,β1,β2,β3
and
β4
are known constants, the functions
ν(t),g1(t),g2(t),F1(t)
and
F2(t)
are assumed to be known functions, and also
Dtν(t)
is the Caputo fractional derivative of order
1<ν(x)≤2.
By taking
g1(t)=g2(t)=1
, the EF model reduces to the Lane-Emden model.
The suggested technique is based on the spectral collocation method with the OMs of variable/integer order differentiation of the SFGPs. To the best of our knowledge, the numerical treatments of VO-EF-DDEs and the system of VO-EF-DDEs have not been established by using SFGPs yet.
Spectral collocation methods are known to be one of the powerful numerical techniques for solving different kinds of differential equations due to their applicability in bounded and unbounded domains. The convergence speed is one of the major advantages of the spectral method. In general, spectral methods are promising candidates for solving fractional differential equations since their global nature fits well with the non-local definition of fractional operators [
29
,
30
]. Spectral collocation methods have been widely used for solving problems with smooth solutions because of their efficiency, exponential rate of convergence, low computational cost, and high order of accuracy achieved using a minimal number of grid points. It is known that the accuracy of spectral methods increases with an increase in grid points but beyond a certain number of grid points, the accuracy rapidly deteriorates. Also, the limitation of spectral methods is that their accuracy deteriorates for complicated domains. OMs were restructured for solving numerous types of fractional differential equations (FDEs). The procedure of numerical methods in combining with OMs of some orthogonal polynomials for solving FDEs created extremely precise solutions for such equations. The essential step in the proposed technique is to represent the numerical solution of the problem as a finite sum using SFGPs as basis functions, then the OMs corresponding to VO/ integer order derivatives are derived. By using these OM in Eq. (
1.5
), and applying the collocation method which requires that the residual of VO-EF-DDEs is vanishing at
N
points that correspond to the Gauss quadrature nodes. When the obtained system of algebraic equations is completed with the initial conditions,
(N+1)
algebraic equations are produced, which can be solved by utilizing any appropriate iterative technique.
The SFGPs have a number of characteristics that make them an ideal candidate for our suggested technique, like:
Due to the unique behavior of fractional differential and/or integral equations, the direct use of spectral methods using traditional orthogonal polynomials such as Legendre, Chebyshev, and Jacobi has low convergence rates. To remedy this difficulty, we employ SFGPs (see [
31
]), in which the variable t in the SCPs is replaced by
tμ,0<μ≤1,t∈[0,L].
One of the benefits of fractional-order polynomials is in dealing with problems with smooth solutions, where any choice of the fractional factor
μ
yields a highly accurate approximate solution. The fractional-order polynomials can lessen the loss of the order of convergence of the non-smoothness problems by suitably chosen of the factor
μ
.
For
μ=1,
the SFGPS reduced to the SGPs. Numerous studies have shown that the Gegenbauer polynomials (GPs) are very effective in solving a wide range of issues [
32
,
33
–
34
].
The parameter
α>-0.5-0.5$$\end{document}]]>
distributes the Gegenbauer polynomials (GPs). Every time this parameter was altered, a new polynomial was formed. As a result, they involve a limitless number of orthogonal polynomials like the Chebyshev polynomials (CPs) of the first kind with the parameter
α=0,
the CPs of the second kind with the parameter
α=1,
and the shifted Legendre polynomials (LPs) with the parameter
α=12
. The same results are assumed for the SFGPs.
The main advantages of the proposed scheme are the low-cost computing, small CPU time, and the simple implementation. Also the present method has the ability to convert the given problem into a system of algebraic equations, which is easy to solve. Many numerical experiments are presented to measure the quality of the presented scheme. The obtained numerical results are compared with other methods in the existing literature.
The following is the contents of this paper: The preliminaries of VO-fractional calculus and SFGPs are presented in Sect.
2
. The various types of shifted fractional Gegenbauer operational matrices (SFGOMs) are deduced in Sect.
3
. In Sect.
4
, the methodology of the proposed technique is explained. In Sect.
5
, some applications of the proposed algorithm to some numerical problems are offered to demonstrate the effectiveness of the proposed algorithm. At last, a conclusion is given in Sect.
6
.
The Caputo fractional derivative of order
ν(t)
for a continuously differentiable function
y:[0,∞]→R
is defined as
Dtν(t)y(t)=1Γ(m-ν(t))∫0t(t-ξ)m-ν(t)-1y(m)(ξ)dξ,
where
ν(t)
is a positive continuous bounded function and
m-1<ν(t)≤m.
Because of the ICs for fractional differential equations with Caputo derivatives are the same as for integer order differential equations; the Caputo type definition is particularly useful in a wide range of applications. The results stated in definition (
2.1
) have the following characteristics:
Shifted fractional Gegenbauer polynomials and some of their properties
Gegenbauer polynomials
The Gegenbauer polynomials
Cj(α)(x),
of degree
j∈Z+,
and associated with the parameter
α>-12\frac{-1}{2}$$\end{document}]]>
are a sequence of real polynomials in the finite domain
[-1,1].
They are a family of orthogonal polynomials which has many applications [
32
].
A suitable standardization of the Gegenbauer polynomials
Cj(α)(x)
of degree j and associated with the parameter
α>-12 \frac{-1}{2}$$\end{document}]]>
dates back to Doha [
34
], where the Gegenbauer polynomials can be represented by
Cj(α)(x)=j!Γα+12Γj+α+12Pjα-12,α-12(x),j=0,1,2,…,
or equivalently
Cj(α)(1)=1,j=0,1,2,⋯,
where
Pj(α-12,α-12)(x)
is the Jacobi polynomial of degree j.
As a result of this standardized
Cj(0)(x)
becomes identical with the Chebyshev polynomials (CPs) of the first kind
Tj(x).Cj(12)(x)
is the Legendre polynomials (LPs)
Lj(x).Cj(1)(x)
is equal to
1(j+1)Uj(x),
where
Uj(x)
is the CPs of the second type [
34
].
The analytical form of the GPs is given by,
Cj(α)(x)=∑k=0j(-1)j-kj!Γ(α+12)Γ(j+k+2α)Γ(k+α+12)Γ(j+2α)(j-k)!k!x+12k.
The Gegenbauer polynomials can be generated using the following useful recurrence equation
(j+2α)Cj+1(α)(x)=2(j+α)xCj(α)(x)-jCj-1(α)(x),j≥1.
The Gegenbauer polynomials satisfy the following orthogonality relation
Ci(α)(x),Cj(α)(x)=∫-11Ci(α)(x)Cj(α)(x)ω(α)(x)dx=λj(α)δi,j,=\int _{-1}^{1}C_{i}^{(\alpha )}(x)C_{j}^{(\alpha )}(x)\omega ^{(\alpha )}(x)dx=\lambda _{j}^{(\alpha )}\delta _{i,j}, \end{aligned}$$\end{document}]]>
where
ω(α)(t)
is the weight function, and it is an even function given by the relation
ω(α)(x)=1-x2α-12,
and
λj(α)=‖Cj(α)(x)‖2=22α-1j!Γ2α+12(j+α)Γj+2α,
is the normalization factor, and
δi,j
is the Kronecker delta function.
We denote the zeros of the Gegenbauer polynomials (Gegenbauer-Gauss nodes)
CN+1(α)(x)
by
xN,j(α),j=0,⋯,N.
Shifted Gegenbauer polynomials
For using Gegenbauer polynomials in the interval [0,
L
], the shifted Gegenbauer polynomials (SGPs) are formed by replacing the variable
x
with
z=2xL-1,0≤z≤L.
So, we can write SGPs as
CS,j(α)(z)=Cj(α)2xL-1.
The explicit analytical form of the SGPs is given as
CS,j(α)(z)=∑k=0j(-1)j-kj!Γα+12Γj+k+2αΓk+α+12Γj+2α(j-k)!k!Lkzk,
These polynomials recover the shifted CPs of the first kind
TS,j(z)≡CS,j(0)(z),
the shifted LPs
LS,j(z)≡CS,j(12)(z),
and the shifted CPs of the second kind
CS,j(1)(z)≡1j+1US,j(z).
The orthogonal relation of SGPs is getting from
CS,i(α)(z),CS,j(α)(z)=∫0LCS,i(α)(z)CS,j(α)(z)ωS(α)(z)dz=λS,j(α)δi,j,=\int _{0}^{L}C_{S,i}^{(\alpha )}(z)C_{S,j}^{(\alpha )}(z)\omega _{S}^{(\alpha )}(z)dz=\lambda _{S,j}^{(\alpha )}\delta _{i,j}, \end{aligned}$$\end{document}]]>
where
ωS(α)(z)
is the weight function, it is an even function given by the relation
ωS(α)(z)=zL-z2α-12,
and
λS,j(α)=L22αλj(α).
We denote the zeros of the shifted Gegenbauer polynomial (shifted Gegenbauer-Gauss polynomial)
CS,N+1(α)(z)
by
zS,N,j(α),j=0,…,N.zS,N,j(α)=L2xN,j(α)+1,j=0,…,N.
The
q
-th derivative of
CS,i(α)(z),
is obtained by
DqCS,i(α)(z)=∑k=0i-q(-1)i-q-ki!Γα+12Γk+i+2α+qΓk+α+q+12Γi+2α(i-q-k)!k!Lq+kzk,DCS,i(α)(0)=(-1)i-1i!Γα+12Γi+2α+1Γα+32Γi+2α(i-1)!Lk+1,D2CS,i(α)(0)=(-1)i-2i!Γα+12Γi+2α+2Γα+52Γi+2α(i-2)!Lk+2.
Shifted fractional Gegenbauer polynomials
The SFGPs can be defined for
t∈[0,L]
by introducing the change of variable
z=tμ
and
μ>00$$\end{document}]]>
on SGPs. Let the SFGPs be denoted by
CS,i(α,μ)(t).
The explicit analytic form of SFGPs is given as
CS,j(α,μ)(t)=∑k=0j(-1)j-kj!Γα+12Γj+k+2αΓk+α+12Γ(j+2α)(j-k)!k!Lktμk,=∑k=0jΩj,ktμk.
Also, the orthogonal relation of SFGPs
∫0LCS,i(α,μ)(t)CS,j(α,μ)(t)ωS(α,μ)(t)dt=λS,j(α)δi,j,
where
ωS(α,μ)(t)=μtμα+12-1(1-tμ)α-12,
and
λS,j(α)
are defined by Eq. (
2.7
).
The SFGPs comprise unlimited number of orthogonal polynomials, among them the shifted fractional-order CPs of the first kind
TS,j(μ)(t)≡CS,j(0,μ)(t),
the shifted fractional-order CPs of the second kind
US,j(μ)(t)≡(j+1)CS,j(1,μ)(t),
and the shifted fractional-order LPs
LS,j(μ)(t)≡CS,j(12,μ)(t).
We denote the zeros of the SFGPs (shifted fractional Gegenbauer-Gauss polynomial)
CS,N+1(α,μ)(t)
by
tS,N,j(μ),j=0,…,N.tS,N,j(μ)=tN,j+121μ.
The fractional Gegenbauer-Gauss quadrature rule. Due to the property of the standard Gegenbauer-Gauss quadrature, it follows that for any
f((t+12)1μ)∈P2N+1(0,L)∫0LwS(α,μ)(t)f(t)dt=L22α∫-11w(α)(t)ft+121μdt,=L22α∑j=0NϖN,j(α)ftN,j+121μ,=∑j=0NϖS,N,j(α,μ)ftS,N,j(μ),
where
ϖS,N,j(α,μ)
are the Christoffel numbers of the SFGPs that are given by
ϖS,N,j(α,μ)=L22αϖN,j(α),j=0,…,N.
The square integrable function
y(t)∈[0,L]
can be approximated by SFGPs as
y(t)=∑j=0∞y~S,jCS,j(α,μ)(t),
where the coefficients
y~S,j
are computed by
y~S,j=(λS,j(α,μ))-1∫0Ly(t)ωS(α,μ)(t)CS,j(α,μ)(t)dt,j=0,1,⋯.
If we approximate
y
(
t
) by the first (N+1)-terms as follows
y(t)=∑j=0Ny~S,jCS,j(α,μ)(t).
The approximation of the function
y
(
t
) can be written in the vector form as
y(t)=YTΦμ(t),
where
YT=[Y0,Y1,...,YN]
is the shifted fractional Gegenbauer (SFG) coefficient vector, and
Φμ(t)=CS,0(α,μ)(t),CS,1(α,μ)(t),…,CS,N(α,μ)(t)T=ATN,μ(t),
is the SFG vector and
A
is a lower triangular matrix of order
(N+1)×(N+1)A=a~S,0,0(α)0⋯0a~S,1,0(α)a~S,1,1(α)⋯0⋮⋮⋱⋮a~S,N,0(α)a~S,N,1(α)⋯a~S,N,N(α),
where
a~S,j,k(α)=(-1)j-kj!Γα+12Γj+k+2αΓk+α+12Γ(j+2α)(j-k)!k!Lk,
and
TN,μ(t)=1tμt2μ⋮tNμ.
Error estimation and convergence
Let
CS,N(α,μ)=Span{CS,0(α,μ)(t),CS,1(α,μ)(t),⋯,CS,N(α,μ)(t)}
and
y
(
t
) be an arbitrary element in
LωS(α,μ)2[0,L].
Since
CS,N(α,μ)
is a finite dimensional vector space,
y
(
t
) has the unique best approximation out of
CS,N(α,μ)
like
yN(t)∈CS,N(α,μ),
such that
∀g(t)∈CS,N(α,μ),‖y(t)-yN(t)‖2≤‖y(t)-g(t)‖2,
where
‖y(t)‖=<y(t),y(t)>.}.$$\end{document}]]>
Theorem 2.1
Suppose that the functionDkμy(t)∈C[0,1]
,
fork=0,1,…,N-1,Re(μ)>0,α>-12, 0,}\quad {\alpha >\frac{-1}{2},}$$\end{document}]]>Re(α+2N+52)μ>0.0.}$$\end{document}]]>IfyN(t)is the best approximation toy
(
t
) from
CS,N(α,μ),then the error bound is presented as‖y(t)-yN(t)‖2≤HμΓ(N+1)μ+1Γα+12Γα+2N+52Γ2α+2N+3,where‖D(N+1)μy(t)‖≤Hμ,t∈[0,1].
Proof
Consider the generalized Taylor formula [
36
]
y(t)=∑i=0NtiμΓ(iμ+1)Diμy(0)+D(N+1)μy(ξ)Γ((N+1)μ+1)t(N+1)μ,
with
ξ∈[0,t],∀t∈[0,1].
Let
PN(t)=∑i=0NtiμΓ(iμ+1)Diμy(0).
Then,
|y(t)-PN(t)|=|D(N+1)μy(ξ)Γ((N+1)μ+1)t(N+1)μ|.
Since
yN(t)
is the best square approximation function of
y
(
t
), then
‖y(t)-yN(t)‖22≤‖y(t)-PN(t)‖22,=∫01ωS(α,μ)(t)(y(t)-pN(t))2dt,=∫01ωS(α,μ)(t)(D(N+1)μy(ξ)Γ((N+1)μ+1)t(N+1)μ)2dt,=D(N+1)μy(ξ)2Γ2(N+1)μ+1∫01μtμα+12-11-tμα-12t2(N+1)μdt,≤Hμ2Γ2(N+1)μ+1Γα+12Γ(α+2N+52)Γ(2α+2N+3).
By taking the square roots, the proof is completed.
□
SFGOMs of VO fractional differentiation
Theorem 3.1
The first order derivative of the SFG vectorΦμ(at+b)is approximated asddtΦμ(at+b)≃ΨΦ1(t),where the OMΨis a square matrix of order(N+1)×(N+1)that takes the following form
,
Ψ≃AΥNΘNA-1,whereAis defined in
(
2.16
),
ΥNandΘNare(N+1)×(N+1)matrices, and their elementsυij,andθij,0≤i,j≤Nare given, respectively, as followsυij=μijajbμi-j,i,j=0,1,…,N,θij=j+1,fori=j+1,j=0,1,…,N,0,o.w..
Proof
From Eq. (
2.15
), we can write
ddtΦμ(at+b)=ddt[CS,0(α,μ)(at+b),CS,1(α,μ)(at+b),…,CS,N(α,μ)(at+b)]T=ddt(ATN,μ(at+b)),=AddtTN,μ(at+b).
Since
ddtTN,μ(at+b)=ddt1(at+b)μ(at+b)2μ⋮(at+b)Nμ.
By using Newton’s generalized binomial theorem, we can write
(at+b)iμ=∑k=0∞μikakbμi-ktk,=∑k=0∞(iμ)(iμ-1)⋯(iμ-k+1)k!akbμi-ktk,i=0,1,…,N,
this series converges for
iμ∈Z+
or
|atb|<1.□
By approximating the series (
3.4
) by using the first (N+1) terms, the relation (
3.3
) can be written as
ddtTN,μ(at+b)=ddt1(at+b)μ(at+b)2μ⋮(at+b)Nμ=ΥN,μddt1tt2⋮tN≃ΥN,μ012t⋮NtN-1≃ΥN,μΘNTN,1(t),
where
ΥN,μ,ΘN
are
(N+1)×(N+1)
matrices, and their elements have the following forms, respectively,
υij=μijajbμi-j,i,j=0,1,…,N,θij=j+1,fori=j+1,j=0,1,…,N,0,o.w..
By substituting Eq. (
3.5
) into Eq. (
3.2
), the following relation is obtained
ddtΦμ(at+b)≃AΥN,μΘNTN,1(t).
Since
A
is invertible, then
ddtΦμ(at+b)≃AΥN,μΘNA-1Φ1(t)≃ΨΦ1(t),
where
Ψ
is a square matrix, and this finishes the proof.
Theorem 3.2
The VO Caputo fractional derivative of the SFG vectorΦμ(at+b)can be approximated asDtν(t)Φμ(at+b)≃ΞΦ1(t),ν(t)>0,0, \end{aligned}$$\end{document}]]>where the OMΞis a square matrix of order(N+1)×(N+1)takes the following form
,
Ξ≃AΥNΛNA-1,whereAis defined in
(
2.16
),
ΥNandΛNare(N+1)×(N+1)matrices, and their elementsυij,andλij,0≤i,j≤N
,
respectively, are given as followsυij=μijajbμi-j,i,j=0,1,…,N,λij=Γ(i+1)Γ(i+1-ν(t))t-ν(t),fori=j,j=n,n+1,…,N,0,o.w..
Proof
From Eq. (
2.15
), we can write
Dtν(t)Φμ(at+b)=Dtν(t)[CS,0(α,μ)(at+b),CS,1(α,μ)(at+b),…,CS,N(α,μ)(at+b)]T=Dtν(t)(ATN,μ(at+b)),=ADtν(t)TN,μ(at+b).
Since
Dtν(t)TN,μ(at+b)=Dtν(t)1(at+b)μ(at+b)2μ⋮(at+b)Nμ.
By using Newton
′
s generalized binomial theorem, we can write
(at+b)iμ=∑k=0∞μikakbμi-ktk,=∑k=0∞(iμ)(iμ-1)⋯(iμ-k+1)k!akbμi-ktk,i=0,1,…,N,
this series converges for
iμ∈Z+
or
|atb|<1.
By approximating the series (
3.10
) by using the first (N+1) terms, the relation (
3.9
) can be written as
Dtν(t)TN,μ(at+b)=Dtν(t)1(at+b)μ(at+b)2μ⋮(at+b)Nμ≃ΥN,μDtν(t)1tt2⋮tN.
By using the definition (2.1) of the VO fractional derivative in the Caputo sense, then Eq. (
3.11
) can be obtained as the below form:
Dtν(t)TN,μ(at+b)≃ΥN,μ0⋮Γ(n+1)Γ(n+1-ν(t))tn-ν(t)Γ(n+2)Γ(n+2-ν(t))tn+1-ν(t)⋮Γ(N+1)Γ(N+1-ν(t))tN-ν(t)≃ΥN,μΛNTN,1(t),
where
ΥN,μ,ΘN
are
(N+1)×(N+1)
matrices, and their elements have the following forms, respectively,
υij=μijajbμi-j,i,j=0,1,…,N,λij=Γ(i+1)Γ(i+1-ν(t))t-ν(t),fori=j,j=n,n+1,…,N,0,o.w..
By substituting from Eq. (
3.12
) into Eq. (
3.8
), the following relation is obtained
Dtν(t)Φμ(at+b)≃AΥN,μΛNTN,1(t).
Since
A
is invertible, then
Dtν(t)Φμ(at+b)≃AΥN,μΛNA-1Φ(t)≃ΞΦ1(t),
where
Ξ
is a square matrix, and this finishes the proof.
□
The methodology
In this section, we will discuss how to use the SFG-OMs to solve the problem Eqs. (
1.5
–
1.6
) numerically. As a result, the problem will be converted into an algebraic system of equations that can be solved easily by any iteration method. By using the properties of the SFGPs which are mentioned in the definition (2.2), the solution can be approximated as
yN(at+b)=∑j=0Ny~S,jCS,j(α,μ)(at+b)=YTΦμ(at+b),
by the same way, the terms
y(ct+d),y(et+f)
are approximated.
Moreover, the Caputo fractional derivative of order
ν(t)
of
yN(at+b)
is estimated as
Dt(ν(t))yN(at+b)=Dt(ν(t))(∑j=0Ny~S,jCS,j(α,μ)(at+b))=YTDt(ν(t))(Φμ(at+b)).
By using Theorem (3.2)
Dt(ν(t))yN(at+b)≃YTΞΦ1(t),
and by using Theorem (3.1)
Dt(1)yN(ct+d)≃YTΨΦ1(t).
Substitute from Eqs. (
4.1
)–(
4.3
) into Eqs. (
1.5
–
1.6
), we get
YTΞΦ1(t)+λtYTΨΦ1(t)+g(t)YTΦμ(t)=F(t),
and the ICs will be
YTΦμ(0)=γ1,YTΨΦμ(0)=γ2.
Suppose the nodes
tk
which are the roots of the SFGPs. By substituting these nodes in Eqs. (
4.4
)–(
4.5
); therefore the collocation scheme can be written as
YTΞΦ1(tk)+λtkYTΨΦ1(tk)+g(tk)YTΦμ(tk)=F(tk),
with the ICs
YTΦμ(0)=γ1,YTΨΦμ(0)=γ2.
This yields a set of (N+1) algebraic equations with the required SFG coefficients
y~S,j,j=0,1,…,N,
which may be solved using an appropriate iterative procedure. As a result, the approximate solution
yN(t)
of the system (
1.7
–
1.8
) can be obtained.
Numerical applications
In this section, we’ll introduce four applications to evaluate the applicability of our approach. By using a Dell laptop with an Intel(R) Core(TM) i3 CPU M 370@ 2.40 GHz and 3.00GB RAM configuration, the software was written in Mathematica version 12. The results obtained by our approach are compared to those obtained in literature using other numerical methods. The results constructed in this paper are measured by means of:
The absolute errors (AEs) given by
E(ti)=|y(ti)-yN(ti)|,1≤i≤N,
where
y(ti)
and
yN(ti)
are the exact solution (ES) and the approximate solution (AS), respectively.
The maximum absolute errors (MAEs)
L∞=max1≤i≤N{E(ti):∀ti∈[0,L]}.
The order of convergence (OC)
OC=log(error(N1)error(N2))log(N2N1),
where error(N) denotes the error corresponding to polynomial degree N.
Applications of EF-DDEs-VO
Application (1)
Consider the following VO-EF-DDEs [
37
],
Dtν(t)y(2t-1)+2tddty(3t)+ty(t+1)=4e2t-1+6te3t+tet+1,
with the ICs
y(0)=1,ddty(0)=1,
where the ES at
ν(t)=2
is
y(t)=et.
Fig. 1
The numerical results of SFGOM technique on application (1) at
N=14,μ=1
and
α=0.5.a
The AEs graph of
y
(
t
).
b
The graphs of the ASs at different values of
ν(t).
In Figure
1
a, the AE function with
N=14,μ=1,α=0.5
displays highly accuracy, and the numerical outcomes at different values of VOs are shown in Fig.
1
b.
Application (2)
Consider the following VO-EF-DDEs- [
37
],
Dtν(t)y(3t-1)+2tddty(2t)+ty(t+1)=F(t),1<ν(t)≤2,
with the ICs
y(0)=1,ddty(0)=0,
where
F(t)=27Γ(4)t3-ν(t)Γ(4-ν(t))-27Γ(3)t2-ν(t)Γ(3-ν(t))+48t+t(1+(t+1)3),
and the ES is
y(t)=t3+1.
Fig. 2
The numerical results of SFGOM technique on application (2) at
N=3,μ=1,α=0.5
and
ν(t)=1+t2.a
The AEs graph of
y
(
t
).
b
The graph of ES and AS in [0,1]
Table 1
The MAEs of application (2) at
α=0.5,ν(t)=1+t2
and different values of
N,μ
and CPU time (second) at
μ=1
N
μ=1
CPU time μ=1
OC (μ=1)
N
μ=12
N
μ=13
2
8.68717×10-1
0.0312002
5
8.31867×10-1
8
5.82913×10-2
3
8.88178×10-16
0.0624004
85.1285
6
3.84462×10-14
9
3.22044×10-12
Table 2
The AEs of application (2) at
N=3,ν(t)=1+t2
and different values of
α
t
α=-0.3
α=0
α=0.5
α=1
0.1
2.22045×10-16
0
8.88178×10-16
4.44089×10-16
0.3
2.22045×10-16
0
4.44089×10-16
4.44089×10-16
0.5
2.22045×10-16
2.22045×10-16
4.44089×10-16
2.22045×10-16
0.7
4.44089×10-16
4.44089×10-16
2.22045×10-16
4.44089×10-16
0.9
6.66134×10-16
8.88178×10-16
2.22045×10-16
0
The numerical results of application (2) are given in Tables
1
and
2
and are plotted in Fig.
2
a–b. Table
1
introduces the computational time at
μ=1,
the MAEs at different values of
N,μ,
and the order of convergence at
μ=1.
Additionally, Table
2
presents the AEs at
N=3,μ=1,
and different values of
α.
In Fig.
2
a the AEs graph is displayed, while Fig.
2
b depicts the agreement between the ES and the AS graphs of
y
(
t
) for
N=3,μ=1,α=0.5
at
ν(t)=1+t2.
Application (3)
Consider the following VO-EF-DDEs [
38
],
12Dtν(t)yt2+3tddtyt2+y2(t)=F(t),1<ν(t)≤2,
with the ICs
y(0)=1,ddty(0)=0,
where
F(t)=Γ(5)t4-ν(t)32Γ(5-ν(t))+3t24+(1+t4)2,
and the ES is
y(t)=t4+1.
Fig. 3
The numerical results of SFGOM technique on application (3) at
N=4,μ=1,α=0.5,
and
ν(t)=1+sin(t)2.a
The AEs graph of
y
(
t
).
b
The graph of ES and AS in [0,1]
Table 3
The MAEs of application (3) at
α=0.5,ν(t)=1+sin(t)2
and different values of
N,μ
and the CPU time (second) at
μ=1
N
μ=1
CPU Time μ=1
OC (μ=1)
N
μ=12
N
μ=13
2
2.74192×10-1
0.0312002
4
1.13406×10-1
8
3.48657×10-3
3
9.92054×10-3
0.0424004
8.18621
6
3.96059×10-3
10
1.55323×10-4
4
2.22045×10-16
0.0624004
109.254
8
6.61693×10-14
12
2.35434×10-10
Table 4
The AEs of application (3) at
N=4,ν(t)=1+sin(t)2
and different values of
α
t
α=-0.3
α=0
α=0.5
α=1
0.1
0
0
2.22045×10-16
0
0.3
0
0
0
2.22045×10-16
0.5
0
0
0
0
0.7
0
2.22045×10-16
0
8.88178×10-16
0.9
0
4.44089×10-16
6.66134×10-16
2.66454×10-15
Table 5
Comparison of the AEs for
μ=1,ν(t)=2
of application (3)
t
N=4
N=6
Our method
Ref. [38]
Our method
Ref. [38]
0.1
0
5.1070×10-14
0
2.2204×10-15
0.3
0
1.6875×10-13
2.22045×10-16
3.9968×10-15
0.5
2.22045×10-16
2.6268×10-13
2.22045×10-16
7.3275×10-15
0.7
4.44089×10-16
1.2870×10-12
6.66134×10-16
1.3411×10-13
0.9
4.44089×10-16
5.3457×10-12
1.33227×10-15
1.2814×10-12
The numerical results of application (3) are given in Tables
3
,
4
, and
5
and are plotted in Fig.
3
a–b. Table
3
introduces the computational time at
μ=1,
the MAEs at
α=0.5,
and different values of
N,μ,
and the order of convergence at
μ=1.
Additionally, Table
4
presents the AEs at
N=4,μ=1,
and different values of
α.
Table
5
displays the efficiency and effectiveness of our method by comparing it to the method in Ref. [
38
] at
ν(t)=2
. In Fig.
3
a we describe the AEs graph, while Fig.
3
b plots the ES and the AS graphs of
y
(
t
) for
N=4,μ=1,α=0.5
at
ν(t)=1+sin(t)2.
Application (4)
Consider the following system of EF-DDEs-VO [
39
],
Dtν(t)y1(2t-1)+3tddty1(3t)+y2(t)y1(t+1)=F1(t),Dtν(t)y2(2t-1)+2tddty2(3t)+y1(t)y2(t+1)=F2(t),1<ν(t)≤2,
with the initial conditions,
y1(1)=2,ddty1(0)=0,y2(1)=0,ddty2(0)=0,
where
F1(t)=8Γ(4)t3-ν(t)Γ(4-ν(t))-12Γ(3)t2-ν(t)Γ(3-ν(t))+243t+(1-t3)(1+(t+1)3),F2(t)=12Γ(3)t2-ν(t)Γ(3-ν(t))-8Γ(4)t3-ν(t)Γ(4-ν(t))-162t+(1+t3)(1-(t+1)3),
and the ESs are
y1(t)=1+t3,y2(t)=1-t3.
Fig. 4
The numerical behavior of the proposed method on application (4).
a
shows the AEs of
y1(t)
at
N=3
and
α=0.5,ν(t)=1+t2.b
shows the graph of the ES and AS
y1(t)
on [0,1].
c
shows the AEs of
y2(t)
at
N=3
and
α=0.5,ν(t)=1+t2.
(
d
) shows the graph of the ES and AS
y2(t)
on [0,1]
Table 6
MAEs of application (4) for
α=0.5,ν(t)=1+t2
and different values of
N,μ
and CPU time (seconds) at
μ=1
N
μ=12
N
μ=13
y1(t)
y2(t)
y1(t)
y2(t)
4
7.60719×10-1
9.06883×10-1
7
1.64448×10-1
2.37767×10-1
5
9.07018×10-2
1.32008×10-1
8
3.35157×10-2
4.10954×10-2
6
4.04156×10-15
4.52355×10-15
9
2.10661×10-12
1.97204×10-12
N
μ=1
CPU time (μ=1)
y1(t)
y2(t)
3
2.22045×10-16
6.66134×10-16
0.0624004
Table 7
The AEs of application (4) at
N=3
and different values of
α
ν(t)=1+t2
t
α=-0.3
α=0
y1(t)
y2(t)
y1(t)
y2(t)
0.1
4.44089×10-16
9.99201×10-16
2.22045×10-16
4.44089×10-16
0.3
6.66134×10-16
7.77156×10-16
2.22045×10-16
1.11022×10-16
0.5
8.88178×10-16
4.44089×10-16
0
6.66134×10-16
0.7
8.88178×10-16
4.44089×10-16
0
7.77156×10-16
0.9
1.11022×10-15
4.44089×10-16
4.44089×10-16
8.88178×10-16
t
α=0.5
α=1
y1(t)
y2(t)
y1(t)
y2(t)
0.1
0
4.44089×10-16
6.66134×10-16
5.55112×10-16
0.3
2.22045×10-16
0
6.66134×10-16
5.55112×10-16
0.5
0
3.33067×10-16
4.44089×10-16
3.33067×10-16
0.7
2.22045×10-16
5.55112×10-16
0
1.11022×10-16
0.9
2.22045×10-16
6.66134×10-16
2.22045×10-16
5.55112×10-17
Table 8
Comparison of the AEs for
N=3,μ=1,ν(t)=2
of application (4)
t
α=-0.3
Ref. [39]
y1(t)
y2(t)
y1(t)
y2(t)
0.1
1.77636×10-15
2.22045×10-16
2.22045×10-16
0
0.3
1.11022×10-15
0
2.22045×10-16
0
0.5
6.66134×10-16
2.22045×10-16
2.22045×10-16
1.11022×10-16
0.7
2.22045×10-16
5.55112×10-16
2.22045×10-16
3.33067×10-16
0.9
2.22045×10-16
7.77156×10-16
2.22045×10-16
4.44089×10-16
The numerical results of application (4) are listed in Tables
6
,
7,
and
8
and graphically illustrated in Fig.
4
a–d. We tabulate the computational time at
μ=1,
and the MAEs in Table
6
using the SFGOMs method with different choices of
μ
at different values of
N
. In addition, Table
7
lists the AEs for various values of
α.
The author of Ref. [
39
] estimated the AS of
y1(t)
and
y2(t)
, and at Table
8,
we calculated the AEs from these results and compared them with the AEs of our technique to demonstrate the effective and precise of our method. To illustrate the convergence of the proposed technique, we plot the AEs between the ES and the AS of
y1(t)
and
y2(t)
with
N=3,μ=1
and
α=0.5
in Fig.
4
a and c, respectively. The agreement between the ES and AS of
y1(t)
and
y2(t)
is depicted in Fig.
4
b and d.
Conclusion
In this article, an effective and accurate numerical technique was used to find the approximate solutions of the VO-EF-DDEs and the system of VO-EF-DDEs. The numerical technique depended on the OMs of variable and integer-order derivatives of SFGPs in conjunction with the spectral collocation method. The used technique has the ability to transform the aforementioned problems into systems of algebraic equations which are easy to solve. Numerical studies for four applications are provided to demonstrate the accuracy, efficiency, and applicability of the used technique. The obtained numerical results indicated that satisfactory results are obtained by using a few terms of the SFGPs, also the efficiency of the estimated technique is increased by using more terms of SFGPs. The SFGPs are simple to compute, as well as they have a fast convergence rate. Another benefit of SFGPs is that for every choice of the SFG parameters
α
a good numerical solution is obtained. Also the computational results of Sect.
5
showed that the value of the parameter
μ
affected on the numerical solutions as follows; when
μ=1
, good solutions are obtained by using few terms of SFGPs while for non-integer values of
μ
more terms of SFGPs are needed to obtain good results. So, we can conclude that the SGPs which are considered as a special case of SFGPs are more flexible and effective than SFGPs in the solved applications. Also, the obtained solutions do not require much CPU time. Finally, the material studied in this paper is novel, and the proposed technique can be extended to solve other fractional delay-differential equations.
Acknowledgements
The authors are very grateful to the referees, for their careful reading of the manuscript and insightful comments, which help to improve the quality of the paper.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
References
Unknown ()
Fowler (1931) Further studies of Emden’s and similar differential equations10.1093/qmath/os-2.1.259
Guirao et al. (2020) Design and numerical solutions of a novel third-order nonlinear Emden-Fowler delay differential model10.1155/2020/7359242
Mall and Chakraverty (2015) Numerical solution of nonlinear singular initial value problems of Emden-Fowler type using Chebyshev Neural Network method10.1016/j.neucom.2014.07.036
Boubaker and Van Gorder (2012) Application of the BPES to Lane-Emden equations governing polytropic and isothermal gas spheres10.1016/j.newast.2012.02.003
Flockerzi and Sundmacher (2011) On coupled Lane-Emden equations arising in dusty fluid models10.1088/1742-6596/268/1/012006
Dehghan and Shakeri (2008) Solution of an integro-differential equation arising in oscillating magnetic fields using He’s homotopy perturbation method10.2528/PIER07090403
Rach et al. (2014) Solving coupled Lane-Emden boundary value problems in catalytic diffusion reactions by the Adomian decomposition method10.1007/s10910-013-0260-6
Taghavi and Pearce (2013) A solution to the Lane-Emden equation in the theory of stellar structure utilizing the Tau method10.1002/mma.2676
Khan et al. (2015) Nature-inspired computing approach for solving non-linear singular Emden-Fowler problem arising in electromagnetic theory10.1080/09540091.2015.1092499
Luo et al. (2016) Nonlinear asymptotic stability of the Lane-Emden solutions for the viscous gaseous star problem with degenerate density dependent viscosities10.1007/s00220-016-2753-1
Ramos (2003) Linearization methods in classical and quantum mechanics 153(2) (pp. 199-208) 10.1016/S0010-4655(03)00226-1
Lane (1870) ART. IX.-on the theoretical temperature of the sun; under the hypothesis of a gaseous mass maintaining its volume by its internal heat, and depending on the laws of gases as known to terrestrial experiment 50(148)
Parand et al. (2010) An approximation algorithm for the solution of the nonlinear Lane-Emden type equations arising in astrophysics using Hermite functions collocation method10.1016/j.cpc.2010.02.018
Unknown ()
Abd-Elhameed et al. (2014) New solutions for singular Lane-Emden equations arising in astrophysics based on shifted ultraspherical operational matrices of derivatives 2(3) (pp. 171-185)
Zhao (1995) Global periodic-solutions for a differential delay system modeling a microbial population in the chemostat10.1006/jmaa.1995.1239
Niculescu (2001) Springer
Erdogan et al. (2020) A finite difference method on layer-adapted mesh for singularly perturbed delay differential equations10.2478/AMNS.2020.1.00040
Brunner et al. (2010) Discontinuous Galerkin methods for delay differential equations of pantograph type10.1137/090771922
Bhrawy and Zaky (2015) Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation10.1007/s11071-014-1854-7
Zhuang et al. (2009) Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term10.1137/080730597
Coimbra (2003) Mechanics with variable-order differential operators10.1002/andp.200310032
Kumar and Chaudhary (2017) Analysis of fractional order control system with performance and stability (pp. 408-416)
Obembe et al. (2017) Variable-order derivative time fractional diffusion model for heterogeneous porous media10.1016/j.petrol.2017.03.015
Cai et al. (2018) A survey on fractional derivative modeling of power-law frequency-dependent viscous dissipative and scattering attenuation in acoustic wave propagation10.1115/1.4040402
Ahmed and Melad (2018) New numerical approach for solving fractional differential-algebraic equations 9(2) (pp. 141-162)
Unknown ()
El-Kalaawy et al. (2018) A computationally efficient method for a class of fractional variational and optimal control problems using fractional Gegenbauer functions 70(2)
El-Gindy et al. (2018) Shifted Gegenbauer operational matrix and its applications for solving fractional differential equations10.21608/JOMES.2018.9463
Ahmed (2019) Gegenbauer collocation algorithm for solving two-dimensional time-space fractional diffusion equations10.7546/CRABS.2019.08.04
Doha (1991) The coefficients of differentiated expansions and derivatives of ultraspherical polynomials10.1016/0898-1221(91)90089-M
Unknown ()
Unknown ()
Sabir et al. (2020) Neuro-swarm intelligent computing to solve the second-order singular functional differential model10.1140/epjp/s13360-020-00440-6
Adel and Sabir (2020) Solving a new design of nonlinear second-order Lane-Emden pantograph delay differential model via Bernoulli collocation method10.1140/epjp/s13360-020-00449-x
Abdelkawy et al. (2020) Numerical investigations of a new singular second-order nonlinear coupled functional Lane-Emden model10.1515/phys-2020-0185
