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.