The Right Formula

In the numerical solution of complex, dynamic systems, you often run into what is called the stiffness problem. If there was ever a sentence designed to lose readers, that may be it. Nevertheless, onward.

Imagine a spring connected to a mass. If you pull on the mass, the spring stretches and exerts a force on the mass. The stronger the force caused by a given extension, the “stiffer” the spring.

Now imagine a lot of different size masses, interconnected by a lot of springs of differing stiffness. Some of the mass-spring combinations will react quickly to any change in the system; other combinations, those with a severe mismatch between mass and spring, will react more slowly. Each combination has what is often called a time constant or a characteristic time for its reaction to change.

When you’re doing a numerical solution of the system (the only option if the system is complex and non-linear), you have to solve the state of the system for one point in time, then for another point somewhat later, a “time step” later, and so on. Some of the mass-spring combinations will require shorter time steps than others, because they have different time constants.

If there are connections between mass-spring combinations of very dissimilar time constants, you have a problem, the stiffness problem in fact. If one part of your system has a time constant of a microsecond, while another has a time constant of an hour, you are going to need billions of time steps to calculate your system for each time step you’d need if you didn’t have that short time step piece to it.

There are some fancy solution algorithms that can deal with the stiffness problem. One of them is called “Gear-Hindemarsh,” and was developed at Lawrence Livermore Labs, originally to facilitate the calculations used in designing thermonuclear weapons. We used it for chemical kinetic calculations in simulating smog chamber experiments. Gear-Hindemarsh is the gold standard for that sort of thing, but it has some problems, especially when you give the system a kick, like turning the lights on or off, or otherwise messing with the boundary conditions in an unsmooth way. Then it becomes pretty inefficient. The first atmospheric smog model produced by Livermore was called LIRAQ, and it spent much of its computing time on the few seconds after sunrise and sunset.

The development group I worked with on the EPA Urban Airshed Model used a different approach to the stiffness problem, called the quasi-steady state approximation. The QSSA makes a few reasonable assumptions, such as the idea that a chemical species that you are treating as being at “steady-state” doesn’t have so much mass that it affects the rest of the system. Imagine an automobile with a bunch of bobble-heads inside. The bobbing of the heads doesn’t affect the behavior of the whole system because their mass is small, relative to the auto itself.

If you react a hydrocarbon with an HO radical, for example, the HO pulls an H off of it to make water plus what is called an acyl radical, a hydrocarbon missing the H. The acyl radical then absorbs an oxygen molecule to form a peroxyacyl radical. This takes a very short time to occur, and we don’t worry about the behavior of the system during the time it takes for the radical to absorb the oxygen. We’ve treated the acyl radical as if it were in steady state.

Of course when there are multiple sources and multiple reaction paths for a QSSA species, the algebra can get more complicated, but it’s not too bad. At least not until the QSSA species begin to react with themselves and each other. In the smog equations, the first place this happened was when we included the reaction of the hydroperoxyl radical, HOO, with itself. That yields oxygen plus hydrogen peroxide (and that’s where the peroxide came from to bleach the guy to death in SunSmoke). The QSSA equations for HOO are quadratic. Fortunately, we have an equation to solve quadratic equations, one we all learned in high school. Unfortunately, it’s the wrong equation.

As you’ll recall (said the character in the pulp novel), the solution to the quadratic equation

aX^2 + bX + c = 0

can be written:

X = -b + [or minus] sqrt(b^2 - 4ac) / 2a

What they don’t usually tell you in high school is what happens when this equation is used on a quadratic that sometimes has the value of “a” as zero. If that happens, we get the dread “divide by zero” condition, and your computer tells you that you’ve just done a Very Bad Thing, and refuses to continue, you naughty person.

It so happens that there is another form of the quadratic equation that they don’t tell you about in high school, or in most colleges, either:

X = 2c / (-b + [or minus] sqrt(b^2 - 4ac))

A little checking tells me that the Wikipedia now gives the alternate formula, no doubt because so many programmers have run into the same problem I did, ‘way back when. I forget exactly where I got the alternate quadratic formula from; it’s penciled into the margins of a handbook I have. Anyway, I used it when I coded the chemistry module for the UAM.

Later, we wound up with more radical-radical cross reactions, and the algebraic QSSA went from quadratic to fifth order. There is no general fifth order solution, so we used a numerical solution called “Newton-Raphson.” I didn’t code the first implementation we did of that, and the program kept blowing up in the QSSA solver. I looked at it and realized that the programmer had used a constant term as the initial value for the Newton-Raphson calculation, and N-R is notoriously sensitive to the initial value. For the best results, you need to start somewhere near the final value. The clever lad that was I realized that if I stripped out all but the HOO quadratic, it was going to be very close to the final value. Using that as the initial value, the N-R calculation usually converged in one or two iterations.

* * *

In 1984, I got very sick. The words “chronic fatigue syndrome” screw up your ability to get health insurance, so I never say that I had CFS on an insurance form, and besides, I was never diagnosed. Nevertheless, I had what was basically a bout of ‘flu that lasted for several years. I was unable to work full time; in 1985, working as a consultant, I averaged maybe 5-10 hours a week working.

I was no longer the go-to guy for working on the kinetics solver, and one person took our module for use in an acid deposition model, Mary, a PhD chemist recently graduated from Cal Tech. She took one look at my “quadratic formula,” saw that it did not conform to what she’d learned in school and replaced it with the “right” version. Of course, it promptly blew up. So she spent the next several weeks putting in all sorts of tests for when the “a” in the formula got too small, switching it over to the linear solution etc. I’m not implying that it took her a lot of effort; she just spent some number of hours over the next few weeks working the bugs out.

When I heard about it, I was, of course, annoyed. It’s one thing to have someone else catch your mistake; it’s quite another when it wasn’t a mistake in the first place.

More recently though, I’ve been working as a technical writer, and I’ve come to understand that I did, in fact, make a mistake. I did not document the tricks I used in the chemical kinetics solver, even the most basic documentation, which is to put comments into the code explaining what was done and why.

I never confronted Mary about the thing in the first place, and she unexpectedly died of a cerebral aneurysm many years ago, so there’s no closure in the cards, unless this essay counts.

This Years Model IV

To recap a few things:

I attended Rensselear Polytechnic Institute from 1968 to 1974, graduating with a Master's degree in Engineering Science. Engineering Science was at that time and remains to this day a funny program, a "roll your own" sort of thing. Engineering Science degrees at RPI require that you convince the curriculum chairman that the course of study that you had designed was an appropriate course of study for an engineer. Then, of course, you actually have to complete the program, which is not as easy as you might think (and certainly not as easy as you thought when you first thought up the idea).

One buddy of mine studied urban planning, transportation, and architecture, so now he designs airports around the country and the world. Another took courses in electrical and biomedical engineering, and he's now a hospital management consultant. All in all, it worked out well, I think, at least sufficiently well that the Engineering Science program at RPI continues to this day.

As for me, my course of study centered on the simulation modeling of urban and environmental systems. These days, it's hard for me to even write that sentence without marveling at youthful hubris. What were we thinking? Well, grand notions were in the air. The Cybernetics movement of the 40s and 50s has flowed into what was called General Systems Theory: the idea that science and engineering has developed a set of analytical tools that were so powerful that they might even be able to handle the social and biological sciences.

I put together a course of study that included linear systems and control theory, voice and image processing, urban analysis, with a solid chunk of statistics and operations research for trying to get the data that most people agreed would be necessary to validate and calibrate these huge models that we were going to build. I'm not sure how coherent the course of study, but I will say that for a while it seemed that every course eventually would up with us trying to invert some damn matrix or another.

Finally, I did my graduate work on rewriting a simple simulation model to compare to the very large model that the Lake George Ecosystem project at RPI had prepared. Then I graduated, moved to California, and began to look for a job. Eventually a dream job fell into my lap (literally, as one of the guys I was living with at the time tossed the phone number into my lap, saying, "We decided I wasn't the right guy for this, but it sounds right up your alley), and I became a smog scientist. Why this was a perfect next step will now require some simple math.

Conceptually, the most general model of a dynamic (changes with time) system is the state variable formulation. It pretty much goes, “Here is a series of variables that describe some phenomenon. Each variable is linked to the other variables such that the state of that variable at an instant in time is a function of the other variables at some previous time.” Often, the “previous time” is the instant immediately preceding (for a differential equation), or a discrete time step back (for a difference equation). Sometimes, however, the functional relationship looks back some period of time, although this can be turned into an instant/single time step formulation just by creating more state variables which then contain a time lag link to previous variables, i. e. “memory variables.”

Having said all that, Keep It Simple, Stupid is a good rule to live by, and it’s a pretty good rule in science and engineering. So let me write a simple equation:

dC/dt = kAB

In case anyone here has math nausea, let me emphasize how simple this is. It just says that the rate that C is changing with time depends on the product of A and B with k as a rate parameter (just multiply A and B and k). A, B, and C represent state variables, while k is a parameter, and t is time.

C can be anything, but it is most interesting when C has an effect on A and/or B. For example, suppose

C = -A-B+D+E or

This is like a chemical reaction, where A and B react to form D and E. Another way of writing it is

A + B => C + D

That’s your basic chemical shorthand.

Or suppose you are dealing with a predator/prey relationship:

Wolf + Deer => (1+∆)Wolf (i. e., a well fed wolf)

Or an aquatic ecosystem:

Phytoplankton + Zooplankton => (1+∆)Zooplankton

Phytoplankton + light +phosphate => (1+∆)Phytoplankton –phosphate

Now why is this so interesting?

The most interesting thing about such equations is how generally applicable they are. As I’ve just shown, you can use them for chemical compounds or ecosystems with equal abandon. Why? Because the setup is similar; the rate of change of the state variables depends upon how often the individuals (molecules, plankton, wolves) come in contact with each other. That is a very common situation.

The equations are non-linear but they can come close to being linear when either A or B is much larger than the other. Sometimes, this is called "well-behaved" which means "I think I sometimes understand how it works." Nevertheless, large, "well-behaved" systems often are "counter-intuitive," which means, "Okay, so I was wrong at first, but this time I'm sure I'm right, maybe."

This is just the chemistry part, of course, and the rest of photochemical modeling also has a lot of physics and mechanics in it. Whitten and I used to joke about how academic lectures of smog modeling would usually begin with somebody writing the diffusion equation on the board, a really daunting looking three-dimensional partial differential equation describing fluid flow, followed by a couple of single letters that represented “emissions” and “chemistry.” The joke was that there were well-established ways of solving the diffusion equation numerically, and the process itself (fluid mechanics) has been pretty well understood for generations. In other words, all the really hard work, preparing the emissions inventory and developing the chemistry, were compressed into two humble little letters.

A while back, in “The Right Formula” (May 7), I wrote a bit on “the stiffness problem” that often occurs when you’re trying to solve dynamic state equations that have widely varying time scales. In that essay, I briefly noted the “Gear-Hindemarsh” routines that we used in simulating smog chamber chemistry, before we hard coded the chemical kinetic mechanisms into urban smog models. When doing the latter, we used various tricks that can only be used if you already know the chemistry you’re dealing with. Obviously this isn’t very good when you’re developing the chemistry, but that’s okay, because we had the Gear routines.

The Gear-Hindemarsh codes were developed at Lawrence Livermore Labs, in order to deal with equations like this:

Li6 + n => He4 + H3

H3 + H2 => He4 + n + 17.2Mev

H2 + H2 => He3 + n

H2 + H2 => H3 + H1

In case you didn't notice, the 17.2 Mev means that if you do this to a substantial amount of Li6 and H2 (deuterium), you get a lot of energy. If the pressure and density of the material is right, you get a very large bang, i.e. a thermonuclear detonation. Hence, Livermore's interest.

Gear-Hindemarsh and similar schemes solve stiffness problem but at a cost: every time the systems see an input with discontinuities (step discontinuities even at the 4th or 5th derivative), the solver drops to a lower order predictor corrector and takes very small steps. This is fine for smog chamber experiments, not so good for urban simulations, where inputs and boundary conditions keep changing by the hour.

As the final bit of something a little like irony, I’ll note that I once revisited my Master’s Degree work and used a Gear-Hindemarsh solver on it instead of DYNAMO, which was a simulation language developed at MIT and used by Jay Forrester for his Urban Dynamics and World Dynamics models. The Gear-Hindemarsh results were substantially different from the DYNAMO results, indicating that the (pretty crude) numerical solver in DYNAMO was still sensitive to step size in my simulations. Oops. I have no idea if Forrester’s results suffered from the same problem, though I don’t actually think it matters that much. Forrester’s work had substantially worse problems than a bad number cruncher.