Like forward Euler, and like backward Euler,
the wolff algorithm is a useful thing to talk
about in a class because it helps people understand
the basic ideas, but no one really uses it anymore.
Although, you will see an occasional exception
to that statement, usually in the literature of a
different field. The primary drawback of the
wolff algorithm is its single point nature - the
fact that you are picking a starting point and
looking for the nearest neighbor at every step.
As you can well imagine, any noise in the data
might really trash that. The Kantz algorithm
and the family of related methods address that
by picking a whole bunch of points, tracking them
in parallel, averaging their distance to some central
point, and watching how that quantity grows.
As in the wolff algorithm, you pick a starting point
in the trajectory, then you draw an epsilon ball
around it, find all the points in the trajectory that
are inside that epsilon ball and measure their
distances to that central point. Finally, you average
all those distances. Now remember, these are all
points in a single trajectory. For each one of them, you
know where it goes next. If the center point,
for example, is the 17th point in the trajectory,
you know where the 18th point is, so what you do
is you figure out where each of those points
goes next...and I'm going to redo this drawing
so it's a little less confusing...and this drawing,
I'm going to draw in green where each point
goes next. This is where the middle one goes next.
This one goes here next.
This one will go here next.
And this one here next, and so on and so forth.
This drawing doesn't exactly match the one on the
left by the way, it's just a schematic. Then what I do
is I compute the distances between the forward image
of that center point - which is where the first one went -
and the forward images of all those other points.
The ratio of the average of all the red distances in
this plot to all the average of all the red distances in
this plot is a measure of how much the dynamics is
stretching out state space around that central point.
It's called the stretching factor and that stretching
factor is exactly what the Lyapunov exponent is
trying to capture. The Kantz algorithm repeats this
calculation for a range of different time lengths
and a bunch of different initial points.
If you plot the results on a log of delta s vs time curve,
you'll see a curve like this. Now what's going on here?
When time is small, that means that you're not giving
things much time to stretch out, so they won't stretch
and the ratio will be 1 and the log will be 0.
If you let them stretch for a really really long time
however, the points in the epsilon ball will spread out all
over the attractor and that's as far apart as they can
get and that's as far as the volume can grow.
That's this upper flat area. In between those
asymptotes, there's a scaling region...
I apologize for the weakly writing by the way,
I downloaded a new copy of the sketchpad software
that I use on the tablet and it's making my writing
look awful... Let's think about what a diagonal
line on a loglinear plot means. What that means
is that the natural log of delta s is proportional
to time and the constant of proportionality is a.
If we take e to the both sides of this, we'll get...
where a is this slope. And state space stretches as
e to the lambda t, so the slope of the scaling region, a,
is the Lyapunov exponent. You saw a similar thing
in the previous segment about calculating the
fractal dimension by the way, although the axes
were different. There, we were after the
relationship between the log of the number of boxes
of size epsilon to cover something and the log of
one over the size of those boxes, so we plotted
log of one over epsilon on the horizontal axis, looked
for a scaling region on that curve and took its
slope as the value of the capacity dimension.
In any calculation like these two things that depends
on the existence of a scaling region, it's critically
important not to impude the existence of one -
that is, to see one there just because you want it there
when it's really not. And that's a really subtle problem
because scaling regions are very hard to define.
You could, for example, fit a line to a chunk of your
curve and insist that the r-squared value of that line
be above 0.9 or something, but that's an arbitrary
threshold. All of this brings up an interesting and
important point about nonlinear time series analysis.
One that's come up before. This is a very very powerful
tool, but it has to be used wisely. And what I mean by
wisely, is by a person in the loop. Someone needs
to look at those plots and see if they have the right
shape before deciding to pick off a fractal dimension
or a Lyapunov exponent. If it doesn't have the right
shape from your perspective, not mine.
Or if that shape doesn't persist as you turn the
knobs on the algorithm, like the size of that
epsilon ball in the Kantz algorithm, you shouldn't
use that curve to pick off those values. And that
means, among other things, that these kinds
of calculations are all but impossible to automate.
By that I mean that if you have 1000 data sets and
you want to compute the Lyapunov exponent of
all 1000 of them, it's not going to work to
simply throw those in the hopper of either of these
algorithms and go off and have a beer.
The answers won't be right. You need to be involved
with the process. That's because the problem is
hard and subtle, and you now know how to do that.
The reason for this is not that these tools are lousy,
but rather that nonlinear time series is a very hard
problem. Linear tools are easy to use and they
always give you an answer, but if you apply them to
a nonlinear system, that answer may well be very
wrong. Remember the lamp post. Knowing how these
algorithms work should also help you understand
and respect the requirements for what kind of data
are needed for them to work properly.
The data set needs to be long enough to cover
the behavior that you are trying to sample. It needs
to be quickly enough sampled to catch all the details
in that, and it needs to not have so much noise
as to obscure those details. And those are tough
requirements sometimes.
The TISEAN toolset, by the way, includes the Kantz
algorithm for computing Lyapunov exponents from
data. All of that was about calculating lambda 1,
the largest positive exponent. There are n-1
other Lyapunov exponents in an n-dimensional
system. How to calculate them? That's much harder.
Instead of following the spreading a 2-D volume, as
sampled by the data, you have to follow higher-
dimensional volumes. The algorithms, and the numerics,
and the data requirements were hard enough for
following the 2-D volumes. The 3-D version occasionally
works, giving you lambda 2, but I've never seen
results beyond that that I'm really happy with.