This paper is fundamental to much of CFD, and that is an understatement.
The paper is acknowledged to have made two essential contributions to the field: a critical theorem, and a basic method. It also has another nice method for shock dissipation that is not widely seen. I will discuss all three and try to put things into perspective.
As I recently said at a conference, “Read the classics, don’t just cite them.” In today’s world of too many research papers, finding a gem is rare. The older mines are full diamonds, and the Lax-Wendroff paper is a massive one. By reading the archival literature you can be sure to read something good, and often the older papers have results in them that are not fully appreciated.
Why would I say this? Textbooks will tell you what is important in a paper, and you either believe what you are told, or just read the relevant portions of the paper. Sometimes developments made later could show the older papers in a new light. I believe this is the case with Lax-Wendroff. By read, I mean really read, cover-to-cover.
1. The theorem has enormous implications on how methods are designed because of the benefits of using discrete conservation form. In this way the numerical method mimics nature and solutions have an intrinsically better change of getting the right answer.
What does the theorem say?
Basically, if a numerical method for a first-order hyperbolic partial differential equation is in discrete conservation form, it must converge to a weak solution (assuming first that its is stable and consistent). Stable and consistent means the solution is well-behaved numerically (doesn’t blow up), and approximated the differential equation.
If one follows the guidance, it is unabashedly good, that is, we want to get valid weak solutions, but there is a catch. The catch is that there are literally infinitely many weak solutions, and most of them are not desirable. What we actually desire is the correct, or physically valid weak solution. This means we want a solution that is an “entropy” solution, which is essentially the limiting solution to the hyperbolic PDE with vanishing dissipation. This leads to numerical dissipation (i.e., things like artificial viscosity) to create the entropy needed to select physically relevant solutions.
The simpliest hyperbolic PDE where this issue comes to a head is Burgers’ equation with an associated upwind discretization assuming the velocity is positive,
which can be rewritten equivalently as
Equations (1) and (2) are identical if and only if the solution is continuously differentiable (and in the limit where error is small), which is not true at a discontinuity. What the conservation form means effectively that one is sure the amount of stuff (u) that leaves on cell on the mesh exactly enters the next cell. This allows all the fluxes on the interior of the mesh to telescopically cancel leaving only the boundary terms, and provides a direct approximation to the weak form of the PDE. Thus the numerical approximation will exactly mimic the dynamics of a weak solution. For this reason we want to solve (1) and not (2) because it is in conservation form. With this theorem, the use of conservation form went from mere suggestion, to an iron clad recommendation.
For a general conservation law the conservation form gives,
which generalizes to higher order methods depending on how the cell-edge fluxes are defined. The last bit of detail is picking the right weak solution since there are an infinite number of them. This depends on numerical dissipation, which mimics the dynamics of vanishingly small physical viscosity. The presence of viscosity produces entropy as demanded by the second law of thermodynamics. The inclusion of sufficient dissipation will choose a unique physically relevant solution to the conservation law.
2. The method introduced in this paper was another centerpiece and icon of the numerical solution of hyperbolic PDEs. This is a stable, second-order, centered method, which was developed because the stable centered method introduced earlier was first-order accurate and quite dissipative (this is the Lax-Friedrichs method from 1954). The Lax-Wendroff method is derived using the Taylor series where the time derivatives are systematically replaced with space derivatives.
For a simple PDE the procedure is quite straightforward,
replace the time derivative with space derivatives,
and discretize with centered differences. This method is second-order accurate in space and time, but will create oscillations near discontinuities. This will provide a motivation for the third result from the paper, a shock dissipation, which is usually overlooked, but needs a new appreciation.
Several years later Richtmyer presented this as a two-step, predictor-corrector method that made for simple and efficient implementation (NCAR Technical Note 63-2 ). The Lax-Wendroff method then rapidly became one of the methods of choice for explicit in time aerodynamic calculations. Perhaps Bob MacCormack at NASA Ames/Stanford developed the most famous of these methods in 1969. This was a variant of Richtmyer’s approach utilizing a sequence of forward, followed by backwards-biased differences.
All of this set the stage for much bigger things in the mid 1970’s, the introduction of genuinely nonlinear schemes where the finite difference stencil depended on the structure of the solution itself.
The procedure of replacing time- with space-derivatives has reappeared in the past few years often as “Lax-Wendroff” time-differencing. It can be very complex for systems of conservation laws, or for higher than second-order accurate differencing of genuinely nonlinear equations. This is a consequence of the explosion of terms arising from the systematic application of the chain rule. Nonetheless it can be an effective alternative to a method-of-lines approach yielding more satisfying CFL stability limits.
3. The dissipation introduced by Lax-Wendroff was needed because of the oscillatory solutions (Godunov’s theorem explains why, something to discuss in detail later). Because of this they introduced a shock dissipation that looked much different at first blush than the Von Neumann-Richtmyer mechanism. We will show briefly that it is basically identical in an approximate mathematical sense. Hopefully this connection will allow this result to be more widely appreciated.
In the paper they introduce a shock dissipation based on the variation of the sound speed (the product of density and speed of sound, “c”, I removed the constant and the mesh spacing from the formula too, from pp. 232 of the paper),
which I will show through the application of the chain rule as being equivalent to,
where the second transformation has used the differential equations themselves (or the Rankine-Hugoniot conditions). This is just like the Von Neumann-Richtmyer approach if we take,
This then implies that the viscosity coefficient actually arises from the equation of state because this quantity is known as the fundamental derivative,
the only proviso is that the fundamental derivative is computed at constant entropy.
Bethe, Zeldovich and Thompson each studied the fundamental derivative, but Menikoff and Plohr’s paper is the best available discussion of its impact on shock wave dynamics (Reviews in Modern Physics, Volume 61, 1989). This is defined by the thermodynamics of the equation of state as opposed to the mindset that the coefficient of the Von Neumann-Richtmyer viscosity is arbitrary. It is not! This coefficient is determined by the state of the material being simulated. Thus the quadratic viscosity can be reinterpreted as
We hope that this observation will allow artificial viscosity to be selected without appealing to arbitrarily selected constants, but rather for the specific material being simulated.
Wendroff still works at Los Alamos (every day as Misha Shashkov notes). The Los Alamos report number is LA-2285 from 1958/1959.
I asked Burt Wendroff whether anyone working on the Lab’s non-conservative codes based on the staggered mesh, artificial viscosity approach of Von Neumann and Richtmyer had inquired about the necessity of conservation form. The answer was no, never. So the theorem that shaped much of CFD outside the Lab had no impact inside the Lab despite being developed there. Los Alamos was literally the birthplace of CFD. Its first baby was artificial viscosity and the staggered mesh methods, and its second baby was driven by Lax and includes Lax-Wendroff. The Lab orphaned the second baby, but found good parents elsewhere where it flourished. Its would be akin to the first child taking over the family business and living a comfortable, but somewhat unfulfilled life, while the second child left home and found greatness.
Why would this be? I believe that after Von Neumann and Richtmyer showed the way to make shock capturing work, the Lab turned its energy almost completely toward the development of the H-Bomb, and the development of improved methods that Lax showed the path toward was not paid attention to. Stunningly this pattern continued for decades. In the past 10-15 years, the work of Lax has started to have an impact there albeit 50 years on from when it should have.
I will note that the term finite volume scheme didn’t come into use until 1973 (in a paper by Rizzi), but the Lax-Wendroff method pushed the field in this direction. Lax used finite volume ideas immediately as published in a Los Alamos report from 1952 (LA-1205). He called it a finite difference method, as was the fashion, but it is clearly a finite volume method, and the conservation form was immediately appreciated by Lax as being special. It is notable that Von Neumann’s method is not in conservation form, and does not conserve in any simple way.
The authors are famous enough to have biographies in Wikipedia.