#C19.txt, April 3, 2014
#Numerical solutions of the Heat Equation
Help:=proc(): print(`Heq1(U,x,h,x1,t1), Heq1f(U,x,h,x1,t1)`):
print(` Heq2(U,x,h,x1,t1) `):
end:
#Heq1(U,x,h,k,x1,t1): Explicit method
#inputs an expression U of the variable x
#two small numbers h and k (the mesh sizes in x and t resp.)
#and a final time t1
#outputs an APPROXIMATION to the value of u(x1,t1)
#with x1 and t1 multiples of h and k.
#u(x,t) is the solution of the 1D-Heat equation
#u_{xx}=u_t, u(x,0)=U(x), u(0,t)=0, u(1,t)=0
Heq1:=proc(U,x,h,k,x1,t1) local T,j,n,N,J,r:
N:=t1/k:
J:=1/h:
r:=k/h^2:
for j from 0 to J do
T[j,0]:=subs(x=j*h,U):
od:
for n from 0 to N-1 do
T[0,n+1]:=0:
T[J,n+1]:=0:
for j from 1 to J-1 do
T[j,n+1]:=(1-2*r)*T[j,n]+r*T[j-1,n]+r*T[j+1,n]:
od:
od:
#T[i,j] is an approximation to u(i*h,k*j)
T[x1/h,t1/k]:
end:
#Heq1f(U,x,h,k,x1,t1): floating version of Heq1(U,x,h,k,x1,t1):
Heq1f:=proc(U,x,h,k,x1,t1) local T,j,n,N,J,r:
N:=t1/k:
J:=1/h:
r:=k/h^2:
for j from 0 to J do
T[j,0]:=evalf(subs(x=j*h,U)):
od:
for n from 0 to N-1 do
T[0,n+1]:=0:
T[J,n+1]:=0:
for j from 1 to J-1 do
T[j,n+1]:=evalf((1-2*r)*T[j,n]+r*T[j-1,n]+r*T[j+1,n]):
od:
od:
#T[i,j] is an approximation to u(i*h,k*j)
T[x1/h,t1/k]:
end:
#Heq2(U,x,h,k,x1,t1): Implicit method same input and output as in Heq1(U,x,h,k,x1,t1):
#inputs an expression U of the variable x
#two small numbers h and k (the mesh sizes in x and t resp.)
#and a final time t1
#outputs a list of lists such that
#L[i][j] is an APPROXIMATION to the value of u(x1,t1)
#of the temperature at x=i*h and t=j*k
#u(x,t) is the solution of the 1D-Heat equation
#u_{xx}=u_t, u(x,0)=U(x), u(0,t)=0, u(1,t)=0
Heq2:=proc(U,x,h,k,x1,t1) local T,j,n,N,J,r,eq,var:
N:=t1/k:
J:=1/h:
r:=k/h^2:
for j from 0 to J do
T[j,0]:=subs(x=j*h,U):
od:
for n from 0 to N-1 do
var:={seq(T[j,n+1], j=0..J)}:
eq:={ T[0,n+1]=0, T[J,n+1]=0}:
eq:= eq union
{seq((1+2*r)*T[j,n+1]-r*T[j-1,n+1]-r*T[j+1,n+1]=T[j,n],j=1..J-1)}:
var:=solve(eq,var):
for j from 0 to J do
T[j,n+1]:=subs(var,T[j,n+1]):
od:
od:
#T[i,j] is an approximation to u(i*h,k*j)
T[x1/h,t1/k]:
end: