Heat Equation 1

Heat Equation (Variable Coeff)

Consider the heat equation 
	(∂_t) 
          u  = α (∂_xx) u		0 ≤ x ≤ L and 0 ≤ t ≤ T
with the initial condition
	u(x, 0)=(Subscript[u, 0])(x) 
and one of the following sets of boundary conditions:
(Dirichlet) 	u(0, t) = 0

      u(L, t) = 0 
or 
(Neumann) 	(∂_x) u(0, t) = 
        0

      (∂_x) u(L, t) = 0.  
or
(Periodic) 	u(0, t) = u(L, t)

Finite Difference Approximations for Derivatives

The fundamental idea is that we are going to split up or discretize space "x" and time "t" into equally spaced points.  The spacing between the points in the x-direction is Δx .  
              The spacing between the points in the t-direction is Δt.   We are going to try to compute a bunch of approximations
	𝒰\(TraditionalForm\`\(\(\_i\%j\)\(≈\)\(u(x\_i, 
          t\_j)\)\(\ \)\)\)
where \(TraditionalForm\`x\_i = \(i\ Δx\ \ and\ t\_j = 
          j\ \(Δt : \ \ we\ started\ both\ x\ and\ t\ at\ zero\ for\ \(\(convenience\)\(.\)\)\)\)\)  There are lots of different ways to approximate the various derivatives in the heat equation.  Just as an example we want to estimate the derivative (∂_t) u((Subscript[x, i]), (Subscript[t, j])) for the left hand side of the heat equation.  We could use either
	\(TraditionalForm\`\((Forward\ Diff)\)\ \ ∂\_t\( u(x\_i, 
                t\_j)\) ≈ \(δ\_t\%+\) 𝒰\_i\%\(j + 1/2\) := \(𝒰\_i\%\(j + 1\) - \ 𝒰\_i\%j\)\/Δt\)
\(TraditionalForm\`\((Backward\ Diff)\)\ \ ∂\_t\( u(x\_i, 
                t\_j)\) ≈ \(δ\_t\%-\) 𝒰\_i\%\(j - 1/2\) := \(𝒰\_i\%j - \ 𝒰\_i\%\(j - 1\)\)\/Δt\)
\(TraditionalForm\`\((Centered\ Diff)\)\ ∂\_t\( u(x\_i, 
                t\_j)\) ≈ \(δ\_t\) 𝒰\_i\%j := \(𝒰\_i\%\(j + 1\) - \ 𝒰\_i\%\(j - 1\)\)\/\(2  Δt\)
\)
or any average of these.

Similarly there are numerous possible approximations for (∂_xx) u but in this case the symmetric centered difference
	\(TraditionalForm\`∂\_xx\( u(x\_i, 
              t\_j)\) ≈ \(δ\_xx\) 𝒰\_i\%j := \(𝒰\_\(i + 1\)\%j\  - \ 2\ 𝒰\_i\%j\  + \ 𝒰\_\(i - 1\)\%j\)\/Δx\^2\)
seems to be the natural choice.  But since we are going to use this to step from j to j + 1it might also be appropriate to use the values at j + 1.

Finite Difference Equation (Explicit)

As a specific choice lets use the centered difference at level j for (∂_xx) u and the forward difference for (∂_t) u.  With this choice the heat equation is approximated by 
	\(TraditionalForm\`\(𝒰\_i\%\(j + 1\) - \ 𝒰\_i\%j\)\/Δt\) =  α(i Δx) \(TraditionalForm\`\(\(\(𝒰\_\(i + 1\)\%j\  - \ 
              2\ 𝒰\_i\%j\  + \ 𝒰\_\(i - 1\)\%j\)\/Δx\^2\)\(.\)\)\)
What we have really done is approximate (∂_t) u and (∂_xx) u and taken the approximate equation we get as a "definition" of the 𝒰's.  Believe it or not this is easy enough to code this but we have to decide what to do about boundary conditions (I am going to do Dirichlet for now) and pick a length L and discretization sizes Δx and Δt.

Dirichlet Boundary conditions

We are going to store the data in a vector which simply lists of all the (Subscript[𝒰, ∘]^j) values.  To enforce the boundary conditions we are going to make the first and last entries zero.

Time Step

A time step is going from (Subscript[𝒰, ∘]^j) to (Subscript[𝒰, ∘]^(j + 1)) i.e. advancing our computation by a time step of Δt.  The discretized equation tells us how to do this.

Clear[VarTimeStep1, 
      UIn, Δt, Δx, α]

VarTimeStep1[
        UIn_List, Δt_, Δx_, α_List] :=  
      Module[
{UOut, i, ν1}, 
 (* 
          UIn is the current u value with disc Δx *) 
 (* Δt is the time step α is a list of conduc. *) 
 (* 
          Initializing UOut *) 
UOut = 
          Table[ 0,  {Length[UIn]}]; 
 (* 
          Premultiplying Δ(t/Δ)(x^2) *) 
ν1 = Δ(t/Δ)(x^2); 
Do[
UOut[[
              i]]  = 
            UIn[[
                i]] + 
α[[i]] ν1 (UIn[[
                      i - 1]]  -  
                    2 UIn[[i]] + 
                    UIn[[i + 
                        1]]), 
{i, 
            2,  Length[UIn] - 
              1}
]; 
UOut
]

Testing the code!

The first thing we need is an initial condition.

Δx = 1;
Data  =  
      Table[ 12.5 - Abs[x - 12.5],  {x,  0,  25,  Δx}]; 
    ListPlot[Data];

Δt = 0.01;
α[x_] = 5  +  x;
αList = 
        Table[ α[(i - 0.5) Δx], {i, 
            Length[Data]}];
NewData  = 
      VarTimeStep1[
        Data, Δt, Δx, αList] 
ListPlot[ NewData];

Putting a number of time steps together.

Lets take 10 time steps and look at what happens.

Clear[NewData, 
        BigData];

Δt = 0.01;
α[x_] = 5  +  x;
αList = 
        Table[ α[(i - 0.5) Δx], {i, 
            Length[Data]}];

NewData  =  Data;
BigData  =  
        Table[
NewData = 
            VarTimeStep1[
              NewData, Δt, Δx, αList], 
{1000}];

\(\(ListContourPlot[\ 
        BigData, 
MeshRange -> {{0, 25}, {0, \ 
              10}}, 
FrameLabel -> {"x", "t"}
];\)\)

We had to take Δt=0.01 because otherwise it went unstable!

Finite Difference Equation (Implicit)

As a specific choice lets use the centered difference at level j + 1 for (∂_xx) u and the forward difference for (∂_t) u.  With this choice the heat equation is approximated by 
	\(TraditionalForm\`𝒰\_i\%\(j + 1\) - \ 𝒰\_i\%j\) =  ν \(TraditionalForm\`\((𝒰\_\(i + 1\)\%\(j + 1\)\  - \ 
          2\ 𝒰\_i\%\(j + 1\)\  + \ 𝒰\_\(i - 1\)\%\(j + 1\))\)\)
where \(TraditionalForm\`ν = \(α(
            i\ Δx)\)\ \ Δt/\(\(Δx\^2\)\(.\)\)\)
Again, we are really done approximatinge (∂_t) u and (∂_xx) u and taking the new approximate equation as a "definition" of the new 𝒰's.  This is still easy enough to code but we again have an issue with bondary conditions. We again pick a length L and discretization sizes Δx and Δt.

Equations

Rewriting the Finite Difference Equation to seperate the (known^j) values form the (unknown^(j + 1)) values gives 
	\(TraditionalForm\`\(-ν\_i\)\ 𝒰\_\(i - 1\)\%\(j + 1\)\  + \((1 + 
                2  ν\_i)\) 𝒰\_i\%\(j + 1\) - ν\_i\ 𝒰\_\(i + 1\)\%\(j + 1\) = 𝒰\_i\%\(j + 1\)\)
for all the non-boundary values of i.  The equations for the boundary points (specified by i = 1and i = N) are simpler since they just say that 	\(TraditionalForm\`𝒰\_0\%\(j + 1\) = \(𝒰\_0\%j = \(… = \(𝒰\_0\%0 = 
              0\)\)\)\)
	\(TraditionalForm\`𝒰\_N\%\(j + 1\) = \(𝒰\_N\%j = \(…\  = \(𝒰\_N\%0 = 
              0\)\)\)\)
Writing this in matrix form we have
	\(TraditionalForm\`A  𝒰\&⇀\^\(j + 1\) = 𝒰\&⇀\^j\)
where \(TraditionalForm\`A ∈ ℝ\^\(\((N + 1)\)×\((N + 1)\)\)\) is the square, tri-diagonal matrix
	A=(({{1, 0, 0, 0, …, 0}, {-ν, 1 + 2  ν, -ν, 0, …, 0}, {0, -ν, 1 + 2  ν, -ν, ⋱, ⋮}, {⋮, ⋱, ⋱, ⋱, ⋱, 0}, {0, …, 0, -ν, 1 + 2  ν, -ν}, {0, …, 0, 0, 0, 1}})).
Note:  
	A) The first and last rows are the Boundary Conditions;
	B) The middle rows are the Finite Difference Equations;
	C) The matrix eqn is simple to solve because it is tri-diagonal.

Time Step

A time step is going from (Subscript[𝒰, ∘]^j) to (Subscript[𝒰, ∘]^(j + 1)) i.e. advancing our computation by a time step of Δt.  The discretized matrix equation tells us how to do this.

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]
]

Testing the code!

The first thing we need is an initial condition.

Δx = N[π/24];
Data  =  Table[ Sin[x],  {x,  0,  π,  Δx}]; 
    ListPlot[Data];

Lets take 10 time steps and look at what happens.

Clear[NewData, 
        BigData];

Δt = 0.1;
α = 1;
ν = α Δt/Δ(x^2);

NewData  =  Data;
BigData  =  
        Table[
NewData = 
            TimeStepImp[
              NewData, ν], 
{10}];

ListPlot3D[ BigData];

Lets take 100 time steps and look at what happens.

Clear[NewData, 
        BigData];

Δt = 
        0.1;

NewData  =  Data;
BigData  =  
        Table[
NewData = 
            TimeStepImp[
              NewData, ν], 
{100}];

ListPlot3D[ 
        BigData, 
PlotRange -> All];

This one works just fine!