# Numerical Solution and Simulation of Second-Order Parabolic PDEs with Sinc-Galerkin Method Using Maple

Authors: Aydin Secer
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

An efficient solution algorithm for sinc-Galerkin method has been presented for obtaining numerical solution of PDEs with Dirichlet-type boundary conditions by using Maple Computer Algebra System. The method is based on Whittaker cardinal function and uses approximating basis functions and their appropriate derivatives. In this work, PDEs have been converted to algebraic equation systems with new accurate explicit approximations of inner products without the need to calculate any numeric integrals. The solution of this system of algebraic equations has been reduced to the solution of a matrix equation system via Maple. The accuracy of the solutions has been compared with the exact solutions of the test problem. Computational results indicate that the technique presented in this study is valid for linear partial differential equations with various types of boundary conditions.


---

## Body

## 1. Introduction

Sinc methods for differential equations were originally introduced by Stenger in [1–3]. The sinc functions were first analyzed in [4, 5] and a detailed research of the method for two-point boundary-value problems can be found in [6, 7]. In [8], parabolic and hyperbolic problems are presented in detail. To solve a problem arising from chemical reactor theory, the properties of the sinc-Galerkin method are used to reduce the computation of nonlinear two-point boundary-value problems to some algebraic equations in [9]. A computer algorithm for sinc method to solve numerically the linear and nonlinear ODEs and their simulations has been presented in [7, 10], respectively. The full sinc-Galerkin method is developed for a family of complex-valued partial differential equations with time-dependent boundary conditions [9]. A study of the performance of the Galerkin method using sinc basis functions for solving Bratu’s problem is presented in [11]. In [12] a numerical algorithm has been presented for recovering the unknown function and obtaining a solution to the inverse ill-posed problem. They have presented a Galerkin method with the sinc basis functions in both space and time domains for solving the direct problem. A sinc-collocation method has been developed for solving linear systems of integrodifferential equations of Fredholm and Volterra type with homogeneous boundary conditions in [13].


## 2. Sinc-Approximation Formula for PDEs

We use the sinc-Galerkin method as mentioned in [1] to derive an approximate solution of the following:
(1)ut-uxx=F(x,t),u(0,t)=u(1,t)=0,0<x<1,u(x,0)=f(x),t>0.
For the equation given above, the sinc-Galerkin scheme can be developed in both space and time directions as follows.In general, approximations can be constructed for infinite, semi-infinite, and infinite intervals and both spatial and time spaces will be introduced. Define the function(2)ϕ(z)=ln(z1-z)
which is a conformal mapping from DE, the eye-shaped domain in the z-plane, onto the infinite strip, DS, where
(3)DE={z=x+iy:|arg(z1-z)|<d≤π2}.
A more general form of sinc basis according to intervals can be given as follows:
(4)S(m,hx)∘ϕ(x)=Sinc(ϕ(x)-mhxhx),m=-Nx,…,Nx,S(k,ht)∘γ(t)=Sinc(γ(t)-khtht),k=-Nt,…,Nt,
where
(5)Sinc(z)={sin(πz)πz,z≠01,z=0,Sinc(k,h)(z)=Sinc(z-khh)={sin(π((z-kh)/h))π((z-kh)/h),z≠kh1,z=kh,k=0,∓1,∓2,∓3,…
and the conformal maps for both directions
(6)ϕ(x)=ln(xl-x),x∈(0,l),γ(t)=ln(t),t∈(0,∞)
are used to define the basis functions on the intervals (0,l) and (0,∞), respectively. hx,ht>0 represents the mesh sizes in the space direction and the time direction, respectively. The sinc nodes xi and tj are chosen so that  xi=ϕ-1(ihx), tj=γ-1(jht).Here the functionx=ϕ-1(x)=ex/(1+ex) is an inverse mapping of ϕ=ϕ(x). We may define the range of ϕ-1 on the real line as
(7)Γ1={ϕ-1(u)∈DE:-∞<u<∞}.
For the evenly spaced nodes {kh}k=-∞∞ on the real line, the image which corresponds to these nodes is denoted by
(8)xk=ϕ-1(kh)=ekh1+ekh,
where 0<xk<1, for all k.The sinc basis functions in (4) do not have a derivative when x tends to 0 or 1. We modify the sinc basis functions as
(9)S(m,hx)∘ϕ(x)ϕ′(x)=Sinc((ϕ(x)-mhx)/hx)ϕ′(x),
where
(10)1ϕ′(x)=x(1-x).
For the temporal space, we construct an approximation by defining the function
(11)w=γ(r)=ln(r)
which is a conformal mapping from DW, the wedge-shaped domain in the r-plane, onto the infinite strip, DS, where
(12)DW={r=t+is:|arg(r)|<d<π2},
derived from composite translated functions
(13)S(k,ht)∘γ(t)=Sinc(γ(t)-khtht),k=-Nt,…,Nt,
for r∈DW.Herew=γ(r) and  γ-1(w)=r=ew. We may define γ-1 on the real line as
(14)Γ2={γ-1(u)∈Dw:-∞<u<∞}.
For the evenly spaced nodes {kh}k=-∞∞ on the real line, the image which corresponds to these nodes is denoted by
(15)tk=γ-1(kh)=ekh,
where 0<tk<∞, for all k.A list of conformal mappings may be found in Table1 [14].Table 1
Conformal mappings and nodes for several subintervals ofR.
(
a
,
b
)
ϕ (z)
z
k
a
b
ln
(
z
-
a
b
-
z
)
a
+
b
e
k
h
1
+
e
k
h
0
1
ln
(
z
1
-
z
)
e
k
h
1
+
e
k
h
0
∞
ln
(
z
)
e
k
h
0
∞
ln
(
sinh
(
z
)
)
ln
(
e
k
h
+
e
2
k
h
+
1
)
-
∞
∞
z
k
h
-
∞
∞
sinh
-
1
(
z
)
k
hDefinition 1.
LetB(DE) be the class of functions F that are analytic in DE and satisfy
(16)∫ψ(L+u)|F(z)|dz⟶0,asu=∓∞,
where
(17)L={iy:|y|<d≤π2}
and on the boundary of DE satisfy
(18)T(F)=∫∂DE|F(z)dz|<∞.
The proof of the following theorems can be found in [1].Theorem 2.
LetΓ be (0,1), F∈B(DE), and then for h>0 sufficiently small
(19)∫ΓF(z)dz-h∑j=-∞∞F(zj)ϕ′(zj)=i2∫∂DF(z)k(ϕ,h)(z)sin(πϕ(z)/h)dz≡IF,
where
(20)|k(ϕ,h)|z∈∂D=|e[(iπϕ(z)/h)sgn(Imϕ(z))]|z∈∂D=e-πd/h.For the sinc-Galerkin method, the infinite quadrature rule must be truncated to a finite sum; the following theorem indicates the conditions under which exponential convergence results.Theorem 3.
If there exist positive constantsα,β, and C such that
(21)|F(x)ϕ′(x)|≤C{e-α|ϕ(x)|,x∈ψ(-∞,∞)e-β|ϕ(x)|,x∈ψ(0,∞),
then the error bound for the quadrature rule (19) is
(22)|∫ΓF(x)dx-h∑j=-NNF(xj)ϕ′(xj)|≤C(e-αNhα+e-βNhβ)+|IF|.
The infinite sum in (19) is truncated with the use of (20) to arrive at (22).
Making the selections(23)h=πdαN,N≡∥αNβ+1∥,
where ∥·∥ is integer part of statement,
(24)∫ΓF(x)dx=h∑j=-NNF(xj)ϕ′(xj)+O(e-(παdN)1/2).Theorems2 and 3 can be used to approximate the integrals that arise in the formulation of the discrete systems corresponding to two-point BVPs.


## 3. Discrete Solutions Scheme for Two-Point BVPs

In ordinary differential equations(25)Lu=f
on Γ1, sinc solution is assumed as an approximate solution um in the form of series with m=2N+1 terms
(26)um(z)=∑j=-NNcjS(j,h)∘ϕ(z).
The coefficients {cj}j=-NN are determined by orthogonalizing the residual Lu-f with respect to the sinc basis functions {Sk}k=-NN where Sk(z)=S(k,h)∘ϕ(z). An inner product for two continuous functions such as f1 and f2 can be given by the following formula
(27)〈f1,f2〉=∫Γf1f2wdz,
where w is the weight function and is chosen depending on boundary conditions. If we implement the above inner product rule in orthogonalization, this yields the discrete sinc-Galerkin system:
(28)∫Γ(Lum-f)(z)S(k,h)∘ϕ(z)·w(z)dz=0,-N≤k≤N.
Now, we are going to derive discrete sinc-Galerkin system for PDEs. Assume umz,mt is the approximate solution of (1).Then, the discrete system takes the following form:
(29)umz,mt(z,t)=∑j=-NN∑k=-NNcjkS(j,h)∘ϕ(z)·S(k,h)∘γ(t).
The coefficients {cjk}j,k=-NN are determined by orthogonalizing the residual Lumz,mt-f with respect to the sinc basis functions {SkSh}k,h=-NN where SjSh(z,t)=S(j,h)∘ϕ(z)S(k,h)∘γ(t) for -N≤j,k≤N. In this case the inner product takes the following form:
(30)〈f1,f2〉=∫Γt∫Γzf1(z,t)f2(z,t)w(z,t)dzdt.
The choice of the weight function w(z,t) in the double integrand depends on the boundary conditions, the domain, and the partial differential equation. Therefore, the discrete Galerkin system is
(31)∫Γt∫Γz(Lumzmt-f)(z,t)·S(j,h)∘ϕ(z)·S(k,h)∘γ(t)·w(z,t)dzdt=0.


## 4. Matrix Representation of the Derivatives of Sinc Basis Functions at Nodal Points

The sinc-Galerkin method actually requires the evaluated derivatives of sinc basis functions at the sinc nodes,z=zj. The rth derivative of Sk(z)=S(k,h)∘ϕ(z) with respect to ϕ, evaluated at the nodal point zj, is denoted by
(32)1hrδpj(r)=drdϕr(S(p,h)∘ϕ(z))|z=zj.
The expressions in (14) for each p and j can be stored in a matrix  I(r)=[δpj(r)]. For r=0,1,2,…(33)I(0)=δjk(0)=[S(j,h)∘ϕ(x)]|x=xk={1,k=j0,k≠j,I(1)=δjk(1)=hddϕ[S(j,h)∘ϕ(x)]|x=xk={0,k=j(-1)k-j(k-j),k≠j,I(2)=δjk(2)=hd2dϕ2[S(j,h)∘ϕ(x)]|x=xk={-π23,k=j-2(-1)k-j(k-j)2,k≠j,
where
(34)Im(0)=[100⋯0010⋯0001⋯0:::⋱:000⋯1]=[δjk(0)],Im(1)=[0-112⋯12N10-1⋯-12N-1-1210⋯12N-2:::⋱:-12N12N-112N-2⋯0]=[δjk(1)],Im(2)=[-π23212-222⋯-2(2N)2212-π23212⋯2(2N-1)2-222212-π23⋯-2(2N-2)2:::⋱:-2(2N)22(2N-1)2-2(2N-2)2⋯-π23]=[δjk(2)].
The chain rule has been used for the z-derivative of product sinc functions. For example, when Sk(z)=S(k,h)∘ϕ(z),
(35)d(Sj(z)w(z))dz=(dSj(z)dϕ(z)·dϕ(z)dz)w(z)+Sj(z)dw(z)dz=dSj(z)dϕϕ′(z)w(z)+Sj(z)w′(z),d2(Sj(z)w(z))dz2=ddz(dSj(z)dϕϕ′(z)w(z)+Sj(z)w′(z))=d2Sj(z)dϕ2(ϕ′(z))2w(z)+dSj(z)dϕϕ′′(z)w(z)+2·dSj(z)dϕϕ′(z)w′(z)+Sj(z)w′′(z).
Now, we are going to develop discrete form for (1). We choose for special case the parameters as follows for the spatial dimension:
(36)ϕ(z)=ln(z1-z),wX(z)=1ϕ′(z),1ϕ′(z)=z(1-z),
and for the temporal space as
(37)γ(t)=ln(t),wT(t)=1γ′(t),1γ′(t)=t.
The discrete form of (1) can be given the following form:
(38)〈Lu-F,Sk·Sl〉=∫Γt∫Γz(Lu-F)S(k,h)∘ϕ(z)·wX(x)S(l,s)∘γ(t)·wT(t)dzdt=∫Γt∫Γz(ut-uxx-F)S(k,h)∘ϕ(z)·wX(x)S(l,s)∘γ(t)·wT(t)dzdt.
We solve this by taking our approximating basis functions to be
(39)Sk(x)=wXS(k,h)∘ϕ(x),wX=1ϕ′(x)=x(1-x),ϕ(x)=ln(x1-x),Sl(t)=wTS(l,s)∘γ(t),wT=1γ′(t)=t,γ(t)=ln(t).
If we apply sinc-quadrature rules with the help of (32)–(37) on the definite integral given (38) by using (39), we can get the following matrix system.LetAm(u) denote a diagonal matrix, whose diagonal elements are u(x-N),u(x-N+1),…,u(xN) and nondiagonal elements are zero. Then (38) reproduces the following matrixes accordingly.Firstly we set the coefficient matrix as follows:(40)C=(c-N,-Nc-N,-N+1c-N,-N+2⋯c-N,Nc-N+1,-Nc-N+1,-N+1c-N+1,-N+2⋯c-N+1,Nc-N+2,-Nc-N+2,-N+1c-N+2,-N+2⋯c-N+2,N:::⋱:cN,-NcN,-N+1cN,-N+2⋯cN,N),B=-2hIm(0)(Am(wX))+Im(1)(Am(wX′))+Im(2)h,G=Am(wT)[sIm(0)-Im(1)],D=hAm(wXϕ′),E=sAm(wTγ′).
Finally, for the right side function F given (1) can be written in the following matrix form:
(41)F=(F-N,-NF-N,-N+1F-N,-N+2⋯F-N,NF-N+1,-NF-N+1,-N+1F-N+1,-N+2⋯F-N+1,NF-N+2,-NF-N+2,-N+1F-N+2,-N+2⋯F-N+2,N:::⋱:FcN,-NFN,-N+1FN,-N+2⋯FN,N).
Using (32)–(37) we arrive at a matrix system given in [1] as follows:
(42)D-1BC+CGE-1=F.
Finally, by using Maple Computer Algebra Software, the matrix system (42) can be solved by using LU or QR decomposition method and unknown coefficients can be found. After calculation of C we get approximate solution as follows:
(43)ux,t=∑j=-NN∑k=-NNcjkS(j,h)∘ϕ(x)·S(k,h)∘γ(t).


## 5. Numerical Simulation

The example in this section will illustrate the sinc method.Example 4.
This problem has been addressed in [1]. The following equation is given in Dirichlet-type boundary condition:
(44)∂u(x,t)∂t-∂2u(x,t)∂x2=0,0<x<1,0<t,u(0,t)=u(1,t)=0,u(x,0)=sin(πx).
The particular solution of (44) can be calculated via separation variables rules and can be given as follows:
(45)u(x,t)=e-π2tsin(πx).
For (44) we choose sinc components here in the following:
(46)h=s=0.75N,xk=ekh1+ekh,tl=esl,ϕ(x)=ln(x1-x),γ(t)=ln(t),wX=1ϕ′(x),wT=1γ′(t).
According to the above parameters, the approximate solution simulation of (44) has been given in Figure 1 and numerical results also can be found in Table 2.Table 2
Numerical results.
t
Exact solution
Sinc-Galerkin solution
Error
N
=
16, x=0.6
0.1
0.35446621870000
0.59102026517764600
0.23655404647764600
0.11
0.00001833412226
−0.08327827240548060
0.08329660652774060
0.21
0.9482992112
×
10
-
9
−0.03068110846399060
0.03068110941228980
0.31
0.4904905672
×
10
-
13
−0.01083286679129120
0.01083286679134030
0.41
0.2536973471
×
10
-
17
−0.000420168677914081
0.000420168677914083
0.51
0.1312203514
×
10
-
21
−0.00108069489771883
0.00108069489771883
0.61
0.6787134677
×
10
-
26
−0.00492511086695295
0.00492511086695295
0.71
0.3510522275
×
10
-
30
0.00331706758955819
0.00331706758955819
0.81
0.1815753976
×
10
34
−0.00303031084399243
0.00303031084399243
0.91
0.9391658013
×
10
-
39
0.00342519240062265
0.00342519240062265Simulation of approximate solution.
(a)
The sinc-Galerkin solutions according to the grid points sizeN
(b)
Intersection of surfaces forN=20. The left-side figure is normal perspective and the second one has medium perspective view
(c)
The exact solution (or particular solution)


## 6. Conclusions

We have developed a Maple algorithm to solve and simulate second-order parabolic PDEs with Dirichlet-type boundary conditions based on sinc-Galerkin approximation on some closed real intervals and the method has been compared with the exact solutions. When compared with other computational approaches, this method turns out to be more efficient in the sense that selection parameters and changing boundary conditions and also giving different problems to the algorithms. The accuracy of the solutions improves by increasing the number of sinc grid pointsN. The method presented here is simple and uses sinc-Galerkin method that gives a numerical solution, which is valid for various boundary conditions. Several PDEs have been solved by using our technique in less than 20 seconds. All computations and graphical representations have been prepared automatically by our algorithm.



---

*Source: 686483-2013-08-13.xml*