next up previous
Next: QT's generalization to arbitrary Up: Shadowing Previous: Shadowing non-hyperbolic systems

Refinement of a numerical trajectory

The refinement procedure of Grebogi et al.  ([28], hereinafter referred to as GHYS) is analagous to Newton's method for finding a zero of a function. GHYS presented their algorithm for the two dimensional case. We will also be looking closely at a generalization by Quinlan and Tremaine ([62, 63], with [63] hereinafter referred to as QT) to arbitrary-dimensional Hamiltonian systems.

The basic idea is as follows. Let tex2html_wrap_inline2546 be a trajectory with S steps that has noise tex2html_wrap_inline2550 , where tex2html_wrap_inline2552 is the machine precision, and tex2html_wrap_inline2554 is some constant significantly greater than 1 that allows room for improvement towards the machine precision. Let tex2html_wrap_inline2282 be the 1-step error at step i+1, where tex2html_wrap_inline2560 for all i. The set of all 1-step errors is represented by tex2html_wrap_inline2564 , and is estimated by a numerical integration technique that has higher accuracy than that used to compute tex2html_wrap_inline2566 . This describes a function, which we call g, which takes as input the entire orbit tex2html_wrap_inline2566 and whose output is the set of 1-step errors tex2html_wrap_inline2572 , i.e.,  tex2html_wrap_inline2574 . Since the 1-step errors are assumed to be small, tex2html_wrap_inline2576 is small. That is, tex2html_wrap_inline2566 is close to a zero of g, if one exists. A zero of the function would represent an orbit in which the 1-step errors are identically zero, i.e.,  a true orbit. This is an ideal situation in which to run Newton's method. If the method converges, then a true orbit lies nearby.

Assume we have a noisy D-dimensional orbit tex2html_wrap_inline2546 , and it has a shadow tex2html_wrap_inline2586 . Then


where tex2html_wrap_inline2590 is an approximation to f with a typical accuracy or noise of tex2html_wrap_inline2248 . Here, tex2html_wrap_inline2248 should be the accuracy typically used by practitioners of the problem being studied. Now suppose we compute the 1-step errors


although numerically we use tex2html_wrap_inline2600 , a better approximation to f than tex2html_wrap_inline2590 . Then tex2html_wrap_inline2606 represents a correction term that perturbs point tex2html_wrap_inline2260 towards the shadow. So,


In the spirit of Newton's method, we ignore the tex2html_wrap_inline2610 in (7), and so one refinement iteration defines the corrections along the entire orbit:


where tex2html_wrap_inline2612 is the linearized perturbation map. For a discrete map, tex2html_wrap_inline2614 is just the Jacobian of the map at step i. For a system of ODEs, tex2html_wrap_inline2614 is the Jacobian of the integral of the ODE from step i to step i+1.gif The computation of the linear maps tex2html_wrap_inline2614 , called resolvents for an ODE system, takes most of the CPU time during a refinement. If we let D be the dimension of the problem, then the resolvent has tex2html_wrap_inline2680 terms in it, and it takes tex2html_wrap_inline2682 time to compute. Hayes and Jackson [33, 34] list several optimizations to the procedure that speed the process considerably, and should be considered if one is interested in using the algorithm in higher dimensions. Presumably, if one is interested in studying simpler high-dimensional systems, a chaotic map would be a better choice than an ODE system, because no variational equation integration is needed.

For simplicity of explanation, we assume D=2 for the remainder of this subsection, deferring discussion of the higher D case to the next subsection. If the problem were not chaotic, the correction terms tex2html_wrap_inline2688 could be computed directly from (8). But since tex2html_wrap_inline2614 will amplify any errors in tex2html_wrap_inline2688 not lying in the stable direction, computing the tex2html_wrap_inline2688 's by iterating (8) forward will amplify errors and typically produce nothing but noise. Therefore, we split the error and correction terms into components in the stable ( tex2html_wrap_inline2696 ) and unstable ( tex2html_wrap_inline2698 ) directions at each timestep:


Since it is not known a priori which direction is unstable at each timestep, the unstable vector tex2html_wrap_inline2700 at time tex2html_wrap_inline2668 is initialized to an arbitrary unit vector. The linearized map is then iterated forward with


Since tex2html_wrap_inline2614 magnifies any component that lies in the unstable direction, and assuming we are not so unlucky to choose a tex2html_wrap_inline2700 that lies too close to the stable direction, then after a few e-folding times tex2html_wrap_inline2698 will point roughly in the unstable direction at tex2html_wrap_inline2262 . Similarly, the stable unit direction vectors tex2html_wrap_inline2696 are computed by initializing tex2html_wrap_inline2716 to an arbitrary unit vector and iterating backward,


Substituting (10) into (8) yields


While tex2html_wrap_inline2614 magnifies errors in the unstable direction, it damps them in the stable direction. Likewise, tex2html_wrap_inline2720 damps errors in the unstable direction and magnifies errors in the stable direction. Thus the tex2html_wrap_inline2722 terms should be computed backward, and the tex2html_wrap_inline2724 terms forward. Taking components of (13) in the unstable direction at step i+1, iterate backward on


and taking components in the stable direction, iterate forward on


The initial choices for tex2html_wrap_inline2728 and tex2html_wrap_inline2730 are arbitrary as long as they are small -- smaller than the maximum shadowing distance -- because (15) damps initial conditions and (14) damps final conditions. GHYS and QT choose them both as 0. This choice is probably as good as any, but it can be seen here that if one shadow exists, there are infinitely many of them.gif Another way of looking at these initial choices for tex2html_wrap_inline2728 and tex2html_wrap_inline2730 is that they ``pinch'' the growing components at the end point, and the backward-growing components at the initial point, to be small. That is, boundary conditions are being forced on the problem so that the exponential divergence is forcibly masked, if possible, making the solution of (8) numerically stable.

Note that, counter-intuitively, these boundary conditions demand that the initial condition for the shadow and noisy orbits differ along the unstable direction. In fact, this must occur if the change in initial conditions is to have any effect. That is, when looking for a shadow, if perturbations were only allowed in the stable direction, those perturbations would die out, leading the ``shadow'' to follow the true orbit that passes through our initial conditions -- the one that is already known to diverge exponentially from the computed orbit.

next up previous
Next: QT's generalization to arbitrary Up: Shadowing Previous: Shadowing non-hyperbolic systems

Wayne Hayes
Fri Dec 27 17:41:39 EST 1996

Access count (updated once a day) since 1 Jan 1997: 9032