0
$\begingroup$

I am reading section 9.3 in Leveque about stability. The explanations are very short and brief so I am asking for some help to understand and elaborations. The result of the section is that for FTCS method stability occures when $$dt/dx^2=k/h^2\leq1/2$$ Can someone please provide a more thorough explaination for his "proof"?

This is the book I am talking about: https://faculty.washington.edu/rjl/fdmbook/

$\endgroup$
5
  • $\begingroup$ Isn't this the same question that you asked before? $\endgroup$ – nicoguaro May 4 at 19:10
  • $\begingroup$ Yes I asked something similar, but deleted it to create this post $\endgroup$ – econmajorr May 4 at 19:53
  • $\begingroup$ But somebody pointed you where to find your answer there. $\endgroup$ – nicoguaro May 4 at 21:17
  • $\begingroup$ As far as I remember it was Von Neuman analysis which is something entirely different than what I am asking about here. In this question I am really all about the his proof $\endgroup$ – econmajorr May 4 at 21:30
  • $\begingroup$ You obtain that inequality using von Neumann analysis. That article presented your differential equation as example. $\endgroup$ – nicoguaro May 4 at 21:45
4
$\begingroup$

The explanation in the book does not use von Neumann analysis at all but the absolute stability regions and the eigenvalues of the discrete Laplacian operator. For the result you specifically mentioned we use the fact that the maximum eigenvalues is $$ \lambda_m \approx -\frac{4}{h^2} $$ from the expression given. We then want this eigenvalue to lie inside of forward Euler's absolute stability region $$ |1 + k \lambda | \leq 1. $$ This stability region implies the inequality expression given in the text of $$ -2 \leq -4 \frac{k}{h^2} \leq 0 $$ from which of course the condition $$ \frac{k}{h^2} \leq \frac{1}{2} $$ comes from. If that inequality looks unfamiliar these are simply the bounds of the forward Euler stability region that encompass the negative real line, which is a bit of the short cut.

$\endgroup$
2
  • $\begingroup$ One quick question: do we want the smallest eigenvalue of A to be within stability region or do we want them all to be in the stability region? :) $\endgroup$ – k.dkhk May 5 at 22:34
  • 1
    $\begingroup$ All of them need to be in the stability region. The eigenvalues for the Laplacian though all lie on the negative real axis so if we can keep the largest eigenvalue in magnitude in the stability region then the rest will be in there as well. $\endgroup$ – Kyle Mandli May 6 at 14:10
3
$\begingroup$

I believe Von Neumann's stability analysis would give you the answer here. Consider the heat transfer equation:

$$\frac{\partial \mathcal{T}}{\partial t} = \alpha \frac{\partial^{2} \mathcal{T}}{\partial x^{2}}$$

By using Forward Euler time integration and central difference in space discretization:

$$\mathcal{T}^{t+\Delta t}_{x} = \mathcal{T}^{t}_{x} + \tilde{\alpha} (\mathcal{T}^{t}_{x+\Delta x} + \mathcal{T}^{t}_{x-\Delta x} - 2 \mathcal{T}^{t}_{x})$$

Where:

$$\tilde{\alpha} = \frac{\alpha \Delta t}{\Delta x^{2}}$$

Now consider the fact that numerical solution $\mathcal{T}_{x}^{t}$ could be expanded by using the discrete Fourier series:

$$\mathcal{T}^{t}_{x} = \sum_{j=-N}^{N} \mathcal{A}_{j}(t) \exp(i k_{j}x)$$

Where:

$$k_{j} = \frac{\pi j}{L}$$

$L$ is the size of the computational domain (e.g. here a 1D line with length $L$) and $N = \frac{L}{\Delta x}$.

Note that $\mathcal{A}_{j}(t)$ shows the temporal amplitude.

Putting the discrete Fourier expansion into the discrete equation:

$$\mathcal{A}_{j}(t+\Delta t) \exp(ik_{j}x) = \mathcal{A}_{j}(t) \exp(ik_{j}x) + \tilde{\alpha} \mathcal{A}_{j}(t) (\exp(i k_{j}(x+\Delta x)) + \exp(i k_{j} (x-\Delta x)) -2 \exp(i k_{j} x))$$

or:

$$\frac{\mathcal{A}_{j}(t+\Delta t)}{\mathcal{A}_{j}(t)} = 1 + \tilde{\alpha}(\exp(i k_{j} \Delta x) + \exp(-i k_{j} \Delta x) -2)$$

We know that:

$$\exp(i k_{j} \Delta x) + \exp(-i k_{j} \Delta x) -2 = -4 \sin^{2}\Big( \frac{i k_{j} \Delta x}{2}\Big)$$

Finally:

$$\frac{\mathcal{A}_{j}(t+\Delta t)}{\mathcal{A}_{j}(t)} = 1 - 4 \tilde{\alpha} \sin^{2}\Big( \frac{i k_{j} \Delta x}{2}\Big)$$

For the numerical solution, in order to remain bounded in time, we need to make sure amplitude would not grow when we advance in time:

$$\frac{\mathcal{A}_{j}(t+\Delta t)}{\mathcal{A}_{j}(t)} \leq 1$$

or:

$$\Bigg |1 - 4 \tilde{\alpha} \sin^{2}\Big( \frac{i k_{j} \Delta x}{2}\Big)\Bigg| \leq 1$$

Due to the fact that $4 \tilde{\alpha} \sin^{2}\Big( \frac{i k_{j} \Delta x}{2}\Big)$ is always positive, we should have:

$$4 \tilde{\alpha} \sin^{2}\Big( \frac{i k_{j} \Delta x}{2}\Big) \leq 2$$

or finally:

$$\tilde {\alpha} \sin^{2}\Big( \frac{i k_{j} \Delta x}{2}\Big) \leq \frac{1}{2}$$

Note that $\sin^{2}\Big( \frac{i k_{j} \Delta x}{2}\Big) \leq 1$, so:

$$\tilde{\alpha} \leq \frac{1}{2}$$

or finally:

$$\frac{\alpha \Delta t}{\Delta x^{2}} \leq \frac{1}{2}$$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.