Taking the last lap around the Rosetta Stone - Embedded.com

Taking the last lap around the Rosetta Stone

We're on the last lap of our study of the equation I call the Rosetta Stone:


With this equation, we can relate events occurring in the world of continuous time, in which we live, to the world of discrete time, in which the digital computer lives.

Together, we've walked a road that began with simple concepts like extrapolating a simple function, moved on through the derivation and understanding of the Taylor series, past Heaviside's notion of operational calculus, to finally arrive at Equation 1. Recall that this equation is not a relationship between variables but between two operators, both operating on things that change with time.

Let's be clear on one thing: it's not my intention to give you formulas that you can use for every possible situation. Rather, I've tried to explain where Equation 1 comes from to such a depth and clarity that you'll be able to manipulate it as needed to solve your own problem. Last time, I picked the easiest problem I could think of, which was to find the derivative of a table of sampled values.

In this current and final installment, I'm going to show you how to solve what's probably the hardest problem I can think of, which is to numerically integrate a function.

If I'm going to do this, I need to be as clear as I can possibly be to make sure you have no doubts how Equation 1 works. Just so there are no misunderstandings or doubts from the get-go, allow me to set the stage. This time, I'll change the notation a bit to make the relationships even more clear than in previous episodes.

Suppose we have some measurable parameter–pressure, perhaps, or temperature, or the RPMs of a shaft. If that function is changing with time, we usually assign it a name something like F (t ). We can also draw a graph of it, as in Figure 1.

Now, this point is critically important:

The digital computer can never know the actual behavior of F (t ).

That's because the digital computer is . . . um . . . digital, so can only do things at discrete time steps. If it's going to measure F (t ) at all, it'll be doing it with an analog-to-digital converter (A/D). Between measurements, the computer has other stuff to do. So the best we can hope for is that we sample F (t ) at discrete times. Equation 1 only works if those samples are taken with a fixed time, h , between samples (some folks prefer t ).

Conceptually, each measurement occurs at a time, and has a value. Together the two values represent a pair, {t , F (t )}. Subsequent measurements create an ordered set of such pairs: {t 0 , F (t 0 )}, {t 1 , F (t 1 )}, {t 2 , F (t 2 )}, {t 3 , F (t 3 )}, and so on. People sometimes actually store data in this fashion. It's called time-stamped data , and we use this approach when the data are sampled at unequal intervals. But we're only interested here in cases where the time interval between samples is constant, so there's no reason to store the t 's; we can always calculate them from h :


We are left, then, with a table–a list–of measurements: F (t 0 ), F (t 1 ), F (t 2 ), and so on. To keep the notation simple, let those measurements be denoted by f 0 , f 1 , f 2 , and so forth. These values are also shown in Figure 1. fn is also a function, just like F (t ), but it's a function of n , not t , and it only exists at the specified time points. In short, F (t ) is a continuous function in continuous time; fn is a discrete function in discrete time. The two functions both represent the same physical parameter–pressure, temperature, or whatever–but they do so in very different ways. Our challenge is to relate the two; in other words, take either function and use it to figure out what the other ought to be doing. That's where the Rosetta Stone comes in.

Remember, the symbols in Equation 1 are operators , which mean that they must operate on something. The operator D represents the time derivative; that is:


The operator z , on the other hand, operates on the tabular values, by advancing one step further into the table.


If z advances one step in the table, its inverse, z -1 , should advance one step backwards, and it does:


While the Rosetta Stone, Equation 1, is a short and very cool way to remember the relationships, it makes more sense if we include the things the operators are operating on. Therefore, a more meaningful form might be:


At first glance, you may be thinking that I should have put the F (t ) as part of the exponent on the right, but that would be wrong. Remember, to get anything useful we must expand the exponential into its power series form, to get:


The powers of D we take to mean higher- and higher-order derivatives of F (t ). Equation 6 is, in fact, a compact way of writing the Taylor series, and Equation 1 is far more compact yet.

In his notion of operational calculus, Oliver Heaviside made the leap of faith that, since Equation 1 looks like a simple relationship between variables, we can treat it as such. Which means we can, for example, revert it to get:


Following the explicit symbolism of Equation 3, we should include the functions the operators operate on, to write:


As before, we don't get anything usable unless we expand the function ln z into a power series. When we do, we get a value for the derivative DF (t ), by applying higher and higher powers of z .

This is the problem we addressed in the last episode. It turns out that there is no convenient series form for ln z , but there is one in the backward difference:


Using this relationship, we were able to write:


Before we go forward, it's important that you understand all the nuances of this equation. Given a table, fn , of sampled values of F (t ), it lets us calculate the first derivative of that function, evaluated at the time tn . In practice the computation is going to be approximate (though remarkably good), because in practice we must truncate the power series to some order. But theoretically, the equation is exact . The only error is incurred in the truncation.

You might think that, if we wanted higher-order derivatives of F (t ), we should raise both sides of Equation 8 to higher powers, like this:


You could, in fact, do just this, but that would be doing things very much the hard way. Multiplying power series, while possible, is not fun. Fortunately, it's much easier to manipulate the relationships in z , or rather, in . Consider that:


In other words, don't try to get a single expression operating on fn . Instead, use Equation 11 to get a new table of derivatives, DF (tn ). Then apply it again to get the second derivative:


With these values, we can generate a new table, which is a table of derivatives. Then, applying Equation 11 gives us a new set of values that are the second derivatives.

I'm sure you can see that we can continue this process indefinitely, for higher- and higher-order derivatives.

There is one catch that complicates things a bit. Suppose you have, say, 100 entries in the table for fn . Remember that the backwards difference operator represents the difference of two adjacent values, so it needs those backwards values to operate. Given a table of, say, 100 values of fn , you can only generate 99 values for fn . Similarly, you can only have 98 values for 2 fn , 97 values for 3 fn , and so on. If we're planning to truncate the series in Equation 11 to fourth order, we can only get 96 values of DF (t ). If we repeat the process to get the second derivative, we're going to be down to 92 values.

This may seem like a real problem, but in practice, it's not. If we have that table with 100 entries, we don't need to store and save all 100; we only need to save the ones we need. That'll mean four values back for each power of D . Storing such a small number of entries is not a problem; we can do it easily, using either a shift register or a circular buffer.

This does mean, however, that we can't get the values of the derivatives we need until we've taken enough steps to “prime the pump.” All the tables need to be fully populated before we can calculate that first value of Dfn . Populating the tables, however, is not usually a problem, either. In practice we're going to be firing up a real-time system that's going to run essentially forever. We'll be taking constant measurements of F (t ), and so after the first few time ticks, the tables will be fully populated.

Integrating a function
So much for the fundamentals. Let's now turn our attention to the problem at hand, which is to find an integral of some differential equation.

Again, let's be very precise here. If I have some function F (t ), the integral of the function is:


This operation is called a quadrature . It's characterized by the fact that F (t ) is just that–a function only of time. Give me a time, and I'll give you F (t ) at that time. Quadrature represents the fundamental problem of integral calculus, and is equivalent to finding the area under the curve of F (t ).

For the record, Simpson's rule solves this problem very nicely and gives a result that's accurate to fourth order. But that's not the problem we want to solve here. We want to solve what is potentially a much more challenging problem, which is to solve an ordinary differential equation (ODE). The standard form of the ODE is:


Here I'm overloading the function f ; it's different from the fn that we used earlier. Sorry about that, but the form in Equation 16 is so standardized that I'd be drummed out of the corps for using any other symbol.

As in Equation 15, the solution of Equation 16 represents the integral:


The difference between Equations 15 and 17 is subtle, but critical. Here, the function f is a function of both t and x . This makes all the difference. Equation 15 represents an area under the curve F (t ). Given an interval {t 0 , tf }, I can evaluate that area once and for all time. The area won't change unless we change the limiting values, t 0 and tf . But the integrand of Equation 17 depends on both t and x . That means two important things:

1. We can't get a solution until we know what the initial value x 0 should be, and

2. If we change x 0 , the solution will change as well

To get a solution, then, we have to add a set of initial conditions to Equation 16:


For obvious reasons, this problem is called the initial value problem, and it's crucial to predicting the motion of everything that moves in physics. This is the equation we use for simulating the motion of everything from rockets and missiles to airplanes in flight simulators, to chemical reactions–to you name it. Virtually everything in nature, if it changes with time at all, can be described by an ODE. Therefore, almost any physics problem is going to need an ODE solver. That's the kind of integration we're talking about here.

The next step is critical. You must very carefully and very precisely write down what you know and what you hope to get as a result.

What we know are the two relationships in Equation 18. We have the initial time and x -value, t 0 and x 0 . Because these two values are known, we can also compute f (t 0 , x 0 ). Call it f 0 .

What we hope to get out is a time history of x (t ), meaning a set of ordered pairs {tn , xn } over whatever time interval we're interested in. From that data, we can draw a graph of the motion, find out what x is at the final time tf , and do myriad other calculations. In the process of solving the problem, we also must, by necessity, compute the values of f (t , x ) at each time tick, so we get the history of the derivative, gratis.

Our problem reduces to this one, then: given values for t 0 , x 0 , and f 0 , calculate t 1 . From there, we can evaluate f 1 , shift everything by one step, and repeat the process. The trick is to have a way to compute that first new value, x 1 .

Euler's method
Is there an easy way to do that? Actually, there is. And it's quite easy to both derive and implement. It's based on the notion that, to a first approximation, we can replace the derivative by a divided difference of two small numbers. Write:


where the operator has its usual connotation:


Of course, the derivative on the left in Equation 19 is f (t , x ). Also, t = h , so we can write:


Now we've achieved our goal, which is to compute the next value of x using only values for the point t = tn . This is called Euler's method , and it works–sort of. You may find it useful in a pinch, but it has one serious shortcoming: it's not very accurate. You can make it more accurate by taking smaller steps, but only to a degree. It's easy to show that, to get computer-quality accuracy, you'd have to take millions–even billions–of integration steps. Aside from the computer time this is likely to take, by that time the solution may well be lost in the noise of floating-point roundoff error. To do serious numerical integration, we need a formula that's more capable.

Multistep methods
What we seek is a more accurate method that lets us estimate x 1 , given only t 0 , x0 , and, of course, f 0 = f (t 0 , x 0 )? Does such a method exist? Sadly, no. The information simply isn't there. Sorry about that.

More accurately, methods do exist to move from one data point to the next without using any values except the current ones. For obvious reasons, these are called single-step methods , of which Runge-Kutta is the best-known example. But these methods cheat a bit, since they probe around for off-mesh values of f (t , x ). They're not usable, in general, for real-time systems because our sensors aren't psychic. We can't ask them what they would give us for hypothetical values of t and x . No fair peeking into the future.

Peering into the past, however, is allowed. If the current values of t , x , and f aren't enough to get an accurate integration, what if we use past values as well? This, it turns out, leads to the multistep methods , and they come right from the Rosetta Stone.

Don't ever think, however, that deriving the formula is going to be easy. It's not. It's going to take some serious thought and foresight. There are many ways to manipulate these operators, but only one takes us where we want to go. And, as you'll see, the process can get very tricky. That's why I keep stressing that you have to focus on what you have, what you want to get, and what result you seek. I can lead you through the process for the specific problem at hand, and you'll see what sorts of tricks we can use to get the outcome we desire. But I won't be there to help you with the next problem, so you need to understand the process, not the details.

Let's begin again by listing what we know. We've already seen that looking only at the most recent values of t , x , and f aren't enough. But what if we retain past values as well? Imagine that we've already been computing for awhile. We've built the time history of x , up through the last value. That is, we have the sequence {x 0 , x 1 , x 2 , …, xn }. We also have the implicit sequence of t 's, which come from:


Since we have both t and x for all n so far, we can also compute a sequence of derivatives, {f 0 , f 1 , f 2 , …, f n }. Indeed, we must already have done all but the last one, else we wouldn't have advanced at all. To complete the picture, calculate:


and we have all the information available.

Next, let's ask what we hope to get out? That's easy: we seek the next value, xn +1 in the sequence. From the definition of z , we have:


For the next step, I've wracked my brains trying to think of a foolproof, guaranteed approach to get the form we need, based only on the things we know. I've not been able to find that approach. I can only think of arm waving, trial-by-error approaches. But hey, trial by error isn't always so bad, and maybe by the error, you'll be able to see how important it is to get the forms right.

Since we want to use both current and past values, we can guess that we need to express z in terms of the backward difference operator, . From Equation 10, we can write:


Although it may not seem obvious, this form can be expressed as a power series in . You can prove it either by synthetic division or simply multiplying the series by 1-. Substituting the power series for z in Equation 24, we can write:


which involves only current and past values of x . So we have the solution to our problem, right? Wrong! Equation 26 is the formula for extrapolation, which may be helpful sometimes, but doesn't let us integrate the ODE. In fact, Equation 26 doesn't contain any f -values at all.

Ok, let's back up and try again. From Equation 18, we can write:


What we need is a relation like Equation 26, but with f 's on the right instead of x 's. To do that, we need a formula having one and only one power of D. We can get this using the trick of writing:


If, again, we replace the term z/D by a power series in , we should have the result we seek.

But we don't. This approach requires us to use only current and past values of f (t , x ). We haven't used the most important data item of all, which is xn . By experience and a little common sense, I know that I'm going to want a formula that's an extension of Euler's method, with the form:


where P () is some polynomial–or rather, a truncated power series–in . Since we need that one value of x , let's move it to the left-hand side for safety. Write:


But, from Equation 25:


Now, if we pull that same D/D trick, we can write:


This, finally, is the form we need. That critical last value of x is safely on the left, hidden in the one forward difference. All the tabular values on the right are values of the derivative function, f . To finish, we need only express D in terms of , and manipulate the form of the result. According to the Rosetta Stone:


The form we get, then, is:

or, “dividing out” the xn :


From here on, it's all algebra. Since we have all the pieces in the right places, we only need to look at the expression in and manipulate it so as to get a power series. Granted, it's tricky algebra, but still just algebra. We've already seen the series for ln(1-). It's:


Multiplying this by 1- gives:


Collecting terms gives a new power series:


We're not done yet, though. Remember, this series goes in the denominator of Equation 35. So we must divide it into – using synthetic division. The process is do-able, but not easy; I'll just show the first few steps here. Write:


Here, I've just divided the first term of the denominator, just as you'd do with ordinary long division. And just as in long division, we must next multiply the first term of the quotient–here just the 1–by the denominator, and subtract.


Divide the first term on the left again, to get:


Now repeat the process. A simple matter of algebra.

At this point, I must confide that the first time I did this algebra, I did it by hand and could only get a few terms. Later, I used some 14-inch, fanfold computer printout paper, and got the expansion to order 25. It took a long, long time.

These days, we have math aids like Mathcad, which will do the symbolic algebra for us. That makes life a lot easier, and for the record, I'm now recommending Mathcad to my readers again.

After the synthetic division is complete (with some horrible algebra), we end up with:


Putting this back into our equation for xn +1 , we can write:


This is called the Adams predictor formula , and it's very good. Granted, the synthetic division gives us some pretty horrible and unpredictable coefficients, but we have what we were seeking: a higher-order formula for computing xn +1 .

Let me remind you that, once we have decided where to truncate the power series, we can always express in terms of z -1 , to get formulas that need only the backwards values of fn . Whether you use or z -1 is mostly a matter of preference. The old astronomers, working without computers, used the difference tables because the math, though tedious, is straightforward. A computer tends to prefer z -1 . Of course, getting the expressions in z -1 takes yet another round of substitutions and collecting of terms. I'll just show three of the formulas here, which should meet all your needs.

First order: (Euler's method)


Second order:


Fourth order:


Since terms of the form z-k fn represent simple delays, they're easy to get using a software or hardware FIFO.

Finish line
I expect you can guess that there are many more things about numerical integration that we can talk about. I'll be having a few things to say about the topic in the future, but for the most part, we're done with the Rosetta Stone–hopefully for good and all. You've seen some easy transformations based on operational calculus. You've also seen about the worst case I've found, which is the numerical integration formula. Next time, I'm going to start off down a completely different road.

See you then.

Jack Crenshaw is a senior software engineer at General Dynamics and the author of Math Toolkit for Real-Time Programming, from CMP Books. He holds a PhD in physics from Auburn University. E-mail him at .

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.