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

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?

-
The so-called Method of Manufactured solutions makes exact solutions available for all problems. It may not be able to generate the sort of problematic solutions you describe, but it's not the case that exact solutions are never available. –  Bill Barth Dec 2 '11 at 13:08
I think this is difficult here since you would need to guess a solution with the kind of discontinuity which is not well-approximated by the solution method. –  Matt Knepley Dec 2 '11 at 13:53
I agree that it's likely difficult to manufacture solutions that excite the problematic modes that Jed mentions. I just wanted to point out that exact solutions are always available for testing. I don't know what happens if you manufacture a solution to the 1D linearized shallow water equations using, say, a mix of trig and exponential functions (typical of MoM exact solutions), turn the crank to get the corresponding source terms, and run them through a 1st-order Godunov scheme. Maybe Jed can give it a shot and report back. –  Bill Barth Dec 2 '11 at 15:30
MoM is a great tool, but in this case, the issue is that diffusion is mis-applied inside a shock. Everywhere else, diffusion converging to zero on each equation equally is acceptable, but diffusion does not converge to zero inside a shock, so applying numerical diffusion to each term equally results in incorrect dynamics. I will write a long answer to this question when I have time, if nobody beats me to it. –  Jed Brown Dec 2 '11 at 15:39
@Jed, shouldn't LET apply to the linearized equations? –  Matt Knepley Dec 2 '11 at 16:29
show 1 more comment

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.

-
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.