I am performing a stability analysis of the 1D heat equation:
$$
\frac{\partial u}{\partial t} = k\frac{\partial ^{2}u}{\partial x^{2}},
$$
Which I have discretised using a forward euler in time and a central difference in space to give the explicit discrete form:
$$
U^{n+1}_{j}=D(U^{n}_{j+1}-2U^{n}_{j}+U^{n}_{j-1})+U^{n}_{j}.
$$
Where $D=k\frac{\Delta t}{\Delta x^{2}}$.
I have started by assuming a solution of the form $U^{n}_{j}=a^{(n)}(\omega)e^{ij\omega \Delta x}$, where $\omega$ is the wave frequency, $i$ is the imaginary unit, and $\Delta x$ is the spatial time step size. Substituting the solution into the discrete equation and then ploughing through some algebra leaves me with the inequality in question:
$$
-2 \leqslant 2D[cos(\omega \Delta x)-1] \leqslant 0 \hspace{3cm} for \hspace{3mm}0 \leqslant \omega \Delta x \leqslant \pi.
$$
To be clearer my question is: How do I solve this inequality to get a condition on $D$?
Best Answer
For ease, split the inequality into two: $$ -2\leqslant 2D[cos(\omega \Delta x)-1] \hspace{1cm}\mathbf{and}\hspace{1cm}2D[cos(\omega \Delta x)-1]\leqslant 0 \hspace{3mm} for \hspace{1mm}\omega \Delta x \in[0,\pi] $$ Now make a substitution $\gamma=cos(\omega \Delta x)-1 \hspace1mm for \hspace{1mm}\omega \Delta x \in[0,\pi]$. The inequality above now becomes: $$ -2\leqslant 2\gamma D \hspace{1cm}\mathbf{and}\hspace{1cm}2\gamma D\leqslant 0 \hspace{3mm} for \hspace{1mm}\gamma \in[-2,0] $$ Now since D is the diffusion coefficient $D=k\frac{\Delta t}{\Delta x^{2}}$ it will always be positive, and for the range we have on $\gamma$ it is always negative. Hence the right hand inequality is always true, which makes the problem much simpler.
$-1 \leqslant \gamma D \hspace{3mm} for \hspace{1mm}\gamma \in[-2,0)$. Note: the range for $\gamma$ can no longer include $0$ otherwise we will have a division by 0.
Rewrite as $D\leqslant \frac{-1}{\gamma} \hspace{3mm} for \hspace{1mm}\gamma \in[-2,0)$. Now we need to think qualitatively about what this represents. Clearly $\frac{-1}{\gamma} \rightarrow \infty \hspace{2mm} as \hspace{2mm} \gamma \rightarrow 0$. So the tighter inequality on $D$ is produced when $\gamma =-2$ which is $D\leqslant \frac{1}{2}$.
In summary, the numerical scheme is stable if and only if $k\frac{\Delta t}{\Delta x^{2}} \leqslant \frac{1}{2}$.