# Φ-Haar Wavelet Operational Matrix Method for Fractional Relaxation-Oscillation Equations Containing Φ-Caputo Fractional Derivative

Authors: Pongsakorn Sunthrayuth; Noufe H. Aljahdaly; Amjid Ali; Rasool Shah; Ibrahim Mahariq; Ayékotan M. J. Tchalla
Journal: Journal of Function Spaces
 (2021)

---

## Abstract

This paper proposes a numerical method for solving fractional relaxation-oscillation equations. A relaxation oscillator is a type of oscillator that is based on how a physical system returns to equilibrium after being disrupted. The primary equation of relaxation and oscillation processes is the relaxation-oscillation equation. The fractional derivatives in the relaxation-oscillation equations under consideration are defined in theΦ-Caputo sense. The numerical method relies on a novel type of operational matrix method, namely, the Φ-Haar wavelet operational matrix method. The operational matrix approach has a lower computational complexity. The proposed scheme simplifies the main problem to a set of linear algebraic equations. Numerical examples demonstrate the validity and applicability of the proposed technique.


---

## Body

## 1. Introduction

The history of fractional or noninteger order differential and integral operators can be traced back to the origins of integer order calculus [1]. In recent years, fractional differential equations have attracted a lot of attention. In fields such as damping laws, diffusion processes, and other physical phenomena, fractional differential equations have proven to be adequate models. Since the majority of fractional differential equations do not have analytical solutions, we must use an approximate method. Many studies have analyzed solution techniques for fractional differential equation such as the collocation method, Adomian decomposition method, variational iteration method, tau method, and operational matrix method [2–12].A relaxation oscillator is a type of oscillator that is based on a physical system’s ability to return to equilibrium after being disrupted. The main equation of relaxation and oscillation processes is the relaxation-oscillation (R-O) equation. A relaxation equation in its standard form is given by(1)u′χ+μuχ=gχ,where μ is a real number and g is a given function. Equation (1) can be used to represent a variety of physical processes, such as the Maxwell model, which uses a spring and a dashpot in succession to explain the characteristics of a viscoelastic material. A simple physical process with a regulated phase shift is described by the standard oscillation equation. The equation that defines the oscillation u of a system corresponding to an external force g has the simplest linear form as
(2)u″χ+μuχ=gχ,where μ is the oscillator’s natural frequency. In the relaxation and oscillation models given in Equations (1) and (2), fractional derivatives are used to depict slow relaxation and damped oscillation (see [13, 14]). The fractional relaxation-oscillation differential equation (FRODE) is given by
(3)Dαuχ+μuχ=gχ,χ>0,0<α<2whereα≠1,having the following initial conditions: u0=u0, if 0<α<1, and u′0=u1, if 1<α<2, where Dα denotes a fractional differential operator of order α.The numerical investigation of FRODEs has received a lot of interest recently. The numerical solution of problem (3) (with gχ=0) was investigated in [15] by taking into account the positive fractional and fractal derivatives. The authors of [16] employed a Taylor matrix method to find the numerical solution of problem (3) by taking into account the Caputo fractional derivative. This approach is based on a fractional version of Taylor’s formula, which was first proposed in [17]. The numerical solution of problem (3) is achieved by the optimal homotopy asymptotic approach in [18], where the fractional derivative is given in the Caputo sense. In [19], a trapezoidal approximation of the fractional integral is used to get the numerical solution of problem (3) with Caputo fractional derivative. To solve problem (3), [20] proposes a generalized wavelet collocation operational matrix approach based on the Haar wavelet (HW), where the fractional derivative is represented in the Caputo sense. Inspired by the above-mentioned studies, this paper focuses on a numerical solution of the fractional differential equation of the form
(4)Dα,Φuχ+μuχ=gχ,χ∈a,b,with initial conditions dΦnua=un,n=0,1,2,⋯,m−1,wheremaxm−1,1/2<α<m, and m is a natural number. Dα,Φ is the Φ-Caputo fractional derivative of order α, and
(5)dΦnuχ=uχ,ifn=0,1Φ′ddχnuχ,ifn=1,2,3,⋯,m−1.These problems are studied in [21] by using an operational matrix of Φ-shifted Legendre polynomials.As far as we know, there is no open literature article dealing with the numerical treatment of FRODEs involving theΦ-Caputo fractional derivative employing HW. Therefore, the basic aim of this paper is to provide a numerical technique for solving Φ-FRODEs that arise in physics. Our method is based on a new type of operational matrix of fractional integration called the Φ-HW operational matrix. We provide a rigorous verification of convergence for the suggested method. Furthermore, numerical experiments are presented to demonstrate the convergence of the procedure by comparing the exact values to numerical approximation.Different types of orthogonal polynomials, including Chebyshev polynomials [22], Legendre polynomials [23], and Laguerre and Hermite polynomials [24], have been utilized with the operational matrix of integer order integration. Many authors then expanded it to the fractional case, as seen in [25–31] and the references therein. Only Caputo or Riemann-Liouville fractional derivatives were examined in all of the above listed papers.The following is a description of this paper’s structure.In Section2, we go through the basic fundamentals of fractional calculus. We introduce HW and function approximation using HW in Section 3. In Section 4, we construct an explicit formula for the Φ-fractional integration of the HW and the Φ-Haar wavelet operational matrix. Section 5 discusses the numerical scheme as well as the method’s convergence.


## 2. Preliminaries

We will go over some definitions ofΦ-fractional integral and differential operators in this section.Let the functiong:a1,a2⟶ℝ be integrable, α is a positive real number, n is a natural number, Φ∈C1a1,a2, and Φ′χ≠0∀χ∈a1,a2, where Φ is increasing.Definition 1 (see [32–34]).
TheΦ-Riemann-Liouvile (Φ-RL) fractional integral operator of order α is defined by
(6)Ja1α,Φgχ=1Γα∫a1χΦ′sΦχ−Φsα−1gsds.TheΦ-RL differential operator of fractional α is defined by
(7)Da1α,Φgχ=1Φ′χddχnJa1n−α,Φgχ=1Γn−α1Φ′χddχn∫a1χΦ′sΦχ−Φsn−α−1gsds,wheren=α+1.Definition 2 (see [21, 35, 36]).
Letα be a positive real number, n a natural number, g,Φ∈Cna1,a2, and Φ′χ≠0∀χ∈a1,a2, where Φ is increasing. The Φ-Caputo differential operator of fractional-order α is defined by
(8)DCa1α,Φgχ=1Γn−α∫a1χΦ′sΦχ−Φsn−α−1Dn,Φgsds,where Dn,Φgχ=1/Φ′χd/dχngχ, n=α+1 for α∉ℕ and n=α when α∈ℕ.


### 2.1. Function Approximation by Haar Wavelet

HW are the simplest wavelets with a compact support among the various wavelet families. These wavelets have been shown to be an effective method for numerical function approximation. The Haar functions, which are orthogonal, contain only one wavelet during some subinterval of time and remain zero elsewhere.Theith HW, hiχ, where χ∈a1,a2 in the HW family is defined as [37]
(9)hiχ=1,whenχ∈a1+a2−a1κm,a1+a2−a12κ+12m,−1,whenχ∈a1+a2−a12κ+12m,a1+a2−a1κ+1m,0,otherwise,where m=2j,j=0,1,2,3,⋯,J, is the dilation parameter and κ=0,1,2,3,⋯,m−1 is the translation parameter. J is the max resolution level of HW. The parameters i, j, and κ are related by the equation i=2j+κ+1; i is called the wavelet number. Equation (9) holds true for i≥3.Fori=1 and i=2, the scaling functions for HW family are defined, respectively, by
(10)h1χ=1,whenχ∈a,b,0,otherwise,h2χ=1,whenχ∈a,a+b2,−1,whenχ∈a+b2,b,0,otherwise.A square integrable function ona,b, uχ can be approximated by HW in the following way:
(11)uχ=∑i=0∞cihiχ,where ci is defined by the inner product of uχ and hiχ and . represents the inner product. The first m terms are employed for function approximation, that is,
(12)uχ≅umχ=∑i=0m−1cihiχ,which can be written in the matrix notation as
(13)uχ≅umχ=CmTHmχ,where C=c0,c1,c2,⋯,cm−1T is the coefficient matrix determined by ci=uχ,hiχ and H=h0χ,h1χ,h2χ,⋯,hm−1χT is the vector of Haar functions.


## 2.1. Function Approximation by Haar Wavelet

HW are the simplest wavelets with a compact support among the various wavelet families. These wavelets have been shown to be an effective method for numerical function approximation. The Haar functions, which are orthogonal, contain only one wavelet during some subinterval of time and remain zero elsewhere.Theith HW, hiχ, where χ∈a1,a2 in the HW family is defined as [37]
(9)hiχ=1,whenχ∈a1+a2−a1κm,a1+a2−a12κ+12m,−1,whenχ∈a1+a2−a12κ+12m,a1+a2−a1κ+1m,0,otherwise,where m=2j,j=0,1,2,3,⋯,J, is the dilation parameter and κ=0,1,2,3,⋯,m−1 is the translation parameter. J is the max resolution level of HW. The parameters i, j, and κ are related by the equation i=2j+κ+1; i is called the wavelet number. Equation (9) holds true for i≥3.Fori=1 and i=2, the scaling functions for HW family are defined, respectively, by
(10)h1χ=1,whenχ∈a,b,0,otherwise,h2χ=1,whenχ∈a,a+b2,−1,whenχ∈a+b2,b,0,otherwise.A square integrable function ona,b, uχ can be approximated by HW in the following way:
(11)uχ=∑i=0∞cihiχ,where ci is defined by the inner product of uχ and hiχ and . represents the inner product. The first m terms are employed for function approximation, that is,
(12)uχ≅umχ=∑i=0m−1cihiχ,which can be written in the matrix notation as
(13)uχ≅umχ=CmTHmχ,where C=c0,c1,c2,⋯,cm−1T is the coefficient matrix determined by ci=uχ,hiχ and H=h0χ,h1χ,h2χ,⋯,hm−1χT is the vector of Haar functions.


## 3.Φ-HW Operational Matrix

The fractional-orderedΦ-RL integral of the HW is defined by
(14)Jα,Φh1χ=1Γα+1Φχ−Φa1α,Plα,Φχ=Jα,Φhlχ=1Γα∫a1χΦ′sΦχ−Φsα−1hisds=1Γα+10,ifχ<ζ1l,Φχ−Φζ1lα,ifχ∈ζ1l,ζ2l,Φχ−Φζ1lα−2Φχ−Φζ2lα,ifχ∈ζ2l,ζ3l,Φχ−Φζ1lα−2Φχ−Φζ2lα+Φχ−Φζ3lα,ifχ>ζ3l,where ζ1l=a1+a2−a1κ/m,ζ2l=a1+a2−a12κ+1/2m,ζ3l=a1+a2−a1κ+1/m.TheΦ-HW operational matrix Pα,Φ is computed in the interval [0, 1] for Φχ=χ3 and α=0.75. The numerical and exact Φ‐RL integration of the function Φχ=χ3 for J=6 and different values of α is plotted in Figure 1.
(15)Pα,Φ=0.3331−0.2641−0.0553−0.2141−0.0115−0.0449−0.0850−0.1295−0.12580.1948−0.05530.2584−0.0115−0.04490.12010.1433−0.06620.03950.0403−0.0179−0.01150.0537−0.0146−0.00530.0386−0.038600.145500−0.10260.1988−0.0130−0.00070.0080−0.00160.0078−0.0035−0.0010−0.0005−0.01750.0251−0.0076−0.008800.0310−0.0081−0.0022−0.00510.005100.0680000.0623−0.01350.0285−0.02850−0.05710000.0999.Figure 1
The exact and numericalΦ-RL integral of gχ=Φχ for J=5 and 1<α≤2 and their absolute error.


## 4. Error Analysis

Caputo-type FDEs have recently been investigated in context of error analysis in [38]. In addition, the convergence analysis of solution of nonlinear Fredholm integral equations by HW is given in [39].Using theΦ-Caputo fractional differential operator, we estimated the max error, demonstrating the efficiency of the Φ-HW approach for Φ-FDEs.Theorem 3.
LetDℓy be continuous on a1,a2 and suppose that M>0 so that Dℓ,Φuχ≤M∀χ∈a1,a2, where a1,a2∈ℝ+, Dℓ,Φuχ=1/Φ′χd/dχℓuχ and DCa1α,Φuχ is approximated by DCa1α,Φuℓχ, then we have
(16)DCa1α,Φuχ−DCa1α,ΦuℓχE≤a2−a1MΦ′a2ℓ−αΓℓ−α+11κℓ−α11−22α−ℓ1/2.Proof.
DCa1α,Φy can be approximated by HW as
(17)DCa1α,Φuχ=∑i=a∞cihiχ.
Here,ci is given by
(18)ci=DCa1α,Φuχ,hiχ=∫a1bDCa1α,Φuχhiχdχ.
Let the approximation ofCDa1α,Φu be DCa1α,Φuℓ which is defined by
(19)DCa1α,Φuℓχ=∑i=0ℓ−1cihiχ,in which ℓ=2β+1,β=1,2,3,⋯.
Therefore,(20)DCa1α,Φuχ−DCa1α,Φuℓχ=∑i=m∞cihiχ=∑i=2β+1∞cihiχ.This gives
(21)DCa1α,ΦuχDCa1α,ΦuℓχE2=∫a1χDCa1α,Φuχ−DCa1α,Φuℓχ2dχ=∑i=2β+1∞∑i′=2β+1∞cici′∫a1χhiχhi′χdχ.
The sequencehmχ being orthogonal, we get ∫a1a2hmχhmχdχ=Im, Im which represents the identity matrix of order m.☐Therefore, from Equation (21) we have
(22)DCa1α,Φuχ−DCa1α,ΦuℓχE2=∑i′=2β+1∞ci2.Equation (18) gives
(23)ci=∫a1bDCCa1α,Φuχhiχdχ=2j2∫a1+a2−a1κ2−ja1+a2−a1κ+122−jDCa1α,Φuχdχ−∫a1+a2−a1κ+122−ja1+a2−a1κ+12−jDCa1α,Φuχdχ.Employing mean value theorem of integration∃χ1,χ2∈a1,a2 where
(24)a1+a2−a1κ2−j<χ1<a1+a2−a1κ+122−j,a1+a2−a1κ+122−j<χ2<a1+a2−a1κ+12−j,so that we arrive at
(25)ci=2j/2a2−a1a1+κ+122−j−a1+κ2−jDCa1α,Φuχ1−a1+κ+12−j−a1+κ+122−jDCa1α,Φuχ2=2j/2a2−a12−j−1DCa1α,Φuχ1−DCa1α,Φuχ2.Therefore,(26)ci2=2−j−2a2−a12DCa1α,Φuχ1−DCa1α,Φuχ22.Applying theΦ-Caputo fractional differential operator along with the facts that Φ is increasing and Dℓ,Φuχ≤M, we get
(27)DCa1α,Φuχ1−DCa1α,Φuχ2=1Γℓ−α∫a1χ1Φ′χΦχ1−Φχℓ−α−1Dℓ,Φuχdχ−∫a1χ2Φ′χΦχ2−Φχℓ−α+1Dℓ,Φuχdχ=1Γℓ−α∫a1χ1Φ′χΦχ1−Φχℓ−α−1Dℓ,Φuχdχ−∫a1χ1Φ′χΦχ2−Φχℓ−α+1Dℓ,Φuχdχ−∫χ1χ2Φ′χΦχ2−Φχℓ−α+1Dℓ,Φuχdχ≤1Γℓ−α∫a1χ1Φ′χΦχ1−Φχℓ−α+1Dℓ,Φuχdχ−∫a1χ1Φ′χΦχ2−Φχℓ−α+1Dℓ,Φuχdχ+∫χ1χ2Φ′χΦχ2−Φχℓ−α+1Dℓ,Φuχdχ=1Γℓ−α∫a1χ1Φ′χΦχ1−Φχℓ−α+1−Φχ2−Φχℓ−α+1Dℓ,Φuχdχ+∫χ1χ2Φ′χΦχ2−Φχℓ−α+1Dℓ,Φuχdχ≤1Γℓ−α∫a1χ1Φ′χΦχ1−Φχℓ−α+1−Φχ2−Φχℓ−α+1Dℓ,Φuχdχ+∫χ1χ2Φ′χΦχ2−Φχℓ−α+1Dℓ,Φuχdχ≤1Γℓ−α∫a1χ1Φ′χΦχ1−Φχℓ−α+1−Φχ2−Φχℓ−α+1Dℓ,Φuχdχ+∫χ1χ2Φ′χΦχ2−Φχℓ−α+1Dℓ,Φuχdχ,whereℓ>1+α≤MΓℓ−α∫a1χ1Φ′χΦχ1−Φχℓ−α+1−Φχ2−Φχℓ−α+1dχ+∫χ1χ2Φ′χΦχ2−Φχℓ−α+1dχ=MΓℓ−α1ℓ−αΦχ1−Φa1ℓ−α+Φχ2−Φχ1ℓ−α−Φχ2−Φa1ℓ−α+Φχ2−Φχ1ℓ−α=MΓℓ−α+1Φχ1−Φa1ℓ−α−Φχ2−Φa1ℓ−α+2Φχ2−Φχ1ℓ−α.Asχ1,χ2>a and χ1<χ2 also Φχ are increasing, so
(28)Φχ1−Φa1ℓ−α−Φχ2−Φa1ℓ−α<0.Therefore,(29)DCa1α,Φyχ1−DCa1α,Φyχ2≤2MΓℓ−α+1Φχ2−Φχ1ℓ−α.According to the mean value theorem,∃ζ∈χ1,χ2⊆a1,a2 such that Φχ2−Φχ1≤χ2−χ1Φ′ζ, we get
(30)DCa1α,Φyχ1−DCa1α,Φyχ2≤2MΓℓ−α+1χ2−χ1Φ′ζℓ−α≤2MΓℓ−α+12jℓ−αΦ′a2ℓ−α,which implies that
(31)DCa1α,Φyχ1−DCa1α,Φyχ22≤4M2Γ2ℓ−α+122jℓ−αΦ′a22ℓ−α.Putting (31) in (26), we get
(32)ci2≤2−j−2a2−a124M2Γ2ℓ−α+122jℓ−αΦ′a22ℓ−α.Putting together equations (22) and (32), we have
(33)DCa1α,Φuχ−DCa1α,ΦuℓχE2=∑i=2β+1∞ci2=∑j=β+1∞∑i=2j2j+1−1ci2≤∑j=β+1∞a2−a12M2Γ2ℓ−α+122jℓ−α+jΦ′a22ℓ−α2j+1−1−2j+1=a2−a12M2Φ′a22ℓ−αΓ2ℓ−α+1∑j=β+1∞122jℓ−α=a2−a12M2Φ′a22ℓ−αΓ2ℓ−α+1122β+1ℓ−α11−22α−ℓ,which implies that
(34)CDa1α,Φuχ−DCa1α,ΦuℓχE≤a2−a1MΦ′a2ℓ−αΓℓ−α+112β+1ℓ−α11−22α−ℓ1/2.Letκ=2β+1, (34) can also be written as
(35)DCa1α,Φuχ−DCa1α,ΦuℓχE≤a2−a1MΦ′a2ℓ−αΓℓ−α+11κℓ−α11−22α−ℓ1/2.To compute the error bound, we need the value ofM.So we will estimateM first. Since Dℓuχ is continuous and bounded on a1,a2, so is Dℓ,Φuχ and it is approximated by
(36)Dℓ,Φuχ≅∑i=0r−1cihiχ=CrTHrχ,where Cr=c0,c1,c2,⋯,cr−1T and Hrχ=h0χ,h1χ,h2χ,⋯,hr−1χT.Integration of (36) gives
(37)Dℓ−1,Φuχ=∫a1χDℓ,Φuχdχ+Dℓ−1,Φua1=∫a1χDℓ,Φuχdχ≅CrTP1,ΦHℓχ.Similarly,(38)Dℓ−2,Φuχ=∫a1χDℓ−1,Φuχdχ+Dℓ−2,Φya1=∫a1χDℓ−1,Φuχdχ≅CℓTP2,ΦHℓχ.Continuing in the same manner, we arrive at(39)DΦuχ≅CℓTPℓ,ΦHℓχ.Takingχj=j−1/2/ℓ,j=0,1,2,⋯,m, and putting it in (39), we have
(40)DΦuχj≅CℓTPℓ,ΦHℓχj.The matrix form of (40) is as
(41)DΦUT≅CℓTPℓ,ΦHℓχjwhereDΦUT=DΦuχ1,DΦuχ2,DΦuχ3,⋯,DΦuχℓT.The linear system in (41) determines the value of the vector CℓT; by putting this value in (36), Dℓ,Φχ can be obtained ∀χ∈a1,a2.Letτi∈a1,a2, then Dℓ,Φuti can be computed for the equidistant points i=1,2,3,⋯,ℓ, then ε+maxDℓuti is the approximation of M.Theorem 4.
Assume thatDCa1α,Φuℓ, computed from Φ-HW is estimated by DCa1α,Φu, then we have
(42)uχ−uℓχE≤MNΓα+1Γℓ−α+11κℓ−α11−22α−ℓ1/2,where N=max∣a2−a1Φa2ℓ−αΦχ−Φ0α∣.Theorem4 can be proven easily by following the procedure of Theorem 3. From Equation (42), we noted that uχ−uℓχE⟶0 as ℓ⟶∞. Thus, the convergence of the Φ-HW method is inferred.


## 5. Numerical Examples

We present several examples of how to get numerical solutions ofΦ-FRODEs using the Φ-HW operational matrix approach.Example 5.
Consider theΦ-FRODE
(43)Dα,Φuχ+auχ=gχ,0<α≤1,χ∈0,1,u0=0.
Fora=2/Γ3−α and gχ=2/Γ3−αΦχ2−α+Φχ2, the actual solution of Equation (43) is uχ=Φχ2. We use the Φ-HW technique to solve problem (43).
Let(44)DCα,Φuχ=CmTHmχ.Integrating Equation (44) with respect to Jaα,Φ and using the initial conditions, we have
(45)uχ=Jα,ΦCmTHmχ=CmTPm×mα,ΦHmχ.Substituting (44) and (45) into (43), we have
(46)CmTHmχ+aPm×mα,ΦHmχ=gχ.Equation (46) has the following matrix form:
(47)CmTHmχ+aPm×mα,ΦHmχ=G,where G is the matrix representation of g at the collocation points.Solving the algebraic system given by Equation (47) for CmT and substituting this value into Equation (45), we will have the required numerical solution. In Table 1, the max absolute error is given for J=6 and φζ=ζ3. Approximate solutions for J=6 and α=0.5 and different choices of the function Φ are plotted in Figure 2. Also, actual and approximate results and the absolute error are given for J=6, α=0.8, and Φχ=1/3χ3−χ2−χ in Figure 2.Example 6.
Consider the compositeΦ-FRODE
(48)Dα,Φuχ+buχ=Γ2α+1Γα+1Φχα1+Φχα,0<α≤1,x∈0,1,(49)u0=0.
Forb=Γ2α+1/Γα+1.Table 1
Max absolute error for various choices ofα and J.
αJ=5J=6J=7J=8J=90.61.0141×10−43.2169×10−51.0296×10−53.3179×10−61.0751×10−60.79.2421×10−52.7180×10−58.0595×10−62.4056×10−67.2213×10−70.87.8481×10−52.1470×10−55.9091×10−61.6349×10−64.5452×10−70.96.3749×10−51.6437×10−54.2464×10−61.0994×10−62.8527×10−71.05.3010×10−51.3250×10−53.3127×10−68.2817×10−72.0704×10−7Figure 2
Approximate and exact results of Equation (43) and their absolute error.The exact solution for the problem (48) is uχ=Φχ2α.For the numerical solution, we employ theΦ-HW technique.Let(50)DCα,Φuχ=CmTHmχ.Integrating Equation (50) with respect to Jaα,Φ and utilizing the initial condition, we have
(51)uχ=Jα,ΦCmTHmx=CmTPm×mα,ΦHmχ.Substituting Equations (50) and (51) in Equation (48), we get
(52)CmTHmχ+Pm×mα,ΦHmχ=gχ,where Γ2α+1/Γα+1Φχα1+Φχα. Equation (52) in matrix form is given as
(53)CmTHmχ+Pm×mα,ΦHmχ=G.where G is the matrix representation of gχ.Required approximate solutions can be obtained by using the value ofCmT from Equation (53) into Equation (51).Themax absolute error is tabulated for Φχ=χ3/5 and various choices of J and α in Table 2, which shows that the Maximum Absolute Error decreases by increasing the values of J. Figure 3 represents approximate solutions for different choices of α. Also, comparison of actual and approximate results and their absolute error is displayed in Figure 3.Example 7.
Consider theΦ-FRODE
(54)Dα,Φuχ+uχ=1−4Φχ+5Φχ2−4Γ2−αΦχ1−α+10Γ3−αΦχ2−α,where 0<α≤1,x∈0,1,andu0=1. It is easy to verify that uχ=1−4Φχ+5Φχ2 is the actual solution of Equation (54). For numerical approximation, we employ the Φ-HW technique.Table 2
Max absolute error for various choices ofα and J.
αJ=5J=6J=7J=8J=90.61.2082×10−43.9804×10−51.3061×10−54.2822×10−61.4045×10−60.79.1203×10−52.8002×10−58.5466×10−62.6042×10−67.9372×10−70.86.7053×10−51.9253×10−55.4799×10−61.5545×10−64.4065×10−70.94.8650×10−51.3151×10−53.5128×10−69.3312×10−72.4729×10−71.03.5205×10−59.0544×10−62.2954×10−65.7785×10−71.4496×10−7Figure 3
Approximate results of Equation (48) for various choices of α, actual, approximate results, and the absolute error.Let(55)DCα,Φuχ=CmTHmχ.Integrating Equation (55) with respect to Jaα,Φ and using the initial conditions, we have
(56)uχ=Jα,ΦCmTHmx+u0=CmTPm×mα,ΦHmχ+1.Substituting (55) and (56) into (54), we have
(57)CmTHmχ+Pm×mα,ΦHmχ=gχ,where gχ=−4Φχ+5Φχ2−4/Γ2−αΦχ1−α+10/Γ3−αΦχ2−α.The matrix form of Equation (57) is
(58)CmTHmχ+Pm×mα,ΦHmχ=G,where G is the matrix representation of gχ.Required approximate solutions can be obtained by using the value ofCmT from Equation (58) in Equation (56). Table 3 shows that the Maximum Absolute Error decreases by increasing the values of J. Approximate solutions are displayed in Figure 4 for various values of Φ. Also, Figure 4 represents approximate and exact solutions and their max absolute error for α=0.75,J=6, and Φχ=χ3/15.Example 8.
Consider theΦ-FRODE
(59)Dα,Φuχ+μuχ=gχ,0<α≤1,x∈0,1u0=0.
Forμ=1 and gχ=Γ2α+1/Γ1+αΦχα+Γ2/Γ2−αΦχ1+α+Φχ2α+Φχ, the exact solution of Equation (59) is Φχ2α+Φχ.Table 3
Max absolute error forΦχ=χ3/15 and various choices of J and α.
αJ=5J=6J=7J=8J=90.69.0149×10−52.9995×10−59.9647×10−63.3058×10−61.0954×10−60.74.3216×10−51.3568×10−54.2468×10−61.3254×10−64.1271×10−70.81.5201×10−54.5357×10−61.3466×10−63.9805×10−71.1721×10−70.92.2599×10−56.0704×10−61.5981×10−64.1680×10−71.0825×10−71.02.2302×10−55.7525×10−61.4605×10−63.6795×10−79.2342×10−8Figure 4
Approximate solutions for different choices ofα and functions Φx.For approximate solutions, we use theΦ-HW technique.Let(60)CDα,Φuχ=CmTHmχ.Integrating Equation (60) with respect to J0α,Φ and using the initial conditions, we have
(61)uχ=J0α,ΦCmTHmx=CmTPm×mα,ΦHmχ.Substituting (60) and (61) into (59), we have
(62)CmTHmχ+μPm×mα,ΦHmχ=gχ.The matrix representation of Equation (62) is
(63)CmTHmχ+μPm×mα,ΦHmχ=G,where G is the matrix representation of gχ at the collocation points.Required approximate solutions can be obtained by using the value ofCmT from Equation (58) in Equation (61). Table 4 shows that the Maximum Absolute Error decreases by increasing the values of J. Also the approximate solutions are displayed in Figure 5 for various values of α.Table 4
Max absolute error forΦχ=χ2/15 and various choices of J and α.
αJ=5J=6J=7J=8J=90.63.8255×10−51.2688×10−54.1987×10−61.3876×10−64.5829×10−70.72.0382×10−56.3431×10−61.9679×10−66.0937×10−71.8843×10−70.89.0481×10−62.6383×10−67.6857×10−72.2372×10−76.5069×10−80.93.0670×10−68.0025×10−72.1048×10−75.5848×10−81.4938×10−81.02.0360×10−65.1692×10−71.3022×10−73.2678×10−88.1852×10−9Figure 5
Numerical solutions forα=1,J=8 and for different functions Φχ.


## 6. Conclusion

This study introduces a numerical approach for solving a class of fractional differential equations with aΦ-Caputo fractional derivative based on a novel type of operational matrix of fractional integration, namely, the Φ-HW operational matrix. The convergence of the proposed method is demonstrated, and the numerical tests reported in Section 5 corroborate the efficacy of our approach.



---

*Source: 7117064-2021-10-11.xml*