Science is the process that takes us from confusion to understanding…

― Brian Greene

stability-in-lifeStability is essential for computation to succeed. Better stability principles can pave the way for greater computational success. We are in dire need of new, expanded concepts for stability that provides paths forward toward uncharted vistas of simulation.

Without stability numerical methods are completely useless. Even a modest amount of instability can completely undermine and destroy the best simulation intentions. Stability became a thing; right after computers became a thing. Early work on ordinary differential equations encountered instability, but the computations being handcrafted was always suspect. The availability of automatic computations via computers ended the speculation, and now it became clear, numerical methods could become unstable. With the proof of a clear issue in hand great minds went to work to put this potential chaos to order. This is the kind of great work we should be asking applied math to be doing today, and sadly are not because of our over reliance on raw computing power.

Von Neumann devised the first technique for stability analysis after encountering ijohn-von-neumann-2nstability at Los Alamos during World War 2. This method is still the gold standard for analysis today in spite of rather profound limitations and applicability. In the early 1950’s Lax came up withrichtmyer_robert_b1the equivalence theorem (interestingly both Von Neumann and Lax worked with Robert Richtmyer,, which only highlighted the importance of stability more boldly. Remarkably ordinary differential equation methods came to stability later than partial differential equations in Dahlquist’s groundbreaking work. He produced a stability theory and equivalence theory that paralleled the work of Von Neumann and Lax for PDEs. All he needed were computers to drive the need for the work. We will note that the PDE theory is all for linear methods and linear equations, while the ODE theory is for
linear methods, but applies to nonlinear ODEs,

Once a theory was established for stability computations could proceed with enough guarantee of solution to progress. For a very long time this stability work was all that was needed. Numerical methods, algorithms and general techniques galore came into being and application covering a broad swath of the physics and engineering World. Gradually, over time, we started to see computation become spoken as a new complementary practice in science that might stand shoulder to shoulder with theory and experiment. These views are a bit on the grandiose side of things where a more balanced perspective might rationally note that numerical methods allow complex nonlinear models to be solved where classical analytical approaches are quite limited. At this point its wise to confront the issue that might be creeping into your thinking, our theory is mostly linear while the utility for computation is almost all nonlinear. We have a massive gap between theory and utility with virtually no emphasis or focus or effort to close it.

This is so super important that I’ve written about it before, doing the basic Von Neumannstability-1methodology using Mathematica, &, in the guise of thoughts about robustness,,

and practical considerations for hyperbolic PDEs, Running headlong through this arc of thought are lessons learned from hyperbolic PDEs.Peter_LaxHyperbolic PDEs have always been at the leading edge of computation because they are important to applications, difficult and this has attracted a lot of real unambiguous genius to solve it. I’ve mentioned a cadre of genius who blazed the trails 60 to 70 years ago (Von Neumann, Lax, Richtmyer, We are in dire need of new geniuses to slay the nonlinear dragons that stand in the way of progress. Unfortunately there is little or no appetite or desire for progress, and the general environment stands squarely in opposition. The status quo is viewed as all we need,, and progress in improving basic capabilities and functionality has disappeared except for utilizing ever more complex and complicated hardware (with ever more vanishing practical returns).

Antifragility is beyond resilience or robustness. The resilient resists shocks and stays the same; the antifragile gets better.

― Nassim Nicholas Taleb

I’ve been thinking about concepts of nonlinear stability for a while wondering if we can move past the simple and time-honored concepts developed so long ago. Recently I’ve taken a look a Taleb’s “anti-fragile” concept and realized that it might have some traction in this arena. In a sense the nonlinear stability concepts developed here to fore are akin to anti-fragile where the methodology is developed and works best in the worst-case scenario. In the case of nonlinear stability for hyperbolic PDEs the worst-case scenario is a linear discontinuity where the solution has a jump and the linear solution is utterly unforgiving of any errors. In this crucible, all the bad things that can happen to a solution arise, either overly diffusive low accuracy solutions, or oscillations with high accuracy producing demonstrably unphysical results.

No structure, even an artificial one, enjoys the process of entropy. It is the ultimate fate of everything, and everything resists it.

― Philip K. Dick

When the discontinuity is associated with a real physical solution these twin and opposing maladies are both unacceptable. The diffusion leads to a complete was of computational effort that dramatically inhibits any practical utility of numerical methods. For more complex systems of equations where turbulent chaotic solutions would naturally arise, the diffusive methods drive all solutions to be laminar and boring (not matching physical reality in essential aspects!). On the plus side of the ledger, the diffusive solution is epically robust and reliable, a monument to stability. On the other high-order methods based on the premise that the solution is smooth and differentiable (i.e., nice and utterly ideal) is the epitome of fragility. The oscillations can easily put the solution into unphysical states that render the solution physically absurd.

Difficulty is what wakes up the genius

― Nassim Nicholas Taleb

Now we get to the absolute genius of nonlinear stability that arose from this challenge. Rather than forcing us to accept one or the other, we introduce a concept by which we can have the best of both, using whatever discretization is most appropriate for the local circumstances. Thus we have a solution adaptive method that chooses the right approximation for the solution locally. Therefore a different method may be used for every place in space and time. The key concept is the rejection of the use of a linear discretization where the same method is applied everywhere in the solution, which caused the entire problem elucidated above. Instead we introduce a mechanism to analyze the solution and introduce an approximation appropriate for the local structure of the solution.

Tmediocritydemotivatorhe desired outcome is to use the high-order solution as much as possible, but without inducing the dangerous oscillations. The key is to build upon the foundation of the very stable, but low accuracy, dissipative method. The theory that can be utilized makes the dissipative structure of the solution a nonlinear relationship. This produces a test of the local structure of the solution, which tells us when it is safe to be high-order, and when the solution is so discontinuous that the low order solution must be used. The result is a solution that is high-order as much as possible, and inherits the stability of the low-order solution gaining purchase on its essential properties (asymptotic dissipation and entropy-principles). These methods are so stable and powerful that one might utilize a completely unstable method as one of the options with very little negative consequence. This class of methods revolutionized computational fluid dynamics, and allowed the relative confidence in the use of methods to solve practical problems.

Instead of building on the lessons learned with these methods, we seem to have entered an era where the belief that no further progress on this front is needed. We have lost sight of the benefits of looking to produce better methods. Better methods like those that are nonlinearly stable open new vistas of simulation that presently are closed to systematic and confident exploration. The methods described above did just this and produced the current (over) confidence in CFD codes.

A good question is why haven’t these sort of ideas spread to a wider swath of the numerical world? The core answer is “good enough” thinking, which is far too pervasive today. Part of the issue is the immediacy of need for hyperbolic PDEs isn’t there in other areas. Take time integration methods where the consensus view is that ODE integration is good enough, thank you very much, and we don’t need effort. To be honest, it’s the same idea as CFD codes, which are they are good enough and we don’t need effort. Other areas of simulation like parabolic and elliptic PDEs might also use such methods, but again thestability-3.hiresneed is far less obvious than for hyperbolic PDEs. The truth is that we have made rather stunning progress in both areas and the breakthroughs have put forth the illusion that the methods today are good enough. We need to recognize that this is awesome for those who developed the status quo, but a very bad thing if there are other breakthroughs ripe for the taking. In my view we are at such a point and missing the opportunity to make the “good enough,” “great” or even “awesome”. Nonlinear stability is deeply associated with adaptivity and ultimately more optimal and appropriate approximations for the problem at hand.

If you’re any good at all, you know you can be better.
― Lindsay Buckingham

So what might a more general principle look like as applied to ODE integration? Let’s explore this along with some ideas of how to extend the analysis to support such paths. The acknowledged bullet-proof method is backwards Euler, u^{n+1} = u^n + h f(u^{n+1}), which is A-stable and the ODE gold standard of robustness. It is also first-order accurate. One might like to use something higher order with equal reliability, but alas this is exactly the issue we have with hyperbolic PDEs. What might a nonlinear approximation look like?

Let’s assume we will stick with the BDF (backwards differentiation formula) methods, and we have an ODE that should produce positive definite (or some sort of sign- or property-preserving character). We will stick with positivity for simplicity’s sake. The fact is that perfectly linearly stable solutions may well produce solutions that lose positivity. Staying with simplicity of exposition, the second-order BDF method is \frac{3}{2} u^{n+1} = 2 u^n - \frac{1}{2} u^{n-1} + h f(u^{n+1}). This can be usefully rearranged to \frac{3}{2} u^{n+1} - h f(u^{n+1}) = 2 u^n -\frac{1}{2} u^{n-1}. If the right hand side of this expression is positive, 2 u^n -\frac{1}{2} u^{n-1} > 0, and the eigenvalue of f(u^{n+1}) < 0, we have faith that u^{n+1} > 0. If the right hand side is negative, it would be wise to switch to the backwards Euler scheme for this time step. We could easily envision taking this approach to higher and higher order.

StabilityFurther extensions of nonlinear stability would be useful for parabolic PDEs. Generally parabolic equations are fantastically forgiving so doing anything more complicated is not prized unless it produces better accuracy at the same time. Accuracy is imminently achievable because parabolic equations generate smooth solutions.  Nonetheless these accurate solutions can still produce unphysical effects that violate other principles. Positivity is rarely threatened although this would be a reasonable property to demand. It is more likely that the solutions will violate some sort of entropy inequality in a mild manner. Instead of producing something demonstrably unphysical, the solution would simply not have enough entropy generated to be physical. As such we can see solutions approaching the right solution, but in a sense from the wrong direction that threatens to produce non-physically admissible solutions. One potential way to think about is might be in an application of heat conductions. One can examine whether or not the flow of heat matches the proper direction of heat flow locally, and if a high-order approximation does not either choose another high-order approximation that does, or limit to a lower order method wit unambiguous satisfaction of the proper direction. The intrinsically unfatal impact of these flaws mean they are not really addressed.

Mediocrity will never do. You are capable of something better.
― Gordon B. Hinckley

Thinking about marvelous magical median ( spurred these thoughts to mind. I’ve been taken with this functions ability to produce accurate approximations from approximations that are lower order, and pondering whether this property can extend to stability either linear or nonlinear (empirical evidence is strongly in favor). Could the same functionality be applied to other schemes or approximation methods? In the case of the median if we take two of the three arguments to have a particular property, the third argument will inherit this property if the other two bound it. This is immensely powerful and you can be off to the races by simply having two schemes that possess a given desirable property. Using this approach more broadly than it has been applied thus far would be an interesting avenue to explore.

If you don't have time to do it right, when will you have the time to do it over?
― John Wooden