I want to make sure I am clear with all the reasons why oscillations and how wild they get. I will summarize what I understood from the books with respect to the application of the Crank-Nicolson to the parabolic PDE. Consider the PDE $$u_t=au_x+bu_{xx}$$ and some initial data $f(x)$. Assume there are some boundary conditions but they are irrelevant for discussion.

Problem #1: If the initial data is rough, step function for example, then the oscillations might occur, but those are "controlled" and do decay, but slowly. So, I will not have an explosive solution but close to the region of the discontinuity I might have wiggles. However, if the time step is fine enough I can prevent them.

Problem #2: if $a>>b$, or even of $b=0$, regardless of the initial condition, using central differences for the $au_x$ term will bring me trouble(implicit or explicit in time). This problem has nothing to do with kinks in the initial data. And this issue is a more serious one as the solution does explode. Regardless how small my time and space grid is I will have oscillations. Would I ever have problems if the initial data is very smooth?

Did I get this two common problems correctly? I just want to make sure they have distinct roots! Thanks!


1 Answer 1


The two problems you mention indeed have different roots. The first one is a consequence of numerical dispersion and the other one of numerical instability. Let me elaborate a little:

Problem #1: The PDE in question is non-dispersive, i.e. waves of any frequency and wavenumber will travel equally fast - at speed $a$. However, after discretisation, this is no longer true. In fact, for the Crank-Nicolson scheme, the wavespeed is reduced somewhat to (roughly) $$ a\left( 1 - \frac{2 + \lambda^2}{12}(\xi \Delta x)^2 \right), $$ where $\lambda = \Delta t / \Delta x$ is the time to space mesh ratio, and $\xi$ is the wavenumber of a propagating wave. This means that waves of higher wavenumber will travel at a slower wavespeed than waves of lower wavenumbers. Now, assume that our initial solution can be expressed in terms of a Fourier transform. Then, if it has sharp gradients, it will typically have contributions from a large range of wavenumbers (i.e. both high and low), and hence, the solution "decomposes" into a train of waves travelling at different speeds. This typically looks like the oscillations you describe. These oscillations are no cause of instability and will not grow in time. In fact, if $b \geq 0$ they are dissipated along with the rest of the solution. Clearly, a small $\Delta x$ and $\Delta t$ reduces the problem.

Problem #2: The issue you describe here is due to numerical instability. This is seen e.g. for the FTCS (Forward in Time, Central in Space) discretisation which is a common text book example. Von Neumann analysis shows that this discretisation is unconditionally unstable, even for very smooth initial data. However, don't dismiss central differences as useless. There are ways to modify the discretisation such that central differences can be used "almost everywhere", except near the domain boundaries, without experiencing any instability.

Update: Below is an example of a stencil that uses central differences almost everywhere, yet is stable for the following advection problem with $a > 0$: $$ u_t + au_x = 0 \\ u(x,0) = f(x) \\ u(0,t) = g(t). $$ Let us discretise in space with a difference matrix $D = P^{-1}Q$ such that the semi-discrete problem becomes $$ \mathbf{u}_t + aP^{-1}Q\mathbf{u} = \sigma P^{-1} (u_0 - g(t)) \mathbf{e}_0. $$ Here, $\mathbf{u}=(u_0, u_1, \dots, u_{n-1}, u_n)^T$ is the semi-discrete solution vector defined on a mesh with grid spacing $\Delta x$. The right-hand side is known as a simultaneous approximation term (SAT). It weakly imposes the boundary condition for $u_0$. Here, $\mathbf{e}_0 = (1,0,0,\dots)^T$ and $\sigma$ is a scalar that we will soon specify. The matrices involved are given by $$ P = \Delta x \, \text{diag}(1/2, 1, 1, \dots, 1, 1, 1/2), $$ $$ Q = \begin{pmatrix} -1/2 & 1/2 & 0 & 0 & 0 & \dots \\ -1/2 & 0 & 1/2 & 0 & 0 & \dots \\ 0 & -1/2 & 0 & 1/2 & 0 & \dots \\ & \ddots & \ddots & \ddots & \ddots & \\ & \dots & 0 & -1/2 & 0 & 1/2 \\ & & 0 & 0 & -1/2 & 1/2 \end{pmatrix}. $$

To show that this discretisation is stable, note first that $P$ is symmetric and positive definite. Note further that $Q$ satisfies the so called summation-by-parts (SBP) property $$ Q + Q^T = \text{diag}(-1, 0, \dots, 0, 1). $$ Now, let us perform a stability analysis:

Define the norm $\| \mathbf{u} \|^2 = \mathbf{u}^T P \mathbf{u}$. Then it follows (I'll leave a few details for you to fill in) that \begin{align} \| \mathbf{u} \|^2_t &= \mathbf{u}_t^T P \mathbf{u} + \mathbf{u}^T P \mathbf{u}_t \\ &= -a \mathbf{u}^T (Q + Q^T) \mathbf{u} + 2 \sigma u_0 (u_0 - g(t)) \\ &= -a u_n^2 + u_0^2 (a + 2 \sigma) - 2 \sigma u_0 g(t). \end{align} Choosing $\sigma = -a$, then adding and subtracting $a g(t)^2$ gives the final result $$ \| \mathbf{u} \|^2_t = a \left( g(t)^2 - u_n^2 - (u_0 - g(t))^2 \right) \leq a g(t)^2. $$ Thus, the norm of the solution is bounded by the boundary data and hence, by definition, the discretisation is stable.

Conclusion: In summary, this discretisation uses the central difference stencil everywhere except at the first and last grid points. It is stable, and will remain stable upon discretisation in time, provided a suitable $\Delta t$ is chosen.

  • $\begingroup$ thank you very much for your answer. I got an idea and look on those topics more now that I have some name terms. Can you provide some references to how central differences still can be used with Pr#2? VN analysis says it is unconditionally unstable, so how could I make it stable? Also, in the second problem, would non smooth initial data make the problem more pronounced? Why in textbooks the examples are usually for rough initial data? $\endgroup$
    – Kamil
    Jul 25, 2016 at 13:59
  • $\begingroup$ I have updated my answer with an example of a stencil that uses central differences almost everywhere while remaining stable. For unstable discretisations, typically the errors are more pronounced if the solution has steep gradients. This is probably why the initial data in your text books is chosen to be "rough", i.e. to simplify illustration. $\endgroup$
    – ekkilop
    Jul 25, 2016 at 18:16
  • $\begingroup$ but for Pr#2, If I used CN there is no numerical instability anymore as amplification factor<=1, but there is still dispersion(due to central differences), so don't the problems come from the same root of the FDM being dispersive and PDE itself is not? In that case I see the same reason for oscillations in a case of a>>b or a<b. So is not Pr#2 a subproblem of Pr#1? $\endgroup$
    – Kamil
    Jul 26, 2016 at 12:37
  • $\begingroup$ and in this case, as above comment, the solution will not explode, however I will observe oscillations, so does it mean they will be more severe in Pr#2? $\endgroup$
    – Kamil
    Jul 26, 2016 at 12:57
  • $\begingroup$ If you have discretised your problem in a stable way and still see oscillation then sure, dispersion may be the issue. This would be the case if the wavenumbers involved are large and the mesh is coarse. However, without specifying things like what boundary conditions you use and how they are imposed, it's hard to tell. Oscillations may be a result of a b.c. imposed at the wrong place, in the wrong way or taking the wrong form, resulting in reflections. In my answer I only considered rough initial data and explosion due to instability as this is what you described in the original question. $\endgroup$
    – ekkilop
    Jul 26, 2016 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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