# Pell Collocation Method for Solving the Nonlinear Time–Fractional Partial Integro–Differential Equation with a Weakly Singular Kernel

Authors: M. Taghipour; H. Aminikhah
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

This article focuses on finding the numerical solution of the nonlinear time–fractional partial integro–differential equation. For this purpose, we use the operational matrices based on Pell polynomials to approximate fractional Caputo derivative, nonlinear, and integro–differential terms; and by collocation points, we transform the problem to a system of nonlinear equations. This nonlinear system can be solved by the fsolve command in Matlab. The method’s stability and convergence have been studied. Also included are five numerical examples to demonstrate the veracity of the suggested strategy.


---

## Body

## 1. Introduction

Nowadays, fractional partial differential equations (FPDEs) have emerged as one of the most crucial issues due to their vast applications in various branches of science, such as medicine [1, 2], control theory [3–5], engineering [6, 7], viscoelasticity [8], mathematical physics [9], geo–hydrology [10], signals [11], stochastic models [12], electrical engineering, [13], and financial economics [14]. Due to the fact that analytical solutions of FPDEs are rarely available, the use of numerical methods is inevitable. Hitherto, a number of numerical methods for FPDEs have been suggested, such as finite difference [15, 16], spectral methods [17–21], homotopy methods [22, 23], and finite element [24, 25]. The nonlinear FPDEs have been extensively analyzed using numerical methods. Dehghan et al. used the homotopy analysis method to construct a scheme for solving the fractional KdV equation [26]. Nikan et al. proposed a meshless technique in order to solve the nonlinear fractional fourth–order diffusion equation [27]. Safari and Azarsa introduced a meshless method based on Muntz polynomials to solve nonlinear and linear space fractional partial differential equations [28]. Yaslan applied the Legender collocation method for solving nonlinear fractional partial differential equations [29].When it comes to solving differential equations, spectral approaches are extremely effective. The solution to the differential equation is sought as a series of basis polynomials using this method. The Galerkin, Tao, and collocation approximations are the most common spectral methods [30, 31]. For example, Samiee et al. [20] designed a Petrov–Galerkin spectral method for distributed–order PDEs. Agarwal et al. [32] suggested a spectral collocation approach for variable–order fractional integro–differential equations. In [33], the authors used the polynomial–sinc collocation method for solving distributed order fractional differential equations. Abbaszadeh et al. proposed a Crank–Nicolson Galerkin spectral method for distributed order weakly singular integro–partial differential equations [34].In the present paper, we offer a numerical technique for solving the nonlinear time–fractional partial integro–differential equation (TFPIDEs) with a weakly singular kernel(1)0CDηαVξ,η+Vξ,ηVξξ,η=∫0ηη−sβ−1Vξξξ,sds+Gξ,η,ξ∈0,L,η∈0,T,with initial and boundary conditions
(2)V0,η=Φ1η,VL,η=Φ2η,0<η≤T,(3)Vξ,0=Ψξ,0<ξ≤L,where0<α,β<1, gξ,η∈C0,L×0,T, and 0CDηα signify the fractional operator. This problem appears in the modeling of heat transfer materials with memory, population dynamics [35], and nuclear reaction theory [36].To the best of the author’s knowledge, little work has been done on problem (1). For example, Guo et al. [37] proposed a numerical technique for solving (1)–(3). In the case of α=1, Zheng et al. [35] described three semi–implicit compact finite difference schemes for problem (1)–(3). This encourages us to suggest a numerical scheme for the problem (1)–(3). The finite difference schemes are the easiest methods for solving these equations. It is, however, difficult to apply the mathematical study of finite difference methods to nonlinear TFPIDEs. Polynomial spectral techniques are effective tools for solving PDEs. To build spectral methods, many polynomials have been developed (see [38–41]). The coefficients of Pell polynomials are integers, and the number of terms increases slowly. This leads to less CPU time and fewer computational errors. Because of this, the Pell polynomials with both of these two characteristics will be employed.In this paper, we will focus on the spectral collocation method based on two–variable Pell polynomials. We use them as the basis polynomials to solve the main problem numerically. With the use of operational matrices, the problem is turned into a system of nonlinear equations in the approach based on these polynomials. The error analysis is presented. Several test problems are provided to illustrate the method’s efficacy.The following is the body of the article: Section2 introduces a number of key themes. To remedy the main problem, we suggest a polynomial spectral technique in Section 3. The error analysis is investigated in Section 4. Section 5 contains the experiments. The conclusion is addressed in Section 6.


## 2. Definitions

Definition 1 (see [15]).
The Riemann–Liouville integral of a functionzξ,η on 0,L×0,T is defined as follows
(4)0RLIηαzξ,η=1Γα∫0ηη−τα−1zξ,τdτ,0<α<1.Definition 2 (see [15]).
The Riemann–Liouville derivative of a functionzξ,η on 0,L×0,T is defined as follows
(5)0RLDηαzξ,η=1Γ1−αddη∫0ηη−τ−αzξ,τdτ,0<α<1.Definition 3 (see [15]).
The Caputo derivative of a functionzξ,η on 0,L×0,T is defined as follows
(6)0CDηαzξ,η=1Γ1−α∫0ηη−τ−αz′ξ,τdτ,0<α<1.With respect to these definitions, we have the following properties(7)0RLIηα0CDηαzξ,η=zξ,η−zξ,0,0<α<1,0.4cm,(8)0CDηαc=0,c=constant,(9)0CDηαηm=0,m=0,Γm+1Γm+1−αηm−α,m=1,2,⋯.The following relation can be used to construct Pell polynomials: [42]:
(10)Pl+2ξ=2ξPl+1ξ+Plξ,P0ξ=0,P1ξ=1.According to [42], Pnξ has the following form
(11)Plξ=∑k=0l−1/2l−k−1k2ξl−2k−1.We can represent a continuous functionVξ via the Pell polynomials as follows:
(12)Vξ≈VKξ=∑i=0Kv¯i+1Pi+1ξ=VTPKξ,where
(13)V=v¯1,v¯2,⋯,v¯K+1T,PKξ=P1ξ,P2ξ,⋯,PK+1ξT.Analogously, a functionVξ,η defined on 0,L×0,T may be described as follows:
(14)Vξ,η≈VKJξ=∑i=0K∑j=0Jv¯i+1,j+1Pi+1ξPj+1η=PKξTWPJη=VTPKJξ,η,andW denotes a matrix of suitable dimensions, as well as V and PKJξ,η are the following:
(15)V=v¯11,v¯12,⋯,v¯1J+1,v¯21,⋯,u¯2J+1,⋯,v¯K1,⋯,v¯K+1J+1T,(16)PKJξ,η=P11ξ,η,⋯,P1J+1ξ,η,P21ξ,η,⋯,P2J+1ξ,η,⋯,PK+11ξ,η,⋯,PK+1J+1ξ,ηT.We may also rephrasePKξ,PJη in the following way:
(17)PKξ=QξTξ,(18)PJη=QηTη,(19)Tξ=1,ξ,⋯,ξKT,Tη=1,η,⋯,ηJT,where
(20)Qξ=a0,0000⋯00a1,100⋯0a2,00a2,20⋯0⋮⋮⋮⋮⋮⋮aK,10aK,3⋯0aK,K,Qη=a0,0000⋯00a1,100⋯0a2,00a2,20⋯0⋮⋮⋮⋮⋮⋮aJ,10aJ,3⋯0aJ,J,with
(21)ai,j=i−i−j2i−j22i−2i−j2,ifi≥j,i,j=oddori,j=even0,otherwise.


## 3. Analysis of the Numerical Method

Here, we find several operational matrices with the help of Pell polynomials, which are useful in developing the suggested technique.To begin, we estimate the fractional operator as follows:(22)0CDηαVξ,η≈0CDηαVKJξ,η=0CDηαPKξTWPJη=PKξTW0CDηαPJη=PKξTW0CDηαQηTη,(23)=PKξTWQη0CDηαTη.Thanks to using relation (9), we obtain
(24)0CDηαVKJξ,η=PKξTWQη0,Γ2Γ2−αη1−α,Γ3Γ3−αη2−α,⋯,ΓJ+1ΓJ+1−αηJ−αT,(25)=PKξTWQη000⋯00Γ2Γ2−αη−α0⋯000Γ3Γ3−αη−α⋯0⋮⋮⋮⋮⋮000⋯ΓJ+1ΓJ+1−αη−α1ηη2⋮ηJ,(26)=PKξTWQηMαQη−1PJη,where
(27)Mα=diag0,Γ2Γ2−αη−α,Γ3Γ3−αη−α,⋯,ΓJ+1ΓJ+1−αη−α.(28)0CDηαVξ,η≈PKξTWQηMαQη−1PJη.Next, we approximate the nonlinear and integro–differential terms in equation (1). First, we compute PK′ξ and PK′′ξ.
(29)PK′ξ=QξT′ξ=Qξ012ξ⋮KξK−1=Qξ0000⋯0001000⋯0000200⋯0000030⋯000⋮⋮⋮⋮⋮⋮⋮⋮0000⋯0N01ξξ2⋮ξK,(30)=QξD′Tξ=QξD′Qξ−1PKξ,where
(31)D′=0000⋯0001000⋯0000200⋯0000030⋯000⋮⋮⋮⋮⋮⋮⋮⋮0000⋯0N0.Similarly, we obtain(32)PK′′ξ=QξT′′ξ=QξD′′Qξ−1PKξ,where D′′=di,j1≤i,j≤K+1.
(33)di,j=i−1i−2,i≥3,j=i−2,0,otherwise.If we use (30), we get
(34)Vξ,ηVξξ,η≈VKJξ,ηVKJxξ,η=PKξTWPJηPKξTQξ−TD′TQξTWPJη.For integro–differential term, using (32), we have
(35)∫0ηη−sβ−1Vξξξ,sds≈∫0ηη−sβ−1PK′′ξTWPJsds=∫0ηη−sβ−1QξD′′Qξ−1PKξTWPJsds,(36)=PKξTQξ−TD′′TQξTW∫0ηη−sβ−1QηTsds,(37)=PKξTQξ−TD′′TQξTWQη∫0η1,s,⋯,sJTη−s1−βds.On the other hand, the following relationship is valid:(38)∫0ηξkη−ξ1−βdξ=Γk+1ΓβΓk+β+1ηk+β,0<β<1,k=0,1,2,⋯.So, by substituting (38) into (37), we have
(39)∫0ηη−sβ−1Vξξξ,ηds≈PKξTQξ−TD′′TQξTWQηΓβΓβ+1ηβ,Γ2ΓβΓβ+2ηβ+1,⋯,ΓJ+1ΓβΓβ+J+1ηβ+JT,(40)=PKξTQξ−TD′′TQξTWQηΓβΓβ+100⋯00Γ2ΓβΓβ+20⋯000Γ3ΓβΓβ+3⋯0⋮⋮⋮⋮⋮000⋯ΓJ+1ΓβΓβ+J+1ηβηβ+1ηβ+2⋮ηβ+J,(41)=PKξTQξ−TD′′TQξTWQηSTβη,where
(42)S=diagΓβΓβ+1,Γ2ΓβΓβ+2,⋯,ΓJ+1Γ1−βΓβ+J+1,(43)Tβη=ηβ,ηβ+1,⋯,ηβ+J.Hence, using relations (28), (34), and (41), as a result, we get
(44)R1ξ,η=PKξTWQηMαQη−1PJη+PKξTWPJηPKξTQξ−TD′TQξTWPJη−PKξTQξ−TD′′TQξTWQηSTβη−Gξ,η≈0,R2ξ=PKξTWPJ0−Ψξ≈0,R3η=PK0TWPJη−Φ1η≈0,R4η=PKLTWPJη−Φ2η≈0.Now, from relation (44), we create the nonlinear system below.
(45)R1ξi,ηj≈00≤i≤K−2,0≤j≤J−1,R2ξi≈0i=0,1,⋯,K,R4ηj≈0R3ηj≈0j=0,1,⋯,J−1,j=0,1,⋯,J−1,whereξi=2i+1/2K+2 and ηj=2j+1/2J+2.By solving this system, the unknown matrixW can be determined. It is worth noting that we have used the fsolve command in Matlab.


## 4. Convergence

Here, we prove that the numerical scheme for solving (1)–(3) is convergent, and we follow references [43, 44]. We assume that
(46)G=spanP1ξ,P2ξ,⋯,PK+1ξ,(47)Q=spanP1η,P2η,⋯,PJ+1η,(48)Gx=spanP1′ξ,P2′ξ,⋯,PK+1′ξ,(49)Gxx=spanP1′′ξ,P2′′ξ,⋯,PK+1′′ξ.(50)Vξ,η−VKJξ,η2≤H1LTL+TK+J+1K+J+1!,(51)Vξξ,η−VKJξξ,η2≤H2LTL+TK+J+1K+J+1!,(52)Vξξξ,η−VKJξξξ,η2≤H3LTL+TK+J+1K+J+1!,Theorem 4.
LetVξ,η∈CK+J+30,L×0,T and VKJξ,η, VKJξξ,η, and VKJξξξ,η be the best approximations of Vξ,η, Vξξ,η, and Vξξξ,η in the spaces G×Q, Gξ×Q, and Gξξ×Q, respectively. The following inequalities are true.whereHi=maxk=0,⋯,J+K+iVkξ,η.
(53)Vξ,η=V0,0+Vξ0,0ξ+Vη0,0η+12!Vξξ0,0ξ2+2Vξη0,0ξη+Vηη0,0η2+⋯,(54)+1n!∑i+j=J+KnjVξiηj0,0ξiηj+RJ+K,0,0ξ,η=PJ+K,0,0ξ,η+RJ+K,0,0ξ,η,where
(55)RJ+K,0,0ξ,η=1J+K!∑i+j=J+K+1J+K+1jξiηj∫01Vξiηjmξ,mη1−mJ+Kdm.Proof.
Using Taylor expansion, we have [45].According to the best approximation theorem(56)Vξ,η−VJKξ,η2≤Vξ,η−PJ+K,0,0ξ,η2,(57)=1J+K!∑i+j=J+K+1J+K+1jξiηj∫01Vξiηjmξ,mη1−mK+Kdm2,(58)=∫0T∫0L1J+K!∑i+j=J+K+1J+K+1jξiηj∫01Vξiηjmξ,mη1−mJ+Kdm2dξdη1/2≤∫0T∫0LH1J+K+1!∑i+j=J+K+1J+K+1jξiηj2dξdη1/2=H1J+K+1!∫0T∫0Lξ+η2K+2J+2dξdη1/2≤H1LTL+TJ+K+1J+K+1!.Similarly, other inequalities can also be proved.Theorem 5.
LetVKJξ,η is the exact solution and V~KJξ,η=V~TPKJξ,η is the approximation solution of equations (45). Then, one has
(59)VNMξ,η−V~NMξ,η2≤V−V~2LT∑r=0K∑s=0J∑k=0r2r−kk2Lr−2k2∑k=0s2s−kk2Ts−2k2.Proof.
We have(60)VKJξ,η−V~KJξ,η22=∫0T∫0LVKJξ,η−V~KJξ,η2dξdη=∫0T∫0L=∑r=0K∑s=0Jv¯r+1,s+1−v¯~r+1,s+1Pr+1ξPs+1η2dξdη,(61)≤∫0T∫0L∑r=0K∑s=0Jv¯r+1,s+1−v¯~r+1,s+12∑r=0K∑s=0JPr+1ξPs+1η2dξdη,(62)=∑r=0K∑s=0Jv¯r+1,s+1−v¯~r+1,s+12∑r=0K∑s=0J∫0T∫0LPr+1ξPs+1η2dξdη,(63)=V−V~22∑r=0K∑s=0J∫0T∫0L∑k=0r2r−kk2r−2kξr−2k2∑k=0s2s−kk2s−2kηs−2k2dξdη≤V−V~22∑r=0K∑s=0J∫0T∫0L∑k=0i2r−kk2r−2kLr−2k2∑k=0s2s−kk2d−2kTs−2k2dξdη=V−V~22LT∑r=0K∑s=0J∑k=0r2r−kk2r−2kLr−2k2∑k=0s2s−kk2s−2kTs−2k2.
So that(64)VKJξ,η−V~KJξ,η2≤V−V~2LT∑r=0K∑s=0J∑k=0r2r−kk2Lr−2k2∑k=0s2s−kk2Ts−2k2.
Now, we prove that the presented numerical method is convergent.Theorem 6.
LetSKJξ,η be the perturbation term and VKJξ,η be the approximate solution to the main problem derived using the proposed approach. Then, the perturbation term tends to zero as K,J⟶∞.Proof.
Thanks to (9), we deduce that
(65)Vξ,η=Ψξ−0RLIηαVξ,ηVξξ,η+0RLIηα∫0ηη−sβ−1Vξξξ,sds+0RLIηαGξ,η.Assume thatVNMξ,η is an approximate solution of the above equation. It means that
(66)VKJξ,η=Ψξ−0RLIηαVKJξ,ηVKJξξ,η+0RLIηα∫0ηη−sβ−1VKJξξξ,sds+0RLIηαGξ,η+SKJξ,η,where SKJ is the perturbation term. From equations (65) and (66), we have
(67)EKJξ,η=−0RLIηαVξ,ηVxξ,η−VKJξ,ηVKJξξ,η+0RLIηα∫0ηη−sβ−1Vξξξ,η−VKJξξξ,sds−SKJξ,η=−0RLIηαVKJξ,ηVξξ,η−VKJξξ,η+Vξξ,ηVξ,η−VKJξ,η+0RLIηα∫0ηη−sβ−1Vξξξ,η−VKJξξξ,sds−SKJξ,η=−0RLIηαVKJξ,ηEKJξξ,η+Vξξ,ηEKJξ,η+0RLIηα∫0ηη−sβ−1EKJξξξ,ηds−SKJξ,η=−0RLIηαVKJξ,ηEKJξξ,η+Vξξ,η−VKJξξ,η+VKJξη,ηEKJξ,η+0RLIηα∫0ηη−sβ−1EKJξξξ,ηds−SKJξ,η=−0RLIηαVKJξ,ηEKJξξ,η+VKJξξ,ηEKJξ,η+EKJξξ,ηEKJξ,η+0RLIηα∫0ηη−sβ−1EKJξξξ,ηds−SKJξ,η,where EKJξ,η=Vξ,η−VKJξ,η. So that
(68)SKJξ,η2≤RL0IηαVKJξ,ηEKJξξ,η+VKJξξ,ηENMξ,η+EKJξξ,ηENMξ,η2+EKJξ,η2+RL0Iηα∫0ηη−sβ−1ENMξξξ,ηds2=I12+ENMξ,η2+I2.Now, we computeI12 and I22.
(69)I222=1Γα2∫0T∫0L∫0η∫0τη−τα−1τ−sβ−1EKJξξξ,sdsdτ2dξdη≤1Γα2∫0T∫0L∫0η∫0τη−τα−1τ−sβ−1eKJξξξ,sdsdτ2dξdη≤1Γα2∫0T∫0L∫0η∫0τη−τ2α−2τ−s2β−2dsdτ∫0η∫0τEKJξξξ,s2dsdτdξdη≤TΓα2∫0T∫0L∫0η∫0τη−τ2α−2τ−s2β−2dsdτ∫0TEKJξ,η2dsdξdη≤TΓα2∫0T∫0LT2β−1×t2α−1∫0TEKJξ,η2dsdξdη=T2βΓα2∫0Tη2α−1dη∫0T∫0LEKJξξξ,s2dξds=T2β+2αΓα22αEKJξξξ,η22.Similarly(70)I122≤T2αΓα22αVKJξ,ηEKJξξ,η+VKJξξ,ηEKJξ,η+EKJξη,ηEKJξ,η22.Therefore, we conclude(71)SKJξ,η2≤T2αΓα22αVKJξ,ηEKJξξ,η+VKJξξ,ηEKJξ,η+EKJξξ,ηEKJξ,η2+EKJξ,η2+T2β+2αΓα22αEKJξξξ,η2≤C1VKJξ,η2EKJξξ,η2+C1VKJξξ,η2EKJξ,η2+C1EKJξξ,η2EKJξ,η2+EKJξ,η2+C2EKJξξξ,η2≤C1H2LTL+Tn+1n+1!VKJξ,η2+C1H1LTL+Tn+1n+1!VKJxξ,η2+C1H1H2L2T2L+T2n+2n+12+H1LTL+Tn+1n+1+C2H3LTL+Tn+1n+1,where n=K+J. Now SKJξ,η⟶0 as K,J⟶∞.Finally, we provide a theorem about the convergence of the series of Pell polynomials. We follow Atta et al. [38, 40], Abd-Elhameed and Youssri [43, 44], and Youssri [46]. According to [47, 48], a square-integrable function Vξ on 0,1 has the following Pell expansion
(72)Vξ=∑m=0∞v¯m+1Pm+1ξ,where
(73)v¯m=m∑s=0∞−1sb2s+m−122s+m−1s+m2s+m−1s,bs=Vs0s!.Lemma 7 (see [46]).
LetIuz denote the modified Bessel function of order u of the first kind. The following identity holds:
(74)∑s=0∞zss!s+m+1!=t−m+1/2Im+12z.Lemma 8 (see [46]).
The modified Bessel function of the first kindIuz satisfies the following inequality:
(75)Iuz≤zucoshz2uΓu+1.Theorem 9.
SupposeVx∈L20,1, Vi0≤Li,i≥0, where L is a positive constant and Vx=∑m=0∞v¯m+1Pm+1x. Then
(76)v¯m+1≤LmcoshL2mm!,and the approximate solution series converges absolutely. Also, ifEKξ=∑m=K+1∞v¯m+1Pm+1ξ, the following error estimation is satisfied
(77)EKξ≤coshLeL/2L2K+1K+1!.Proof.
From (73), Lemma (7), and Lemma (8), we have
(78)v¯m+1=m+1∑s=0∞−1sV2s+m022s+ms+m+12s+m2s+ms=m+1∑s=0∞−1sV2s+m022s+ms+m+1ss+m=m+1∑s=0∞−1sV2s+m022s+mss+m+1≤m+1∑s=0∞L2s+m22s+mss+m+1=s+1×L2m∑s=0∞L2/22sss+m+1≤m+1×L2m×L222−m+1/2×Im+12L2=2m+1LIm+1L≤2m+1L×Lm+1coshL2m+1Γm+2=LkcoshL2mm.
As a result, the first portion of the theorem is established.We now go on to the second section.(79)v¯m+1Pm+1ξ≤LmcoshL2mmPm+11,additionally
(80)∑m=0∞Lm2mmPm+11=∑m=0∞Lm2mm×1+2m+1−1−2m+122=∑m=0∞1+222L1+2/2mm−∑m=0∞1−222L1−2/2m!m=1+222eL1+2/2−1−222eL1−2/2.As a result, we can conclude that the series is absolutely convergent using the comparison test.For the third part, we have(81)EKξ≤coshL∑m=K+1∞Lm2mm=coshL∑m=K+1∞L/2mm=coshLγK+1,L/2ΓK+1eL/2,where γK+1,L/2 is the lower incomplete gamma function [49], then
(82)EKξ≤coshLeL/2ΓK+1∫0L/2ηKe−ηdη≤coshLeL/2L2k+1ΓK+1.


## 5. Numerical Experiments

Five test problems are offered in this part to demonstrate the correctness and validity of the presented method. On a Windows 10 (64 bit) Intel(R) Core(TM) i7-7500U CPU operating at 2.70 GHz with 8.0 GB of RAM, all computations are done with Matlab R2020b software. In all examples, we use theL∞ error norm and L2 error norm
(83)e∞=max0≤r≤K0≤s≤JVξr,ηs−VKJξr,ηs,e2=h∑s=0nesm21/2,Example 1.
Consider the following nonlinear time–fractional partial integro–differential equation on0,1×0,1 with the exact solution Vξ,η=η3sinπξ(84)0CDηαVξ,η+Vξ,ηVξξ,η=∫0ηη−sβ−1Vξξξ,sds+Gξ,η,with conditions (2), (3), and
(85)Gξ,η=6η3−αΓ4−α+6π2Γβη3+βΓ4+βπη6cosπξsinπξ.We solved this problem numerically using the polynomial spectral scheme provided in this paper. We employed fsolve in Matlab to solve a nonlinear system of equations (45). Table 1 shows the absolute errors for α=0.5 and various β values. We can observe from this table that the recommended strategy is effective. Also, we portrayed the numerical solution and absolute error surfaces in Figure 1. Furthermore, the norm of errors and CPU times is reported in Table 2. Table 1, Table 2, and Figure 1 show that the numerical method provides acceptable results.Example 2.
Consider the following equation(86)0CDηαVξ,η+Vξ,ηVξξ,η=∫0ηη−sβ−1Vξξξ,sds+Gξ,η,with conditions (2) and (3). The source term is taken as
(87)Gξ,η=2ξ1−ξη2−αΓ3−α+ξ1−2ξ1−ξη4+4ηβ+2β3+3β2+2β,and Vξ,η=η2ξ1−ξ.Table 1
Numerical reports in Example1.
ξi,ηiα=0.5,K=11α=0.5,K=11α=0.5,K=9α=0.5,K=9β=0.1,J=4β=0.3,J=4β=0.7,J=4β=0.9,J=40.1,0.15.2665e−117.5026e−119.4157e−101.5935e−090.2,0.28.8713e−107.0142e−101.7084e−094.5505e−090.3,0.34.6062e−093.6713e−095.8173e−081.3764e−090.4,0.41.4766e−081.2353e−083.6620e−071.2902e−070.5,0.53.6426e−083.1853e−081.2736e−066.7104e−070.6,0.67.6243e−086.9345e−083.4065e−062.2596e−060.7,0.71.4258e−071.3432e−077.7630e−066.0393e−060.8,0.82.4555e−072.3833e−071.5689e−051.3673e−050.9,0.93.8165e−073.7801e−072.5814e−052.4145e−051,17.3344e−151.1944e−161.6793e−152.1690e−14Figure 1
Pictorial results in Example1 with α=0.5,β=0.1, and N=11,M=4.Table 2
Norm of errors forα=β=0.5 and CPU time in Example 1.
J=4,K=6J=4,K=7J=4,K=8J=4,K=9J=4,K=10e∞e∞e∞e∞e∞1.7427e−031.7093e−033.7497e−053.7046e−055.5452e−07CPU0.7787s1.3037s1.4066s1.5607s1.9541sThe absolute errors for equal values ofα,β, and K=J=3 are illustrated in Table 3. This table shows quite revealingly that the expressed method has good precision. In addition, we can observe from the table that only a small number of basis functions have produced the necessary outcomes. The CPU times and the norm of errors are provided in Table 4. Figure 2 shows a visualization of the approximate solution as well as absolute errors.Example 3.
Consider the following equation on0,1×0,1:
(88)0CDηαVξ,η+Vξ,ηVξξ,η=∫0ηη−sβ−1Vξξξ,sds+Gξ,η,with conditions (2) and (3). The source term is taken as
(89)Gξ,η=600ξ31−ξη3−αΓ4−α−24η4−αΓ5−α+10000η61−η23−4ξξ2ξ3−ξ4−600ξ1−2ξ6Γβηβ+3Γ4+β−24Γβηβ+4Γ5+β.Table 3
Numerical results in Example2.
ξi,ηiα=0.5α=0.7α=0.9α=1β=0.5β=0.5β=0.5β=0.50,01.3878e−171.3878e−171.3878e−176.9389e−170.1,0.12.9677e−152.4715e−152.1732e−152.4820e−150.2,0.21.3601e−141.0819e−148.4386e−158.1749e−150.3,0.34.2369e−143.3196e−142.4869e−142.1178e−140.4,0.41.0748e−138.5140e−146.4282e−145.2874e−140.5,0.52.2699e−131.8235e−131.4061e−131.1730e−130.6,0.64.0491e−133.2892e−132.5840e−132.2102e−130.7,0.76.1083e−135.0086e−133.9917e−133.4946e−130.8,0.87.5499e−136.2282e−135.0181e−134.4802e−130.9,0.96.5503e−135.4298e−134.4097e−134.0014e−131,17.5493e−163.3867e−161.0991e−151.0437e−15Table 4
Norm of errors forα=β=0.5 and CPU time in Example 2.
J=3,K=3J=3,K=4J=3,K=5J=3,K=6J=3,K=7e∞e∞e∞e∞e∞1.9604e−124.9638e−128.0466e−111.3686e−131.4568e−12CPU0.7492s0.9354s1.1262s1.2330s1.2589sFigure 2
Pictorial results in Example2 for α=0.5,β=0.5 and N=M=3.The exact solution isVξ,η=100ξ3η31−ξ1−η.The numerical results are reported in Table5. We have chosen K=4,J=5, and for equal α and β, the obtained results are fruitful. This table confirms that the presented method has high performance and produces accurate results. For α=β=0.5 and different J and K, the norm of errors and the CPU times are provided in Table 3. The absolute error functions for equal α and β are sketched in Figure 3. These figures show that the numerical and exact solutions are almost identical (Table 6).Example 4.
Consider the following equation:(90)0CDηαVξ,η+Vξ,ηVξξ,η=∫0ηη−sβ−1Vξξξ,sds+Gξ,η,with conditions (2), (3), and
(91)Gξ,η=6eξη3−αΓ4−α+η6e2ξ−6eξΓβηβ+3Γ4+β.Table 5
Numerical results in Example3.
ξi,ηiα=0.1α=0.4α=0.6α=0.8β=0.1β=0.4β=0.6β=0.80,04.7184e−143.4972e−155.5511e−153.9191e−140.1,0.13.0197e−141.8513e−179.9751e−135.3904e−140.2,0.29.1656e−141.8249e−153.8079e−121.1759e−120.3,0.31.9482e−132.8203e−158.9047e−123.4628e−120.4,0.44.2674e−137.0777e−151.8213e−116.6157e−120.5,0.59.4080e−131.0436e−143.9364e−111.0878e−110.6,0.61.9982e−123.2196e−148.7362e−111.6063e−110.7,0.73.9228e−124.3965e−141.7459e−101.9690e−110.8,0.86.4684e−128.2379e−142.8359e−101.6850e−100.9,0.97.9272e−126.1506e−143.2058e−105.0651e−121,15.7618e−128.9421e−145.1744e−113.9714e−12Figure 3
Pictorial results in Example3 for J=4,K=5: (a) (α=0.1,β=0.1), (b) (α=0.4,β=0.4), (c) (α=0.6,β=0.6), and (d) (α=0.8,β=0.8).Table 6
Norm of errors forα=β=0.5 and CPU time in Example 3.
J=5,K=4J=5,K=5J=5,K=6J=5,K=7J=5,K=8e∞e∞e∞e∞e∞1.7474e−134.6830e−105.1812e−102.2715e−093.7225e−08CPU0.6941s1.4075s1.4990s1.7516s2.1537sThe exact solution isVξ,η=η3eξ.The absolute errors forα=β=0.5 are presented in Table 7. Table 8 illustrates the norm of errors and CPU times for α=β=0.5 with various J and K values. Table 9 also contains data on L2 errors. Numerical solutions and pointwise error graphs are demonstrated in Figure 4. This figure shows the behavior of the numerical solution and the error function. Numerical results are in good settlement with theoretical results.Example 5.
Finally, we investigate the following equation on0,1×0,1(92)0CDηαVξ,η+Vξ,ηVξξ,η=∫0ηη−sβ−1Vξξξ,sds+Gξ,η,with conditions (2) and (3). The source term is taken as
(93)Gξ,η=21−ξcosξη2−αΓ3−α+η41−ξcosξ−sinξ−cosξ+ξsinξ−2sinξ+ξ−1cosξ2Γβηβ+2Γ3+β,and Vξ,η=η21−ξcosξ.Table 7
Numerical reports in Example4.
ξi,ηiJ=4J=4J=4J=4K=6K=7K=8K=90,03.9549e−201.9230e−226.7320e−201.1194e−180.1,0.11.0519e−105.3038e−124.6608e−132.8112e−120.2,0.22.7128e−091.4373e−107.5190e−128.1029e−120.3,0.32.6072e−081.4112e−096.8470e−111.6060e−110.4,0.41.0383e−075.6159e−092.7027e−103.5035e−110.5,0.52.8045e−071.5208e−087.3186e−108.5052e−110.6,0.66.1575e−073.3347e−081.6119e−091.9935e−100.7,0.71.2047e−066.5140e−083.1553e−094.1556e−100.8,0.82.1149e−061.1777e−075.7928e−097.3660e−100.9,0.92.8168e−061.7066e−078.9641e−099.7501e−101,11.3234e−142.6645e−151.0658e−143.5527e−15Table 8
Norm of errors forα=β=0.5 and CPU time in Example 4.
J=4,K=5J=4,K=6J=4,K=7J=4,K=8J=4,K=9e∞e∞e∞e∞|e∞5.4599e−053.6853e−062.2397e−071.1841e−081.2327e−09CPU0.6648s1.3185s1.2689s1.3789s1.5100sTable 9
Norm of errors forα=0.5,β=0.5 in Example 4.
JKe2453.0728e−05462.0621e−06471.1742e−07486.0069e−09497.7425e−10Figure 4
Pictorial results in Example4 for α=0.5,β=0.5 and K=9,J=4.In Table10, the L2 error is computed for α=β=0.5 and different N and M. Table 11 also shows the norm of errors and CPU times. Figure 5 shows the numerical solution and absolute error plots. This figure shows that for K=9 and J=4, the numerical solution is close to the exact solution. Table 10, Table 11, and Figure 5 affirm the validity and efficacy of the presented method.Table 10
Norm of errors forα=0.5,β=0.5 in Example 5.
α=β=0.3α=β=0.5α=β=0.7α=β=0.9JKe2e2e2e2357.2576e−057.1011e−057.1448e−056.0514e−05361.5363e−051.5000e−051.4978e−051.3806e−05373.8656e−073.7824e−073.7609e−073.4637e−07385.1965e−085.0757e−085.0432e−083.9225e−08391.1091e−092.4799e−091.6532e−086.7070e−08Table 11
Norm of errors forα=β=0.5 and CPU time in Example 5.
J=3,K=5J=3,K=6J=3,K=7J=3,K=8J=3,K=9e∞e∞e∞e∞|e∞1.1228e−042.3945e−056.0186e−078.0155e−084.0315e−09CPU0.5944s1.0808s1.1061s1.1594s1.2277sFigure 5
Pictorial results in Example5 for α=0.5,β=0.5 and K=9,J=3.


## 6. Conclusions

The purpose of this study was to suggest a collocation approach for solving a nonlinear TFPIDE based on Pell polynomials. In the Caputo sense, the fractional derivative is considered. The equation’s solution was expressed as a series of Pell polynomials with two variables. An algebraic system of nonlinear equations is obtained using the numerical technique. We proved that the method is convergent. Five test problems are provided to show that the method is efficacious. In numerical results, a small number of basis Pell polynomials is used to obtain good accuracy. In all examples, the CPU time was about one second. All of the tables and graphs demonstrated that the strategy is effective.



---

*Source: 8063888-2022-05-23.xml*