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,(-Pi2)/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-> t2; |
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 |
Ai,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((Pi2-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:LinearAlgebraGenerateMatrix(evalf(SYS),vars): |
c:=linsolve(A,b): |
CoeffMatrix=Matrix(SIZE): |
cnt:=1; |
for i to SIZE do |
for j to SIZE do |
CoeffMatrixj,i:=ccnt: |
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(-Pi2*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; |