Solution of Unsteady Heat Conduction in an Insulated Rod using the Orthogonal Collocation Method
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Orthogonal Collocation:
First and Second Derivative Matrices (A and B) using Shifted Legendre Polynimilas (Finlayson, The method of Weighted Residuals and Variational Principles, Academic Press, 1972)

Seven internal collocation points (N=7)

x[1] = 1.0 10^-20; x[2] = 0.0254; x[3] = 0.1292; x[4] = 0.2971; 
    x[5] = 0.5000; x[6] = 0.7029; x[7] = 0.8708; x[8] = 0.9746; 
    x[9] = 1. ;

Matrix Q and (Q^(-1))

Q = 
      Chop[Table[Table[x[j]^(i - 1), {i, 1, 9}], {j, 1, 9}]]

IQ = Inverse[Q]

Matrix C

Ci = 
      Chop[Table[
          Table[(i - 1) x[j]^(i - 2), {i, 1, 9}], {j, 1, 
            9}]]

Matrix A

A=Ci.IQ

Ci . IQ // MatrixForm

Matrix D

Di = 
      Chop[Table[
          Table[(i - 1) (i - 2) x[j]^(i - 3), {i, 1, 9}], {j, 
            1, 9}]]

Matrix B

B=Di.IQ

Di . IQ // MatrixForm

Unsteady Heat Conduction in an Insulated Rod:

Must solve following PDE:  		\(TraditionalForm\`\(\(			\)\(∂\_t\ 
          u = ∂\_\(x, x\)u\ \ \ \ \ 

      	\ \ \ \ \ with\ \ \(u(0, t)\) = \(u(1, t) = \(0\ and\ \(u(x, 0)\) = 
            1\)\)\)\)\)

(Ramirez, Computational Methods for Process Simulation, 2nd Edition, Butterworth-Heinemann,1997).

Tempeartures at internal collocation points are (Subscript[T, i]) [t]    for  i = 2.   .8

DT=Table[D[\!\(T\_i\)[t],{t,1}],{i,1,9}]

Te=Table[\!\(T\_i\)[t],{i,1,9}]

Writting the system of ODEs to be solved

f[a_,b_]:=a==b

eq=MapThread[f,{DT,B.Te}]

eq7=Take[eq,{2,8}]

Initial conditions. Initial temperature is equal to one.

init=Table[\!\(T\_i\)[0]==1,{i,1,9}]

init7=Take[init,{2,8}]

Tempeartures at boundary of rod are constant and they are equal to zero.

\!\(T\_1\)[t_]:=0

\!\(T\_9\)[t_]:=0

Equations and initial conditions

EQ=Join[eq7,init7]

Solving the set of ODEs

sol=NDSolve[EQ,Table[\!\(T\_i\)[t],{i,2,8}],{t,0.03}]

Tempeartures at collocation points at t=0.03

tbl1=Table[{x[i],\!\(T\_i\)[t]/.sol/.t->0.03},{i,1,9}]

tbl2=Map[Flatten,tbl1]

Tempeartures profile in the rod at t=0.03

plt1=ListPlot[tbl2,PlotStyle->{RGBColor[1,0,0],PointSize[0.02]}]

Solving the PDE using Mathematica

Boundary conditions must be set to (e^(-1000 t)) at x=0 and x=1
in order for NDSolve to operate properly. In fact, when the boundary conditions and initial conditions in a partial differential equation give different values for the solution at a single point, one must replace the discontinuous condition with a similar condition that changes quickly but smoothly.

solution = 
    NDSolve[{D[u[x,t],{t,1}]==D[u[x, t],{x,2}], 
        u[x, 0]==1, u[0, t]==Exp[-1000 t], u[1, t]==Exp[-1000 t]}, 
      u, {x, 0, 1}, {t,0,0.03}]

plt2=Plot[u[x,t]/.solution/.t->0.03,{x,0,1}]

Solutions using Mathematica NDSolve and orthogonal collocation method are identical

Show[plt1,plt2]