2 D Heat Equation

We are going to examine a number of the different 2D Heat equation solvers. We need an intial condition with some detail:  the book has an M I am going to do an A.

Initial Conditions

We need an initial condition so lets make an image of the "A".

\(\(Export["A.bmp", 
Show[
Graphics[
            Text["A", {0, 
                0}, 
TextStyle -> {FontSize -> 
                    72}
]
]
], 
ImageSize -> {47, 
            66}
];\)\)

Now if we reimport it we will get a 0-1 matrix of the A.

IC  = (255 - Import["A.bmp"][[1, 1, All, 
                All, 1]])/
          255. ;
ListDensityPlot[
        IC, 
Mesh -> False];

Explicit Const-Coeff Solver

The explicit solver is simple.

ExplicitStep[UIn_, νx_, νy_] :=  
      Module[
{Dimx, Dimy, UOut, i, 
          j}, 
{Dimx, Dimy} = 
          Dimensions[UIn]; 
UOut  =  
          Table[0, {Dimx}, {Dimy}]; 
Do[
UOut[[i, 
              j]]  =  
            UIn[[i, 
                j]] 
 +  νx (UIn[[i - 1, j]] - 
                    2 UIn[[i, j]]  +  
                    UIn[[i + 1, 
                      j]]) 
 + νy (UIn[[i, j - 1]] - 
                    2 UIn[[i, j]]  +  
                    UIn[[i, 
                      j + 1]]) , 
{i, 
            2, Dimx - 1}, {j, 2, 
            Dimy - 1}]; 
UOut
]

It seems to take about 0.06 seconds per step.

Timing[ExplicitStep[
          IC, νx, νy];]
Dimensions[IC]

The stability criteria we developed for the explicit scheme was 
	νx + νy < 1/2.
Lets see what happens with some values just below the stability cutt-off.

{νx, νy}  =  
        Dimensions[IC]* 0.001;
νx + νy
UCurrent = IC;
UData  = 
        Table[
UCurrent  =  
            ExplicitStep[
              UCurrent, νx, νy]; 
ListPlot3D[ 
            UCurrent, 
Mesh -> 
              False, 
PlotRange -> {0, 
                1}, 
PlotLabel -> 
              n], 
{n, 1, 40}];

With a slightly larger values for the νx and νy we get the following.

{νx, νy}  =  
        Dimensions[IC]* 0.005;
νx + νy
UCurrent = IC;
UData  = 
        Table[
UCurrent  =  
            ExplicitStep[
              UCurrent, νx, νy]; 
ListDensityPlot[ 
            UCurrent, 
Mesh -> 
              False, 
ColorFunctionScaling -> 
              False, 
PlotLabel -> 
              n], 
{n, 1, 20}];

Implicit Const-Coeff Solver
Version 1 - Sparse Matrix

The implicit solver is also simple once you get the idea.  The simplest way to list of the {i, j} positions in our rectangular example is by flattening out the matrix.  As discussed in class this gives the labelling scheme
	{i, j} ⟶ n = i + (j - 1) dimx
or equivalently
	n ⟶ {i, j} where i = Mod[n, dimx, 1] and j = 1 + (n - i)/dimx.

Once we have decided that this is how we are going to count we just need to build the sparse implicit matrix and repeatedly use the sparse solver to take timesteps.  Building the sparse matrix is OK once we realize that each {i, j} has an equation.  If {i, j} is a boundary point the equation (for our Dirichlet boundary condition) is given by the boundary condition
	(Subscript[𝒰, i, j]^(n + 1)) = 0.
The right hand side "b" entry is a zero and the only entry in the matrix is a 1 on the diagonal.  If {i, j} is an internal point the equation is given by 
	\(TraditionalForm\`𝒰\_\(i, j\)\%\(n + 1\) - νx(𝒰\_\(i - 1, j\)\%\(n + 1\) - 
              2  𝒰\_\(i, j\)\%\(n + 1\) + 𝒰\_\(i + 1, j\)\%\(n + 1\)) - νy(𝒰\_\(i, j - 1\)\%\(n + 1\) - 2  𝒰\_\(i, j\)\%\(n + 1\) + 𝒰\_\(i, j + 1\)\%\(n + 1\)) = 𝒰\_\(i, j\)\%n\)
The right hand side "b" entry is (Subscript[𝒰, i, j]^n) and the matrix entries are 1 + 2  νx + 2  νy, -νx (twice), and -νy (twice).

Clear[SparseImplicitA, νx, νy, Dimx, Dimy, i, 
      j]
SparseImplicitA[{Dimx_, Dimy_}, νx_, νy_] :=  
      Module[
{Eqs, n, n1, i, j}, 
 (* 
          Setting up the Counting *) 
n[i_, j_]  = 
          i  +  (j - 1) Dimx;  
Eqs  =  
          Flatten[
Table[
n1  =  
                n[i, j]; 
If[
i == 
                    1  ||  j == 1 || i == Dimx || 
                  j == Dimy, 
{n1, n1} -> 
                  1.0, 
{
{n1, 
                      n1} -> 
                    1 + 2  νx + 2  νy, 
{n1,  
                      n[i - 1, 
                        j]} -> -νx, 
{n1,  
                      n[i + 1, 
                        j]} -> -νx, 
{n1,  
                      n[i, j - 
                          1]} -> -νy, 
{n1,  
                      n[i, j + 
                          1]} -> -νy
}
], 
{i, Dimx}, 
{j, 
                Dimy}]
]
]

We need to be able to convert between the matrix and vector representation of our data.  First converting a matrix to a vector.

Clear[MatrixToVec, Mat];
MatrixToVec[Mat_] :=  
      Module[
{Dimx, Dimy, Vec, n, i, 
          j}, 
{Dimx, Dimy} = 
          Dimensions[Mat]; 
 (* 
          Setting up the Counting *) 
n[i_, j_]  = 
          i  +  (j - 1) Dimx;  
Vec  =  
          Table[0, {Dimx×
                Dimy}]; 
Do[
Vec[[n[i, j]]] = 
            Mat[[i, 
              j]], 
{i, 
            Dimx}, 
{j, 
            Dimy}
]; 
Vec
]

Next converting a vector to a matrix.

Clear[VecToMatrix, Vec,  Dimx, Dimy]
VecToMatrix[Vec_, {Dimx_, Dimy_}] :=  
      Module[
{Mat, n, i, j}, 
 (* 
          Initializing a matrix *) 
Mat  =  
          Table[ 0, {Dimx}, {Dimy}]; 
 (* 
          Setting up the counting *) 
i[n_]  = 
          Mod[n, Dimx, 1]; 
j[n_] = 
          1 + (n - i[n])/
              Dimx; 
Do[
Mat[[i[n], j[n]]]  =  
            Vec[[
              n]], 
{n, 
            Dimy×
              Dimx}
]; 
Mat
]

After putting this all together its is easy to run the implicit scheme.  Here it is in action with a step-size ten times larger than we can take with the explicit scheme.

{νx, νy}  =  
        Dimensions[IC]* 0.05;
νx + νy
SparseMatA  =  
        SparseImplicitA[
          Dimensions[IC], νx, νy];
VecIC  =  MatrixToVec[IC];

With the current data sizes the solver takes roughly 1.18 seconds per step.  This compares to the speedy 0.06 seconds per step for the explicit scheme!  It is 20 times slower!

Timing[Hmm  = 
          Developer`SparseLinearSolve[SparseMatA, 
            VecIC];]
Dimensions[IC]

Here is the code doing a bunch of steps.

UCurrent = VecIC;
UData  = 
        Table[
UCurrent  =  
            Developer`SparseLinearSolve[SparseMatA, 
              UCurrent]; 
ListDensityPlot[
VecToMatrix[UCurrent, 
              Dimensions[IC]], 
Mesh -> 
              False, 
ColorFunctionScaling -> 
              False, 
PlotLabel -> 
              n
], 
{n, 
            10}];

Implicit Const-Coeff Solver
Version 2 - Normal Matrix

The implicit solver is also simple once you get the idea.  The simplest way to list of the {i, j} positions in our rectangular example is by flattening out the matrix.  As discussed in class this gives the labelling scheme
	{i, j} ⟶ n = i + (j - 1) dimx
or equivalently
	n ⟶ {i, j} where i = Mod[n, dimx, 1] and j = 1 + (n - i)/dimx.

Once we have decided that this is how we are going to count we just need to build the sparse implicit matrix and repeatedly use the sparse solver to take timesteps.  Building the sparse matrix is OK once we realize that each {i, j} has an equation.  If {i, j} is a boundary point the equation (for our Dirichlet boundary condition) is given by the boundary condition
	(Subscript[𝒰, i, j]^(n + 1)) = 0.
The right hand side "b" entry is a zero and the only entry in the matrix is a 1 on the diagonal.  If {i, j} is an internal point the equation is given by 
	\(TraditionalForm\`𝒰\_\(i, j\)\%\(n + 1\) - νx(𝒰\_\(i - 1, j\)\%\(n + 1\) - 
              2  𝒰\_\(i, j\)\%\(n + 1\) + 𝒰\_\(i + 1, j\)\%\(n + 1\)) - νy(𝒰\_\(i, j - 1\)\%\(n + 1\) - 2  𝒰\_\(i, j\)\%\(n + 1\) + 𝒰\_\(i, j + 1\)\%\(n + 1\)) = 𝒰\_\(i, j\)\%n\)
The right hand side "b" entry is (Subscript[𝒰, i, j]^n) and the matrix entries are 1 + 2  νx + 2  νy, -νx (twice), and -νy (twice).

Clear[ImplicitA, νx, νy, Dimx, Dimy, i, 
      j]
ImplicitA[{Dimx_, Dimy_}, νx_, νy_] :=  
      Module[
{Mat, n, n1, i, j}, 
 (* 
          Setting up the Counting *) 
n[i_, j_]  = 
          i  +  (j - 1) Dimx; 
 (* 
          Initializing the Matrix *) 
Mat  =  
          Table[ 0,  {Dimx Dimy},  {Dimx Dimy}];   
Do[
n1  =  
            n[i, j]; 
If[
i == 
                1  ||  j == 1 || i == Dimx || 
              j == Dimy, 
Mat[[n1, 
                n1]] = 
              1.0, 
{
Mat[[n1, n1]] = 
                1 + 2  νx + 
                  2  νy, 
Mat[[n1,  
                  n[i - 1, 
                    j]]] = -νx, 
Mat[[n1,  
                  n[i + 1, 
                    j]]] = -νx, 
Mat[[n1,  
                  n[i, j - 
                      1]]] = -νy, 
Mat[[n1,  
                  n[i, j + 
                      1]]] = -νy
}
], 
{i, 
            Dimx}, 
{j, 
            Dimy}]; 
Mat
]

After putting this all together its is easy to run the implicit scheme with the non-sparse storage.  Here it is in action with a step-size ten times larger than we can take with the explicit scheme.

{νx, νy}  =  
        Dimensions[IC]* 0.05;
νx + νy
MatA =  
        ImplicitA[
          Dimensions[IC], νx, νy];
VecIC  =  MatrixToVec[IC];

With the current data sizes the non-sparse solver takes roughly 6.20 seconds per step.  This compares to the speedy 0.06 seconds per step for the explicit scheme!  It is 103 times slower!

Timing[Hmm  = 
          LinearSolve[MatA, VecIC];]
Dimensions[IC]

66 47

The non-sparse algorithm is the same as the sparse algorithm.  The sparse code is just 6 times more efficient and uses a lot less memory.  It is really the memory that will kill us on larger problems.  The non sparse matrix occupies 132 times as much memory!

N[ByteCount[MatA]/ByteCount[SparseMatA]]

ADI Const-Coeff Solver

The Alternating Direction Implicit scheme with constant coefficients is very easy to implement.  We already have code that performs explicit and implicit steps in one dimension.  We just need to combine these two pieces of code into a single step.  It is pretty easy once you have realized that our 1-D codes will work with vector input.

1-D codes

Here are the 1-D codes snitched from the earlier notebooks.

Clear[TimeStepExp, 
      UIn, ν]

TimeStepExp[UIn_List, ν_] :=  
      Module[
{UOut, i}, 
 (* 
          UIn is the current u value with disc Δx *) 
 (* ν is the ν value *) 
 (* Initializing UOut *) 
UOut = 
          UIn; 
 (* 
          Setting the BC *) 
UOut[[1]] = 
          0×
            UOut[[1]]; 
UOut[[-1]] = 
          0×
            UOut[[-1]]; 
 (* 
          Stepping through the Data *) 
Do[
UOut[[i]]  = 
            UIn[[
                i]] + 
ν (UIn[[i - 1]]  -  
                    2 UIn[[i]] + 
                    UIn[[i + 
                        1]]), 
{i, 
            2,  Length[UIn] - 
              1}
]; 
UOut
]

Clear[TimeStepImp, UIn, ν]
Needs["LinearAlgebra`Tridiagonal`"]

TimeStepImp[UIn_List, ν_] :=  
      Module[
{Dia, LowDia, 
          HighDia}, 
Dia =  
          Join[ {1}, 
Table[
              1 + 2  ν, {Length[UIn] - 
                  2}], 
{1}]; 
LowDia = 
          Join[
-ν Table[
                1, {Length[UIn] - 
                    2}], 
{0}]; 
HighDia = Join[{0}, -ν Table[
                1, {Length[UIn] - 2}]]; 
TridiagonalSolve[
          LowDia, Dia, HighDia, UIn]
]

TimeStepθ[UIn_List, ν_, θ_] := 
      Module[
{b}, 
b =  
          TimeStepExp[
            UIn, (1 - θ) ν]; 
TimeStepImp[
          b, θ ν]
]

ADI

One ADI scheme is
	\(TraditionalForm\`𝒰\_\(i, j\)\%\(n + 1/2\) - \(1\/2\) \(νx(𝒰\_\(i - 1, j\)\%\(n + 1/2\) - 
                2  𝒰\_\(i, j\)\%\(n + 1/2\) + 𝒰\_\(i + 1, j\)\%\(n + 1/2\))\) = 𝒰\_\(i, j\)\%n + \(1\/2\) \(νy(𝒰\_\(i, j - 1\)\%n - 
                2  𝒰\_\(i, j\)\%n + 𝒰\_\(i, j + 1\)\%n)\)\)
	\(TraditionalForm\`𝒰\_\(i, j\)\%\(n + 1\)\ \ \  - \(1\/2\) \(νy(𝒰\_\(i, j - 1\)\%\(n + 1\) - 
                2  𝒰\_\(i, j\)\%\(n + 1\) + 𝒰\_\(i, j + 1\)\%\(n + 1\))\) = 𝒰\_\(i, j\)\%\(n + 1/2\) + \(1\/2\) \(νx(𝒰\_\(i - 1, j\)\%\(n + 1/2\) - 
                2  𝒰\_\(i, j\)\%\(n + 1/2\) + 𝒰\_\(i + 1, j\)\%\(n + 1/2\))\)\)
i.e. a 1/2 of an explicit step in y followed by 1/2 of an implicit step in x followed by 1/2 an explicit step in x with a final 1/2 an implicit step in y.  This looks complicated but it it really pretty easy to code using our 1-D code.

The (Subscript[A, 1]^2)

The quick

ADI1[UIn_List, νx_, νy_]  :=  
      Module[
{ννx, ννy, Q1Step, Q2Step, 
          Q3Step, Q4Step}, 
ννx = 
          0.5 νx;  ννy = 0.5 νy; 
  
Q1Step = 
          TimeStepExp[UIn, ννy]; 
  
  
Q1Step = Transpose[Q1Step]; 
  
Q2Step = 
          TimeStepImp[Q1Step, ννx]; 
Q3Step = 
          TimeStepExp[Q2Step, ννx]; 
  
Q3Step = Transpose[Q3Step]; 
  
TimeStepImp[
          Q3Step, ννy]
]

Testing

ListDensityPlot[IC]

Timing[Hmm  = ADI1[IC, 0.5, 0.5];]

1.2/0.06

Dimensions[Hmm]

ListDensityPlot[ Hmm]

The non-sparse algorithm is the same as the sparse algorithm.  The sparse code is just 6 times more efficient and uses a lot less memory.  It is really the memory that will kill us on larger problems.  The non sparse matrix occupies 132 times as much memory!

N[ByteCount[MatA]/ByteCount[SparseMatA]]