How to determine if a numerical solution to a PDE is converging to a continuum solution?

Question

The Lax equivalence theorem states that consistency and stability of a numerical scheme for a linear initial value problem is a necessary and sufficient condition for convergence. But for nonlinear problems, numerical methods can converge very plausibly to incorrect results, despite being consistent and stable. For example, this paper shows how a first order Godunov method applied to the 1D linearized shallow water equations converges to an incorrect solution.

Evidently self-convergence under mesh and time step refinement is not sufficient, but exact solutions are generally not available for nonlinear PDEs, so how can one determine if a numerical method is converging to a genuine solution?

Answer 1

There are two main classes of solutions to be discussed in this regard.

"Sufficiently" Smooth Solutions

In Strang's classical paper it is shown that the Lax equivalence theorem (i.e., the idea that consistency plus stability implies convergence) extends to nonlinear PDE solutions if they have a certain number of continuous derivatives. Note that that paper is focused on hyperbolic problems, but the result carries over to parabolic problems. The number of derivatives needed is a technical point, but this approach is usually applicable to solutions that satisfy the PDE in a strong sense.

Discontinuous solutions

At the other extreme, we have PDE "solutions" with discontinuities, which typically arise from nonlinear hyperbolic conservation laws. In this situation, of course, the solution cannot be said to satisfy the PDE in the strong sense, as it is non-differentiable at one or more points. Instead, a notion of weak solution must be introduced, which essentially amounts to requiring that the solution satisfy an integral conservation law.

Proving convergence of a sequence of solutions is also more difficult in this case, as $L_p$-stability is not enough; usually the sequence must be shown to lie in a compact space, such as the set of $L_\infty$ functions with some finite maximum total variation.

If the sequence can be shown to converge to something, and if the method is conservative, then the Lax-Wendroff theorem guarantees that it will converge to a weak solution of the conservation law. However, such solutions are not unique. Determining which weak solution is "correct" requires information that is not contained in the hyperbolic PDE. Generally, hyperbolic PDEs are obtained by neglecting parabolic terms in a continuum model, and the correct weak solution can depend on exactly what parabolic terms were discarded (this last point is the focus of the paper linked to in the question above).

This is a rich and involved topic, and the mathematical theory is far from complete. Most convergence proofs are for 1D problems and rely on specialized techniques. Thus nearly all of the actual computational solutions of hyperbolic conservation laws in practice cannot be proved convergent with existing tools. For a practical discussion from a computational standpoint, see LeVeque's book (chapters 8, 12, and 15); for a more rigorous and detailed treatment I'd suggest Dafermos.

Answer 2

I have little to contribute here other than to point out that whenever numerical methods have trouble with hyperbolic equations (and converge to the wrong solution), it isn't usually because of shocks. Rather, the areas they have difficulty with it are rarefaction waves -- where the solution is smooth.

Another example of something that appears to be difficult is for a simple equation of the kind $$ u_t + \beta \cdot \nabla F(u) = g $$ you get into trouble in those places where $F'(u)=0$. This is called the "sonic point" of an equation. (Note that the equation can be rewritten as $u_t + \beta F'(u) \cdot \nabla u = g$, which explains why $F'=0$ is a problem). What happens here is that the equation degenerates: it's just an ODE wherever $F'=0$ but the solution is of course affected by surrounding areas where this is not the case. Consequently, if you happen to have $F'=0$ in a whole part of the domain at the initial time, say $\omega\subset\Omega$ where $|\omega|>0$, then this part of the domain may grow or shrink. Because the solution is typically smooth in these areas (in fact, it is typically constant), numerical methods don't add artificial diffusion in these areas and consequently do not converge to a viscosity solution.

As a final note: An example for a flux $F(u)$ from multiphase flow in porous media is $$ F(u) = \frac{u^4}{u^4 + (1-u)^2 \cdot (1-u^2)} $$ where $u$ is the saturation. You have $F'(u)=0$ for $u=0$, i.e., in areas where the saturation is zero.