Abstract and Applied Analysis

Research Article

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

Algorithm 1

restart:
with(linalg):
with(LinearAlgebra):
N:=16:
d:=Pi/2:
h:=0.75/sqrt(N);
s:=.5/sqrt(N);
SIZE:=2*N+1;
delta0:= i,j ->piecewise(j=i,1,j<>i,0):
delta1:= i,j ->piecewise(i=j,0,i<>j,((-1) (i-j))/(i-j)):
delta2:= i,j ->piecewise(i=j,(-Pi 2)/3,i<>j,-2*(-1) (i-j)/(i-j) 2):
I_0:=Matrix(SIZE,delta0):
I_1:=Matrix(SIZE,delta1):
I_2:=Matrix(SIZE,delta2):
xk:=1/2+1/2*tanh(k*h/2);
xk_func:=unapply(xk,k);
phi:=unapply(log((x)/(1-x)),x);
Dphi:=unapply(simplify(diff(phi(x),x)),x);
g:=unapply(simplify(1/diff(phi(x),x)),x);
Dg:=unapply(diff(g(x),x),x);
gk:=unapply(subs(x=xk,g(x)),k);
Dgk:=unapply(subs(x=xk,Dg(x)),k);
g_Div_Dphi:=unapply(g(x)/Dphi(x),x);
g_Div_Dphi_k:=unapply(subs(x=xk,g_Div_Dphi(x)),k);
#Temporal  Spaces;
tl:=exp(l*s);
tl_func:=unapply(tl,l);
gamm:=unapply(log(t),t);
Dgam:=unapply(diff(gamm(t),t),t);
g_gamm:=unapply(1/diff(gamm(t),t),t);
g_gamm_l:=unapply(subs(t=tl,g_gamm(t)),l);
gamm_Div_Dgam:=unapply(g_gamm(t)/Dgam(t),t);
gamm_Div_Dgam:=  t->  t 2;
gamm_Div_Dgam_l:=unapply(subs(t=tl,gamm_Div_Dgam(t)),l);
GenerateDiagonalAm  :=  proc(  x  )
            local  i:=1:
              local  A:=Matrix(SIZE):
               for  i  from  1  by  1  to  SIZE
                do
                A i,i :=evalf(x(-N+i-1)):
                end  do:
               return  A;
          end  proc;
B:=  -2*h*MatrixMatrixMultiply(I_0,GenerateDiagonalAm(gk))+
   MatrixMatrixMultiply(I_1,GenerateDiagonalAm(Dgk))+1/h*I_2:
DD:=h*GenerateDiagonalAm(g_Div_Dphi_k):
C:=MatrixMatrixMultiply(GenerateDiagonalAm(g_gamm_l),s*(I_0-I_1)):
E:=s*GenerateDiagonalAm(gamm_Div_Dgam_l):
Fkl:=unapply((Pi 2-4)*sin(Pi*xk_func(k-N-1))*exp(-tl_func(l-N-1)),k,l);
F:=evalf(Matrix(SIZE,Fkl)):
V:=Matrix(SIZE,v):
EQN_SYS:=evalf(
 MatrixMatrixMultiply(
        MatrixMatrixMultiply(
          Matrix(inverse(DD)),B
                ),V
       )
     )
  +evalf(
    MatrixMatrixMultiply(
      MatrixMatrixMultiply(
          V,C),Matrix(inverse(E)
                 )
               )
     ):
SYS:= :
   for  i  from  1  by  1  to  SIZE
      do
      for  j  from  1  by  1  to  SIZE
        do
           SYS:= op(SYS),EQN_SYS(i,j)=F(i,j) ;
        end  do:
      end  do:
vars:=seq(seq(v(i,j),i=1..2*N+1),j=1..2*N+1):
A,b:LinearAlgebra GenerateMatrix (evalf(SYS), vars ):
c:=linsolve(A,b):
CoeffMatrix=Matrix(SIZE):
cnt:=1;
for  i  to  SIZE  do
  for  j  to  SIZE  do
   CoeffMatrix j,i :=c cnt :
   cnt:=cnt+1
  end  do;
end  do;
CoeffMatrix:=Matrix(CoeffMatrix,SIZE):
ApproximateSol:=unapply(
            evalf(
              sum(
               sum("CoeffMatrix"[m+N+1,n+N+1]
          *sin(Pi*(phi(x)-m*h)/h)/(Pi*(phi(x)-m*h)/h)
         *sin(Pi*(gamm(t)-n*s)/s)/(Pi*(gamm(t)-n*s)/s)
             ,m=-N..N),
            n=-N...N)
             )
               +exp(-4*t)*sin(Pi*x)
           ,x,t):
plot3d(ApproximateSol(x,t),x=0..1,t=0..1):
Exact:=unapply(exp(-Pi 2*t)*sin(Pi*x),x,t);
plot3d(Exact(x,t),x=0..1,t=0..1);
plot3d({Exact(x,t),ApproximateSol(x,t)},x=0..1,t=0..1);
XX:=.6;
#Numerical  Comparision  EXACT
for  j  from  0.1  to  10  by  1
   do
  evalf(Exact(XX,j)):
   od;
#Numerical  Comparision  APPROX
for  j  from  0.1  to  10  by  1
   do
  evalf(ApproximateSol(XX,j)):
   od;
#Numerical  Comparision  ERROR
for  j  from  0.1  to  10  by  1
   do
  abs(evalf(evalf(ApproximateSol(XX,j)-Exact(XX,j)))):
   od;