Solution of an Unsteady Diffusion 
Problem in Spherical 
Coordinates 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

Acknowledgement :

This problem was presented in the excellent book by W. F. Ramirez, Computational Methods for Process Simulation, Second Edition, 1997, Butterworth Heinemann, Oxford.

Solution using NDSolve

sole = 
      NDSolve[{D[c[x, t], {t, 1}] ==  
            1/x^2 D[x^2 D[c[x, t], {x, 1}], {x, 1}], 
c[
              x, 0] ==  Exp[-10^8 (1 - x)],  
          c[1, t] ==  
            1,  (D[c[x, t], {x, 1}] /. 
                x ->  10^-20) ==  0}, {c[x, t]}, {x, 
          10^-20, 1}, {t, 0, 1}]

Plotting the mass flux versus time at x=1.

<< Graphics`Graphics`

plt1 = 
      LogLinearPlot[
        Evaluate[D[c[x, t] /. sole, {x, 1}] /. x ->  1], {t, 0.001, 
          0.01}, PlotStyle ->  RGBColor[1, 0, 0], 
        PlotRange -> All, PlotPoints ->  10000, 
        DisplayFunction ->  Identity]

plt2 = 
      LogLinearPlot[
        Evaluate[D[c[x, t] /. sole, {x, 1}] /. x ->  1], {t, 0.01, 1}, 
        PlotStyle ->  RGBColor[1, 0, 0], PlotRange -> All, 
        PlotPoints ->  1000, 
        DisplayFunction ->  Identity]

Show[plt1, plt2, PlotRange -> All, AspectRatio ->  0.7, 
      DisplayFunction ->  $DisplayFunction, 
      AxesOrigin ->  {-3, 0}, Frame ->  True]

Plotting the concentrations versus time at x=0.91760 and x=0.44849.

plt3 = 
      Plot[c[x, t] /. sole /. x ->  0.91760, {t, 0, 0.4}, 
        PlotRange ->  All, 
        PlotStyle -> RGBColor[1, 0, 0]]

plt4 = 
      Plot[c[x, t] /. sole /. x ->  0.44849, {t, 0, 0.4}, 
        PlotRange ->  All, 
        PlotStyle -> RGBColor[1, 0, 0]]

Solution using Orthogonal Collocation

Six internal collocation points (N=6)

x[1] = 0.23046; x[2] = 0.44849; x[3] = 0.64235; x[4] = 0.80156; 
    x[5] = 0.91760; x[6] = 0.98418; x[7] = 1. ;

Matrix Q and (Q^(-1))

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

IQ = Inverse[Q]

Matrix C

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

Matrix A

A=Ci.IQ

Ci . IQ // MatrixForm

Matrix D

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

Matrix B

B=Di.IQ

Di . IQ // MatrixForm

DC = Table[D[(Subscript[c, i])[t], {t, 1}], {i, 1, 7}]

Ce = Table[(Subscript[c, i])[t], {i, 1, 7}]

Writting the system of ODEs to be solved

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

eq=MapThread[f,{DC,B.Ce}]

(Subscript[c, 7])[t_] := 1

init = Table[(Subscript[c, i])[0] ==  0, {i, 1, 6}]

EQ = Delete[Join[eq, init], -7]

Length[EQ]

sol = NDSolve[EQ, Table[(Subscript[c, i])[t], {i, 1, 6}], {t, 100}]

Plotting the mass flux versus time at x=1.

tbl = Table [A[[7, i]], {i, 1, 7}]

J[t_] := 
      tbl .  {(Subscript[c, 1])[t], (Subscript[c, 2])[t], (Subscript[c, 3])[t], (Subscript[c, 4])[t], (Subscript[c, 5])[t], (Subscript[c, 6])[t], 
            (Subscript[c, 7])[t]} /. sol

J[t]

plt5 = 
      LogLinearPlot[First[J[t]], {t, 0.005, 0.1}, 
        PlotStyle -> RGBColor[0, 1, 0], PlotRange -> All, 
        PlotPoints ->  10000, 
        DisplayFunction ->  Identity]

plt6 = 
      LogLinearPlot[First[J[t]], {t, 0.1, 1}, 
        PlotStyle -> RGBColor[0, 1, 0], PlotRange -> All, 
        PlotPoints ->  10000, 
        DisplayFunction ->  Identity]

Show[plt6, plt5, PlotRange -> All, AspectRatio ->  0.7, 
      DisplayFunction ->  $DisplayFunction, 
      AxesOrigin ->  {-3, 0}, Frame ->  True]

Plotting the concentrations versus time at x=0.91760 and x=0.44849.

plt7 = 
      Plot[(Subscript[c, 5])[t] /. sol, {t, 0, 0.4}, PlotRange ->  All, 
        PlotStyle -> RGBColor[0, 0, 1]]

plt8 = 
      Plot[(Subscript[c, 2])[t] /. sol, {t, 0, 0.4}, PlotRange ->  All, 
        PlotStyle -> RGBColor[0, 0, 1]]

Comparing both solutions. We get perfect agreement betwen the solution using NDSolve and the orthogonal collocation method. In both cases, the mass flux is equal to 0.7842 at t=0.1 and x=1.

J[t] /. sol /. t ->  0.1

D[c[x, t] /. sole, {x, 1}] /. x ->  1 /. 
      t -> 0.1

Show[plt5, plt6, plt1, plt2, AspectRatio ->  1.5, 
      DisplayFunction ->  $DisplayFunction]

Show[plt7, plt3]

Show[plt8, plt4]