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;
Algorithm 1