Wayne B. Hayes

Robust and reliable defect control for Runge-Kutta methods

The quest for reliable integration of initial value problems (IVPs) for ordinary differential equations (ODEs) is a long-standing problem in numerical analysis. At one end of the reliability spectrum are fixed stepsize methods implemented using standard floating point, where the onus lies entirely with the user to ensure the stepsize chosen is adequate for the desired accuracy. At the other end of the reliability spectrum are rigorous interval-based methods, that can provide provably correct bounds on the error of a numerical solution. This rigour comes at a price, however: interval methods are generally two to three orders of magnitude more expensive than fixed stepsize floating-point methods. Along the spectrum between these two extremes lie various methods of different expense that estimate and control some measure of the local errors and adjust the stepsize accordingly. In this article, we continue previous investigations into a class of interpolants for use in Runge-Kutta methods that have a defect function whose qualitative behavior is asymptotically independent of the problem being integrated. In particular the point, in a step, where the maximum defect occurs as h → 0 is known a priori. This property allows the defect to be monitored and controlled in an efficient and robust manner even for modestly large stepsizes. Our interpolants also have a defect with the highest possible order given the constraints imposed by the order of the underlying discrete formula. We demonstrate the approach on three Runge-Kutta methods of orders 5, 6, and 8, and provide Fortran and preliminary Matlab interfaces to these three new integrators. We also consider how sensitive such methods are to roundoff errors. Numerical results for four problems on a range of accuracy requests are presented.

Optimal Network Alignment with Graphlet Degree Vectors

Important biological information is encoded in the topology of biological networks. Comparative analyses of biological networks are proving to be valuable, as they can lead to transfer of knowledge between species and give deeper insights into biological function, disease, and evolution. We introduce a new method that uses the Hungarian algorithm to produce optimal global alignment between two networks using any cost function. We design a cost function based solely on network topology and use it in our network alignment. Our method can be applied to any two networks, not just biological ones, since it is based only on network topology. We use our new method to align protein-protein interaction networks of two eukaryotic species and demonstrate that our alignment exposes large and topologically complex regions of network similarity. At the same time, our alignment is biologically valid, since many of the aligned protein pairs perform the same biological function. From the alignment, we predict function of yet unannotated proteins, many of which we validate in the literature. Also, we apply our method to find topological similarities between metabolic networks of different species and build phylogenetic trees based on our network alignment score. The phylogenetic trees obtained in this way bear a striking resemblance to the ones obtained by sequence alignments. Our method detects topologically similar regions in large networks that are statistically significant. It does this independent of protein sequence or any other information external to network topology.

Algorithm 908: Online Exact Summation of Floating-Point Streams

We present a novel, online algorithm for exact summation of a stream of floating-point numbers. By “online” we mean that the algorithm needs to see only one input at a time, and can take an arbitrary length input stream of such inputs while requiring only constant memory. By “exact” we mean that the sum of the internal array of our algorithm is exactly equal to the sum of all the inputs, and the returned result is the correctly-rounded sum. The proof of correctness is valid for all inputs (including nonnormalized numbers but modulo intermediate overflow), and is independent of the number of summands or the condition number of the sum. The algorithm asymptotically needs only 5 FLOPs per summand, and due to instruction-level parallelism runs only about 2--3 times slower than the obvious, fast-but-dumb “ordinary recursive summation” loop when the number of summands is greater than 10,000. Thus, to our knowledge, it is the fastest, most accurate, and most memory efficient among known algorithms. Indeed, it is difficult to see how a faster algorithm or one requiring significantly fewer FLOPs could exist without hardware improvements. An application for a large number of summands is provided.

Fitting Selected Random Planetary Systems to Titius–Bode Laws☆

Abstract Simple “solar systems” are generated with planetary orbital radii r distributed uniformly random in log r between 0.2 and 50 AU, with masses and order identical to our own Solar System. A conservative stability criterion is imposed by requiring that adjacent planets are separated by a minimum distance of k times the sum of their Hill radii for values of k ranging from 0 to 8. Least-squares fits of these systems to generalized Bode laws are performed and compared to the fit of our own Solar System. We find that this stability criterion and other “radius-exclusion” laws generally produce approximately geometrically spaced planets that fit a Titius–Bode law about as well as our own Solar System. We then allow the random systems the same exceptions that have historically been applied to our own Solar System. Namely, one gap may be inserted, similar to the gap between Mars and Jupiter, and up to 3 planets may be “ignored,” similar to how some forms of Bodes law ignore Mercury, Neptune, and Pluto. With these particular exceptions, we find that our Solar System fits significantly better than the random ones. However, we believe that this choice of exceptions, designed specifically to give our own Solar System a better fit, gives it an unfair advantage that would be lost if other exception rules were used. We compare our results to previous work that uses a “law of increasing differences” as a basis for judging the significance of Bodes law. We note that the law of increasing differences is not physically based and is probably too stringent a constraint for judging the significance of Bodes law. We conclude that the significance of Bodes law is simply that stable planetary systems tend to be regularly spaced and conjecture that this conclusion could be strengthened by the use of more rigorous methods of rejecting unstable planetary systems, such as long-term orbit integrations.

Is the Outer Solar System Chaotic

One-sentence summary: Current observational uncertainty in the positions of the Jovian planets precludes deciding whether or not the outer Solar System is chaotic. 100 word technical summary: The existence of chaos in the system of Jovian planets has been in question for the past 15 years. Various investigators have found Lyapunov times ranging from about 5 millions years upwards to infinity, with no clear reason for the discrepancy. In this paper, we resolve the issue. The position of the outer planets is known to only a few parts in 10 million. We show that, within that observational uncertainty, there exist Lyapunov timescales in the full range listed above. Thus, the “true” Lyapunov timescale of the outer Solar System cannot be resolved using current observations. 100 word summary for general public: The orbits of the inner planets (Mercury, Venus, Earth, and Mars) are practically stable in the sense that none of them will collide or be ejected from the Solar System for the next few billion years. However, their orbits are chaotic in the sense that we cannot predict their angular positions within those stable orbits for more than about 20 million years. The picture is less clear for the outer planets (Jupiter, Saturn, Uranus and Neptune). Again their orbits are practically stable, but it is not known for how long we can accurately predict their positions within those orbits.

Surfing on the edge: chaos versus near-integrability in the system of Jovian planets

We demonstrate that the system of Sun and Jovian planets, integrated for 200 Myr as an isolated five-body system using many sets of initial conditions all within the uncertainty bounds of their currently known positions, can display both chaos and near-integrability. The conclusion is consistent across four different integrators, including several comparisons against integrations utilizing quadruple precision. We demonstrate that the Wisdom‐Holman symplectic map using simple symplectic correctors as implemented in MERCURY 6.2 gives a reliable characterization of the existence of chaos for a particular initial condition only with time-steps less than about 10 d, corresponding to about 400 steps per orbit. We also integrate the canonical DE405 initial condition out to 5 Gyr, and show that it has a Lyapunov time of 200‐400 Myr, opening the remote possibility of accurate prediction of the Jovian planetary positions for 5 Gyr.

Tinkerbell Is Chaotic

Shadowing is a method of backward error analysis that plays a important role in hyperbolic dynamics. In this paper, the shadowing by containment framework is revisited, including a new shadowing theorem. This new theorem has several advantages with respect to existing shadowing theorems: It does not require injectivity or differentiability, and its hypothesis can be easily verified using interval arithmetic. As an application of this new theorem, shadowing by containment is shown to be applicable to infinite length orbits and is used to provide a computer assisted proof of the presence of chaos in the well-known noninjective Tinkerbell map.

Rigorous Shadowing of Numerical Solutions of Ordinary Differential Equations by Containment

An exact trajectory of a dynamical system lying close to a numerical trajectory is called a shadow. We present a general-purpose method for proving the existence of finite-time shadows of numerical ODE integrations of arbitrary dimension in which some measure of hyperbolicity is present and there are either 0 or 1 expanding modes, or 0 or 1 contracting modes. Much of the rigor is provided automatically by interval arithmetic and validated ODE integration software that is freely available. The method is a generalization of a previously published containment process that was applicable only to two-dimensional maps. We extend it to handle maps of arbitrary dimension with the above restrictions, and finally to ODEs. The method involves building n-cubes around each point of the discrete numerical trajectory through which the shadow is guaranteed to pass at appropriate times. The proof consists of two steps: first, the rigorous computational verification of a simple geometric property, which we call the inductive containment property, and second, a simple geometric argument showing that this property implies the existence of a shadow. The computational step is almost entirely automated and easily adaptable to any ODE problem. The method allows for the rescaling of time, which is a necessary ingredient for successfully shadowing ODEs. Finally, the method is local, in the sense that it builds the shadow inductively, requiring information only from the most recent integration step, rather than more global information typical of several other methods. The method produces shadows of comparable length and distance to all currently published results. Finally, we conjecture that the inductive containment property implies the existence of a shadow without restriction on the number of expanding and contracting modes, although proof currently eludes us.

Correct Rounding and a Hybrid Approach to Exact Floating-Point Summation

We present two algorithms for computing correctly rounded sums of arrays of floating-point numbers. First, iFastSum improves upon our previous FastSum by requiring no additional space beyond the original array, which is destroyed. It runs about 20% faster than FastSum in the general case and two times faster when extremely ill-conditioned data are used. The second algorithm is HybridSum, which combines three summation ideas together: splitting the mantissa, radix sorting, and using iFastSum. The result is that when the number of summands is greater than about

SANA: simulated annealing far outperforms many other search algorithms for biological network alignment

10^4