Today I wanted to share as series of interesting finite difference schema that is much less about the implementation and more about the mathematical underpinnings of the system. This is intended for someone who may be new to finite difference models yet found feel confident with the material I presented last time (Finite Difference Schema).
The family of schema presented here is a set of explicit approximations that represent an evolution of the Forward in Time Centered in Space scheme we saw before. Just as before, mathematical rigour will be sacrificed for the sake of conceptual understanding; but the conceptual underpinnings will hopefully reach much deeper this time around. Without further ado, let’s get started.
As far as notation goes, I’ll be using the standard super/sub script for description of two dimensions (space and time) which works quite well here. Nevertheless there are a variety of notational standards used elsewhere–especially when three and four dimensional extensions are formulated–so be cognizant that what you may see here might look quite different from elsewhere. Here’s an example of the notation to be used here.
[latex]U^{timestep}_{grid point} = U^{n}_{i} = f(n\cdot\Delta t,i\cdot\Delta x)[/latex]
CTCS
Looking back at CTCS, we’ll need to define a few more things this time around to help us in understanding and comparing the other schema. Remember that the CTCS scheme was stable so long as the CFL criteria was maintained (i.e. when [latex]\lambda = \frac{c\Delta t}{\Delta x}\leq 1[/latex]). The standard CTCS scheme can be written as follows where [latex]\lambda=c\cdot\Delta t/\Delta x[/latex].
[latex]U^{n+1}_i = U^{n-1}_i + \frac{\lambda}{2}(U^n_{i+1}-U^n_{i-1})[/latex]
We have already seen how this scheme produces an unstable approximation whereby the amplitude grows exponentially and the error (approximation solution – exact solution) grows rapidly in time. The way we will formally measure this trait is through the Gain (G). The gain is simply the ratio of how the signal grows (or shrinks over time). Since we cannot compare the signal height at different times directly due to phase error, we first employ a Fourier Decomposition to the scheme whereby we can compare signals directly.
[latex]U^n_i = \sum_kA_k^ne^{ikx_i} \rightarrow A^{n+1}_k=A^{n-1}_k-A^n_k[2i\lambda\sin(k\Delta x)][/latex]
Don’t let the maths get you down, what’s important here is that we can now compare apples to apples. If we now let [latex]a=2i\lambda\sin(k\Delta a)[/latex] we can look at our gain.
[latex]\left(\begin{array}{1} A^{n+1}_k \\ A^n_k \\ \end{array}\right) =\left(\begin{array}{cc} -a & 1 \\ 1 & 0\\ \end{array}\right)\left(\begin{array}{1} A^n_k \\ A^{n-1}_k \\ \end{array}\right) [/latex]
The matrix here is known as the amplitude matrix and it provides all of the information about the behaviour of the scheme. If we look at the eigenvalues of this matrix and then solve for [latex]\lambda[/latex] perhaps we can understand the stability of this scheme.
[latex]|G-\mu I| = \left|\begin{array}{cc} -a-\mu & 1 \\ 1 & -mu\\ \end{array}\right|=\mu^2+\mu a+1 = 0[/latex]
Now that we know that our eigenvalues [latex]\mu_{1,2} = \frac{1}{2}(-a\pm\sqrt{a^2+4}[/latex] we can carry out the standard analysis of the different cases (negative vs positive discriminant). This will give us several relations between [latex]\mu[/latex] and [latex]\lambda[/latex], and while I won’t show the math here, we end up with this:
[latex]\text{If }\lambda^2\sin^2(k\Delta x) \leq 1 \text{ then } |\mu_{1,2}|^2=1[/latex]
In other words our scheme is stable if [latex]\lambda \leq 1[/latex] (since a gain of greater than 1 is inherently unstable!). This requirement is precisely what we expected to find for CTCS since our model became unstable for values of [latex]\lambda > 1[/latex] (CFL Criteria).
FTCS
[latex]|G|^2=1+\lambda^2\sin^2(k\Delta x)[/latex]
Diffusive Scheme
Also know as the Friedrich or Euler Diffusive scheme, this scheme takes the FTCS scheme and adds a diffusive term in order to stabilize the model.
[latex]U^{n+1}_i = U^n_i-\frac{\lambda}{2}(U^n_{i+1}-U^n_{i-1}) + \frac{1}{2}(U^n_{i+1}+U^n_{i-1}-2U^n_i)[/latex]
The first term you’ll recognize as the FTCS scheme, and the last term is simply diffusion [latex]\left(\frac{\partial^2 f}{\partial x^2}\right)[/latex] which dampens and smooths out the instability inherent in the FTCS scheme. If we run the same sort of analysis on the amplification matrix we find:
[latex]G=\cos(k\Delta x) – i\lambda\sin(k\Delta x)[/latex]
Therefore the gain is
[latex]|G|^2=\cos^2(k\Delta x)+\lambda^2\sin^2(k\Delta x)[/latex]
which is less than 1 as long as [latex]|\lambda|< 1[/latex]. This is a great improvement over FTCS as you can see in the figure below. The gain of the Diffusive scheme with typical values, [latex]\lambda=0.9, k=0.5, \Delta x = 0.2[/latex], comes out to 0.9990 which is pretty good compared to the gain of the FTCS scheme which comes out to 1.0081. Remember that a gain of 1 is ideal and anything over 1 is strictly unstable.
Even with this, there is still room for improvement. Since [latex]\lambda[/latex] must remain less than 1, the small gap in the gain from 1.0 means that the signal will eventually dampen out and become flat over time. While this is typically better than blowing up, hopefully we can find a scheme which can minimize this dampening as much as possible while still remaining stable.
Upstream Differencing
Just as the last scheme this one is an extension of FTCS wherein diffusion is used to control the instabilities. Upstream Differencing, also know as the Donor Cell scheme, attempts to reduce the displacement of the gain from 1 as we shall see. The standard form is quite similar to the Diffusive Scheme except for the coefficient of the last term.
[latex]U^{n+1}_i = U^n_i-\frac{\lambda}{2}(U^n_{i+1}-U^n_{i-1}) + \frac{|\lambda|}{2}(U^n_{i+1}+U^n_{i-1}-2U^n_i)[/latex]
Looking at the gain again, we do see an improvement. Recall that the Diffusive Scheme gave us a “typical” gain of 0.9990, and now we are seeing a gain of 0.9991 with the Upstream Differencing. Well it might not be something to write home about, it’s a step in the right direction.
Lax-Wendroff
Sensing a trend in the two previous scheme, you might be include to ask whether there is simply a scheme that uses the best possible gain (i.e. the gain that minimizes [latex]1-|G|^2[/latex]), and that would be a good question. The answer is, naturally, yes and it is the Lax-Weidroff scheme. Once again the standard form replicates the priors and again tweaks the coefficient on the final term.
[latex]U^{n+1}_i = U^n_i-\frac{\lambda}{2}(U^n_{i+1}-U^n_{i-1}) + \frac{\lambda^2}{2}(U^n_{i+1}+U^n_{i-1}-2U^n_i)[/latex]
The remarkable aspect of the Lax-Wendroff scheme is that it is second order accurate in both time and space. Unlike all the previous methods that we’ve inspected which are first order accurate in time, this method will maintain it’s accuracy in phase for a much longer period of time.
Let’s see how are gain is improved this time around. Using the same values as before we now see a gain of 0.9999 and that is definitely an improvement! Since this method is mathematically the best we can do (at least for this lineage of development and improvement), we shouldn’t be surprised to see this pushing against the 1.0 limit. For ease of comparison, checkout the chart below.
Putting it all together
At least until I can get my image magick installation working again, here is the best way to understand the differences between the five schema we’ve just looked at. The FTCS signal is not displayed since it far outstripes the bounds and runs off towards [latex]\infty[/latex] incredibly fast. Notice the improvement of Lax-Weidroff over the other two familial schemes (Upstream Difference and Diffusive).
One thought on “Computational Fluid Dynamics: Improving FTCS”
Comments are closed.