# A Numerical Scheme Based on the Chebyshev Functions to Find Approximate Solutions of the Coupled Nonlinear Sine-Gordon Equations with Fractional Variable Orders

Authors: MohammadHossein Derakhshan
Journal: Abstract and Applied Analysis
 (2021)

---

## Abstract

In this article, a numerical method based on the shifted Chebyshev functions for the numerical approximation of the coupled nonlinear variable-order fractional sine-Gordon equations is shown. The variable-order fractional derivative is considered in the sense of Caputo-Prabhakar. To solve the problem, first, we obtain the operational matrix of the Caputo-Prabhakar fractional derivative of shifted Chebyshev polynomials. Then, this matrix and collocation method are used to reduce the solution of the nonlinear coupled variable-order fractional sine-Gordon equations to a system of algebraic equations which is technically simpler for handling. Convergence and error analysis are examined. Finally, some examples are given to test the proposed numerical method to illustrate the accuracy and efficiency of the proposed method.


---

## Body

## 1. Introduction

The basic and important applications of differential equations with fractional variable order which are considered as a generalization of ordinary differential equations with integer order in various fields of computational science [1–3], engineering [4–6], physics [7, 8], and chemistry [9, 10] have attracted significant attention in the literature. One of the important discussions in differential equations is the study of nonlinear partial and nonlinear differential equations that are used in physics and applied mathematics such as the sine-Gordon equation which is a significant nonlinear integrable evolution partial differential system in space-time coordinates [11]. The coupled nonlinear partial differential sine-Gordon equations of integer order on the interval p,q∈a,b×0,T are defined by [12]
(1)zqqp,q−zppp,q+θ2sinzp,q−vp,q=h1p,q,(2)vqqp,q−ξ2vppp,q−sinzp,q−vp,q=h2p,q,where symbols θ and ξ2are the ratios of the acoustic velocities of the functionsz,v, and the ratios of masses of particles in the “lower” and “upper” parts of the crystal given in the generalized Frenkel-Kontorova dislocation equation [12–14] are the defined parameters θ and ξ2, respectively, which are known, and the functions zp,q,vp,q and h1p,q,h2p,q are unknown and known, respectively. If in Equation (2), ξ=1 is considered, then Equation (2) explains the open states in deoxyribonucleic acid (DNA) [15], and if in Equation (2), ξ2≠1 is considered, then Equation (2) shows a gap in the spectrum of velocities of solitary waves [12]. This paper focuses on finding approximate solutions based on the Chebyshev function scheme for fractional differential equations of variable orders which are called coupled nonlinear sine-Gordon equations, and these type of equations of fractional variable orders α1p,q,α2p,q∈1,2 are defined as follows:
(3)DC−Pqα1p,qzp,q−DC−Pqα1p,qzp,q+θ2sinz−v=h1p,q,(4)DC−Pqα2p,qvp,q−ξ2DC−Pqα2p,qvp,q−sinz−v=h2p,q,where DC−P is the Caputo-Prabhakar operator. Also, the initial and boundary conditions for problems (3) are considered as follows:
(5)zp,0=l0p,zqp,0=l1p,(6)vp,0=l2p,vqp,0=l3p,(7)za,q=f0q,zb,q=f1q,(8)va,q=f2q,vb,q=f3q,where functions l0p,l1p,l2p,l3p and f0q,f1q,f2q,f3q are continuous. The Caputo-Prabhakar operator considered in this article is an extension of the Caputo derivative which is introduced in [16], and it is defined by
(9)D0+μ1zx=I0+n−μ1dndxnzx,where n−1<μ1≤n,z∈L10,b,b∈ℝ and I is the Riemann–Liouville integral (RLI). So by putting a function containing the Prabhakar function which is studied in [17], in Equation (9), by developing the RLI kernel with its function, this function is given by [17]
(10)eϱ,vυx=xν−1Eϱ,vvωxϱ,Rϱ>0,Rv>0,x∈ℝ,ϱ,v,υ,ω∈ℂ,where Eϱ,vυωxϱ is the Prabhakar function and is defined as follows:
(11)Eϱ,vυωxϱ=∑l=0∞Γv+lxl/l!ΓvΓϱl+v.Some applications of the Prabhakar function can be seen in mathematics and physics as a fractional Poisson process [16], Havriliak–Negami relaxation functions [18, 19], irregular case of the dielectric relaxation responses [20], a model of anomalous relaxation in dielectrics of fractional order [21], fractional thermoelasticity [10], telegraph equations [22], thermodynamics [23], and fractal time random [24]. By placing α1p,q=2,α2p,q=2 in Equation (3), the coupled nonlinear sine-Gordon equations of fractional variable orders given in (3) change into the classical coupled nonlinear sine-Gordon equations which are defined by Equation (2), and the classical coupled nonlinear sine-Gordon equations have many applications in physics as nonlinear models [25, 26], plasma [27], quantum [28], optics [29], and mathematics [13, 30, 31]. Getting analytic solutions to fractional differential equations in general are not easy; therefore, numerical methods are used to obtain the solutions of this type of equations. For instance, [32] defined a new method based on the linearization formula to find the fuzzy approximate solution of fractional differential equations under uncertainty, a new spectral tau method for solving equations of Kelvin–Voigt equations is proposed in [33], and [34] studied a numerical method based on the Bernstein functions to solve the linear cable equation of variable order. An implicit RBF meshless method [35], spectral method [36], homotopy analysis scheme [37], meshless method [38], collocation and finite difference-collocation methods [39], Ritz–Galerkin method [40], new collocation scheme for solving fractional partial differential equations [41], Sinc-Chebyshev collocation scheme [42], numerical method based on the shifted Chebyshev functions [43], numerical method based on the shifted Legendre functions [44], piecewise integroquadratic spline interpolation and finite difference method [45], wavelet scheme [46], Legendre wavelet method [47], Bernoulli wavelet method [48], and Lagrange multiplier scheme [49] and other methods [50–54]. Several methods are presented for solving variable-order fractional equations. For example, Bhrawy and Zaky [55] used the shifted Jacobi polynomials to obtain solution variable-order fractional Schrödinger equations, Bhrawy and Zaky [56] studied the Jacobi–Gauss–Lobatto collocation method to obtain solution variable-order fractional Schrödinger equations, Mahmoud et al. [57] proposed the Jacobi wavelet collocation method to obtain solution variable-order fractional equations, and Zaky et al. [58] proposed the shifted Chebyshev polynomials to obtain solution variable-order fractional equations; in [59], a proper discrete form of fractional Grönwall-type inequality is introduced and other methods such as shifted Jacobi collocation method [60] and shifted Jacobi method [61]. Since this mathematical system which is given in Equation (3), due to the variable-order fractional operators and nonlinearity, is very complex, we need to reduce it by using a highly accurate and efficient expansion scheme. Thus, for this aim, we first extract an operational matrix of variable-order fractional derivatives for the shifted Chebyshev functions, then we use them for extending the unknown solution. In other words, the problem shown in Equation (3) is converted into an algebraic system of equations by exploiting the operational matrix of variable-order derivative. In this paper, we consider a class of coupled nonlinear variable-order fractional sine-Gordon equations. For solving the given equations, operational matrices based on the shifted Chebyshev functions are applied. First, we approximate the unknown function and its derivatives in terms of the shifted Chebyshev functions. Then, by substituting these approximations into the equation and applying the properties of the shifted Chebyshev functions together with the collocation points, the main problem is reduced to a set of nonlinear algebraic equations. By solving this system, the approximate solution is calculated. This article is divided into the following sections. In Section 2, we state the essential definitions and lemmas about fractional derivatives and integrals of the variable order, also in this section are introduced the Chebyshev polynomials and the shifted Chebyshev polynomials. The approximate function is obtained in Section 3. In Section 4, an explanation of the suggested method for obtaining approximate solutions of the equations introduced in (3) is provided and four numerical examples are shown in Section 5.


## 2. Preliminaries of Fractional Calculus

This section discusses the important topics of fractional derivative and integral of the variable order such as the Prabhakar integral and the Caputo-Prabhakar derivative. The Chebyshev polynomials and the shifted Chebyshev polynomials are introduced in this section.Definition 1.
The Prabhakar integral of variable orderRμp,q>0,q≥0 is defined by [17, 62]
(12)ε1ϱ,μp,q,ω,0v+zp,q=∫0qq−τμp,q−1Eϱ,μp,qvωq−τϱzp,τdτ,ε2ϱ,μp,q,ω,0v+zp,q=∫0pp−τμp,q−1Eϱ,μp,qvωp−τϱzτ,qdτ,where f∈L10,∞. Also, for υ=0, the Prabhakar integral coincides with the Riemann-Liouville fractional integral of variable order Rμp,q>0,q≥0.Definition 2.
The Caputo-Prabhakar derivative of variable-orderα1p,q,n−1<α1p,q≤n,n∈ℕ,q≥0 is defined by [17, 62, 63]
(13)DC−Pqα1p,qzp,q=ε1ϱ,n−α1p,q,ω,0+−vdndqnzp,q,DC−Pqα1p,qzp,q=ε2ϱ,n−α1p,q,ω,0+−vdndpnzp,q,where Rϱ>0,Rα1p,q>0,ϱ,α1p,q,v,ω∈ℂ. Here, we consider n=2.Lemma 3.
Letk>1. Then, the following formula based on Definition 1 is established [63]:
(14)ε1ϱ,μp,q,ω,0+vqk−1=Γkqk+μp,q−1eϱ,k+μp,q+1vq,ε2ϱ,μp,q,ω,0+vpk−1=Γkpk+μp,q−1eϱ,k+μp,q+1vp.Lemma 4.
Letk>2. Then the following formula based on the Definition 2 and the Lemma 3 for n=2 are established:
(15)DC−Pqα1p,qqk−1=Γkqk+μp,q−1eϱ,k+μp,q+1−vq,DC−Ppα1p,qpk−1=Γkpk+μp,q−1eϱ,k+μp,q+1−vp.Definition 5.
The Chebyshev polynomials on−1,1 are defined as follows [64]:
(16)T0p=1,T1p=p,Tm+1p=2pTmp−Tm−1p,where Tmp=m∑j=0m/2−1j2m−2j−1m−j−1!/j!m−2j!pm−2j. The orthogonality manner for Chebyshev polynomials is given by
(17)∫−11ΩpTipTjpdp=π,i=j=0,π2,i=j≠0,0,i≠j,where Ωp=1/1−p2 is a weight function for Chebyshev polynomials. Also, the shifted Chebyshev polynomials on 0,1 are defined as follows:
(18)T^0q=T02p−1=1,T^1q=T12p−1=q,Tm+12p−1⏟T^m+1q=22p−1Tm2p−1−Tm−12p−1⏟2qT^mq−T^m−1q,where T^mq=m∑j=0m−1m−j22jm+j−1!/2j!m−j!qj. The orthogonality conditions for shifted Chebyshev polynomials are given by
(19)∫01Ω∗qT^iqT^jqdq=π2,i=j=0,π4,i=j≠0,0,i≠j,where Ω∗q=1/1−q2 is a weight function.Lemma 6.
LetΨq=T∧0q,T∧1q,⋯,T∧mqT. Then, the function Ψq has a matrix display as follows [64]:
(20)Ψq=AHmq,where A and Hq are defined by
(21)A=100⋯0−120⋯02−121!2!2−11222!2!2−10243!4!⋯0⋮⋮⋮⋱⋮m−1mm−1!m!m−1m−122m!2!m−1!m−1m−224m+1!4!m−2!⋯m−1022m2m−1!2m!,Hmq=1q⋮qm.
Also, Lemma6 is established for the function Ψp.


## 3. Approximate Function

The approximation functions aszp,q,vp,q on L20,1×0,1 can be developed in shifted Chebyshev polynomials as follows:
(22)zp,q≃zmp,q=∑i=0m∑j=0mzijT^ipT^jq=ΨTpZΨq,(23)vp,q≃vmp,q=∑i=0m∑j=0mvijT^ipT^jq=ΨTpVΨq,where Z=zijm+1×m+1, zij=T^ip,zp,q,T^jq and V=vijm+1×m+1, vij=T^ip,vp,q,T^jq. By giving effect Caputo-Prabhakar derivative C−PDqα1p,q on Equation (22), we obtain
(24)C−PDqα1p,qzp,q≃C−PDqα1p,qzmp,q=ΨTpZC−PDqα1p,qΨq,using Lemma 6, then from Equation (24), we have
(25)C−PDqα1p,qzmp,q=ΨTpZC−PDqα1p,qΨq=ΨTpZAC−PDqα1p,qHmp=ΨTpZAC−PDqα1p,q1C−PDqα1p,qq⋮C−PDqα1p,qqm,applying Lemma 4 on Equation (25), we have
(26)C−PDqα1p,qzmp,q=ΨTpZAC−PDqα1p,q1C−PDqα1p,qq⋮C−PDqα1p,qqm=ΨTpZA00⋮Γm+1qm−α1p,qeϱ,m−α1p,q+1q−ν=ΨTpZA00⋯00Γ2q1−α1p,qeϱ,2−α1p,qq−v⋯0⋮⋮⋱⋮000Γm+1qm−α1p,qeϱ,m−α1p,q+1q−v×1q⋮qm=ΨTpZAM1Hmq=ΨTpZAM1A−1Ψq,C−PDqα1p,qzp,q≅C−PDqα1p,qzmp,q=ΨTpZAM1A−1Ψq,where
(27)M1=00⋯00Γ2q1−α1p,qeρ,2−α1p,q−υq⋯0⋮⋮⋱⋮000Γm+1qm−α1p,qeρ,m−α1p,q+1−υq.With a similar process as Equation (26), these results are easily achieved:
(28)C−PDqα2p,qvmp,q=ΨTpVAM2A−1Ψq,where
(29)M2=00⋯00Γ2q1−α2p,qeρ,2−α2p,q−υq⋯0⋮⋮⋱⋮000Γm+1qm−α2p,qeρ,m−α2p,q+1−υq,C−PDpα1p,qzmp,q=C−PDpα1p,q1,p,⋯,pmATZΨq=0,0,⋯,Γm+1pm−α1p,qeρ,m−α1p,q+1−υpATZΨq=1,p,⋯,pm00⋯00Γ2p1−α1p,qeρ,2−α1p,q−υp⋯0⋮⋮⋱⋮000Γm+1pm−α1p,qeρ,m−α1p,q+1−υp=HmTpM3ATVΨq=ΨTpA−1TM3ATVΨq,where(30)M3=00⋯00Γ2p1−α1p,qeρ,2−α1p,q−υp⋯0⋮⋮⋱⋮000Γm+1pm−α1p,qeρ,m−α1p,q+1−υp,(31)C−PDpα2p,qvmp,q=C−PDpα2p,q1,p,⋯,pmATVΨq=0,0,⋯,Γm+1pm−α1p,qeρ,m−α1p,q+1−υpATVΨq=1,p,⋯,pm00⋯00Γ2p1−α2p,qeρ,2−α2p,q−υp⋯0⋮⋮⋱⋮000Γm+1pm−α2p,qeρ,m−α2p,q+1−υp=HmTpM4ATVΨq=ΨTpA−1TM4ATVΨq,where
(32)M4=00⋯00Γ2p1−α2p,qeρ,2−α2p,q−υp⋯0⋮⋮⋱⋮000Γm+1pm−α2p,qeρ,m−α2p,q+1−υp.


### 3.1. Description of Convergence

Letzp,q,vp,q be the exact solutions of Equations (3) and zmp,q,vmp,q be the numerical solutions of Equations (3) which are introduced in (22)–(23), then we show that asm⟶∞, we have∥z−zm∥2⟶0and∥v−vm∥2⟶0. For this purpose, we show that the sequences zmp,q,vmp,q are a Cauchy sequence in Hilbert space L2a,b×0,T. Then, for m^>m, we obtain
(33)∥zm^−zm∥2=∥∑i=0m^∑j=0m^zijT^ipT^jq−∑i=0m∑j=0mzijT^ipT^jq∥2=∥∑i=m+1m^∑j=m+1m^zijT^ipT^jq∥2=∑i=m+1m^∑j=m+1m^zijT^ipT^jq,∑l=m+1m^∑k=m+1m^zlkT^lpT^kq=∑i=m+1m^∑j=m+1m^∑l=m+1m^∑k=m+1m^zijz¯lkT^ipT^jq,T^xpT^tq=π216∑i=m+1m^∑j=m+1m^zij2.From Bessel’s inequality, we obtain that the series is presented in (33) asm⟶∞andm^⟶∞is convergent, then∥zm^−zm∥2⟶0.∥zm^−zm∥2⟶0 show that the sequence zm is a Cauchy sequence in Hilbert space L2a,b×0,T and it converges to say z1p,q. Now, we show z=z1p,q, then
(34)z−z1p,q,T^ipT^jq=z,T^ipT^jq−z1p,q,T^ipT^jq=limm⟶∞zm^,T^ipT^jq−zij=zij−zij=0,where zij=∫0T∫abzp,qT^ipT^jqω∗pω∗qdpdq. Then, ∑i=0m∑j=0mzijT^ipT^jq converges to zp,q,vp,q as m⟶∞. The same process can be considered for the sequence vm.


### 3.2. Error Analysis

Suppose thatzp,q,vp,q∈CΩ,Ω=0,1×0,1 be the exact solutions, respectively. If zmp,q=Λ′×Λ,vmp,q=Λ′×Λ be the best numerical solution tozp,q that Λ=SpanT^0q,T^1q,⋯,T^mq,Λ′∈SpanT^0p,T^1p,⋯,T^mp. We define Em which is the error as follows:
(35)C−PDqα1p,qzp,q−C−PDpα1p,qzp,q+θ2sinz−v−h1p,q=C−PDqα1p,qzmp,q−C−PDpα1p,qzmp,q+θ2sinzm−vm−h1p,q+Em.Then, from Eq. (35), we have:
(36)∥Em∥L2Ω=∥C−PDqα1p,qzm−z−C−PDpα1p,qzm−z+θ2sinzm−vm−sinz−v∥L2Ω≤∥C−PDqα1p,qzm−z∥L2Ω+∥C−PDpα1p,qzm−z∥L2Ω+∥θ2sinzm−vm−sinz−v∥L2Ω≤∥C−PDqα1p,qzm−z∥L2Ω+∥C−PDpα1p,qzm−z∥L2Ω+θ2∥sinzm−vm∥L2Ω+θ2∥sinz−v∥L2Ω.To get Equation (36), we define
(37)Δmαp,q,α′p,qp,q=∑i=0m∑j=0mpiαp,qqjα′p,qΓiαp,q+1Γjα′p,q+1C−PDpiαp,qDqjα′p,qC−Pzp,q0,0.From the multivariable Taylor definition and the Generalized Taylor’s formula in [65, 66], we have
(38)zp,q−Δmαp,q,α′p,qp,q≤pm+1αp,qqm+1α′p,qΓm+1αp,q+1Γm+1α′p,q+1Ymαp,q,α′p,q,where Ymαp,q,α′p,q=supp,q∈ΩC−PDpiαp,qC−PDqjα′p,qzp,q∣. By using Lemma 4, we have
(39)∣C−PDqα1p,qzp,q−C−PDqα1p,qΔmαp,q,α′p,qp,q∣≤pm+1αp,qqm+1α′p,q−α1p,qeρ,m+1α′p,q−α1p,q+1−υqΓm+1αp,q+1×Ymαp,q,α′p,q,∣C−PDpα1p,qzp,q−C−PDpα1p,qΔmαp,q,α′p,qp,q∣≤qm+1α′p,qpm+1α′p,q−α1p,qeρ,m+1α′p,q−α1p,q+1−υpΓm+1α′p,q+1×Ymαp,q,α′p,q.Then,(40)∥C−PDqα1p,qzm−z∥L2Ω2≤∥C−PDqα1p,qzp,q−C−PDqα1p,qΔmαp,q,α′p,qp,q∥L2Ω2≤∫01∫01xm+1αp,qtm+1α′p,q−α1p,qeρ,m+1α′p,q−α1p,q+1−υtΓm+1αp,q+1Ymαp,q,α′p,q2dtdx.According to the definitioneρ,m+1α′p,q−α1p,q+1−υt, we have
(41)eρ,m+1α′p,q−α1p,q+1−υt=tm+1α′p,q−α1p,qEρ,m+1α′p,q−α1p,q+1υωtρ.From [16, 67, 68], the Prabhakar function is convergent and the convergent boundary for its M is considered, then we obtain
(42)∥C−PDqα1p,qzm−z∥L2Ω2≤M2Ymαp,q,α′p,q22m+1αp,q+14m+1α′p,q−α1p,q+1Γm+1αp,q+12,also
(43)∥C−PDpα1p,qzm−z∥L2Ω2≤∥C−PDpα1p,qzp,q−C−PDpα1p,qΔmαp,q,α′p,qp,q∥L2Ω2≤∫01∫01tm+1α′p,qxm+1α′p,q−α1p,qeρ,m+1α′p,q−α1p,q+1−υxΓm+1α′p,q+1Ymαp,q,α′p,q2dtdx,from
(44)eρ,m+1α′p,q−α1p,q+1−υx=xm+1α′p,q−α1p,qEρ,m+1α′p,q−α1p,q+1υωxρ,∣Eρ,m+1α′p,q−α1p,q+1υωxρ∣≤M1,we have
(45)∥C−PDpα1p,qzm−z∥L2Ω2≤M12Ymαp,q,α′p,q22m+1α′p,q+14m+1α′p,q−α1p,q+1Γm+1α′p,q+12.So by placing Equations (42) and (45) into (36), we have
(46)∥Em∥L2Ω≤MYmαp,q,α′p,q2m+1αp,q+14m+1α′p,q−α1p,q+1Γm+1αp,q+1+M1Ymαp,q,α′p,q2m+1α′p,q+14m+1α′p,q−α1p,q+1Γm+1α′p,q+1+θ2R+R1,where R=supp,q∈Ω∣sinz−v∣ and R1=supp,q∈Ω∣sinzm−vm∣. With a similar process, the following inequality can be considered for the function E′m that is defined as follows:
(47)C−PDqα2p,qvp,q−ξ2Dpα2p,qC−Pvp,q−sinz−v−h1p,q=C−PDqα2p,qvmp,q−ξ2Dpα2p,qC−Pvmp,q−sinzm−vm−h1p,q+E′m,∥E′m∥L2Ω≤MYmαp,q,α′p,q2m+1αp,q+14m+1α′p,q−α2p,q+1Γm+1αp,q+1+ξ2M1Ymαp,q,α′p,q2m+1α′p,q+14m+1α′p,q−α2p,q+1Γm+1α′p,q+1+R+R1.


## 3.1. Description of Convergence

Letzp,q,vp,q be the exact solutions of Equations (3) and zmp,q,vmp,q be the numerical solutions of Equations (3) which are introduced in (22)–(23), then we show that asm⟶∞, we have∥z−zm∥2⟶0and∥v−vm∥2⟶0. For this purpose, we show that the sequences zmp,q,vmp,q are a Cauchy sequence in Hilbert space L2a,b×0,T. Then, for m^>m, we obtain
(33)∥zm^−zm∥2=∥∑i=0m^∑j=0m^zijT^ipT^jq−∑i=0m∑j=0mzijT^ipT^jq∥2=∥∑i=m+1m^∑j=m+1m^zijT^ipT^jq∥2=∑i=m+1m^∑j=m+1m^zijT^ipT^jq,∑l=m+1m^∑k=m+1m^zlkT^lpT^kq=∑i=m+1m^∑j=m+1m^∑l=m+1m^∑k=m+1m^zijz¯lkT^ipT^jq,T^xpT^tq=π216∑i=m+1m^∑j=m+1m^zij2.From Bessel’s inequality, we obtain that the series is presented in (33) asm⟶∞andm^⟶∞is convergent, then∥zm^−zm∥2⟶0.∥zm^−zm∥2⟶0 show that the sequence zm is a Cauchy sequence in Hilbert space L2a,b×0,T and it converges to say z1p,q. Now, we show z=z1p,q, then
(34)z−z1p,q,T^ipT^jq=z,T^ipT^jq−z1p,q,T^ipT^jq=limm⟶∞zm^,T^ipT^jq−zij=zij−zij=0,where zij=∫0T∫abzp,qT^ipT^jqω∗pω∗qdpdq. Then, ∑i=0m∑j=0mzijT^ipT^jq converges to zp,q,vp,q as m⟶∞. The same process can be considered for the sequence vm.


## 3.2. Error Analysis

Suppose thatzp,q,vp,q∈CΩ,Ω=0,1×0,1 be the exact solutions, respectively. If zmp,q=Λ′×Λ,vmp,q=Λ′×Λ be the best numerical solution tozp,q that Λ=SpanT^0q,T^1q,⋯,T^mq,Λ′∈SpanT^0p,T^1p,⋯,T^mp. We define Em which is the error as follows:
(35)C−PDqα1p,qzp,q−C−PDpα1p,qzp,q+θ2sinz−v−h1p,q=C−PDqα1p,qzmp,q−C−PDpα1p,qzmp,q+θ2sinzm−vm−h1p,q+Em.Then, from Eq. (35), we have:
(36)∥Em∥L2Ω=∥C−PDqα1p,qzm−z−C−PDpα1p,qzm−z+θ2sinzm−vm−sinz−v∥L2Ω≤∥C−PDqα1p,qzm−z∥L2Ω+∥C−PDpα1p,qzm−z∥L2Ω+∥θ2sinzm−vm−sinz−v∥L2Ω≤∥C−PDqα1p,qzm−z∥L2Ω+∥C−PDpα1p,qzm−z∥L2Ω+θ2∥sinzm−vm∥L2Ω+θ2∥sinz−v∥L2Ω.To get Equation (36), we define
(37)Δmαp,q,α′p,qp,q=∑i=0m∑j=0mpiαp,qqjα′p,qΓiαp,q+1Γjα′p,q+1C−PDpiαp,qDqjα′p,qC−Pzp,q0,0.From the multivariable Taylor definition and the Generalized Taylor’s formula in [65, 66], we have
(38)zp,q−Δmαp,q,α′p,qp,q≤pm+1αp,qqm+1α′p,qΓm+1αp,q+1Γm+1α′p,q+1Ymαp,q,α′p,q,where Ymαp,q,α′p,q=supp,q∈ΩC−PDpiαp,qC−PDqjα′p,qzp,q∣. By using Lemma 4, we have
(39)∣C−PDqα1p,qzp,q−C−PDqα1p,qΔmαp,q,α′p,qp,q∣≤pm+1αp,qqm+1α′p,q−α1p,qeρ,m+1α′p,q−α1p,q+1−υqΓm+1αp,q+1×Ymαp,q,α′p,q,∣C−PDpα1p,qzp,q−C−PDpα1p,qΔmαp,q,α′p,qp,q∣≤qm+1α′p,qpm+1α′p,q−α1p,qeρ,m+1α′p,q−α1p,q+1−υpΓm+1α′p,q+1×Ymαp,q,α′p,q.Then,(40)∥C−PDqα1p,qzm−z∥L2Ω2≤∥C−PDqα1p,qzp,q−C−PDqα1p,qΔmαp,q,α′p,qp,q∥L2Ω2≤∫01∫01xm+1αp,qtm+1α′p,q−α1p,qeρ,m+1α′p,q−α1p,q+1−υtΓm+1αp,q+1Ymαp,q,α′p,q2dtdx.According to the definitioneρ,m+1α′p,q−α1p,q+1−υt, we have
(41)eρ,m+1α′p,q−α1p,q+1−υt=tm+1α′p,q−α1p,qEρ,m+1α′p,q−α1p,q+1υωtρ.From [16, 67, 68], the Prabhakar function is convergent and the convergent boundary for its M is considered, then we obtain
(42)∥C−PDqα1p,qzm−z∥L2Ω2≤M2Ymαp,q,α′p,q22m+1αp,q+14m+1α′p,q−α1p,q+1Γm+1αp,q+12,also
(43)∥C−PDpα1p,qzm−z∥L2Ω2≤∥C−PDpα1p,qzp,q−C−PDpα1p,qΔmαp,q,α′p,qp,q∥L2Ω2≤∫01∫01tm+1α′p,qxm+1α′p,q−α1p,qeρ,m+1α′p,q−α1p,q+1−υxΓm+1α′p,q+1Ymαp,q,α′p,q2dtdx,from
(44)eρ,m+1α′p,q−α1p,q+1−υx=xm+1α′p,q−α1p,qEρ,m+1α′p,q−α1p,q+1υωxρ,∣Eρ,m+1α′p,q−α1p,q+1υωxρ∣≤M1,we have
(45)∥C−PDpα1p,qzm−z∥L2Ω2≤M12Ymαp,q,α′p,q22m+1α′p,q+14m+1α′p,q−α1p,q+1Γm+1α′p,q+12.So by placing Equations (42) and (45) into (36), we have
(46)∥Em∥L2Ω≤MYmαp,q,α′p,q2m+1αp,q+14m+1α′p,q−α1p,q+1Γm+1αp,q+1+M1Ymαp,q,α′p,q2m+1α′p,q+14m+1α′p,q−α1p,q+1Γm+1α′p,q+1+θ2R+R1,where R=supp,q∈Ω∣sinz−v∣ and R1=supp,q∈Ω∣sinzm−vm∣. With a similar process, the following inequality can be considered for the function E′m that is defined as follows:
(47)C−PDqα2p,qvp,q−ξ2Dpα2p,qC−Pvp,q−sinz−v−h1p,q=C−PDqα2p,qvmp,q−ξ2Dpα2p,qC−Pvmp,q−sinzm−vm−h1p,q+E′m,∥E′m∥L2Ω≤MYmαp,q,α′p,q2m+1αp,q+14m+1α′p,q−α2p,q+1Γm+1αp,q+1+ξ2M1Ymαp,q,α′p,q2m+1α′p,q+14m+1α′p,q−α2p,q+1Γm+1α′p,q+1+R+R1.


## 4. Description of the Suggested Method

This section describes the suggested method to obtain solutions of the equations considered in (3) with the conditions (6) and (8). Then, we put relations (26)–(30) into Equation (3), and the following results are obtained:
(48)ΨTpZAM1A−1Ψq−ΨTpA−1TM3ATVΨq+θ2ΨTpsinZ−VΨq−h1p,q≃0,(49)ΨTpZAM1A−1−A−1TM3ATV+θ2sinZ−VΨq−h1p,q≃0,(50)ΨTpVAM2A−1Ψq−ξ2ΨTpA−1TM4ATVΨq−ΨTpsinZ−VΨq−h2p,q≃0,ΨTpVAM2A−1−ξ2A−1TM4ATV−sinZ−VΨq−h2p,q≃0,where sinz−v≃sinzm−vm=∑i=0m∑j=0msinzij−vijT^ipT^jq=ΨTpsinZ−VΨq. Considering the formulas lkp≃ΨTpLk,k=0,1,2,3 and fkq≃FkTΨq,k=0,1,2,3, Fk,Lk are known. Then, by substituting functions lkp,fkq,k=0,1,2,3, and Equations (22) and (23) in the conditions introduced in Equations (6) and (4), we obtain
(51)ΨTpZΨ0−ΨTpL0≃0,ΨTpZAM1A−1Ψ0−ΨTpL1≃0,(52)ΨTpVΨ0−ΨTpL2≃0,ΨTpVAM2A−1Ψ0−ΨTpL3≃0,(53)ΨTaZΨq−F0TΨq≃0,ΨTbZΨq−F1TΨq≃0,(54)ΨTaVΨq−F2TΨq≃0,ΨTbVΨq−F3TΨq≃0.To find solutions of problems (3) which are z,v, we solve the Equations (49)–(50) corresponding to (51)–(54) in collocation points to obtain functions Z,V.


## 5. Numerical Examples

In this section, we show the numerical results of the suggested method on four examples to illustrate the efficiency and accuracy of our scheme. Here, obtained numerical results of some test problems are shown to illustrate the computational efficiency and applicability of the reported method. All numerical examples are done by applying MATLAB (R2020b) software. Also, the convergence order (C-order) of the depicted method is calculated by the following relation:
(55)C‐order=∣logE1/E2∣logm2+1m2+1/m1+1m1+1.Here,E1 and E2 show the first and second values of the maximum absolute error (MAE) provided by the proposed method, respectively. Moreover, mi+1mi+1 for i=1,2 are the number of the basis functions used in the first and second implementations, respectively.Example 7.
We consider the following coupled nonlinear variable-order time-fractional sine-Gordon equations:(56)C−PDqα1p,qzp,q−C−PDpα1p,qzp,q+14sinz−v=h1p,q,C−PDqα2p,qvp,q−Dpα2p,qC−Pvp,q−sinz−v=h2p,q,where
(57)h1p,q=q2−α1p,qE1,3−α1p,q−q+e−qsinp+14sine−qsinp−cosp,h2p,q=q2−α2p,qE1,3−α2p,q−q+e−qcosp−sine−qsinp−cosp,where p,q∈0,1×0,1 and Eζ,ηq=∑k=0∞qk/Γkζ+η. For this example, the exact solutions to the case where υ=0,α1p,q=α2p,q=2 are zp,q=sinpe−q,vp,q=cospe−q. The numerical results obtained using the proposed method for α1p,q=1.5+0.25sinpq,α2p,q=1.75+0.25sinpq are shown in Figure 1. An image of the absolute error for value m=6 is shown in Figure 2, and numerical results for absolute error for different values m are shown in Tables 1 and 2. Table 3 shows comparison between the maximum absolute error (MAE) given for different values of m based on the present and the Chebyshev cardinal function (CCF) technique [69]. In other words, the comparison between our numerical method result and the numerical method result in [69] is done using two different operators, i.e., the Caputo-Prabhakar derivative and the Caputo derivative. From the obtained numerical results in Table 3, we see that the proposed method is an effective and good tool for solving this problem.Example 8.
We consider the following coupled nonlinear variable-order time-fractional sine-Gordon equations:(58)C−PDqα1p,qzp,q−C−PDpα1p,qzp,q+sinz−v=h1p,q,C−PDqα2p,qvp,q−Dpα2p,qC−Pvp,q−sinz−v=h2p,q,where
(59)h1p,q=−q3−α1p,qE2,4−α1p,q−q2+sinqcoshp+sinsinqcoshp−sinhp,h2p,q=−q3−α2p,qE2,4−α2p,q−q2+sinqsinhp−sinsinqcoshp−sinhp,where p,q∈0,1×0,2. For this example, the exact solutions to the case where υ=0,α1p,q=α2p,q=2 are zp,q=sinqcoshp,vp,q=sinqsinhp. The numerical results obtained using the proposed method for α1p,q=1.5+0.25p2q,α2p,q=1.75+0.25sinpq are shown in Figure 3. An image of the absolute error for value m=6 is shown in Figure 4, and numerical results for the absolute error for different values m are shown in Tables 4 and 5. The values of the MAE, C-order, and CPU time of the presented method are listed in Table 6.Example 9.
We consider the following coupled nonlinear variable-order time-fractional sine-Gordon equations:(60)C−PDqα1p,qzp,q−C−PDpα1p,qzp,q+sinz−v=h1p,q,C−PDqα2p,qvp,q−14Dpα2p,qC−Pvp,q−sinz−v=h2p,q,where
(61)h1p,q=p3p−2q2−α1p,qE2,3−α1p,q−q2+6p2p−2cosq+sinp32−pcosq−sinq,h2p,q=p3p−2q3−α2p,qE2,4−α1p,q−q2+32p2p−2sinq−sinp32−pcosq−sinq,where p,q∈0,2×0,1. For this example, the exact solutions to the case where υ=0,α1p,q=α2p,q=2 are zp,q=p32−pcosq,vp,q=p32−psinq. The numerical results obtained using the proposed method for α1p,q=1.5+0.25sinpq,α2p,q=1.75+0.25sinpq are shown in Figure 5. An image of the absolute error for value m=6 is shown in Figure 6, and numerical results for absolute error for different values m are shown in Tables 7 and 8. The values of the MAE, C-order, and CPU time of the presented method are listed in Table 9.Example 10.
We consider the following coupled nonlinear variable-order time-fractional sine-Gordon equations:(62)C−PDqα1p,qzp,q−C−PDpα1p,qzp,q+4sinz−v=h1p,q,C−PDqα2p,qvp,q−9Dpα2p,qC−Pvp,q−sinz−v=h2p,q,where
(63)h1p,q=Γ7/2Γ7/2−α1p,qp2q5/2−α1p,q−2q5/2+4sinp2−p3q5/2,h2p,q=Γ7/2Γ7/2−α2p,qp3q5/2−α1p,q−54pq5/2−sinp2−p3q5/2,where p,q∈0,1×0,1. For this example, the exact solutions to the case where υ=0,α1p,q=α2p,q=2 are zp,q=p2q5/2,vp,q=p3q5/2. The numerical results obtained using the proposed method for α1p,q=1.65+0.25sinpq,α2p,q=1.55+0.35sinpq are shown in Figure 7. An image of the absolute error for value m=6 is shown in Figure 8, and numerical results for absolute error for different values m are shown in Tables 10–12. Table 12 shows the comparison between the maximum absolute error (MAE) given for different values of m based on the present and the Chebyshev cardinal function (CCF) technique [69]. In other words, the comparison between our numerical method result and the numerical method result in [69] is done with using two different operators, i.e., the Caputo-Prabhakar derivative and the Caputo derivative.Figure 1
Numerical solution ofzp,q,vp,q with ρ=ω=υ=1 for Example 7 when m=6.Figure 2
Charts of absolute errorz (a) and v (b) with ρ=ω=υ=1 for Example 7 when m=6.Table 1
The absolute error for different valuesm for Example 7 at various points p,q.
p,q∣zm−z∣m=6m=7m=8m=90.1,0.18.9329e−176.2530e−167.1463e−168.0396e−160.2,0.21.6084e−161.1259e−151.2867e−151.4476e−150.3,0.32.1649e−161.5154e−151.7319e−151.9484e−150.4,0.42.5813e−161.8069e−152.0650e−152.3232e−150.5,0.52.8755e−162.0128e−152.3004e−152.5880e−150.6,0.63.0643e−162.1450e−152.4515e−152.7579e−150.7,0.73.1635e−162.2144e−152.5308e−152.8471e−150.8,0.83.1874e−162.2312e−152.5499e−152.8687e−150.9,0.93.1493e−162.2045e−152.5195e−152.8344e−151,13.0612e−162.1428e−152.4489e−152.7550e−15Table 2
The absolute error for different valuesm for Example 7 at various points p,q.
p,q∣vm−v∣m=6m=7m=8m=90.1,0.18.9031e−166.2322e−157.1225e−158.0128e−150.2,0.27.9349e−165.5544e−156.3479e−157.1414e−150.3,0.36.9986e−164.8990e−155.5989e−156.2988e−150.4,0.46.1054e−164.2738e−154.8843e−155.4949e−150.5,0.55.2636e−163.6845e−154.2109e−154.7373e−150.6,0.64.4792e−163.1354e−153.5833e−154.0312e−150.7,0.73.7558e−162.6291e−153.0047e−153.3803e−150.8,0.83.09572e−162.1670e−152.4765e−152.7861e−150.9,0.92.4991e−161.7494e−151.9993e−152.2492e−151,11.9655e−161.3759e−151.5724e−151.7690e−15Table 3
Comparison of the numerical results for Example7.
mα1p,q=1.5+0.25sinpq,α2p,q=1.75+0.25sinpqMAEzMAEvOur method62.1048e−161.328e−1674.2097e−182.6561e−1786.3146e−193.9842e−1996.3146e−213.9842e−20MAuMAwSecond-kind CCF method [69]61.4597e−073.4301e−0772.7615e−091.2573e−0881.5434e−105.3747e−1094.0590e−122.0088e−11Figure 3
Numerical solution ofzp,q,vp,q with ρ=ω=υ=1 for Example 8 when m=6.Figure 4
Charts of absolute errorz (a) and v (b) with ρ=ω=υ=1 for Example 8 when m=6.Table 4
The absolute error for different valuesm for Example 8 at various points p,q.
p,q∣zm−z∣m=6m=7m=8m=90.2,0.21.9744e−161.3821e−151.5795e−151.7770e−150.4,0.43.9281e−162.7497e−153.1425e−153.5353e−150.6,0.65.8368e−164.0857e−154.6694e−155.2531e−150.8,0.87.6689e−165.3682e−156.1351e−156.9020e−151,19.3832e−166.5682e−157.5065e−158.4449e−150.2,1.21.0926e−157.6483e−158.7409e−159.8336e−150.4,1.41.2231e−158.5621e−159.7853e−151.1008e−140.6,1.61.3220e−159.2540e−151.0576e−141.1898e−140.8,1.81.3801e−159.6607e−151.1040e−141.2420e−141,21.3875e−159.7127e−151.1100e−141.2487e−14Table 5
The absolute error for different valuesm for Example 8 at various points p,q.
p,q∣vm−v∣m=6m=7m=8m=90.2,0.21.9678e−171.3775e−161.5743e−161.7711e−160.4,0.47.7532e−175.4272e−166.2026e−166.9779e−160.6,0.61.7003e−161.1902e−151.3602e−151.5303e−150.8,0.82.9138e−162.0396e−152.3310e−152.6224e−151,14.3361e−163.0353e−153.4689e−153.9025e−150.2,1.25.8679e−164.1075e−154.6943e−155.2811e−150.4,1.47.3924e−165.1746e−155.9139e−156.6531e−150.6,1.68.7786e−166.1450e−157.0229e−157.9007e−150.8,1.89.8856e−166.9199e−157.9085e−158.8970e−151,21.0567e−157.3971e−158.4538e−159.5106e−15Table 6
The MA-error,C-order, and CPU time of the presented method in Example 8.
mzp,qvp,qCPU‐time (s)MAEC‐orderMAEC‐order61.1801e−15—1.0547e−15—6.1577.5127e−165.70337.3971e−1607.672618.1481.0100e−1615.02606.4538e−1611.314644.1391.0487e−1716.95709.5106e−1711.402587.14Figure 5
Numerical solution ofzp,q,vp,q with ρ=ω=υ=1 for Example 9 when m=6.Figure 6
Charts of absolute errorz (a) and v (b) with ρ=ω=υ=1 for Example 9 when m=6.Table 7
The absolute error for different valuesm for Example 9 at various points p,q.
p,q∣zm−z∣m=6m=7m=8m=90.2,0.21.4168e−179.9182e−171.1335e−161.2751e−160.4,0.49.9243e−176.9470e−167.9395e−168.9319e−160.6,0.62.8568e−161.9997e−152.2854e−152.5711e−150.8,0.85.5961e−163.9172e−154.4769e−155.0365e−151,18.6783e−166.0748e−156.9426e−157.8104e−151.2,0.21.1282e−157.8978e−159.0261e−151.0154e−141.4,0.41.2452e−158.7167e−159.9619e−151.1207e−141.6,0.61.1288e−157.9016e−159.0304e−151.0159e−141.8,0.87.1699e−165.0189e−155.7359e−156.4529e−15Table 8
The absolute error for different valuesm for Example 9 at various points p,q.
p,q∣vm−v∣m=6m=7m=8m=90.2,0.21.4216e−189.9514e−181.1373e−171.2794e−170.4,0.42.0117e−171.4082e−161.6094e−161.8105e−160.6,0.68.8372e−176.1860e−167.0697e−167.9535e−160.8,0.82.3660e−161.6562e−151.8928e−152.1294e−151,14.7409e−163.3186e−153.7927e−154.2668e−151.2,0.27.7188e−165.4032e−156.1751e−156.9470e−151.4,0.41.0488e−157.3419e−158.3908e−159.4397e−151.6,0.61.1622e−158.1358e−159.2980e−151.0460e−141.8,0.89.0352e−166.3246e−157.2281e−158.1316e−15Table 9
The MAE,C-order, and CPU time of the presented method in Example 9.
mzp,qvp,qCPU‐time (s)MAEC‐orderMAEC‐order67.2352e−15—1.1522e−15—8.1774.7067e−152.70238.1258e−166.670620.0481.9519e−1612.02507.2880e−1610.315646.1096.8004e−1713.95609.4297e−1710.401589.04Figure 7
Numerical solution ofzp,q,vp,q with ρ=ω=υ=1 for Example 10 when m=6.Figure 8
Charts of absolute errorz (a) and v (b) with ρ=ω=υ=1 for Example 10 when m=6.Table 10
The absolute error for different valuesm for Example 10 at various points p,q.
p,q∣zm−z∣m=6m=7m=8m=90.1,0.13.1271e−202.1890e−192.50179e−192.8144e−190.3,0.34.3872e−183.0710e−173.5098e−173.9485e−170.5,0.54.3703e−173.0592e−163.4962e−163.9332e−160.7,0.71.9865e−161.3905e−151.5892e−151.7878e−150.9,0.96.1551e−164.3086e−154.9241e−155.5396e−15Table 11
The absolute error for different valuesm for Example 10 at various points p,q.
p,q∣vm−v∣m=6m=7m=8m=90.1,0.13.1271e−212.1890e−202.5017e−202.8144e−200.3,0.31.3161e−189.2132e−181.0529e−171.1845e−170.5,0.52.1851e−171.5296e−161.7481e−161.9666e−160.7,0.71.3905e−169.7338e−161.1124e−151.2514e−150.9,0.95.5396e−163.8777e−154.4317e−154.9856e−15Table 12
Comparison of the numerical results for Example10.
mα1p,q=1.65+0.25sinpq,α2p,q=1.55+0.35sinpqMAEzMAEvOur method63.1131e−182.1440e−1971.3161e−187.2122e−2082.1941e−191.4596e−2091.2605e−198.5238e−21MAEuMAEwSecond-kind CCF method [69]66.5157e−041.7826e−0574.1287e−049.8284e−0682.7169e−045.8208e−0692.0393e−042.9727e−06


## 6. Conclusion

This paper is focused on a numerical method based on the orthogonal polynomials known as the shifted Chebyshev polynomials of degreem. The fractional sine-Gordon equations of variable order which includes a fractional derivative in terms of the Caputo-Prabhakar sense in this paper are studied. Using this numerical method, we obtain the numerical solutions of this type of equations. In this article, the convergence and error analysis for numerical solutions are studied. Some numerical examples are shown for the application of numerical methods.



---

*Source: 8830727-2021-03-09.xml*