> <\body> ||<\author-affiliation> Laboratoire d'informatique, UMR 7161 CNRS Campus de l'École polytechnique 1, rue Honoré d'Estienne d'Orves Bâtiment Alan Turing, CS35003 91120 Palaiseau |>>|>|> This paper concerns the reliable integration of dynamical systems with a focus on the computation of one specific trajectory for a given initial condition at high precision. We describe several algorithmic tricks which allow for faster parallel computations and better error estimates. We also introduce so called \PLagrange models\Q. These serve a similar purpose as the more classical Taylor models, but we will show that they allow for larger step sizes, especially when the truncation orders get large. ||> Let \\,\,F|]>> and consider the dynamical system <\eqnarray*> >||.>>>> Given an initial condition =I\\> at \>, a target point u> such that is analytic on>, the topic of this paper is to compute >. On actual computers, this problem can only be solved at finite precisions, although the user might request the precision to be as large as needed. One high level way to formalize this is to assume that numbers in the input ( the coefficients of > and , as well as and ) are and to request > again to be avector of computable complex numbers. From a more practical point of view, it is customary to perform the computations using . In our complex setting, we prefer to use a variant, called or . In our main problem, this means that we replace our input coefficients by complex balls, and that the output again to be a vector of balls. Throughout this paper, we assume that the reader is familiar with interval and ball arithmetic. We refer to for basic details on ball arithmetic. The website provides a lot of information on interval analysis. It will be convenient to denote balls using a bold font,

>>. The explicit compact ball with center and radius will be denoted by

>>. Vector notation will also be used systematically. For instance, if \> and

>|)>> with >=\:x\0|}>>, then =,r|)>,\,\,r|)>|)>>>. Sometimes, it is useful to obtain further information about the dependence of the value > on the initial conditions; this means that we are interested in the >, which satisfies the same differential equation() and the initial condition =I>. In particular, the f/\ I> is an important quantity, since it measures the sensitivity of the output on the initial conditions. If > denotes the condition number of , then it will typically be necessary to compute with a precision of at least \> bits in order to obtain any useful output. There is a vast literature on the reliable integration of dynamical systems. Until recently, most work focussed on low precision, allowing for efficient implementations using machine arithmetic. For practical purposes, it is also useful to have higher order information about the flow. are currently the most efficient device for performing this kind of computations. In this paper, we are interested in the time complexity of reliable integration of dynamical systems. We take a more theoretical perspective in which the working precision might become high. We are interested in certifying one particular trajectory, so we do not request any information about the flow beyond the first variation. From this complexity point of view it is important to stress that there is a tradeoff between efficiency and quality: faster algorithms can typically designed if we allow for larger radii in the output. Whenever one of these radii becomes infinite, then we say that the integration method : a ball with infinite radius no longer provides any useful information. Now some \Pradius swell\Q occurs structurally, as soon as the condition number of becomes large. But high quality integration methods should limit all other sources of precision loss. The outline of this paper is as follows: <\enumerate> For problems from reliable analysis it is usually best to perform certifications at the outermost level. In our case, this means that we first compute the entire numeric trajectory with asufficient precision, and only perform the certification at a second stage. We will see that this numeric computation is the only part of the method which is essentially sequential. The problem of certifying a complete trajectory contains a global and a local part. From the global point of view, we need to cut the trajectory in smaller pieces that can be certified by local means, and then devise a method to recombine the local certificates into a global one. For the local certification, we will introduce . As in the case of Taylor models, this approach is based on Taylor series expansions, but the more precise error estimates allow for larger time steps. The first idea has been applied to many problems from reliable computation (it is for instance known as Hansen's method in the case of matrix inversion). Nevertheless, we think that progress is often possible by applying this simple idea even more systematically. In Section

, we will briefly recall some facts about the efficient numeric integration of(). In Section

we discuss the way in which the problem of certying a global trajectory can be reduced to more local problems. It is interesting to compare this approach with the more classical stepwise certification scheme, along with the numerical integration itself. The problem with stepwise schemes is that they are essentially sequential and thereby give rise to a linear precision loss in the number of steps. The global approach reduces this to a logarithmic precision loss only. The global strategy already pays off in the linear case and it is possible to reorganize the computations in such a way that it can be re-incorporated into an interative scheme. Our current presentation has the merit of conceptual simplicity and ease of implementation. The main contribution of this paper concerns the introduction of \PLagrange models\Q and the way that such models allow for larger time steps. Classical Taylor models approximate analytic functions on the compact unit disk (say) by a polynomial \> of degree n> and an error \0> with the property that -P|\|>\\> for all \1>. The idea behind Lagrange models is to give amore precise meaning to the \Pbig Oh\Q in =f+\+f*z+O|)>>. More precisely, it consists of a polynomial \

>> of degree n> with ball coefficients and an \0> such that \\+\|)>*z>. The advantage comes from the fact that the integration operator has norm for general analytic functions on the unit disk but only norm > for analytic functions that are divisible by >. Although Lagrange model arithmetic can be a constant times more expensive than Taylor model arithmetic, the more precise error estimates allow us to increase the time step. An earlier (non refereed) version of the material from Section

appeared in the lecture notes. The algorithm from Section

has been implemented in the package of the system, both for Taylor models and Lagrange models. For high precision computations, we indeed observed improved step sizes in the case of Lagrange models.

Since the main focus of this paper is on certification, we will content ourselves with a quick survey of the most significant facts about numerical integration schemes.

From a high level perspective, the integration problem between two times z> can be decomposed into two parts: finding suitable intermediate points \z\\\z\z=z> and the actual computation of |)>> as afunction of |)>> for ,s> (or as a function of |)>,\,f|)>> for some schemes). The optimal choice of intermediate points is determined by the distance to the closest singularities in the complex plane as well as the efficiency of the scheme that computes |)>> as afunction of |)>>. For >, let > be the convergence radius of the function at . High order integration schemes will enable us to take -z|\|>\c*\

|)>> for some fixed constant 0>. Lower order schemes may force us to take smaller steps, but perform individual steps more efficiently. In some cases ( when admits many singularities just above the real axis, but none below), it may also be interesting to allow for intermediate points ,\,z> in > that keep alarger distance with the singularities of . For small working precisions, Runge-Kutta methods provide the most efficient schemes for numerical integration. For instance, the best Runge-Kutta method of order requires evaluations of for each step. For somewhat larger working precisions ( quadruple precision), higher order methods may be required in order to produce accurate results. One first alternative is to use relaxed power series computations which involve an overhead > for small orders and n> for large orders. For very large orders, a power series analogue of Newton's method provides the most efficient numerical integration method . This method actually computes the first variation of the solution along with the solution itself, which is very useful for the purpose of this paper.

Another interesting question concerns the amount of computations that can be done in parallel. In principle, the integration process is essentially sequential (apart from some parallelism which may be possible at each individual step). Nevertheless, given a full numeric solution at asufficiently large precision , we claim that a solution at a roughly doubled solution can be computed in parallel. More precisely, for each , and , let > be the solution of() with =I>, and denote = f/\I|)>>. We regard > as the \Ptransitional flow\Q between and, assuming the initial condition at . Notice that =Id> and, for v\z>, <\eqnarray*> >|||)>>>|>|||)>*V.>>>> Now assume that we are given \f|)>> for ,s> and at precision . Then we may compute \V,z,f|)>|)>> in parallel at precision . Using adichotomic procedure of depth|log s|\>>, we will compute \f|)>> at precision in parallel for ,s>, together with \V,z,f|)>|)>> at precision . Assume that 2> and let |s/2|\>>. We start with the recursive computation of \f|)>> and \f,z,f|)>> at precision for ,m> ( ,s-m>), together with \V,z,f|)>|)>> and \V,z,f|)>|)>> at precision . Setting =f-f>, we have <\eqnarray*> |)>>|>|,z,f+\|)>>>||>|+V*\>>>> at precision (for ,s-m>) and <\eqnarray*> ,z,f|)>|)>>|>|*V>>>> at precision . It thus suffices to take \f+V*\> and \V*V>. The above algorithm suggests an interesting practical strategy for the integration of dynamical systems on massively parallel computers: the fastest processor(s) in the system plays the rôle of a\Pspearhead\Q and performs a low precision integration at top speed. The remaining processors are used for enhancing the precision as soon as a rough initial guess of the trajectory is known. The spearhead occasionally may have to redo some computations whenever the initial guess drifts too far away from the actual solution. The remaining processors might also compute other types of \Penhancements\Q, such as the first and higher order variations, or certifications of the trajectory. Nevertheless, the main bottleneck on a massively parallel computer seems to be the spearhead.

A of the dynamical system() can be defined to be a ball function <\equation*> \:|)>\\

|)> with the property that \\

|)>> for any z> and \>. An additionally requires a ball function <\equation*> \:|)>\\

|)> with the property that \\

|)>> for any z> and \>. A l of() is a special kind of certified integrator which only produces meaningful results if and are sufficiently close (and in particular \\>). In other words, we allow the radii of the entries of

|)>> to become infinite whenever this is not the case. Extended local certified integrators are defined similarly. One interesting problem is how to produce global (extended) certified integrators out of local ones. The most naive strategy for doing this goes as follows. Assume that we are given a local certified integrator >, as well as z> and >. If the radii of the entries of

|)>> are \Psufficiently small\Q (finite, for instance, but we might require more precise answers), then we define

|)>\\

|)>>. Otherwise, we take /2> and define