# Double Laplace Decomposition Method and Finite Difference Method of Time-fractional Schrödinger Pseudoparabolic Partial Differential Equation with Caputo Derivative

Authors: Mahmut Modanli; Bushra Bajjah
Journal: Journal of Mathematics
 (2021)

---

## Abstract

In this paper, an initial-boundary value problem for a one-dimensional linear time-dependent fractional Schrödinger pseudoparabolic partial differential equation with Caputo derivative of orderα∈0,1 is being considered. Two strong numerical methods are employed to acquire the solutions to the problem. The first method used is the double Laplace decomposition method where closed-form solutions are obtained for any α∈0,1. As the second method, the implicit finite difference scheme is applied to obtain the approximate solutions. To clarify the performance of these two methods, numerical results are presented. The stability of the problem is also investigated.


---

## Body

## 1. Introduction

The initial-boundary value problem for one-dimensional pseudoparabolic partial differential equation (PPPDE) is considered as follows:(1)utx,t=λuxxtx,t+kuxxx,t,with the same initial-boundary conditions mentioned for one-dimensional parabolic partial differential equation [1]:(2)utx,t=kuxxx,t.Equation (1) has been established considering the hydrostatic overflow pressure within a portion of mud during unification, where the parameter λ is a composite soil property with the dimensions of viscosity [2]. With regard to the assumption, if the resistance to compression is elastic (proportional to the rate of compression), then equation (1) will result with λ>0. On the contrary, equation (2) is led by the classic Terzaghi hypothesis, which presumes that any increase in hydrostatic overflow pressure is proportional to an increase in the ratio of pore volume to solids’ volume in mud [3].PPPDEs appear in several fields of mechanics and physics. They are used to investigate homogeneous fluid flow in fractured rocks, quasi-stationary paths in semiconductors, thermodynamic and transportation phenomena, and different physical systems (see, e.g., [4–6]).In the next step, the initial-boundary value problem for one-dimensional linear time-dependent Schrödinger pseudoparabolic partial differential equation (SPPPDE) in the domain0,L×0,T is considered as follows:(3)iutx,t=λuxxtx,t+uxxx,t+fx,t,x,t∈0,L×0,T,ux,0=ϕx,x∈0,L,u0,t=uL,t=0,t∈0,T,where i=−1 is an imaginary number, λ>0, and f and ϕ are known sufficiently smooth functions that may be complex-valued.Since Schrödinger found his famous equation in 1926, researchers have examined the counterintuitive laws of nature for microscopic systems [7]. The time-dependent Schrödinger partial differential equation (SPDE) is originated from the vector wave equation for the electric field which rules the propagation of electromagnetic waves in an inhomogeneous medium [8]. It arises in many areas such as in quantum physics, chemistry, electromagnetic wave diffusion, underwater acoustics, and design of certain optoelectronic devices [9, 10].Fractional analysis, which contains derivatives and integrals of the arbitrary order, is the generalized form of classical analysis. Recently, its applications in modeling real world phenomena have developed greatly in different fields of science and engineering, inclosing biology, physics, chemistry, medicine, finance, control theory, nanotechnology, viscoelasticity, anomalous transport, etc. The main reason for this is that the fractional modeling can be more accurate and effective, than the classical versions.In recent years, great effort has been made to obtain analytic and approximate solutions to fractional partial differential equations (FPDEs), and numerous effective methods have been proposed, for instance, the double Laplace transform method [11], Laplace Adomian decomposition method [12], fractional iteration algorithm method [13], fractional variational iteration method [14], finite difference method [15], finite element method [16], fractional subequation method [17], modified G′/G-expansion method [18], expψζ-expansion method [19], and q-homotopy analysis method [20].In this work, the initial-boundary value problem for one-dimensional linear time-dependent fractional Schrödinger pseudoparabolic partial differential equation (FSPPPDE) with Caputo derivative in the domain0,L×0,T is presented as follows:(4)iD0ctαux,t=λuxxtx,t+uxxx,t+fx,t,0<α≤1,λ>0,x,t∈0,L×0,T,ux,0=ϕx,x∈0,L,u0,t=uL,t=0,t∈0,T,where i=−1 is an imaginary number, f and ϕ are given sufficiently smooth and may possibly be complex-valued functions, and the Caputo derivative of order α is given as follows:(5)D0ctαux,t≔1Γm−α∫0tt−τm−α−1∂αux,τ∂ταdτ,m−1<α≤m,∂mux,t∂tm,α=m,where m=α, u∈Cm, and Γ denotes Gamma function.The double Laplace decomposition method will be constructed to solve problem (4) semianalytically. It considers one of the best integral transform utilized by many researchers. Kaabar et al. [21] defined a novel generalized double Laplace transform associated with the Adomian decomposition method and applied it to solve a reformulated nonlinear SPDE with spatiotemporal dispersion. A modified double Laplace transform decomposition method has been employed to solve the nonlinear coupled Hirota equations by Khan et al. [22]. Karatas Akgül et al. [23] utilized Laplace transform to solve the economic models based on the market balance in the sense of constant proportional Caputo derivative. Jarad and Abdeljawad [24] generalized the Laplace transform to make it viable to generalized fractional integrals and derivatives.This work is arranged as follows. In Section2, some necessary definitions and theorems on the double Laplace transform with an application are provided. In Section 3, the implicit finite difference scheme (FDS) for problem (4) is evoked then the stability is proved. In Section 4, we show an application with figures to illustrate the sufficiency as well as the accuracy of the approximate results. In the end, Section 5 is devoted to our conclusions.


## 2. Double Laplace Decomposition Method

In this section, we give some basic definitions and theorems of the double Laplace transform.Definition 1.
The double Laplace transform of the functionux,t is defined by the following double improper integral:(6)LxLtux,t≔8<Up,s≔∫0∞∫0∞e−px+stdtdx≔limm,n⟶∞∫0m∫0me−px+stdtdx,where p and s are complex numbers and ℜep>0 and ℜes>0.Definition 2.
If a piecewise continuous functionux,t on the intervals 0,X and 0,T is of the exponential order, that is,(7)ux,t≤Mee1x+e2t,e1,e2>0,for a positive constant M, then the double Laplace transform of ux,t exists for all p and s provided ℜep>e1 and ℜes>e2.Theorem 1.
Letu∈Cmℝ+×ℝ+ and m=maxn,k; there exist M,e1,e2>0 such that(8)∂i+jux,t∂xi∂tj≤Mee1x+e2t,i=0,1,…,n,j=0,1,…,k.
Then, the following formulas hold:(9)LxLt∂nux,t∂xn=pnUp,s−∑i=0n−1pn−1−iLt∂iu0,t∂xi,LxLt∂kux,t∂tk=skUp,s−∑j=0k−1sk−1−jLx∂jux,0∂tj,LxLt∂n+kux,t∂xn∂tk=pnskUp,s−∑i=0n−1p−1−iLt∂iu0,t∂xi−∑j=0k−1s−1−jLx∂jux,0∂tj+∑i=0n−1∑j=0k−1p−1−is−1−j∂i+ju0,0∂xi∂tj.Theorem 2.
If the Laplace transforms ofu and ∂ju/∂tj,j=0,1,…,k, exist, then the double Laplace transform of the partial fractional Caputo derivative with respect to the time variable is defined as(10)LxLtD0ctαux,t=sαUp,s−∑j=0k−1s−j−1Lx∂jux,0∂tj.


### 2.1. Application of the Method

We consider the following one-dimensional linear time-dependent FSPPPDE in terms of Caputo derivative:(11)iD0ctαux,t=λuxxtx,t+uxxx,t+fx,t,0<α≤1,λ>0,x>0,t>0,where(12)fx,t=1+tα+λαtα−1+iΓα+1sinx,subject to the initial condition,(13)ux,0=sinx,and boundary conditions,(14)u0,t=uπ,t=0.Applying the double Laplace transform with respect tox and t to equation (11) yields(15)LxLtux,t=1s+1isα+1+Γα+1sα+1+Γα+1is2α+1+λΓα+1is2α1p2+1+1isαLxLtλutxxx,t+uxxx,t.Taking the inverse double Laplace transform to both sides of equation (15), it results in(16)ux,t=1+1iΓα+1tα+tα+Γα+1iΓ2α+1t2α+λΓα+1iΓ2αt2α−1sinx+Lx−1Lt−11isαLxLtλutxxx,t+uxxx,t.Using the decomposition series forux,t into equation (16) leads to(17)∑n=0∞unx,t=1+1iΓα+1tα+tα+Γα+1iΓ2α+1t2α+λΓα+1iΓ2αt2α−1sinx+Lx−1Lt−11isαLxLtλ∂3∂t∂x2∑n=0∞unx,t+∂2∂x2∑n=0∞unx,t.From equation (17), the following recursive relationship can be obtained:(18)u0x,t=1−iΓα+1tα+tα−iΓα+1Γ2α+1t2α−iλΓα+1Γ2αt2α−1sinx,(19)un+1x,t=Lx−1Lt−11isαLxLtλ∂3∂t∂x2unx,t+∂2∂x2unx,t,n≥0.The following formula is reached from equation (19):(20)u1x,t=λ1+iΓα+1Γ2αt2α−1+Γα+1Γ3αt3α−1+iλΓα+1Γ3α−1t3α−2+iΓα+1tα+1+iΓα+1Γ2α+1t2α+iΓα+1Γ3α+1t3α+λΓα+1Γ3αt3α−1sinx,u2x,t=λ2i−Γα+1Γ3α−1t3α−2+iΓα+1Γ4α−1t4α−2−λΓα+1Γ4α−2t4α−3+λ−1Γ2αt2α−1+2i−Γα+1Γ3αt3α−1+i−1Γα+1Γ4αt4α−1+λi−1Γα+1Γ4α−1t4α−2+−1Γ2α+1t2α+i−Γα+1Γ3α+1t3α−Γα+1Γ4α+1t4α+iλΓα+1Γ4αt4α−1sinx,and so on. Hence, the required analytical solution is(21)ux,t=1+tαsinx,where the other terms vanish in the limit.


## 2.1. Application of the Method

We consider the following one-dimensional linear time-dependent FSPPPDE in terms of Caputo derivative:(11)iD0ctαux,t=λuxxtx,t+uxxx,t+fx,t,0<α≤1,λ>0,x>0,t>0,where(12)fx,t=1+tα+λαtα−1+iΓα+1sinx,subject to the initial condition,(13)ux,0=sinx,and boundary conditions,(14)u0,t=uπ,t=0.Applying the double Laplace transform with respect tox and t to equation (11) yields(15)LxLtux,t=1s+1isα+1+Γα+1sα+1+Γα+1is2α+1+λΓα+1is2α1p2+1+1isαLxLtλutxxx,t+uxxx,t.Taking the inverse double Laplace transform to both sides of equation (15), it results in(16)ux,t=1+1iΓα+1tα+tα+Γα+1iΓ2α+1t2α+λΓα+1iΓ2αt2α−1sinx+Lx−1Lt−11isαLxLtλutxxx,t+uxxx,t.Using the decomposition series forux,t into equation (16) leads to(17)∑n=0∞unx,t=1+1iΓα+1tα+tα+Γα+1iΓ2α+1t2α+λΓα+1iΓ2αt2α−1sinx+Lx−1Lt−11isαLxLtλ∂3∂t∂x2∑n=0∞unx,t+∂2∂x2∑n=0∞unx,t.From equation (17), the following recursive relationship can be obtained:(18)u0x,t=1−iΓα+1tα+tα−iΓα+1Γ2α+1t2α−iλΓα+1Γ2αt2α−1sinx,(19)un+1x,t=Lx−1Lt−11isαLxLtλ∂3∂t∂x2unx,t+∂2∂x2unx,t,n≥0.The following formula is reached from equation (19):(20)u1x,t=λ1+iΓα+1Γ2αt2α−1+Γα+1Γ3αt3α−1+iλΓα+1Γ3α−1t3α−2+iΓα+1tα+1+iΓα+1Γ2α+1t2α+iΓα+1Γ3α+1t3α+λΓα+1Γ3αt3α−1sinx,u2x,t=λ2i−Γα+1Γ3α−1t3α−2+iΓα+1Γ4α−1t4α−2−λΓα+1Γ4α−2t4α−3+λ−1Γ2αt2α−1+2i−Γα+1Γ3αt3α−1+i−1Γα+1Γ4αt4α−1+λi−1Γα+1Γ4α−1t4α−2+−1Γ2α+1t2α+i−Γα+1Γ3α+1t3α−Γα+1Γ4α+1t4α+iλΓα+1Γ4αt4α−1sinx,and so on. Hence, the required analytical solution is(21)ux,t=1+tαsinx,where the other terms vanish in the limit.


## 3. Finite Difference Scheme

The usual notationsxn=nh,0≤n≤M, tk=kδt,0≤k≤N, and uxn,tk≡unk are used, where unk denotes the numerical solution of the analytical solution of ux,t at the grid points xn,tk.The numerical approximation based on the Caputo derivative is defined as(22)D0ctαunk≈σα,δt∑j=0kdjαunk−j+1−unk−j+Oδt2−α,where σα,δt=δt−α/Γ2−α and djα=j+11−α−j1−α.Since the solutionu∈Cmandm∈ℕ, then the following three-layer implicit FDS can be written at any point xn,tk as follows:(23)uxxtnk=utxxnk≈1h2utn+1k−2utnk+utn−1k+Oh2≈1δtun+1k+1−2unk+1+un−1k+1h2−un+1k−2unk+un−1kh2+Oδth2+h2.By neglecting the following discretization errorsOδt2−α, Oδt/h2+h2, and Oh2, the implicit FDS construction of problem (4) can be written as(24)iσα,δt∑j=0kdjαunk+1−j−unk−j≈λh2δtun+1k+1−2unk+1+un−1k+1−un+1k−2unk+un−1k+1h2un+1k+1−2unk+1+un−1k+1+fnk,0<n<M,0<k≤N,0<α≤1,λ>0,fnk=fxn,tk,un0=ϕxn,0≤n≤M,u0k=uMk=0,0≤k≤N,and arrange system (24) as(25)iσα,δtunk+1−unk+iσα,δt∑j=1kdjαunk+1−j−unk−j≈λh2δtun+1k+1−2unk+1+un−1k+1−un+1k−2unk+un−1k+1h2un+1k+1−2unk+1+un−1k+1+fnk,0<n<M,0<k≤N,0<α≤1,λ>0,fnk=fxn,tk,un0=ϕxn,0≤n≤M,u0k=uMk=0,0≤k≤N.System (25) can be rewritten as follows:(26)λh2δtun+1k−1h2+λh2δtun+1k+1−2λh2δt+iσα,δtunk+2h2+2λh2δt+iσα,δtunk+1+λh2δtun−1k−1h2+λh2δtun−1k+1+iσα,δt∑j=1kdjαunk+1−j−unk−j≈fnk,0<n<M,0<k≤N,0<α≤1,λ>0,fnk=fxn,tk,un0=ϕxn,0≤n≤M,u0k=uMk=0,0≤k≤N.Writing system (26) in a matrix form is demonstrated as(27)Qun+1+Gun+Qun−1=ηn,where un=un0,un1,…,unNT, ηn=ηn0,ηn1,…,ηnNT, and ηnk=fnk=fxn,tk, 0<n<M,0<k<N. Also, Q and G are two matrices of the size N+1×N+1 and have the following form:(28)Q=000⋯⋯⋯00ab0⋯⋯000ab0⋯0⋮⋮⋱⋱⋮⋮⋮⋱⋱⋮00⋯⋯0ab000⋯⋯⋯0,G=100⋯⋯⋯0−ξ1c+γ1d0⋯⋯0−ξ2γ2c+γ1d⋯⋯0⋮⋮⋱⋱⋮⋮⋮⋱⋱⋮−ξk−1γk−1⋯⋯γ2c+γ1d−ξkγk⋯⋯⋯γ2c+γ1,where a=λ/h2δt, b=−1/h2+λ/h2δt, c=−2λ/h2δt, d=2/h2+2λ/h2δt+iσα,δt, ξn=iσα,δtdnα, and γn=ξn−ξn−1, 1≤n≤k.


### 3.1. Stability and Convergence of the FDM

To mutate the two-step FDS (27) into one-step, the modified Gauss elimination method is used to achieve(29)un=θn+1un+1+βn+1,n=1,2,…,M−1.By using the given Dirichlet boundary conditions in equation (29), we obtain(30)θ1u1+β1=u0.Thus, it gives(31)θ1=0N+1×N+1,β1=0N+1×1.To determine the matricesθn+1 and βn+1, we substitute equation (29) into equation (27), and we obtain(32)Q+Gθn+1+Qθnθn+1un+1+Gβn+1+Qβn+Qθnβn+1=ηn.Hence, equation (32) gives(33)Q+Gθn+1+Qβnθn+1=0,Gβn+1+Qβn+Qαnβn+1=ηn.Finally, the following equalities are acquired:(34)θn+1=−G+Qθn−1Q,βn+1=G+Qαn−1ηn−Qβn,where n=1,2,…,M−1.Stability can be acquired by applying the method of analyzing the eigenvalues of the iteration matrices of the schemes. The maximum of the absolute value of the eigenvalues of the matrixQ has been represented by ρQ, which refers to the spectral radius of the given matrix.Theorem 3.
FDS (26) is conditionally stable subject to the following space step limitation:(35)h>Γ2−αd1α−d0α2δt+4λδt1−2α4.Proof.
Using the analytical methods used in [25], we have to prove that ρθn<1,n=1,2,…,M, under condition (35).(i)
Obviously,ρθ1=0<1.(ii)
(36)ρθ2=ρ−G−1Q≤−G−1Q≤−G−1Q=Qmin1≤k≤N−1akk−∑m≠km=1N−1akm=λ/h2δt+−1/h2+λ/h2δt−2λ/h2δt+iδt−α/Γ2−αd1α−d0α=2λ/h2δt+1/h22λ/h2δt2+δt−α/Γ2−αd1α−d0α2<1⇔h>Γ2−αd1α−d0α2δt+4λδt1−2α4.(iii)
Assume thatρθn<1 subject to condition (35) and we need to find the value of ρθn+1<1. We notice that(37)θn+1=00ξ∗γ∗+aθn2,20∗ξ∗γ∗+aθn3,3⋮⋮⋱0∗⋯⋯ξ∗γ∗+aθnN+1,N+1,where ξ∗=2λ/h2δt+1/h2, γ∗=2λ/h2δt2+δt−α/Γ2−αd1α−d0α2, and a=λ/h2δt.
As we know thatθni,i=ρθn for 2≤i≤N+1 and from the induction hypothesis, we conclude that ρθn+1<1 subject to condition (35). Thus, we acquire the required result by induction.


## 3.1. Stability and Convergence of the FDM

To mutate the two-step FDS (27) into one-step, the modified Gauss elimination method is used to achieve(29)un=θn+1un+1+βn+1,n=1,2,…,M−1.By using the given Dirichlet boundary conditions in equation (29), we obtain(30)θ1u1+β1=u0.Thus, it gives(31)θ1=0N+1×N+1,β1=0N+1×1.To determine the matricesθn+1 and βn+1, we substitute equation (29) into equation (27), and we obtain(32)Q+Gθn+1+Qθnθn+1un+1+Gβn+1+Qβn+Qθnβn+1=ηn.Hence, equation (32) gives(33)Q+Gθn+1+Qβnθn+1=0,Gβn+1+Qβn+Qαnβn+1=ηn.Finally, the following equalities are acquired:(34)θn+1=−G+Qθn−1Q,βn+1=G+Qαn−1ηn−Qβn,where n=1,2,…,M−1.Stability can be acquired by applying the method of analyzing the eigenvalues of the iteration matrices of the schemes. The maximum of the absolute value of the eigenvalues of the matrixQ has been represented by ρQ, which refers to the spectral radius of the given matrix.Theorem 3.
FDS (26) is conditionally stable subject to the following space step limitation:(35)h>Γ2−αd1α−d0α2δt+4λδt1−2α4.Proof.
Using the analytical methods used in [25], we have to prove that ρθn<1,n=1,2,…,M, under condition (35).(i)
Obviously,ρθ1=0<1.(ii)
(36)ρθ2=ρ−G−1Q≤−G−1Q≤−G−1Q=Qmin1≤k≤N−1akk−∑m≠km=1N−1akm=λ/h2δt+−1/h2+λ/h2δt−2λ/h2δt+iδt−α/Γ2−αd1α−d0α=2λ/h2δt+1/h22λ/h2δt2+δt−α/Γ2−αd1α−d0α2<1⇔h>Γ2−αd1α−d0α2δt+4λδt1−2α4.(iii)
Assume thatρθn<1 subject to condition (35) and we need to find the value of ρθn+1<1. We notice that(37)θn+1=00ξ∗γ∗+aθn2,20∗ξ∗γ∗+aθn3,3⋮⋮⋱0∗⋯⋯ξ∗γ∗+aθnN+1,N+1,where ξ∗=2λ/h2δt+1/h2, γ∗=2λ/h2δt2+δt−α/Γ2−αd1α−d0α2, and a=λ/h2δt.
As we know thatθni,i=ρθn for 2≤i≤N+1 and from the induction hypothesis, we conclude that ρθn+1<1 subject to condition (35). Thus, we acquire the required result by induction.


## 4. Numerical Applications

Considering the following initial-boundary value problem for one-dimensional linear time-dependent FSPPPDE in the sense of Caputo derivative,(38)iD0ctαux,t=λutxxx,t+uxxx,t+fx,t,fx,t=1+tα+λαtα−1+iΓα+1sinx,0<α≤1,λ>0,0<x<π,0<t≤1,u0,t=uπ,t=0,0≤t≤1,ux,0=sinx,0≤x≤π.We have shown that the analytical solution of problem (38) is(39)ux,t=1+tαsinx.The maximum norm of the error of the numerical solution can be reached with(40)ε=maxn=0,1,…,Mk=0,1,…,Nux,t−uxn,tk,where ux,t and uxn,tk indicate the analytical and numerical solution, respectively.Figures1 and 2 illustrate the comparison between the profile evolution of the analytical and numerical solutions of problem (38). The numerical solution was approximated by using implicit FDS scheme (26), whereas the analytical solution was approximated by using the double Laplace decomposition method.Figure 1
Profile evolution solution of problem (38) at N=210, M=25, α=0.25, and λ=2.480. (a) Numerical simulation. (b) Analytical solution.
(a)(b)Figure 2
Profile evolution solution of problem (38) at N=100, M=29, α=0.90, and λ=2.412. (a) Numerical simulation. (b) Analytical solution.
(a)(b)We observe from Figures1 and 2 that the numerical solution improves as the value of α approaches 1 (which is the case where the FPDE becomes a classical PDE). This also illustrates how conditional stability (35) of our scheme (26) is affected by the value of α.To quantitatively show the accuracy and efficiency of our scheme, we ran several simulations with different input parameters to compare the numerical with the analytical solutions. Table1 shows the error ε between the two solutions.Table 1
Comparison between the numerical and analytical solutions for problem (38) with h=π/M and δt=1/N, where the chosen values of the initial data satisfy stability condition (35).
NMαλNumerical solutionAnalytical solutionε37150.101.981.98901.98830.0008210250.252.4801.99791.99610.0018263160.332.6601.99722.00000.00281080.502.1502.00032.00000.000315150.802.171.98901.98500.0040100290.902.4121.99811.99710.0010


## 5. Conclusion

In this manuscript, the initial-boundary value problem for a one-dimensional linear time-dependent FSPPPDE in the sense of Caputo derivative of orderα∈0,1 is presented. Double Laplace decomposition and implicit FDS methods were investigated to obtain the closed-form and numerical solutions of problem (4), respectively. We have proved that implicit FDS (32) of problem (4) is conditionally stable. From the abovementioned table, it can be concluded that the obtained numerical solution has a good agreement with the analytic solution of various α∈0,1. The results obtained in this study may be useful in many branches of quantum physics.In the future work, we hope to extend this model to multidimensional nonlinear time-dependent coupled SPPPDEs with Caputo fractional derivatives and solve them by the finite element method.



---

*Source: 7113205-2021-10-23.xml*