more on orbit problems, ellipsoids, and naive interval analysis

[David Hough]

After claiming that naive interval analysis was no panacea, and citing a
very old lecture, I wondered if I could prove what I claimed or at least
come up with some plausible details.
 
It took me all day to figure out what I really wanted to show, and what was
important, which is this:  conventional interval analysis can be appropriate
when the sup-norm is, and can be a bad deal when another norm is appropriate,
such as the two-norm.
 
The example is simpler than I thought at first:
a simple 1-dimensional orbit problem
will suffice (the attractor at the origin has a hole in it for the
satellite to pass through).  The differential equation is
 
        x'' = -x
        x(0) = x0
        x'(0) = x1
 
The solution is
 
        x(t) = x0 * cos(t) + x1 * sin(t)

Suppose that you can solve the differential equation perfectly (no analytical
errors and no rounding errors) so the only cause of dispersion is
error in the initial values x0 and x1.  Setting up as a conventional
linear system

        X(t) = ( x(t) ) = A * X(0) = ( c s) (x0)
               ( x'(t))              (-s c) (x1)

where c=cos(t) and s=sin(t),
we find that any error D(0) in the initial conditions X(0) maps linearly
into the error D(t) = A * D(0) in the evolved solution X(t).
How should we measure the evolved error D(t)?

We could use a p-norm

        Np(D) ** p = sum (|wi * di|**p)

which in one case of interest is

        Ninf(D) = max{ |w0 * d0|, |w1 * d1| }

Ninf(D) corresponds to naive interval analysis: the "size" of a set
is the largest dimension of an enclosing box whose sides are parallel to
the coordinate axes.

The other case of interest is

        N2(D) = sqrt ((w0 * d0)**2 + (w1 * d1)**2)
              = sqrt (D*.W.D)
where
                W = (w0**2 0    )

                    (0     w1**2)

N2(D) corresponds to measuring the size of a set by the largest
dimension of an enclosing ellipsoid.  Note that because A is
orthogonal and W is diagonal,

        N2(D(t)) = sqrt(D(0)*.W.D(0))

independent of t!  There is no growth in the size of the error bound.
That reflects the fact that this is a very stable system.

Ninf(D(t)) doesn't reflect that stability, since it's really

        Ninf(D(t)) = max { w0 * | c * d0 + s * d1 |,
                           w1 * |-s * d0 + c * d1 | }


If Ninf(D(0)) == N0, looking at t=pi/2 and considering the extremes
of the starting enclosing box, where |w0 * d0| == N0 and |w1 * d1| == N0,
we see that

        Ninf(D(pi/2)) = max { w0 * |d1|, w1 * |d0| }
                      = max {w0/w1,w1/w0} * N0

which means that the size of the bounding set must grow by an exponential
factor with time unless w0 == w1; but in that case, consider t=pi/4
whence

        Ninf(D(pi/4)) = w0 * sqrt(0.5) * max {|d0+d1|,|d0-d1|}
                      = w0 * sqrt(0.5) * (|d0|+|d1|)
                      = w0 * sqrt(0.5) * (N0/w0 + N0/w1)
                      = sqrt(2) * N0


so the size of the bounding set must grow by an exponential factor in
this case as well.

What's really going on here is that the actual reachable set is shaped
much more like an ellipsoid than a bounding box, even if it starts
out as a box; bounding that set by
a box incorporates much a much larger set than the reachable set;
and the reachable set from that box will be somewhat larger than from
the ellipsoid; compounding this over many steps, you get exponentially
larger growth in the size of the bounding box than in the bounding
ellipsoid.

People who teach or work with this sort of thing on a daily basis will no doubt
point out the errors in the foregoing, but hopefully the main idea is
correct and instructive.

Realistic orbit problems involve three x dimensions and a corresponding
set of x' dimensions, and the attractor is a lumpy planet with an atmosphere
rather than a point mass, but I think the fundamental issue is the same.