B4 - discretising heat equation PDEs

lax equivalence theorem

a consistent finite difference method applied to a well-posed linear initial value problem is convergent if and only if it is stable

for simple, homogeneous PDEs, stability is quantified by von neumann stability criterion, equivalent to checking if the amplitude of a fourier mode is stable
considering a heat equation
an approximation for the second derivative is needed
try fitting a quadratic function through three neighbouring points, or take the finite difference derivative of the first derivative:

\begin{matrix} \frac{\partial y}{\partial t} = α \frac{\partial^{2} y}{\partial x^{2}} \to \frac{y_{k}^{j + 1} - y_{k}^{j}}{Δ t} = α \frac{y_{k + 1}^{j} + y_{k - 1}^{j} - 2 y_{k}^{j}}{(Δ x)^{2}} \\ y_{k}^{j + 1} = y_{k}^{j} + \frac{α Δ t}{(Δ x)^{2}} [y_{k + 1}^{j} + y_{k - 1}^{j} - 2 y_{k}^{j}] \end{matrix}

considering an infinitely long rod with an initial temperature profile, $T (x, t_{0})$
taking the fourier transform at each numerical step, and checking that there are no growing modes, which would imply instability
a homogeneous PDE always has a plain wave solution
a fourier decomposition in space of a numerical solution is:

y_{k}^{j} = \sum_{m = 0}^{M - 1} \exp (i \frac{2 π Δ x}{L} m k) {\bar{y}}_{m}^{j}

where, ${\bar{y}}_{m}^{j}$ are fourier amplitudes of each mode, $m$

note: $k$ is an arbitrary grid point, and $m$ is a wave number of a mode
to find the time evolution for the fourier amplitude, considering just one mode, ie. one term in the sum:

y_{k}^{j} = \exp (i m q_{0} k) {\bar{y}}_{m}^{j}

where, $q_{0} = \frac{2 π Δ x}{L}$

\exp (i m q_{0} k) {\bar{y}}_{m}^{j + 1} = \exp (i m q_{0} k) {\bar{y}}_{m}^{j} + \frac{α Δ t}{(Δ x)^{2}} [\exp (i m q_{0} (k + 1)) {\bar{y}}_{m}^{j + 1} + \exp (i m q_{0} (k - 1)) {\bar{y}}_{m}^{j} - 2 \exp (i m q_{0} k) {\bar{y}}_{m}^{j}]

{\bar{y}}_{m}^{j + 1} = {\bar{y}}_{m}^{j} [1 + \frac{α Δ t}{(Δ x)^{2}} (e^{i m q_{0}} + e^{- i m q_{0}} - 2)] = {\bar{y}}_{m}^{j} [1 + \frac{2 α Δ t}{(Δ x)^{2}} (\cos (m q_{0}) - 1)]

the term made up of constants in $[\dots]$ is the amplification factor ( $g$ ), which determines the stability

g = 1 + \frac{2 α Δ t}{(Δ x)^{2}} [\cos (m q_{0}) - 1] = 1 - \frac{4 α Δ t}{(Δ x)^{2}} \sin^{2} (\frac{m q_{0}}{2})

for stability: $| g | \leq 1 \forall m, k$
the 'worst case scenario' is when $\sin^{2} (x) = 1$ , ie. $m q_{0} = π$
for this upper threshold, if $| g | \leq 1$ , it will remain so for all $k$ and $m$
the only non-trivial inequality is $| g | \leq 1 \to g \geq - 1$ , when the stability condition is: $\frac{α Δ t}{(Δ x)^{2}} \leq 1 / 2$
now, finding the most unstable mode, which has $m$ that maximises $g$
so, it can be found by solving $d_{m} g = 0$ such that $m q_{0} \sim π$ (excluding $m = 0$ )

m = \frac{L}{2 Δ x}

\begin{matrix} y_{k}^{j} = \exp (i m q_{0} k) {\bar{y}}_{m}^{j} = i π k \\ y_{k + 1}^{j} = \exp [i m (k + 1) q_{0}] {\bar{y}}_{m}^{j} = \exp [i π (k + 1)] \end{matrix}

so the phase changes by $π$
this is called grid-scale instability and it is purely numerical (as the heat equation is stable)
the short wavelength fluctuations must be resolved to make this method stable
as $Δ x \to 0$ , this method has poor performance, as $Δ t$ needs to be decreased to resolve fast time scales of short wavelengths
it does not converge due to the numerical instability on small scale
this is similar to what was found previously: for $\dot{y} = - γ y$ , $Δ t \leq 2 γ^{- 1}$
LTE analysis indicated that the explicit finite difference scheme is consistent and of order one