Orbits Worth Betting On

Rob  Knapp
Wolfram Research, Inc.
100 Trade Center Drive, Champaign, IL 61820
rknapp@wri.com

and

Stan Wagon
Department of Mathematics and Computer Science
Macalester College, St. PaAul, MN 55105
wagon@macalstr.edu

1. The Problem

The numerical approach to solving a differential equation is by no means infallible. Interesting equations are often very sensitive to initial conditions, and therefore great care must be taken when using numerical software to attack such equations. Things may be fine for small values of t, but it often happens that, even for moderately sized t-values, machine precision (typically 16 digits) and default error goals (typically, 6 digits) are insufficient to obtain reliable results. Nevertheless, many “solutions” to such equations have been published in the belief that they are correct. In this note we will discuss in some detail the steps that one can take to ensure, or at least to put the odds in one’s favor, that a solution to a sensitive equation is correct.

The Duffing equation is a good example. The Spring 1994 issue of CODEE [C, Figure 2] shows a Duffing orbit with t running from 0 to 500. But, as we shall see, for this Duffing equation machine precision is inappropriate for values of t beyond 70. The published and actual trajectories are quite different. A similar situation occurs in [BCB, p. 298] where the same equation is solved for t up to 100. Again, the orbit is wrong. In fact, our interest in this topic arose because of our inability to duplicate published orbits. Both [BCB] and [C] show Poincaré sections with t running up to 10,000 or more; indeed, such an image appears on the CODEE poster. Another example: the Duffing solution given by Edwards and Penney [EP, p. 575, Fig. 7.6.4] is incorrect. Both [BCB] and [EP] give as a student exercise the task of reproducing their figures. Instructive as such attempts might be, they are essentially doomed, since the results will be strongly machine and software dependent. Indeed, one of the authors of this note has produced a computer laboratory project [K] where students are asked to use code that generates solutions (inaccurate in this same sense for long time intervals) for the Duffing equation,
but are not asked to consider the validity of the solutions.

It is important to keep in mind that in a real study of the Duffing equation, the information needed is generally not the solution to some given accuracy. In fact, a tremendous amount of information can be accurately obtained through clever use of machine-precision calculations. Further, since the system is believed to have a strange attractor, numerical approximations should be drawn towards this attracting set, if not toward the actual solution. On the other hand, it is important to distinguish between an inaccurate approximate solution that may be used to draw some information about the attractor and an accurate solution that can be used in a numerical proof of properties of the system. In fact, due to this very distinction, Guckenheimer and Holmes [GH, p. 91] state about this system that, “There is a substantial theoretical question as to whether this ‘strange attractor’ of the Duffing equation is an artifact of the noise and is absent in the ideal deterministic system.”

Additionally, there is the so-called “shadowing theory”, which can be used in many cases to show that even if an approximate numerical solution with a particular initial condition may not remain close to the actual trajectory for that initial condition, the approximate solution does lie close to the actual trajectory for some nearby initial condition [GH, p. 251]. This theory can be used to justify significant results, but it is not sufficient in many cases, for example in the proof of the existence of a strange attractor. Finally, although we focus on the Duffing equation here, our goal is to present tools that can be applied to just about any numerical computation, whether or not one suspects that shadowing theory or the theory of strange attractors applies.

We use Mathematica throughout (version 3.0; version 2.2 behaves somewhat differently), but have not included much of the code here. A Mathematica notebook containing this text and the code that generated the output is available from MathSource, which can be accessed via the Wolfram Research web page, http://www.wri.com/.

2. Detecting a Bad Bet

Although there are various seat-of-the-pants methods to detect whether an orbit might be wrong (redo the computation with higher precision, or on a different machine), the proper approach requires a rudimentary understanding of what the numerical algorithm is doing, not a bad idea in any case. Numerical schemes typically proceed step-by-step, but in an adaptive way, attempting to make each step satisfy a certain error tolerance. Mathematica’s NDSolve, for example, tries for 6-digit precision at each step. This aspect of the error is called truncation error (as opposed to roundoff error or propagation error). Now, since the result of the first step is only accurate to 6 digits, this inaccuracy might propagate unacceptably through the perhaps thousands of subsequent steps. Whether or not it does so has to do with the propagation error for the equation, discussed more fully in section 4. An easy way to check for this is to perturb the initial conditions by 10–6.

Duffing = (x''[t] + (15/100) x'[t] - x[t] +(x[t])^3 == (3/10) Cos[t]);

lowPrecSoln = NDSolve[{Duffing, x[0] == -1, x'[0] == 1}, 
   x, {t, 0, 100}, MaxSteps -> 2000][[1]];

perturbed   = NDSolve[{Duffing, x[0] == -1 + 10^-6, x'[0] == 1}, 
   x, {t, 0, 100}, MaxSteps -> 2000][[1]];

x[100] /. {lowPrecSoln, perturbed}

The serious disagreement indicates that a 10–6 difference can lead to a radically different solution. Of course, we are not looking at any true solutions in these computations; rather we are observing that small changes of the sort that will occur in the computation cause major changes in the computed result, and therefore we can place no confidence in the results. A view of the graphs (Figure 1) provides some evidence that the solution is all right up to t = 70.

Code for Figure 1

Plot[Evaluate[(x[t] /. {lowPrecSoln, perturbed})], {t, 0, 100},
  PlotStyle ->{{AbsoluteThickness[4], GrayLevel[0.6]},
               {AbsoluteThickness[1.5]}}];

Figure 1.  Two Duffing orbits starting from nearby initial conditions. They agree only up to about t = 70. The unperturturbed solution is shown with thick light gray lines and the solution with perturbation in black.

Code for 10^-10 perturbation

x[100] /. NDSolve[{Duffing, 
   x[0] == -1 + 10^-10, x'[0] == 1}, x, {t, 0, 100}, MaxSteps->2000][[1]]

A perturbation of 10–10 leads to x(100) = –0.42444, and so has only a small effect on the computed result. This is because the perturbation is negligible compared to the local truncation error, so it is effectively lost. For these tests to be meaningful, the perturbation must be of the same order of magnitude or larger than the local truncation error.

In Mathematica we can change the local error of the numerical algorithm by increasing the accuracy goal to, say, 10. (In Mathematica, “precision” refers to relative error and “accuracy” refers to absolute error; we will use the latter in this article since these Duffing solutions have modest scales.) When we increase the accuracy goal to 10, then a 10-10 perturbation leads to the results 0.23 and –0.52 at t = 100. Thus, as before, we are still getting inaccurate solutions. But the increased accuracy goal (which causes the numerical algorithm to work harder, as will be discussed in section 3) has improved things in that the perturbed and unperturbed solutions now agree out to about t = 80 (Figure 2).

Code for Figure 2

lowPrecSoln = 
   NDSolve[{Duffing, x[0] == -1, x'[0] == 1}, x, {t, 0, 100}, 
           MaxSteps->2000, AccuracyGoal->10][[1]];

perturbed6 = NDSolve[{Duffing, x[0] == -1 + 10^-6, x'[0] == 1}, x,
  {t, 0, 100}, MaxSteps->2000, AccuracyGoal->10][[1]];
  
perturbed10 = NDSolve[{Duffing, x[0] == -1 + 10^-10, x'[0] == 1}, x,
  {t, 0, 100}, MaxSteps->2000, AccuracyGoal->10][[1]];
  
x[100] /. {lowPrecSoln, perturbed6, perturbed10}

Plot[Evaluate[(x[t] /. {lowPrecSoln, perturbed6, perturbed10})],
  {t, 0, 100},
  PlotStyle->{
  	{AbsoluteThickness[4], GrayLevel[.8]},
  	{AbsoluteThickness[2], GrayLevel[.6]},
  	{AbsoluteThickness[1.5]},
  }];

Figure 2.  A higher accuracy goal causes the perturbed solution to agree with the unperturbed case out to t = 80. The unperturturbed solution is shown with thick, light gray lines, the solution with perturbation 10–6 in gray, and the one with perturbation 10–10 in black.

With the more accurate solution, it is clear that the divergence of solutions with initial conditions of order 10–6 is not simply an artifact of the numerical method or the choice of size of the local truncation error. On the contrary, the divergence is a property of the differential equation. However, it is true that the smaller the perturbation (or truncation error for that matter), the longer the time interval for which solutions match before diverging. This will be quantified in Section 4.

One must not yield to the temptation to simply increase the accuracy goal to, say, 20, in the hope that this would yield good results past t = 100. That would be very misleading because the inherent imprecision of hardware arithmetic (typically 16 digits) makes such an accuracy goal unattainable.

To see that this technique does indeed discriminate, consider a different initial condition: (0.6, 1.3). Then, going back to the default accuracy goal, a 10–6 perturbation leads to results that are very close: –0.0147109 and –0.0147116. And Figure 3 shows that the solutions are essentially identical and asymptotically periodic.

Code for Figure 3

orbit = NDSolve[{Duffing, x[0] == 0.6, x'[0] == 1.3}, x, 
          {t, 0, 100}, MaxSteps -> Infinity][[1]];

porbit = NDSolve[{Duffing, x[0] == 0.6 + 10^-6, x'[0] == 1.3}, x, 
          {t, 0, 100}, MaxSteps -> Infinity][[1]];

x[100] /. {orbit, porbit}

Plot[Evaluate[x[t] /. {orbit, porbit}], {t, 0, 100},
  PlotStyle -> {{AbsoluteThickness[4], GrayLevel[.6]},
                {AbsoluteThickness[1.5]}}];

Figure 3.  For a different set of initial conditions, the perturbed and unperturbed solutions agree well.

3. Putting the Odds in Your Favor

The easiest way to obtain a solution worth betting on is to increase the precision (number of digits) of the numbers used to do the calculation and to correspondingly decrease the local truncation error for the approximation to the differential equation.

In Mathematica, the precision of numbers used to do the calculation can be set using the WorkingPrecision option. The error tolerances used to limit the local truncation error are by default set to be 10 digits less than the working precision, which is fine for our purposes. Note that it is crucial that the constants in the equation be symbolic — 3/10 as opposed to 0.3 — because 0.3 is parsed to be a machine number, which doesn’t have the number of digits needed at values of working precision higher than 16 (machine-precision). These computations take a long time, but if accuracy is essential it is worth the time to be careful.

Code for Table 1

Do[wpsol[wp] = NDSolve[{Duffing, x[0] == -1,x'[0] == 1}, x,
    {t, 0, 100}, MaxSteps -> Infinity,
      WorkingPrecision -> wp][[1, 1, 2]],
  {wp, 16, 32, 4}];

We omit the details of the computation that produced Table 1. Briefly, we computed solutions using 16, 20, 24, 28, and 32 digits of working precision. The local truncation errors for these solutions should be 10–6, 10–10, 10–14, 10–18, and 10–22, respectively. Assuming there is convergence, the most accurate solution will be the one computed with 32-digit numbers, so we compare the other solutions to that one. Table 1 shows the results of such a comparison, where the data in the table are base-10 logarithms of the differences. The sequence of solutions with more and more stringent error tolerances is in effect telling us how solutions with perturbations of different size converge to a solution at t = 100.

Code for nicely formatted Table 1

ppdiff[a_,b_][t_] := Sqrt[(a[t] - b[t])^2 + (a'[t] - b'[t])^2];

Table[Join[{wp},
	      Table[Round[Log[10, ppdiff[wpsol[32],wpsol[wp]][t]]],{t,20,100,20}],
      	{Length[wpsol[wp][[1,3,1]]]}],
	{wp, 16, 28, 4}]

TableForm[{
{WP,"t=20", 40, 60, 80, 100, Steps},
{"--","----","--", "--", "--", "--", "-----"},
{16, -4, -2, -1, 0, 0, 1067}, {20, -8, -6, -5, -2, 0, 2279}, 
 
   {24, -12, -10, -9, -7, -5, 5125}, {28, -16, -13, -13, -10, -8, 10521},
   {32, "", "", "", "", "", 19536}},
TableSpacing -> {1, 5},
TableAlignments->{Automatic, Right}]

Table 1.   Logarithms of differences show the convergence to the 32-digit solution.

The first column gives the working precision, the central columns give the base-10 logarithms of the “error” as t equals 20, 40, 60, 80, and 100, and the last column shows the number of steps that NDSolve used. It is, of course, comforting to see the number of steps increase, as we expect because of the increased accuracy demands. More important, we can see that the solutions are happily converging towards the one obtained using 32-digit precision. The results for different t-values show that there is a loss of accuracy as t increases, but that convergence towards the 32-digit solution is consistent.

Based on this data, we would be willing to bet that the solution calculated with 28-digit numbers, illustrated in Figure 4, is accurate to six digits at t = 100.

Code for Figure 4

asol = Compile[{x},Evaluate[wpsol[28][x]]];
(* we will call the accurate solution asol for convenience *)

ParametricPlot[{asol[t],asol'[t]},{t,0,100}, AxesLabel->{"x","x'"}];

Figure 4.  The Duffing trajectory accurate to six digits up to t = 100 starting from (–1, 1). This required 28 digits of precision and 10,521 steps in the numerical algorithm.

If we perturb the initial condition by a perturbation on the order of the local truncation error used to compute the accurate solution (10–18), we find that the solutions are the same in the first nine digits even at t = 100. The difference is 2·10–10, which is smaller than the 10–8 we might have expected from Table 1.  This is further evidence that we have computed an accurate solution. On the other hand, if the initial condition is perturbed by 10–6, this allows us to compute accurately how two solutions really diverge. This important idea is pursued in the next section.

Code for 10^-18 perturbation

pasol = NDSolve[
{Duffing, x[0] == -1 + 10^-18, x'[0] == 1}, x,
    {t, 0, 100}, MaxSteps -> Infinity,
      WorkingPrecision -> 28][[1, 1, 2]];

{wpsol[28][100], pasol[100], 
 wpsol[28][100] - pasol[100]}// TableForm

4. The Lyapunov Characteristic Exponent

Why does this equation give the differential equation solver such fits? The reason is that the system is sensitive to initial conditions in the sense that the trajectories for initial conditions near (–1, 1) tend to diverge exponentially fast from the trajectory starting at (–1, 1). To see if this is indeed the case, we compare the 28-digit solution from (–1, 1) with a perturbed one (by 10–6). Because we are using the high-precision solution, we have some confidence that it is the true solution. Thus we can plot the phase-plane distance between the two solutions as t increases. We do so on a semi-log plot.

Code for Figure 5

psol = NDSolve[{Duffing, x[0] == -1 + 10^-6, x'[0] == 1}, x, {t, 0, 100},
	MaxSteps -> Infinity, WorkingPrecision -> 28, 
	InterpolationPrecision->16][[1, 1, 2]];

ppdiff[a_,b_][t_] := Sqrt[(a[t] - b[t])^2 + (a'[t] - b'[t])^2];

<<Graphics`Graphics`;
LinearLogPlot[ppdiff[asol, psol][t], {t, 0, 100},
	FormatType->AccountingForm];

Figure 5.  Divergence of trajectories with a 10-6 perturbation in the initial conditions.

It appears that, at least initially, the difference is roughly linear on this plot, indicating an exponential relationship. To find the exponent we need to find a line that fits the logarithm of the data. It is appropriate to use only the data up to the point where the difference is order one. This is because, although a perturbation causes exponential divergence locally, solutions near this initial condition are believed to be attracted to something called a “strange attractor”, which is a bounded set with zero area. Since this set is bounded, the divergence cannot continue indefinitely. A regression on the data up to t = 60 gets us a reasonable exponential function to model the divergence: 0.792 e0.178 t.

Code for exponential fit

ComplexExpand[Exp[
  Fit[Table[{t, Log[ppdiff[asol, psol][t]]}, {t, 0, 60, 1}],
      {1, t}, t]]] / 10.^(-6)

A mathematical way of expressing this is to say that, for the initial condition (–1, 1), the equation has a positive Lyapunov characteristic exponent (LCE). Below we will give a relatively nontechnical description of the definition of the LCE’s of this dynamical system. A more precise definition of this important concept can be found in many texts on dynamical systems theory. We use the form of the definition given in [BGGS I].

Consider the linearized equation, y¢¢(t) + (15/100) y¢(t) + (3 x(t)2 – 1) y(t) = 0, where x(t) is the solution to Duffing’s equation starting at X0 = (x(0), x¢(0)). Let at be the linear map such that atv = (y¢(t), y(t)), where (y(0), y¢(0)) = v. Then an LCE is defined by

         c(X0, v) = lim suptÆ∞  (1/t) ln ||at v||.

Notice that the LCE c depends on the initial condition as well as the vector v. The linearization is necessary for the definition because without this, the limit would eventually be zero since solutions are attracted to a bounded set as mentioned above. The linearization allows for an asymptotic result that describes the divergence (or convergence of trajectories).

For the initial condition (–1, 1) it can be shown that there are two different values of the LCE, approximately 0.18 and –0.34. There is a one-dimensional subspace of vectors v that have c((–1,1), v) = –0.34. Any vector not in this subspace will have have some component that has the positive value. The component with positive value will magnify and asymptotically this will dominate. Thus, all vectors u not in the subspace of vectors v with c((–1, 1), v) = –0.34 will have c((–1, 1), u) = 0.18. This is why there are only two values. (For dynamical systems with more than two degrees of freedom, there will be more values.) The model found above gives a crude approximation to one of the LCEs for the Duffing equation for this initial condition. The negative LCE is difficult to compute because any numerical representation of the vector v will have some error that acts like a component outside of the subspace and will eventually dominate. To compute the more accurate values, 0.18 and –0.34, we used the method outlined in [BGGS I, II], which takes into account machine-number errors.

Code for computing the LCEs

This program is based on the algorithm in [BGGS II] for computing all of the LCEs for a dynamical system

Clear[LCE];
LCE[n_,s_,X0_,w0_] := Module[{t,tau,X,solX,
w,solw,alpha,sum = 0,
w1,solw1,alpha1,sum1 = 0},
X = X0;
w = w0/Sqrt[w0.w0];
w1 = {-w[[2]],w[[1]]};
deq := (x''[tau] + (15/100) x'[tau] - x[tau] + (x[tau])^3 == (3/10) Cos[t]);
leq := (x''[tau] + (15/100) x'[tau] + x[tau](3(solX[tau])^2-1) == 0);
Table[
	t = tau + Mod[(k-1)*s,2 Pi];
	solX = NDSolve[
		{deq,x[0] == X[[1]],x'[0] == X[[2]]},
		x,{tau,0,s}][[1,1,2]];
	solw = NDSolve[
		{leq,x[0] == w[[1]],x'[0] == w[[2]]},
		x,{tau,0,s}][[1,1,2]];
	solw1 = NDSolve[
		{leq,x[0] == w1[[1]],x'[0] == w1[[2]]},
		x,{tau,0,s}][[1,1,2]];
	X = {solX[s],solX'[s]};
	w = {solw[s],solw'[s]};
	w1 = {solw1[s],solw1'[s]};
(* Normalize w--the amount of stretching is alpha *)
	alpha = Sqrt[w.w];
	w = w/alpha;
(* Make w1 orthogonal to w -- alpha1 is the size of the 
   compnent in this direction - -this is the shrinking
   from the negative LCE *)
	w1 = w1 - (w.w1) w ;
	alpha1 = Sqrt[w1.w1];
	w1 = w1/alpha1;
(* The LCE is the sum of Log of the stretching.shrinking
   factors divided by (ks) (like 1/t in def'n *)
	sum += Log[alpha];
	sum1 += Log[alpha1];
	{sum,sum1}/(k s),
	{k,1,n}]]

seq = LCE[10000,1,{-1,1},{Random[],Random[]}];
tail = Take[seq,250];

This should give the larger (positive) LCE.

SequenceLimit[First /@ tail]

This should give an estimate for the smaller (negative) LCE.

SequenceLimit[Last /@ tail]

The presence of both a positive and a negative Lyapunov exponent indicates that there are distinct directions along which stretching and shrinking occurs. Under the map induced by the differential equation, a circle of initial conditions ends up being stretched on one axis and shrunk on another. The map also induces some folding. The repeated stretching, shrinking, and folding is effectively producing the sensitivity to initial conditions. The stretching and shrinking can be seen particularly well by observing how a circle of initial conditions evolves around its center. In Figure 6 the evolution of a circle of radius 0.1 centered at (–1, 1) is shown.

Code for Figure 6

sol0 = With[{wp = 16},Part[
	NDSolve[{x''[t] + (15/100) x'[t] - x[t] +(x[t])^3== (3/10) Cos[t],
		x[0] == -1,x'[0] == 1},
		x,
		{t,0,7},
		MaxSteps->100000,
		WorkingPrecision->wp,
		PrecisionGoal->wp/2,
		AccuracyGoal->Infinity],1,1,2]];
r = 10^(-1);
n = 1000;
Do[solc[j] = With[{wp = 16},Part[
	NDSolve[{x''[t] + (15/100) x'[t] - x[t] +(x[t])^3== (3/10) Cos[t],
		x[0] == -1+r*Cos[j],x'[0] == 1+r*Sin[j]},
		x,
		{t,0,7},
		MaxSteps->100000,
		WorkingPrecision->wp,
		PrecisionGoal->wp/2,
		AccuracyGoal->Infinity],1,1,2]],
	{j,2*Pi/n,2*Pi,2*Pi/n}];

tab = Table[ListPlot[
		Table[{solc[j][t],solc[j]'[t]},
			{j,2*Pi/n,2*Pi,2*Pi/n}],
		PlotJoined->True],
	{t,0,7}];

Show[tab,
	ParametricPlot[{sol0[t],sol0'[t]},{t,0,7},
		PlotStyle->{{AbsoluteThickness[2], GrayLevel[0.6]}},
		DisplayFunction->Identity],
	DisplayFunction->$DisplayFunction];

Figure 6.  Evolution of a circle of initial conditions. The trajectory shown in gray is the solution with initial condition (–1, 1). The circle has radius 0.1. The evolution is shown up to t = 7 at intervals Dt = 1.


Notice that the area within the evolving circle seems to shrink. This is indeed the case because the shrinking LCE has larger magnitude than the stretching one. This is an important observation because it tells us that the strange attractor (if indeed it does exist) will have area 0. 

Now, how does this analysis affect the loss of digits? Let’s take a look at how the rate of loss of digits relates to the rate of divergence of solutions. For the solution done with 28 digits of precision, Mathematica chooses the step size so that the error in each step is about 10–18. Over the time interval [0, 100], roughly 10 digits of accuracy are lost (18 – 8; from Table 1). This means that we expect that a small perturbation will diverge on the order of 1010 in time 100. We can now look for a model for the loss of accuracy having the form 10digits lost up to t = ert. A computation using the data for t = 100 gives 1010 = e100r, or r = 0.23, which is not far from the computed LCE of 0.18. Note that this estimate is bound to be a bit off since it does not take into account the number of steps, which will tend to increase the error. Nonetheless, the Lyapunov exponent is a good way to estimate how many digits might be necessary to compute an accurate solution from a given initial condition out to a certain time. For example, if we require the solution at t = 10,000 with three digits of accuracy, we could expect the error to grow by exp(0.18 · 10000), or about 10800. Thus we would have to carry out the computation using at least 803 digits of precision at each step.

Now, 800-digit numbers are themselves not impossible to deal with. But the big problem here is the number of steps that would be necessary. An extrapolation using the data in Table 1 allows us to predict that about 1066 time steps would be necessary!

Note that the definition of the Lyapunov exponent includes a dependence on the initial condition. That this is indeed real can be observed by considering that initial conditions that have trajectories attracted to the stable orbit shown in section 2 have two negative Lyapunov exponents. In this case, solutions tend to do the opposite of diverging — they converge.

Moral: In the immortal words of Sergeant Phil Esterhaus of Hill Street Blues: “Be careful out there!”

References

[BGGS I, II] G. Bennettin, L. Galgani, A. Giorgilli, and J. Strelcyn, Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part I: Theory, Part II: Numerical application, Meccanica 15 (1990) 9–30.

[BCB]  R. L. Borrelli, C. Coleman, W. E. Boyce, Differential Equations Laboratory Workbook, Wiley, New York, 1992.

[C]  C. Coleman, Duffing’s donut, CODEE Newsletter, Spring 1994, 11–13.

[EP]  C. H. Edwards, Jr., and D. E. Penney, Elementary Differential Equations with Applications, Prentice-Hall, Englewood Cliffs, N.J., 1994.

[GH] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 1983.

[K] R. Knapp, Duffing’s equation, a computer laboratory project (Mathematica notebook) available via anonymous ftp from ftp.bowdoin.edu.