---

Numerical Algorithms for Differential Equations and Dynamical Systems / Algorithmes numériques pour les équations différentielles et les systèmes dynamiques
(Org: Tony Humphries, McGill University)

---

DAVID COTTRELL, McGill University
A backward analysis of simple collisions

Molecular dynamics simulations often involve the numerical integration of pair-wise particle interactions via a constant step-size method. Of primary concern in these simulations is the introduction of errors in velocity statistics. We consider the simple example of symplectic Euler applied to two-particle collisions in 1-d governed by linear restoring force and use backward error analysis to predict these errors. For nearly all choices of system and method parameters, the post-collision energy is not conserved and depends upon the initial conditions of the particles and the step-size of the method. The analysis of individual collisions is extended to predict energy growth in systems of particles in 1-d. Generalization of the analysis to predict evolution of statistics in systems with non-linear interaction forces in higher dimensions will be discussed briefly.

EUSEBIUS DOEDEL, Concordia University, 1455 boulevard de Maisonneuve O., Montreal, QC H3G 1M8
Bifurcation of periodic orbits in the Circular Restricted 3-Body Problem

The Circular Restricted 3-Body Problem describes the motion of a body of negligible mass (a "satellite") in the presence of an Earth-Moon like system (the "primaries"). This problem has been well-studied in the literature, with many important contributions, in view of its application to space-mission design. In this talk I will describe some recent results on the classification of certain types of solutions, namely "elemental" periodic orbits, for all values of the mass-ratio of the primaries. Some of these orbits are of practical interest; for example the Genesis mission, which is now on its way back to Earth for a September arrival, was in a so-called Halo orbit. I will also describe some of the methods used in our numerical work.

MARTIN GANDER, McGill University, 805 Sherbrooke Street West, Montreal, CA
Moving Mesh Methods and Energy Conservation

Moving mesh methods are based on an elegant formulation of partial differential equations coupled with mesh equations. In 1987, a Nobel Prize winner in physics, Tsung Dao Lee, proposed a discrete variant of continuous models from continuum and quantum mechanics, which, instead of using a classical discretization on a fixed mesh, leaves the location of the mesh nodes as unknowns and uses additional equations, like energy conservation, to determine both the grid location and the solution on the corresponding grid. We show for a model problem and its discrete form proposed in the original publication by Lee that imposing energy conservation leads to a method for which the discrete trajectory can cease to exist at perfectly regular points of the underlying continuous physical model. The results presented are joint work with Bob Russel from Simon Fraser University.

WAYNE HAYES, IPST Building 224, University of Maryland, College Park, MD 20742, USA
A Practical Shadowing-based Timestep Criterion for Galaxy Simulations

The "Butterfly Effect", more formally known as "sensitive dependence on initial conditions", is exhibited by many nonlinear dynamical systems from integrated circuits to galaxies. This sensitivity causes small numerical errors to become exponentially magnified, leading some to believe that trajectories of such simulations are the result of nothing but magnified noise. To justify the reliability of such simulations, we turn to the study of "shadowing". A shadow is an exact trajectory that stays close to a numerical trajectory for a long time, even in the face of sensitive dependence. From the standpoint of physics, a numerical trajectory that has a shadow can be viewed as an experimental observation of that shadow, which means that the dynamics observed in the simulation are real. This is a very strong statement of simulation reliability. However, verifying the existence of a shadow formally takes time O(N³), where N is the number of components in the system. In this talk I will outline how I demonstrated the existence of shadows of galaxy simulations in which N=10⁸, without having to perform O(N³) work.

TONY HUMPHRIES, McGill University
Numerics and Dynamics

In a dynamical system defined by an ordinary differential equation it is often the asymptotic behaviour of general initial conditions for large time that is of interest. However traditional numerical analysis typically focuses on the solution of a given initial value problem over a finite time interval, and moreover usually provides error bounds which grow exponentially in time, and so are not useful in the dynamical systems context. Over the last two decades there has been an explosion of work in the "numerics of dynamics" to provide techniques for numerically studying dynamical systems and for providing rigorous meanings for some of the pretty pictures which are beyond the scope of traditional numerical analysis. In this session we will see some of these techniques including direct methods for computing special trajectories (e.g. periodic orbits), backward error analysis, stiffness and adaptive time-stepping. In this first general talk we set the scene by introducing some of the issues that arise, and the techniques required to tackle them.

Along the way we will show that the backward Euler method is a very bad method, and also find one or two other surprises.

NED NEDIALKOV, McMaster University
Solving Differential-Algebraic Equations by Taylor Series

We present a method for solving numerically an initial-value problem differential algebraic equation (DAE). The DAE can be of high-index, fully implicit, and contain derivatives of order higher than one.

We do not reduce a DAE to a first-order, lower-index form: we solve it directly by expanding its solution in Taylor series. To compute Taylor coefficients, we employ the structural analysis of J. Pryce and automatic differentiation.

Generally, our approach succeeds for any DAE whose sparsity structure correctly represents its mathematical structure. We show that a failure occurs if and only if the "system Jacobian" of the DAE is structurally singular up to roundoff, a situation recognizable in practice.

This method has been implemented in a C++ code by N. Nedialkov. Numerical results on several standard test problems show that it is both efficient and accurate.

This is a joint work with J. Pryce, Royal Military College of Science, England.

RAYMOND SPITERI, Dalhousie University, Halifax, NS
A Comparison of Stiffness Detection Methods for Initial-Value ODEs

Stiffness is one of the most enigmatic yet pervasive concepts in the numerical solution of initial-value ODEs. This is especially true for ODEs arising from method-of-lines discretizations of PDEs. The goal of stiffness detection methods is to allow software to automatically choose a stiff or non-stiff integrator, where appropriate. Making the proper choice of integrator often has a profound effect on the overall efficiency of the software. However, the many faces of stiffness have resulted in many definitions, many disagreements, and many methods for its detection. In this talk, I will compare three methods for stiffness detection and show their performance on a representative sample of test problems.