# A Note on the Inverse Problem for a Fractional Parabolic Equation

Authors: Abdullah Said Erdogan; Hulya Uygun
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

For a fractional inverse problem with an unknown time-dependent source term, stability estimates are obtained by using operator theory approach. For the approximate solutions of the problem, the stable difference schemes which have first and second orders of accuracy are presented. The algorithm is tested in a one-dimensional fractional inverse problem.


---

## Body

## 1. Introduction

Inverse problems arise in many fields of science and engineering such as ion transport problems, chromatography, and heat determination problems with an unknown internal energy source. Different typed of inverse problems have been investigated, and the main results obtained in this field of research were given by many researchers (see [1–10]). More than three centuries the theory of fractional derivatives developed mainly as a pure theoretical field of mathematics. Fractional integrals and derivatives appear in the theory of control of dynamical systems, when the controlled system or/and the controller is described by a fractional differential equation (see [11]). Recently, many application areas such as bioengineering applications, image and signal processing are also related to fractional calculus. Methods of solutions of problems and theory of fractional calculus have been studied by many researchers [11–28]. Among them finite difference method is used for solving several fractional differential equations (see [20, 22, 23, 27] and the references therein).


### 1.1. Statement of the Problem

Many scientists and researchers are trying to enhance mathematical models of real-life cases for investigating and understanding the behavior of them. Therefore, some phenomena have been modeled and investigated as fractional inverse problems (see [29–33] and the references therein). In this paper, we consider the fractional parabolic inverse problem with the Dirichlet condition
(1.1)∂u(t,x)∂t-a∂2u(t,x)∂x2-Dt1/2u(t,x)+σu(t,x)=p(t)q(x)+f(t,x),0<x<π,0<t≤T,u(t,0)=u(t,π)=0,0≤t≤T,u(0,x)=φ(x),0≤x≤π,u(t,x*)=ρ(t),0<x*<π,
where u(t,x) and p(t) are unknown functions, a(x)≥a>0, and σ>0 is a sufficiently large number. Here, Dt1/2=D0+1/2 is the standard Riemann-Liouville’s derivative of order 1/2.Theorems on the stability of problem (1.1) are analyzed by assuming that q(x) is a sufficiently smooth function, q(0)=q(π)=0 and q(x*)≠0.


## 1.1. Statement of the Problem

Many scientists and researchers are trying to enhance mathematical models of real-life cases for investigating and understanding the behavior of them. Therefore, some phenomena have been modeled and investigated as fractional inverse problems (see [29–33] and the references therein). In this paper, we consider the fractional parabolic inverse problem with the Dirichlet condition
(1.1)∂u(t,x)∂t-a∂2u(t,x)∂x2-Dt1/2u(t,x)+σu(t,x)=p(t)q(x)+f(t,x),0<x<π,0<t≤T,u(t,0)=u(t,π)=0,0≤t≤T,u(0,x)=φ(x),0≤x≤π,u(t,x*)=ρ(t),0<x*<π,
where u(t,x) and p(t) are unknown functions, a(x)≥a>0, and σ>0 is a sufficiently large number. Here, Dt1/2=D0+1/2 is the standard Riemann-Liouville’s derivative of order 1/2.Theorems on the stability of problem (1.1) are analyzed by assuming that q(x) is a sufficiently smooth function, q(0)=q(π)=0 and q(x*)≠0.


## 2. Main Results

In this section, stability estimates for the solution of (1.1) are investigated. For the mathematical substantiation, we introduce the Banach space C∘α[0,π],α∈(0,1), of all continuous functions ϕ(x) defined on [0,π] with ϕ(0)=ϕ(π)=0 satisfying a Hölder condition for which the following norm is finite
(2.1)‖ϕ‖C∘α[0,π]=‖ϕ‖C[0,π]+sup⁡0<x<x+h<π|ϕ(x+h)-ϕ(x)|hα,
where C[0,π] is the space of all continuous function ϕ(x) defined on [0,π] with the norm
(2.2)‖ϕ‖C[0,π]=max⁡0≤x≤π|ϕ(x)|.
With the help of a positive operator A, we introduce the fractional spaces Eα,0<α<1, consisting of all v in a Banach space E for which the following norm is finite:
(2.3)‖v‖Eα=‖v‖E+sup⁡λ>0λ1-α‖Aexp⁡{-λA}v‖E.
Throughout the paper, positive constants will be indicated by Mi(α,β,…). Here variables are used to focus on the fact that the constant depends only on α,β,… and the subindex i is used to indicate a different constant.Theorem 2.1.
Letφ∈C∘2α+2[0,π],F∈C([0,T],C∘2α[0,π]), and ρ′∈C[0,T]. Then for the solution of problem (1.1), the following coercive stability estimates
(2.4)‖ut‖C([0,T],C∘2α[0,π])+‖u‖C([0,T],C∘2α+2[0,π])≤M(x*,q)‖ρ′‖C[0,T]+M(a,δ,σ,α,x*,q,T)×(‖φ‖C∘2α+2[0,π]+‖F‖C([0,T],C∘2α[0,π])+‖ρ‖C[0,T]),‖p‖C[0,T]≤M(x*,q)‖ρ′‖C[0,T]+M(a,δ,σ,α,x*,q,T)×[‖φ‖C∘2α+2[0,π]+‖F‖C([0,T],C∘2α[0,π])+‖ρ‖C[0,T]]
hold.Proof.
Let us search for the solution of inverse problem (1.1) in the following form (see [8]):
(2.5)u(t,x)=η(t)q(x)+w(t,x),
where
(2.6)η(t)=∫0t‍p(s)ds.
Using the overdetermined condition, we get
(2.7)η(t)=ρ(t)-w(t,x*)q(x*),(2.8)p(t)=ρ′(t)-wt(t,x*)q(x*).
Using identity (2.8) and the triangle inequality, it follows that
(2.9)|p(t)|=|ρ′(t)-wt(t,x*)q(x*)|≤M(x*,q)(|ρ′(t)|+|wt(t,x*)|)≤M(x*,q)(max⁡0≤t≤T|ρ′(t)|+max⁡0≤t≤Tmax⁡0≤x≤π|wt(t,x)|)≤M(x*,q)(max⁡0≤t≤T|ρ′(t)|+max⁡0≤t≤T‖wt(t,⋅)‖C∘2α[0,π])
for any t,t∈[0,T]. Here, w(t,x) is the solution of the following problem:
(2.10)∂w(t,x)∂t-a∂2w(t,x)∂x2-aρ(t)-w(t,x*)q(x*)d2q(x)dx2-Dt1/2w(t,x)-Dt1/2ρ(t)-Dt1/2w(t,x*)q(x*)q(x)+σρ(t)-w(t,x*)q(x*)q(x)+σw(t,x)=f(t,x),0<x<π,0<t≤T,w(t,0)=w(t,π)=0,0≤t≤T,w(0,x)=φ(x),0≤x≤π.
For simplicity, we assign
(2.11)F(t,x)=aρ(t)q(x*)d2q(x)dx2-σρ(t)q(x*)q(x)+Dt1/2ρ(t)q(x*)q(x)+f(t,x),G(t,x)=Q1(q,ρ,x,x*,t)w(t,x*)+Q2(q,x,x*)Dt1/2w(t,x*)+Dt1/2w(t,x),
where
(2.12)Q1(q,ρ,x,x*,t)=1q(x*)(-ad2q(x)dx2+σρ(t)),Q2(q,x,x*)=-q(x)q(x*).
Note that functions F(t,x),Q1(q,ρ,x,x*,t) and Q2(q,x,x*) only contain given functions. Then, we can rewrite problem (2.10) as
(2.13)∂w(t,x)∂t-a∂2w(t,x)∂x2+σw(t,x)=F(t,x)+G(t,x),0<x<π,0<t≤T,w(t,0)=w(t,π)=0,0≤t≤T,w(0,x)=φ(x),0≤x≤π.
So, the end of proof of Theorem 2.1 is based on estimate (2.9) and the following theorem.Theorem 2.2.
For the solution of problem (2.10), the following coercive stability estimate
(2.14)‖wt‖C∘2α[0,π]≤M(a,δ,σ,α,x*,q,T)×(‖φ‖C∘2α+2[0,π]+‖F‖C([0,T],C∘2α[0,π])+‖ρ‖C[0,T])
holds.Proof.
In a Banach spaceE=C∘[0,π], with the help of the positive operator A defined by
(2.15)Au=-a(x)∂2u(t,x)∂x2+σu,
with
(2.16)D(A)={u(x):u,u′,u′′∈C[0,π],u(0)=u(π)=0},
where σ is a positive constant, problem (2.10) can be written in the abstract form as an initial-value problem
(2.17)wt+Aw=F(t)+G(t),0<t≤T,w(0)=φ.
By the Cauchy formula, the solution can be written as
(2.18)w(t)=e-tAφ-∫0te-(t-s)A(F(s)+G(s))ds.
Applying the formula
(2.19)Dt1/2u(t)=∫0tu′(ζ)dξπ(t-ξ)1/2,
we get the following presentation of the solution of abstract problem (2.17):
(2.20)D1/2w(t)=-∫0tAe-ξAφπ(t-ξ)1/2dξ-∫0tF(ξ)π(t-ξ)1/2dξ-∫0tG(ξ)π(t-ξ)1/2dξ+∫0t∫0ξAe-(ξ-s)AF(s)π(t-ξ)1/2dsdξ+∫0t∫0ξAe-(ξ-s)AG(s)π(t-ξ)1/2dsdξ.
Changing the order of integration, we obtain that
(2.21)D1/2w(t)=-∫0tAe-ξAφπ(t-ξ)1/2dξ-∫0tF(ξ)π(t-ξ)1/2dξ+∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξF(s)ds-∫0tG(ξ)π(t-ξ)1/2dξ+∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξG(s)ds=∑k=15‍Jk,
where
(2.22)J1(t)=-∫0tAe-ξAφπ(t-p)1/2dξ,J2(t)=-∫0tF(ξ)π(t-ξ)1/2dξ,J3(t)=∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξF(s)ds,J4(t)=-∫0tG(ξ)π(t-ξ)1/2dξ,J5(t)=∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξG(s)ds.
Now, we estimate Jk(t),k=1,2,3,4,5 separately. It is known that [13]
(2.23)‖Aαe-tA‖E→E≤M,0≤α≤1.
Since operators A and exp(-tA) commute,
(2.24)‖Ae-tAφ‖Eα≤‖e-tA‖Eα→Eα‖Aφ‖Eα≤‖e-tA‖E→E‖Aφ‖Eα.
Applying the definition of norm of the spaces Eα and (2.23) and (2.24), we get
(2.25)‖J1(t)‖Eα=‖∫0tAe-ξAφπ(t-p)1/2dξ‖Eα≤M1‖Aφ‖Eα
for any t,t∈[0,T]. Estimation of J2(t) is as follows:
(2.26)‖J2(t)‖Eα=‖∫0tF(ξ)π(t-ξ)1/2dξ‖Eα≤‖F(t)‖C(Eα)∫0t1π(t-ξ)1/2dξ≤M2‖F‖C(Eα).
Let us estimate J3(t):
(2.27)‖J3(t)‖Eα=‖∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξF(s)ds‖Eα≤‖∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξds‖Eα→Eα‖F‖C(Eα).
It is proven that (see [28])
(2.28)‖∫stAe-(ξ-s)Aπ(t-ξ)1/2dξ‖E→E≤Mt-s.
Using the definition of norm of the spaces Eα, we can obtain that
(2.29)‖∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξds‖Eα→Eα=∫0t‖∫stAe-(t-s)Aπ(t-ξ)1/2dξ‖E→Eds+sup⁡λ>0∫0t‖∫stλ1-αAe-λAAe-(t-s)Aπ(t-ξ)1/2dξ‖E→Eds.
Using estimates (2.23) and (2.28), we get
(2.30)‖J3(t)‖Eα≤‖∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξds‖Eα→Eα‖F‖C(Eα)≤M3‖F‖C(Eα).
Expanding G(s), estimation of J4(t) is as follows:
(2.31)‖J4(t)‖Eα≤∫0t‖Q1(q,ρ,x,x*,t)w(ξ,x*)π(t-ξ)1/2‖Eαdξ+∫0t‖Q2(q,x,x*)Dt1/2w(ξ,x*)π(t-ξ)1/2‖Eαdξ+∫0t‖Dt1/2w(ξ,x)π(t-ξ)1/2‖Eαdξ.
It is known that (see [34])
(2.32)‖w‖Eα≤M‖Dt1/2w‖Eα.
Since Q1(q,ρ,x,x*,t) and Q2(q,x,x*) are known functions, it is easy to see that
(2.33)‖J4(t)‖Eα≤M4(q,ρ,x,x*,T)∫0t1π(t-ξ)1/2‖Dt1/2w(ξ)‖Eαdξ.
Estimation of J5(t) can be given similar to the estimation of J4(t). By (2.23) and (2.32),
(2.34)‖J5(t)‖Eα≤∫0t∫stAe-(ξ-s)Aπ(t-ξ)1/2dξG(s)ds≤M5(q,ρ,x,x*,T)∫0t‖Dt1/2w(s)‖Eαds.
Finally combining estimates (2.25), (2.26), (2.30), (2.33), and (2.34), we get
(2.35)‖Dt1/2w‖Eα≤M1‖Aφ‖Eα+(M2+M3)‖F‖C(Eα)+∫0t(M4π(t-ξ)1/2+M5)‖Dt1/2w(s)‖Eαds.
Using the Gronwall’s inequality, we can write
(2.36)‖Dt1/2w‖Eα≤eM6(M1‖Aφ‖Eα+M7‖F‖C(Eα)).
From the last estimate, we can obtain the estimate for wt(t) by using problem (2.17) and well-posedness of the Cauchy problem in C(Eα) (see [35]). So the following theorem finishes the proof of Theorem 2.2.Theorem 2.3 (see, [36]).
For0<α<1/2, the norms of the spaces Eα(C[0,π],A) and C2α[0,π] are equivalent.


## 3. Numerical Results

We have not been able to obtain a sharp estimate for the constants figuring in the stability inequalities. So we will provide the following results of numerical experiments of the following problem:(3.1)∂u(t,x)∂t=∂2u(t,x)∂x2-u(t,x)+Dt1/2u(t,x)+p(t)sin⁡x+f(t,x),f(t,x)=(-3t-1πt-1/2+2πt1/2)sin⁡x,x∈(0,π),t∈(0,1],u(0,x)=sin⁡x,x∈[0,π],u(t,0)=u(t,π)=0,t∈[0,1],u(t,π2)=1-t.
The exact solution of the given problem is u(t,x)=(1-t)sinx and for the control parameter p(t) is 1+t.


### 3.1. The First Order of Accuracy Difference Scheme

For the approximate solution of the problem (3.1), the Rothe difference scheme
(3.2)unk-unk-1τ=un+1k-2unk+un-1kh2-unk+Dτ1/2unk+pkqn+f(tk,xn),f(tk,xn)=(-3tk-1πtk-1/2+2πtk1/2)sin⁡xn,pk=p(tk),qn=sin⁡xn,xn=nh,tk=kτ,1≤k≤N,1≤n≤M-1,Mh=π,Nτ=1,un0=sin⁡xn,0≤n≤M,u0k=uMk=0,0≤k≤N,usk=ρ(tk),ρ(tk)=1-tk,0≤k≤N,s=⌊π2h⌋,
where ⌊x⌋ denotes greatest integer less than x is constructed. Throughout the paper, let us denote
(3.3)ρ(tk)=1-tk,qn=sin⁡xn,tk={tk=kτ,0≤k≤N,Nτ=1},xn={xn=nh,0≤n≤M-1,Mh=π},f(tk,xn)=(-3tk-1πtk-1/2+2πtk1/2)sin⁡xn,F(tk,xn)=ρ(tk)sin⁡⁡(xs)(sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)h2-sin⁡⁡(xn))-1π∑m=1k‍Γ(k-m+(1/2))(k-m)!τ1/2sin⁡⁡(xs)sin⁡(xn)+f(tk,xn).
We search the solution of (3.2) in the following form:
(3.4)unk=ηkqn+wnk,
where
(3.5)ηk=∑i=1k‍piτ,1≤k≤N,η0=0.
Moreover for the interior grid point usk, we have that
(3.6)usk=ηkqs+wsk=ρ(tk).
From (3.4), (3.5), and the condition usk=ρ(tk), it follows that
(3.7)ηk=ρ(tk)-wskqs,(3.8)pk=ηk-ηk-1τ,1≤k≤N,(3.9)unk=ρ(tk)-wskqsqn+wnk,0≤k≤N,0≤n≤M,
where wnk,0≤k≤N,0≤n≤M is the solution of the difference scheme
(3.10)wnk-wnk-1τ=wn+1k-2wnk+wn-1kh2-wnk-wsksin⁡⁡(xs)(sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)h2-sin⁡⁡(xn))-1π∑m=1k‍Γ(k-m+(1/2))(k-m)!(wsm-wsm-1sin⁡⁡(xs)τ1/2sin⁡⁡(xn)-wnm-wnm-1τ1/2)+F(tk,xn),1≤k≤N,1≤n≤M-1,w0k=wMk=0,0≤k≤N,wn0=sin⁡(xn),0≤n≤M.
First, applying the first order of accuracy difference scheme (3.10), we obtain (M+1)×(M+1) system of linear equations and we write them in the matrix form
(3.11)Awk+∑j=0k-1‍Bjwj=Dφk,1≤k≤N,w0={sin⁡(xn)}n=0M,
where
(3.12)A=[1000.0.000xyx0.z1.0000xyx.z2.000..........0000.zM-1.xyx0000.0.001](M+1)×(M+1),B0=[000.0.000a0.f1.0000a.f2.00........000.fs+a.00........000.fM-1.0a000.0.00](M+1)×(M+1),Bj=[000.0.000c0.d1.0000c.d2.00........000.ds+c.00........000.dM-1.0c000.0.00](M+1)×(M+1),
for any j=1,2,…,k-2, and
(3.13)Bk-1=[000.0.000v0.c1.0000v.c2.00........000.cs+v.00........000.cM-1.0v000.0.00](M+1)×(M+1).
Here, for any n=1,2,…,M-1,
(3.14)x=-1h2,y=1τ+2h2+1+1τ,zn=sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)sin⁡(xs)h2-sin⁡(xn)sin⁡(xs)-sin⁡(xn)sin⁡(xs)τ,in(s+1)thcolumn,a=Γ(k-1/2)τπ(k-1)!,fn=-sin⁡(xn)Γ(k-1/2)τπsin⁡(xs)(k-1)!,c=1τπ(Γ(k-j-1/2)(k-j-1)!-Γ(k-j+1/2)(k-j)!),dn=sin⁡(xn)τπsin⁡(xs)(Γ(k-j-1/2)(k-j-1)!-Γ(k-j+1/2)(k-j)!),v=-1τ-1τ,cn=sin⁡(xn)τsin⁡(xs),wr=[w0r⋮wMr](M+1)×1foranyr=0,1,…,k,φk=[0ϕ1k⋮ϕM-1k0](M+1)×1,ϕnk=(ρ(tk)sin⁡(xs)sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)h2-ρ(tk)sin⁡(xs)sin⁡(xn))-1π∑m=1k‍Γ(k-m+1/2)(k-m)!τ1/2sin⁡(xs)sin⁡(xn)+f(tk,xn),
and D is (M+1)×(M+1) identity matrix. Using (3.11), we can obtain that
(3.15)wk=A-1(Dφk-∑j=0k-1‍Bjwj),k=1,2,…,N,w0={sin⁡xn}n=0M.
To solve the resulting difference equations, we apply the method given in (3.15) step by step for k=1,2,…,N. For the evaluation of wr,r=2,3,…,N,wr-1 is needed. It is obtained in the previous step. Then, the solution pairs (u,p) are obtained by using the last formulas (3.9) and (3.8).


### 3.2. The Second Order of Accuracy Difference Scheme

For the approximate solution of the problem (3.1), the Crank-Nicholson difference scheme
(3.16)unk-unk-1τ=un+1k-2unk+un-1k2h2+un+1k-1-2unk-1+un-1k-12h2-unk+unk-12+pk+pk-12qn+Dτ1/2u(tk-τ2,xn)+f(tk-τ2,xn),pk=p(tk),1≤k≤N,1≤n≤M-1un0=sin⁡(xn),0≤n≤M,u0k=uMk=0,0≤k≤N,usk+us+1k-uskh(x*-sh)=ρ(tk),0≤k≤N,s=⌊π2h⌋
is constructed.Here,(3.17)Γ(k-r+12)=∫0∞tk-r+1/2e-tdt.
Moreover, applying the second order of approximation formula for
(3.18)Dt1/2u(tk-τ2)=1Γ(1/2)∫0tk-τ/2‍(tk-τ2-s)-1/2u′(s)ds,
it is obtained (see [27])
(3.19)Dτ1/2u={-d23u0+d23u1+dτ32sin⁡(xn),k=1,-2d65u0+d65u1+d65u2-dτ610sin⁡(xn),k=2,d∑m=2k-1{[(k-m)b1+b2]um-2+[(2m-2k-1)b1-2b2]um-1+[(k-m+1)b1+b2]um}+d62[-uk-2-4uk-1+5uk],3≤k≤N.
Here and throughout the paper,
(3.20)Dt1/2u=Dτ1/2u(tk-τ2,xn),d=2πτ,b1=k-m+12-k-m-12,b2=-13((k-m+12)3/2-(k-m-12)3/2).
We search the solution of (3.16) in the following form:
(3.21)unk=ηkqn+wnk,
where
(3.22)ηk=∑i=1k‍pi+pi-12τ,1≤k≤N,η0=0.
We have that
(3.23)usk+us+1k-uskh(x*-sh)=ηk((1-x*-shh)qs+x*-shhqs+1)+(1-x*-shh)wsk+x*-shhws+1k=ρ(tk).
Let us denote
(3.24)y=x*-shh=x*h-⌊x*h⌋h,
where 0≤y<1. Then, one can write
(3.25)ηk=ρ(tk)-(1-y)wsk-yws+1k(1-y)qs+yqs+1.
So the values of (p(tk)+p(tk-1))/2,1≤k≤N can be obtained by the following formula:
(3.26)pk+pk-12=(ρ(tk)-ρ(tk-1))/τ-(1-y)((wsk-wsk-1)/τ)-y((ws+1k-ws+1k-1)/τ)(1-y)qs+yqs+1.
Let wr denote
(3.27)wr=[w0r⋮wMr](M+1)×1forr=0,1,…,N.
For k=1, one can show that w1 is the solution of the difference scheme
(3.28)wn1-wn0τ=wn+11-2wn1+wn-112h2+wn+10-2wn0+wn-102h2-wn1+wn02+(qn+1-2qn+qn-12h2-qn2)×(ρ(t1)-(1-y)ws1-yws+11(1-y)qs+yqs+1+ρ(t0)-(1-y)ws0-yws+10(1-y)qs+yqs+1)+d23(ρ(t1)-(1-y)ws1-yws+11(1-y)qs+yqs+1q(n)+wn1)-d23(ρ(t0)-(1-y)ws0-yws+10(1-y)qs+yqs+1q(n)+wn0)+dτ32q(n)+f(t1-τ2,xn),1≤n≤M-1,w01=wM1=0,wn0=sin⁡(xn),0≤n≤M.
We have the system of linear equations and we write them in the matrix form
(3.29)A1w1+B1w0=Dφ1,
where
(3.30)A1=[1000⋅0⋅⋅000ay1a0⋅l1c1⋅0000ay1a⋅l2c2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅000a⋅ls+acs⋅0000000⋅ls+1+y1cs+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅lM-1cM-1⋅ay1a0000⋅00⋅001](M+1)×(M+1),B1=[0000⋅0⋅⋅000av1a0⋅d1e1⋅0000av1a⋅d2e2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅000a⋅ds+aes⋅0000000⋅ds+1+v1es+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅dM-1eM-1⋅av1a0000⋅00⋅000](M+1)×(M+1).
Here, for any n=1,2,…,M-1,
(3.31)a=(-12h2),y1=(1τ+1h2+12-d23),v1=(-1τ+1h2+12+d23),ln=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)+d23qn,cn=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)+d23qn,dn=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)-d23qn,en=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)-d23qn,φ1=[0ϕ11⋮ϕM-110](M+1)×1,ϕn1=(sin⁡⁡(xn+1)-2sin⁡⁡(xn)+sin⁡⁡(xn-1)2h2-sin⁡⁡(xn)2)ρ(t1)+ρ(t0)(1-y)qs+yqs+1+d2qn3ρ(t1)-ρ(t0)(1-y)qs+yqs+1+dτ32qn+f(t1-τ2,xn)
and D is (M+1)×(M+1) identity matrix. Using (3.29), we can obtain that
(3.32)w1=A1-1(Dφ1-B1w0),w0={sin⁡xn}n=0M.
For k=2,  w2 is the solution of the difference scheme
(3.33)wn2-wn1τ=wn+12-2wn2+wn-122h2-wn2+wn12+wn+11-2wn1+wn-112h2+(qn+1-2qn+qn-12h2-qn2)(ρ(t2)-(1-y)ws2-yws+12(1-y)qs+yqs+1+ρ(t1)-(1-y)ws1-yws+11(1-y)qs+yqs+1)+d65(ρ(t2)-(1-y)ws2-yws+12(1-y)qs+yqs+1q(n)+wn2)+d65(ρ(t1)-(1-y)ws1-yws+11(1-y)qs+yqs+1q(n)+wn1)-2d65(ρ(t0)-(1-y)ws0-yws+10(1-y)qs+yqs+1q(n)+wn0)-dτ610q(n)+f(t2-τ2,xn),1≤n≤M-1,w02=wM2=0,wn0=sin⁡(xn),0≤n≤M.
The system of linear equations given above can be written in the matrix form
(3.34)A2w2+B2w1+C2w0=Dφ2,
where
(3.35)A2=[1000⋅0⋅⋅000ay2a0⋅g1h1⋅0000ay2a⋅g2h2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅gs+ahs⋅0000000⋅gs+1+y2hs+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅gM-1hM-1⋅ay2a0000⋅00⋅001](M+1)×(M+1),B2=[0000⋅0⋅⋅000av2a0⋅g1h1⋅0000av2a⋅g2h2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅gs+ahs⋅0000000⋅gs+1+v2hs+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅gM-1hM-1⋅av2a0000⋅00⋅000](M+1)×(M+1),C2=[0000⋅0⋅⋅0000z00⋅i1j1⋅00000z0⋅i2j2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅is+zjs⋅0000000⋅is+1js+1+z⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅iM-1jM-1⋅0z00000⋅00⋅000](M+1)×(M+1).
Here, for any n=1,2,…,M-1,
(3.36)a=-12h2,y2=1τ+1h2+12-d65,v2=-1τ+1h2+12-d65,gn=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)+d6qn(1-y)5((1-y)qs+yqs+1),in(s+1)thcolumn,hn=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)+d6qny5((1-y)qs+yqs+1),in(s+2)thcolumn,in=-2d6qn(1-y)5((1-y)qs+yqs+1),jn=-26qny5((1-y)qs+yqs+1),z=2d65,φ2=[0ϕ12⋮ϕM-120](M+1)×1,ϕn2=(sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)2h2-sin⁡(xn)2+d6qn5)×ρ(t2)+ρ(t1)(1-y)qs+yqs+1-2d6qn5ρ(t0)(1-y)qs+yqs+1-dτ610qn+f(t2-τ2,xn).
Using (3.34), we can obtain that
(3.37)w2=A2-1(Dφ2-B2w1-C2w0),w0={sin⁡xn}n=0M.
For 3≤k≤N, we can obtain the following difference scheme:
(3.38)wnk-wnk-1τ=wn+1k-2wnk+wn-1k2h2+wn+1k-1-2wnk-1+wn-1k-12h2-wnk+wnk-12+(qn+1-2qn+qn-12h2-qn2)×(ρ(tk)-(1-y)wsk-yws+1k(1-y)qs+yqs+1)+(qn+1-2qn+qn-12h2-qn2)×(ρ(tk-1)-(1-y)wsk-1-yws+1k-1(1-y)qs+yqs+1)+d∑m=2k-1‍{((2m-2k-1)b1-2b2)wnm-1((k-m)b1+b2)×(ρ(tm-2)-(1-y)wsm-2-yws+1m-2(1-y)qs+yqs+1q(n))+((2m-2k-1)b1-2b2)×(ρ(tm-1)-(1-y)wsm-1-yws+1m-1(1-y)qs+yqs+1q(n))+((2m-2k-1)b1-2b2)wnm-1×(ρ(tm)-(1-y)wsm-yws+1m(1-y)qs+yqs+1q(n))+((k-m)b1+b2)wnm-2+((k-m-1)b1+b2)wnm+((2m-2k-1)b1-2b2)wnm-1}-d62(ρ(tk-2)-(1-y)wsk-2-yws+1k-2(1-y)qs+yqs+1q(n)+wnk-2)-4d62(ρ(tk-1)-(1-y)wsk-1-yws+1k-1(1-y)qs+yqs+1q(n)+wnk-1)+5d62(ρ(tk)-(1-y)wsk-yws+1k(1-y)qs+yqs+1q(n)+wnk)+f(tk-τ2,xn),1≤n≤M-1,w0k=wMk=0,3≤k≤N,wn0=sin⁡(xn),0≤n≤M.
This system can be written in matrix form as
(3.39)A3wk+B3wk-1+C3wk-2+∑j=0k-3‍Ejwj=Dφk,
where
(3.40)A3=[1000⋅0⋅⋅000ay3a0⋅r1m1⋅0000ay3a⋅r2m2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅rs+ams⋅0000000⋅rs+1+y3ms+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅rM-1sM-1⋅ay3a0000⋅00⋅001](M+1)×(M+1),B3=[0000⋅0⋅⋅000aXa0⋅Y1Z1⋅0000aXa⋅Y2Z2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Ys+aZs⋅0000000⋅Ys+1+XZs+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅YM-1ZM-1⋅aXa0000⋅0⋅⋅000](M+1)×(M+1),C3=[0000⋅0⋅⋅0000S00⋅R1T1⋅00000S0⋅R2T2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Rs+STs⋅0000000⋅Rs+1Ts+1+S⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅RM-1TM-1⋅0S00000⋅0⋅⋅000](M+1)×(M+1),E0=[0000⋅0⋅⋅0000G00⋅H1I1⋅00000G0⋅H2I2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Hs+GIs⋅0000000⋅Hs+1Is+1+G⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅HM-1IM-1⋅0G00000⋅0⋅⋅000](M+1)×(M+1),E1=[0000⋅0⋅⋅0000G100⋅J1K1⋅00000G10⋅J2K2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Js+G1Ks⋅0000000⋅Js+1Ks+1+G1⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅JM-1KM-1⋅0G100000⋅0⋅⋅000](M+1)×(M+1),Ej=[0000⋅0⋅⋅0000G200⋅P1L1⋅00000G20⋅P2L2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Ps+G2Ls⋅0000000⋅Ps+1Ls+1+G2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅PM-1LM-1⋅0G200000⋅0⋅⋅000](M+1)×(M+1),
for j=2,3,…k-3.Here, for anyn=1,2,…,M-1,
(3.41)a=-12h2,d=2πτ,y3=1τ+1h2+12-5d62,v3=-1τ+1h2+12+4d62,αn=dq(n)(1-y)62((1-y)qs+yqs+1),βn=dq(n)y62((1-y)qs+yqs+1),rn=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)+5dq(n)(1-y)62((1-y)qs+yqs+1),sn=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)+5dq(n)y62((1-y)qs+yqs+1),X=v3-pp×d,Yn-1=γn-1+62pp×αn-1,Zn-1=δn-1+62pp×βn-1,S=d(162-pp-nn),Rn-1=-αn-1+62(nn+pp)αn-1,Tn-1=-βn-1+62(nn+pp)βn-1,γn=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)-4dq(n)(1-y)62((1-y)qs+yqs+1),δn=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)-4dq(n)y62((1-y)qs+yqs+1),mm=(k-m)b1+b2,nn=(2m-2k-1)b1-2b2,pp=(k-m+1)b1-2b2,nn1=(2(m+1)-2k-1)b1-2b2,pp1=(k-(m+2)+1)b1+b2,G=-d×mm,G1=-d×(mm+nn1),G2=-d×(mm+nn1+pp1),Hn-1=62×α1×mm,In-1=62×β1×mm,Jn-1=62×αn-1×(mm+nn1),Kn-1=62×βn-1×(mm+nn1),Pn-1=62×αn-1×(mm+nn1+pp1),Ln-1=62×βn-1×(mm+nn1+pp1),φk=[0ϕ1k⋮ϕM-1k0](M+1)×1,ϕnk=(sin⁡⁡(xn+1)-2sin⁡⁡(xn)+sin⁡⁡(xn-1)2h2-sin⁡⁡(xn)2)ρ(tk)+ρ(tk-1)(1-y)qs+yqs+1+5dqnρ(tk)62((1-y)qs+yqs+1)-4dqnρ(tk-1)62((1-y)qs+yqs+1)-dqnρ(tk-2)62((1-y)qs+yqs+1)+f(tk-τ2,xn)+dqn(1-y)qs+yqs+1×∑m=2k-1‍{((k-m)b1+b2)ρ(tm-2)+((2m-2k-1)b1-2b2)ρ(tm-1)+((k-m+1)b1+b2)ρ(tm)}.
Finally, from (3.39), it follows that
(3.42)wk=(A3)-1(Dφk-B3wk-1-C3wk-2-∑j=0k-3‍Ejwj),3≤k≤N.
Applying the last formula step by step, we can reach wk. Then, using (3.21), (3.25), and (3.26), we reach the approximate solutions of u(t,x) and (p(tk)+p(tk-1))/2.


### 3.3. Error Analysis

In this part, the results of the numerical analysis is given. The numerical solutions are recorded for different values ofN and M and unk represents the approximate solution of u(t,x) at grid points (tk,xn). Table 1 gives the error analysis between the exact solution and the solutions derived by difference schemes. Table 1 is constructed for N=M=15,45, and 75, respectively. For their comparison, the errors are computed by
(3.43)E=max⁡1≤k≤N1≤n≤M|u(tk,xn)-unk|.Table 1
Comparison of exact solution and approximate solutions.
Method
N = M = 15
N = M = 45
N = M = 75
1st order of accuracy
0.1190
0.0126
0.0045
2nd order of accuracy
0.0055
6.0917 × 10−4
2.1932 × 10−4Thus, the second order of accuracy difference scheme is more accurate comparing with the first order of accuracy difference scheme.


## 3.1. The First Order of Accuracy Difference Scheme

For the approximate solution of the problem (3.1), the Rothe difference scheme
(3.2)unk-unk-1τ=un+1k-2unk+un-1kh2-unk+Dτ1/2unk+pkqn+f(tk,xn),f(tk,xn)=(-3tk-1πtk-1/2+2πtk1/2)sin⁡xn,pk=p(tk),qn=sin⁡xn,xn=nh,tk=kτ,1≤k≤N,1≤n≤M-1,Mh=π,Nτ=1,un0=sin⁡xn,0≤n≤M,u0k=uMk=0,0≤k≤N,usk=ρ(tk),ρ(tk)=1-tk,0≤k≤N,s=⌊π2h⌋,
where ⌊x⌋ denotes greatest integer less than x is constructed. Throughout the paper, let us denote
(3.3)ρ(tk)=1-tk,qn=sin⁡xn,tk={tk=kτ,0≤k≤N,Nτ=1},xn={xn=nh,0≤n≤M-1,Mh=π},f(tk,xn)=(-3tk-1πtk-1/2+2πtk1/2)sin⁡xn,F(tk,xn)=ρ(tk)sin⁡⁡(xs)(sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)h2-sin⁡⁡(xn))-1π∑m=1k‍Γ(k-m+(1/2))(k-m)!τ1/2sin⁡⁡(xs)sin⁡(xn)+f(tk,xn).
We search the solution of (3.2) in the following form:
(3.4)unk=ηkqn+wnk,
where
(3.5)ηk=∑i=1k‍piτ,1≤k≤N,η0=0.
Moreover for the interior grid point usk, we have that
(3.6)usk=ηkqs+wsk=ρ(tk).
From (3.4), (3.5), and the condition usk=ρ(tk), it follows that
(3.7)ηk=ρ(tk)-wskqs,(3.8)pk=ηk-ηk-1τ,1≤k≤N,(3.9)unk=ρ(tk)-wskqsqn+wnk,0≤k≤N,0≤n≤M,
where wnk,0≤k≤N,0≤n≤M is the solution of the difference scheme
(3.10)wnk-wnk-1τ=wn+1k-2wnk+wn-1kh2-wnk-wsksin⁡⁡(xs)(sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)h2-sin⁡⁡(xn))-1π∑m=1k‍Γ(k-m+(1/2))(k-m)!(wsm-wsm-1sin⁡⁡(xs)τ1/2sin⁡⁡(xn)-wnm-wnm-1τ1/2)+F(tk,xn),1≤k≤N,1≤n≤M-1,w0k=wMk=0,0≤k≤N,wn0=sin⁡(xn),0≤n≤M.
First, applying the first order of accuracy difference scheme (3.10), we obtain (M+1)×(M+1) system of linear equations and we write them in the matrix form
(3.11)Awk+∑j=0k-1‍Bjwj=Dφk,1≤k≤N,w0={sin⁡(xn)}n=0M,
where
(3.12)A=[1000.0.000xyx0.z1.0000xyx.z2.000..........0000.zM-1.xyx0000.0.001](M+1)×(M+1),B0=[000.0.000a0.f1.0000a.f2.00........000.fs+a.00........000.fM-1.0a000.0.00](M+1)×(M+1),Bj=[000.0.000c0.d1.0000c.d2.00........000.ds+c.00........000.dM-1.0c000.0.00](M+1)×(M+1),
for any j=1,2,…,k-2, and
(3.13)Bk-1=[000.0.000v0.c1.0000v.c2.00........000.cs+v.00........000.cM-1.0v000.0.00](M+1)×(M+1).
Here, for any n=1,2,…,M-1,
(3.14)x=-1h2,y=1τ+2h2+1+1τ,zn=sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)sin⁡(xs)h2-sin⁡(xn)sin⁡(xs)-sin⁡(xn)sin⁡(xs)τ,in(s+1)thcolumn,a=Γ(k-1/2)τπ(k-1)!,fn=-sin⁡(xn)Γ(k-1/2)τπsin⁡(xs)(k-1)!,c=1τπ(Γ(k-j-1/2)(k-j-1)!-Γ(k-j+1/2)(k-j)!),dn=sin⁡(xn)τπsin⁡(xs)(Γ(k-j-1/2)(k-j-1)!-Γ(k-j+1/2)(k-j)!),v=-1τ-1τ,cn=sin⁡(xn)τsin⁡(xs),wr=[w0r⋮wMr](M+1)×1foranyr=0,1,…,k,φk=[0ϕ1k⋮ϕM-1k0](M+1)×1,ϕnk=(ρ(tk)sin⁡(xs)sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)h2-ρ(tk)sin⁡(xs)sin⁡(xn))-1π∑m=1k‍Γ(k-m+1/2)(k-m)!τ1/2sin⁡(xs)sin⁡(xn)+f(tk,xn),
and D is (M+1)×(M+1) identity matrix. Using (3.11), we can obtain that
(3.15)wk=A-1(Dφk-∑j=0k-1‍Bjwj),k=1,2,…,N,w0={sin⁡xn}n=0M.
To solve the resulting difference equations, we apply the method given in (3.15) step by step for k=1,2,…,N. For the evaluation of wr,r=2,3,…,N,wr-1 is needed. It is obtained in the previous step. Then, the solution pairs (u,p) are obtained by using the last formulas (3.9) and (3.8).


## 3.2. The Second Order of Accuracy Difference Scheme

For the approximate solution of the problem (3.1), the Crank-Nicholson difference scheme
(3.16)unk-unk-1τ=un+1k-2unk+un-1k2h2+un+1k-1-2unk-1+un-1k-12h2-unk+unk-12+pk+pk-12qn+Dτ1/2u(tk-τ2,xn)+f(tk-τ2,xn),pk=p(tk),1≤k≤N,1≤n≤M-1un0=sin⁡(xn),0≤n≤M,u0k=uMk=0,0≤k≤N,usk+us+1k-uskh(x*-sh)=ρ(tk),0≤k≤N,s=⌊π2h⌋
is constructed.Here,(3.17)Γ(k-r+12)=∫0∞tk-r+1/2e-tdt.
Moreover, applying the second order of approximation formula for
(3.18)Dt1/2u(tk-τ2)=1Γ(1/2)∫0tk-τ/2‍(tk-τ2-s)-1/2u′(s)ds,
it is obtained (see [27])
(3.19)Dτ1/2u={-d23u0+d23u1+dτ32sin⁡(xn),k=1,-2d65u0+d65u1+d65u2-dτ610sin⁡(xn),k=2,d∑m=2k-1{[(k-m)b1+b2]um-2+[(2m-2k-1)b1-2b2]um-1+[(k-m+1)b1+b2]um}+d62[-uk-2-4uk-1+5uk],3≤k≤N.
Here and throughout the paper,
(3.20)Dt1/2u=Dτ1/2u(tk-τ2,xn),d=2πτ,b1=k-m+12-k-m-12,b2=-13((k-m+12)3/2-(k-m-12)3/2).
We search the solution of (3.16) in the following form:
(3.21)unk=ηkqn+wnk,
where
(3.22)ηk=∑i=1k‍pi+pi-12τ,1≤k≤N,η0=0.
We have that
(3.23)usk+us+1k-uskh(x*-sh)=ηk((1-x*-shh)qs+x*-shhqs+1)+(1-x*-shh)wsk+x*-shhws+1k=ρ(tk).
Let us denote
(3.24)y=x*-shh=x*h-⌊x*h⌋h,
where 0≤y<1. Then, one can write
(3.25)ηk=ρ(tk)-(1-y)wsk-yws+1k(1-y)qs+yqs+1.
So the values of (p(tk)+p(tk-1))/2,1≤k≤N can be obtained by the following formula:
(3.26)pk+pk-12=(ρ(tk)-ρ(tk-1))/τ-(1-y)((wsk-wsk-1)/τ)-y((ws+1k-ws+1k-1)/τ)(1-y)qs+yqs+1.
Let wr denote
(3.27)wr=[w0r⋮wMr](M+1)×1forr=0,1,…,N.
For k=1, one can show that w1 is the solution of the difference scheme
(3.28)wn1-wn0τ=wn+11-2wn1+wn-112h2+wn+10-2wn0+wn-102h2-wn1+wn02+(qn+1-2qn+qn-12h2-qn2)×(ρ(t1)-(1-y)ws1-yws+11(1-y)qs+yqs+1+ρ(t0)-(1-y)ws0-yws+10(1-y)qs+yqs+1)+d23(ρ(t1)-(1-y)ws1-yws+11(1-y)qs+yqs+1q(n)+wn1)-d23(ρ(t0)-(1-y)ws0-yws+10(1-y)qs+yqs+1q(n)+wn0)+dτ32q(n)+f(t1-τ2,xn),1≤n≤M-1,w01=wM1=0,wn0=sin⁡(xn),0≤n≤M.
We have the system of linear equations and we write them in the matrix form
(3.29)A1w1+B1w0=Dφ1,
where
(3.30)A1=[1000⋅0⋅⋅000ay1a0⋅l1c1⋅0000ay1a⋅l2c2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅000a⋅ls+acs⋅0000000⋅ls+1+y1cs+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅lM-1cM-1⋅ay1a0000⋅00⋅001](M+1)×(M+1),B1=[0000⋅0⋅⋅000av1a0⋅d1e1⋅0000av1a⋅d2e2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅000a⋅ds+aes⋅0000000⋅ds+1+v1es+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅dM-1eM-1⋅av1a0000⋅00⋅000](M+1)×(M+1).
Here, for any n=1,2,…,M-1,
(3.31)a=(-12h2),y1=(1τ+1h2+12-d23),v1=(-1τ+1h2+12+d23),ln=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)+d23qn,cn=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)+d23qn,dn=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)-d23qn,en=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)-d23qn,φ1=[0ϕ11⋮ϕM-110](M+1)×1,ϕn1=(sin⁡⁡(xn+1)-2sin⁡⁡(xn)+sin⁡⁡(xn-1)2h2-sin⁡⁡(xn)2)ρ(t1)+ρ(t0)(1-y)qs+yqs+1+d2qn3ρ(t1)-ρ(t0)(1-y)qs+yqs+1+dτ32qn+f(t1-τ2,xn)
and D is (M+1)×(M+1) identity matrix. Using (3.29), we can obtain that
(3.32)w1=A1-1(Dφ1-B1w0),w0={sin⁡xn}n=0M.
For k=2,  w2 is the solution of the difference scheme
(3.33)wn2-wn1τ=wn+12-2wn2+wn-122h2-wn2+wn12+wn+11-2wn1+wn-112h2+(qn+1-2qn+qn-12h2-qn2)(ρ(t2)-(1-y)ws2-yws+12(1-y)qs+yqs+1+ρ(t1)-(1-y)ws1-yws+11(1-y)qs+yqs+1)+d65(ρ(t2)-(1-y)ws2-yws+12(1-y)qs+yqs+1q(n)+wn2)+d65(ρ(t1)-(1-y)ws1-yws+11(1-y)qs+yqs+1q(n)+wn1)-2d65(ρ(t0)-(1-y)ws0-yws+10(1-y)qs+yqs+1q(n)+wn0)-dτ610q(n)+f(t2-τ2,xn),1≤n≤M-1,w02=wM2=0,wn0=sin⁡(xn),0≤n≤M.
The system of linear equations given above can be written in the matrix form
(3.34)A2w2+B2w1+C2w0=Dφ2,
where
(3.35)A2=[1000⋅0⋅⋅000ay2a0⋅g1h1⋅0000ay2a⋅g2h2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅gs+ahs⋅0000000⋅gs+1+y2hs+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅gM-1hM-1⋅ay2a0000⋅00⋅001](M+1)×(M+1),B2=[0000⋅0⋅⋅000av2a0⋅g1h1⋅0000av2a⋅g2h2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅gs+ahs⋅0000000⋅gs+1+v2hs+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅gM-1hM-1⋅av2a0000⋅00⋅000](M+1)×(M+1),C2=[0000⋅0⋅⋅0000z00⋅i1j1⋅00000z0⋅i2j2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅is+zjs⋅0000000⋅is+1js+1+z⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅iM-1jM-1⋅0z00000⋅00⋅000](M+1)×(M+1).
Here, for any n=1,2,…,M-1,
(3.36)a=-12h2,y2=1τ+1h2+12-d65,v2=-1τ+1h2+12-d65,gn=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)+d6qn(1-y)5((1-y)qs+yqs+1),in(s+1)thcolumn,hn=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)+d6qny5((1-y)qs+yqs+1),in(s+2)thcolumn,in=-2d6qn(1-y)5((1-y)qs+yqs+1),jn=-26qny5((1-y)qs+yqs+1),z=2d65,φ2=[0ϕ12⋮ϕM-120](M+1)×1,ϕn2=(sin⁡(xn+1)-2sin⁡(xn)+sin⁡(xn-1)2h2-sin⁡(xn)2+d6qn5)×ρ(t2)+ρ(t1)(1-y)qs+yqs+1-2d6qn5ρ(t0)(1-y)qs+yqs+1-dτ610qn+f(t2-τ2,xn).
Using (3.34), we can obtain that
(3.37)w2=A2-1(Dφ2-B2w1-C2w0),w0={sin⁡xn}n=0M.
For 3≤k≤N, we can obtain the following difference scheme:
(3.38)wnk-wnk-1τ=wn+1k-2wnk+wn-1k2h2+wn+1k-1-2wnk-1+wn-1k-12h2-wnk+wnk-12+(qn+1-2qn+qn-12h2-qn2)×(ρ(tk)-(1-y)wsk-yws+1k(1-y)qs+yqs+1)+(qn+1-2qn+qn-12h2-qn2)×(ρ(tk-1)-(1-y)wsk-1-yws+1k-1(1-y)qs+yqs+1)+d∑m=2k-1‍{((2m-2k-1)b1-2b2)wnm-1((k-m)b1+b2)×(ρ(tm-2)-(1-y)wsm-2-yws+1m-2(1-y)qs+yqs+1q(n))+((2m-2k-1)b1-2b2)×(ρ(tm-1)-(1-y)wsm-1-yws+1m-1(1-y)qs+yqs+1q(n))+((2m-2k-1)b1-2b2)wnm-1×(ρ(tm)-(1-y)wsm-yws+1m(1-y)qs+yqs+1q(n))+((k-m)b1+b2)wnm-2+((k-m-1)b1+b2)wnm+((2m-2k-1)b1-2b2)wnm-1}-d62(ρ(tk-2)-(1-y)wsk-2-yws+1k-2(1-y)qs+yqs+1q(n)+wnk-2)-4d62(ρ(tk-1)-(1-y)wsk-1-yws+1k-1(1-y)qs+yqs+1q(n)+wnk-1)+5d62(ρ(tk)-(1-y)wsk-yws+1k(1-y)qs+yqs+1q(n)+wnk)+f(tk-τ2,xn),1≤n≤M-1,w0k=wMk=0,3≤k≤N,wn0=sin⁡(xn),0≤n≤M.
This system can be written in matrix form as
(3.39)A3wk+B3wk-1+C3wk-2+∑j=0k-3‍Ejwj=Dφk,
where
(3.40)A3=[1000⋅0⋅⋅000ay3a0⋅r1m1⋅0000ay3a⋅r2m2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅rs+ams⋅0000000⋅rs+1+y3ms+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅rM-1sM-1⋅ay3a0000⋅00⋅001](M+1)×(M+1),B3=[0000⋅0⋅⋅000aXa0⋅Y1Z1⋅0000aXa⋅Y2Z2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Ys+aZs⋅0000000⋅Ys+1+XZs+1+a⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅YM-1ZM-1⋅aXa0000⋅0⋅⋅000](M+1)×(M+1),C3=[0000⋅0⋅⋅0000S00⋅R1T1⋅00000S0⋅R2T2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Rs+STs⋅0000000⋅Rs+1Ts+1+S⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅RM-1TM-1⋅0S00000⋅0⋅⋅000](M+1)×(M+1),E0=[0000⋅0⋅⋅0000G00⋅H1I1⋅00000G0⋅H2I2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Hs+GIs⋅0000000⋅Hs+1Is+1+G⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅HM-1IM-1⋅0G00000⋅0⋅⋅000](M+1)×(M+1),E1=[0000⋅0⋅⋅0000G100⋅J1K1⋅00000G10⋅J2K2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Js+G1Ks⋅0000000⋅Js+1Ks+1+G1⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅JM-1KM-1⋅0G100000⋅0⋅⋅000](M+1)×(M+1),Ej=[0000⋅0⋅⋅0000G200⋅P1L1⋅00000G20⋅P2L2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅Ps+G2Ls⋅0000000⋅Ps+1Ls+1+G2⋅000⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅0000⋅PM-1LM-1⋅0G200000⋅0⋅⋅000](M+1)×(M+1),
for j=2,3,…k-3.Here, for anyn=1,2,…,M-1,
(3.41)a=-12h2,d=2πτ,y3=1τ+1h2+12-5d62,v3=-1τ+1h2+12+4d62,αn=dq(n)(1-y)62((1-y)qs+yqs+1),βn=dq(n)y62((1-y)qs+yqs+1),rn=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)+5dq(n)(1-y)62((1-y)qs+yqs+1),sn=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)+5dq(n)y62((1-y)qs+yqs+1),X=v3-pp×d,Yn-1=γn-1+62pp×αn-1,Zn-1=δn-1+62pp×βn-1,S=d(162-pp-nn),Rn-1=-αn-1+62(nn+pp)αn-1,Tn-1=-βn-1+62(nn+pp)βn-1,γn=(qn+1-2qn+qn-1)(1-y)2h2((1-y)qs+yqs+1)+qn(1-y)2((1-y)qs+yqs+1)-4dq(n)(1-y)62((1-y)qs+yqs+1),δn=(qn+1-2qn+qn-1)y2h2((1-y)qs+yqs+1)+qny2((1-y)qs+yqs+1)-4dq(n)y62((1-y)qs+yqs+1),mm=(k-m)b1+b2,nn=(2m-2k-1)b1-2b2,pp=(k-m+1)b1-2b2,nn1=(2(m+1)-2k-1)b1-2b2,pp1=(k-(m+2)+1)b1+b2,G=-d×mm,G1=-d×(mm+nn1),G2=-d×(mm+nn1+pp1),Hn-1=62×α1×mm,In-1=62×β1×mm,Jn-1=62×αn-1×(mm+nn1),Kn-1=62×βn-1×(mm+nn1),Pn-1=62×αn-1×(mm+nn1+pp1),Ln-1=62×βn-1×(mm+nn1+pp1),φk=[0ϕ1k⋮ϕM-1k0](M+1)×1,ϕnk=(sin⁡⁡(xn+1)-2sin⁡⁡(xn)+sin⁡⁡(xn-1)2h2-sin⁡⁡(xn)2)ρ(tk)+ρ(tk-1)(1-y)qs+yqs+1+5dqnρ(tk)62((1-y)qs+yqs+1)-4dqnρ(tk-1)62((1-y)qs+yqs+1)-dqnρ(tk-2)62((1-y)qs+yqs+1)+f(tk-τ2,xn)+dqn(1-y)qs+yqs+1×∑m=2k-1‍{((k-m)b1+b2)ρ(tm-2)+((2m-2k-1)b1-2b2)ρ(tm-1)+((k-m+1)b1+b2)ρ(tm)}.
Finally, from (3.39), it follows that
(3.42)wk=(A3)-1(Dφk-B3wk-1-C3wk-2-∑j=0k-3‍Ejwj),3≤k≤N.
Applying the last formula step by step, we can reach wk. Then, using (3.21), (3.25), and (3.26), we reach the approximate solutions of u(t,x) and (p(tk)+p(tk-1))/2.


## 3.3. Error Analysis

In this part, the results of the numerical analysis is given. The numerical solutions are recorded for different values ofN and M and unk represents the approximate solution of u(t,x) at grid points (tk,xn). Table 1 gives the error analysis between the exact solution and the solutions derived by difference schemes. Table 1 is constructed for N=M=15,45, and 75, respectively. For their comparison, the errors are computed by
(3.43)E=max⁡1≤k≤N1≤n≤M|u(tk,xn)-unk|.Table 1
Comparison of exact solution and approximate solutions.
Method
N = M = 15
N = M = 45
N = M = 75
1st order of accuracy
0.1190
0.0126
0.0045
2nd order of accuracy
0.0055
6.0917 × 10−4
2.1932 × 10−4Thus, the second order of accuracy difference scheme is more accurate comparing with the first order of accuracy difference scheme.



---

*Source: 276080-2012-09-10.xml*