Stability of the heat equation FD approximation¶
Calcuation of the stability condtions¶
The explicit (forward difference in time) formulation of the FD heat equation has a maximum time step that can be used in the calculations. Exceeding this time step will lead to the instability of the calculations, and typically, misbehaviour or crashing of the computational code (and, certainly, erroneous results).
Consider the simplified heat equation
\begin{equation} \frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2} \end{equation}
where \(\kappa = \frac{k}{\rho C_p}\). The FD approximation is
\begin{equation} \frac{T_i^{q+1}-T_i^q}{\Delta t} = \kappa \frac{T_{i+1}^{q} - 2T_{i}^{q} + T_{i-1}^{q}}{(\Delta x)^2} \end{equation}
and, with \(r=\frac{\kappa\Delta t}{\Delta x^2}\), can be written
\begin{equation} T_i^{q+1} = T_i^q + r (T_{i+1}^{q} - 2T_{i}^{q} + T_{i-1}^{q}) \end{equation}
If we assume that the numerical solution of the temperature, \(T\) is the sum of the correct temperature values \(\theta\) and the error \(\epsilon\) at each grid point (\(T_i = \theta_i + \epsilon_i\)), the difference equations of temperatures and errors can be separated, leading to
\begin{equation} \epsilon_i^{q+1} = \epsilon_i^q + r (\epsilon_{i+1}^{q} - 2\epsilon_{i}^{q} + \epsilon_{i-1}^{q}) \end{equation}
That is, the error diffuses from grid point to next grid point, similar to the actual temperature. We want to suppress the error, that is, we want
\begin{equation} |\epsilon_i^{q+1}| \leq |\epsilon_i^q| \Rightarrow \Big|\frac{\epsilon_i^{q+1}}{\epsilon_i^q}\Big| \leq 1, \end{equation}
since otherwise the error would grow from one time step to the next.
Dividing the error diffusion expression by \(\epsilon_i^q\) and applying the inequality
\begin{equation} \Big|1 + r (\epsilon_{i+1}^{q} - 2\epsilon_{i}^{q} + \epsilon_{i-1}^{q}) / \epsilon_i^q\Big| \leq 1 \end{equation}
Let’s consider the worst case scenario, where \(\epsilon_{i+1} = \epsilon_{i-1} = a\) and \(\epsilon_i = -a\). We get
\begin{equation} \Big|1 - 4r\Big| \leq 1 \Rightarrow 0 \leq r \leq \frac{1}{2} \end{equation}
The first inequality is naturally true (\(r\) cannot be negative), the second leads to
\begin{equation} \frac{\kappa \Delta t}{\Delta x^2} \leq \frac{1}{2} \end{equation}
This inequality can be used to determine a proper time step if grid spacing has been chosen, or proper grid spacing if time step has been chosen. Typically, it is safer to enforce an even stricter condition, e.g., \(\frac{\kappa \Delta t}{\Delta x^2} \leq \frac{1}{4}\).
Note: The stability of the method does not guarantee its correctness!
A physical interpretation exists for the stability criterion: The diffusive time scale \(\sqrt{\kappa\Delta t}\) should not exceed the grid spacing, i.e. heat should not diffuse in significant amounts over more than one grid point within one time step.