LIBRARY Mlchlgan State University This is to certify that the dissertation entitled Rigorous Nunerical Analysis With High—Order Taylor Models presented by Jens Hoefkens has been accepted towards fulfillment of the requirements for Ph . D degree in Physics / Mathemat ics it Mu }W Major professor M. Berz Date jl'll ! Ol MSU is an Affirmative Action/Equal Opportunity Institution 0- 12771 PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 6/01 cJCIRC/DateOuepes-pJS RIGOROUS NUMERICAL ANALYSIS WITH HIGH-ORDER TAYLOR MODELS By Jens Hoefkens AN ABSTRACT OF A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics and Department of Physics and Astronomy 2001 Professor Martin Berz IV ed. '9 v: ’ I" 'r “‘4 ' ll 1." ‘n L, _ .. . ._ 1‘ l‘ V ‘ ALI r ' n .. I". ABSTRACT RIGOROUS NUMERICAL ANALYSIS WITH HIGH-ORDER TAYLOR MODELS By Jens Hoefltens Interval techniques have been utilized for rigorous numerical analysis since the 19608. However, their use has been limited by the dimensionality curse and the dependency problem. The recently developed Taylor model approach alleviates these problems and allows the use of validated numerics in a wide range of applications. To broaden the applicability of the Taylor model method, we introduce new algorithms for the inversion of functional relations and the integration of differential algebraic equations. First we present a new method for computing verified enclosures of the inverses of given functions over large domains. The approach utilizes Taylor models and the sharpness of the enclosures scales with a high order of the domain. An integral part of the new method is the rigorous determination of invertibility of high dimensional functions over large domains, which is reduced to a verified linear algebra problem in- volving only first derivatives of the function of interest. Several examples highlighting various aspects of the methods are discussed. Differential algebraic equations (DAEs) describe important problems in mechani- cal and chemical engineering. Existing algorithms for the integration of DAE initial Value problems have traditionally been restricted to low-index systems and until re- cently, no practical scheme for the verified integration of DAEs existed. Recognizing the antiderivation as a natural operation on Taylor models yields a method that treats DAEs within a fully differential algebraic context as implicit equations made of con- ventional functions and the antiderivation. The resulting integration scheme can be applied to high-index problems and allows the computation of guaranteed enclosures of final coordinates from large initial regions. To demonstrate the general applicability of the Taylor model approach, we present results from verified asteroid orbit integrations and the theory of Hamiltonian systems. We show that the newly developed methods are practical and can indeed outperform conventional interval methods in a wide class of problems. Finally, we discuss some details of implementing interval libraries on general purpose computers and present the concept of language independent software development, which has been used for the design and implementation of the C++ and Fortran 90 interfaces to COSY Infinity. To $00, for her love. iv ..,, LIT, 1 A ~.") - «Wu. 54.. q. . l l A ACKNOWLEDGMENTS Most of all, I would like to thank my parents and my sister for their unconditional love. I am especially thankful to them for supporting me in my decision to leave Germany for the United States. I want them to know that I understand and appreciate their efforts and sacrifices. I also would like to thank Soo Chang for her never ending love, faith, and support over the last four years. Whenever life seems dull, she opens my eyes to the beauty of the world. The field of interval analysis is small but has a close and friendly community of great scientists. I would like to thank the following colleagues for helpful advice and fruitful discussions over the years: Christian Bischof, George Corliss, Paul Hovland, Ken Jackson, Vladik Kreinovich, Rudolf Lohner, Ned Nedialkov, John Pryce, and Bill Walster. Special thanks goes to Ramon Moore for his ground breaking work on intervals. Working with him on the asteroid problem was one of the best experiences I had during my time in graduate school. However, my deepest thanks and gratitude 80130 my thesis advisor Martin Berz. He is a constant source of guidance and wisdom. He 0Pened my eyes to the rigorous analysis of Taylor models and deserves my thanks for giVing me the opportunity to study at MSU. I also would like to thank those professors who have kindly served on my guid- ance committee: John Masterson, Jerry Nolen, Bradley Sherrill, Daniel Stump, and Clifford Weil. I have also enjoyed many discussions with friends and fellow students over the course of the last five years. In particular I would like to mention Jens von Bergmann, Sen Cheng, Jennifer Church, Lars Diening, Bela Erdélyi, Holger Har- reis, Michael Lindemann, Kyoko Makino, Declan Mulhall, Cristian Opazo, Khodr Shamseddine, and Ralf Tiinjes. For financial support of my research I owe thanks to the Studienstiftung des deutschen Volkes, the US Department of Energy, the Alfred P. Sloan Foundation, and the National Science Foundation. I am grateful for their continuous support of young researchers and basic science. Finally, special thanks goes to my father for making me stay in school when I was ready to graduate with a Masters degree. Without him, this dissertation would have never been written. vi C01 llil l lEll Contents LIST OF TABLES x LIST OF FIGURES xii LIST OF ALGORITHMS xv 1 Introduction 1 2 Background Information 4 2.1 The Differential Algebra "D1, ....................... 5 2.1.1 Elementary Operations ...................... 6 2.1.2 Contracting Operators and Fixed Points ............ 7 2.1.3 Functions and Antiderivation .................. 8 2.1.4 Summary ............................. 10 2.2 Interval Arithmetic ............................ 11 2.2.1 Elementary Operations and Functions ............. 13 2.2.2 Interval Arithmetic and Set Theory ............... 16 2.2.3 Applications of Interval Arithmetic ............... 18 2.2.4 Summary ............................. 22 2.3 Taylor Models ............................... 22 2.3.1 Arithmetic and Operations .................... 26 2.3.2 Accuracy of Taylor Models .................... 29 2.3.3 Applications of Taylor Models .................. 32 2.3.4 Summary ............................. 34 3 Verified High-Order Inversion 35 3.1 Rigorous Invertibility Analysis ...................... 36 3.1.1 lnvertibility from First Derivatives ............... 37 3.1.2 Verifying Invertibility with Taylor Models ........... 41 vii 3.1.3 Comparison of Invertibility Criteria ............... 46 3.2 Guaranteed Enclosures of Inverse Functions .............. 54 3.2.1 Polynomial Inverses ........................ 54 3.2.2 . Inverse Taylor Models ...................... 56 3.3 Examples of Taylor Model Inversion ................... 59 3.3.1 One-Dimensional Function .................... 60 3.3.2 Six-Dimensional Function .................... 61 3.4 A Superconvergent Newton Method ................... 63 3.4.1 Examples ............................. 66 3.5 Implicit Equations ............................ 71 3.5.1 Curves, Surfaces, Manifolds ................... 73 3.6 Analysis of Inversion Algorithms ..................... 75 3.6.1 Space and Time Requirements .................. 76 3.6.2 Order Dependence ........................ 77 3.7 Summary ................................. 78 Differential Algebraic Equations 81 4.1 Background and Motivation ....................... 82 4.1.1 First Order Differential Algebraic Equations .......... 83 4.1.2 General Differential Equations .................. 87 4.2 Verified Integration of Implicit ODEs .................. 88 4.2.1 High Order Taylor Model Solutions ............... 90 4.2.2 Mathematical Background .................... 92 4.2.3 Discussion of the Algorithm ................... 97 4.2.4 A Second Order Example .................... 101 4.3 Verified Integration of DAEs ....................... 103 4.3.1 Structural Analysis ........................ 105 4.3.2 Verified Index Analysis ...................... 109 4.4 Examples and Applications ....................... 110 4.4.1 Kirchhoff’s Law and Electric Circuits .............. 110 4.4.2 Constrained Mechanical Systems ................ 113 4.4.3 Optimal Prescribed Path Control ................ 118 4.5 Summary ................................. 120 5 Applications of High Order Taylor Model Methods 122 5.1 Solar System Dynamics and Asteroid Orbits .............. 123 5.1.1 Physical Background ....................... 124 viii Drb. Judi 5.1.2 5.1.3 5.1.4 5.1.5 5.1.6 Geometric Description of Kepler Orbits ............. Verified Orbit Integration .................... Results and Discussion ...................... Comparison with AWA ...................... Summary ............................. 5.2 Existence of Generating Functions .................... 5.2.1 5.2.2 5.2.3 5.2.4 5.2.5 5.2.6 Introduction ............................ Theory of Extended Generating Functions ........... Enclosing Derivatives of Flows .................. Examples ............................. Rigorous Analysis of Symplectification Errors ......... Summary ............................. 6 Implementational Details 6.1 Implementation of Portable Interval Libraries ............. 6.1.1 6.1.2 6.1.3 6.1.4 6.1.5 Floating Point Numbers ..................... Directed Rounding ........................ Implementation of Interval Operations ............. Benchmarks and Results ..................... Summary ............................. 6.2 Language Independent Programming .................. 6.2.1 6.2.2 6.2.3 6.2.4 6.2.5 The Least Common Denominator Approach .......... The C++ Interface to COSY Infinity .............. The Fortran 90/95 Interface to COSY Infinity ......... Performance Analysis ....................... Summary ............................. A Orbital Elements of Major Planets B The COSY C++ Interface B.1 Available COSY Functions ........................ B.2 Available COSY Procedures ....................... C The COSY Fortran 90 / 95 Interface BIBLIOGRAPHY 131 136 143 153 158 158 159 160 162 164 175 176 176 179 185 187 189 190 190 193 202 210 211 213 215 215 216 219 223 I‘HW ‘ 1.1.4 All.) in“. List of Tables 3.1 3.2 3.3 3.4 3.5 3.6 4.1 5.1 5.2 5.3 5.4 5.5 5.6 Taylor model of one-dimensional sine function computed with COSY Infinity. Shown are Taylor coefficients, reference point, domain infor- mation, and remainder bound ....................... Left-inverse Taylor model for the Taylor model shown in Table 3.1. Shown are Taylor coeflicients, reference point, domain information, and remainder bound. Remainder Bounds of the left-inverse Taylor model for six-dimensional exponential function given by Eqns. (3.31) and (3.33) .......... Comparison of different Newton methods for the determination of the fixed point do of Eqn. (3.35) ........................ Enclosures of 7r, obtained by the interval Newton and Taylor model methods ................................... Enclosures of the zero for the six-dimensional exponential function, computed with an interval Newton method and the high order Taylor model approach ............................... Taylor model for the solution of the implicit second order ODE initial value problem given by Eqn. (4.44) .................... SI units for length, time, and mass. ................... Astronomical units and their conversions to SI units. Units of the derived quantities velocity, acceleration, and force. . . . . Masses, mean distances, and periods of the major planets. Masses are given as fractions of the Sun mass M, mean distances in astronomical units, and times in Julian years [125] ................... Geometric shape of conic sections, depending on the eccentricity 5. Dates and minimal distances of the predicted closest Earth approach distances for the near-Earth asteroid 1997 XFll between the years 2000 and 2040 [73]. ............................ 5.7 Taylor model describing final 2: positions as a function of the initial conditions 9:0, 00, yo, b0 for the accelerator cell of Sec. 5.2.4. . 6.1 Average execution times for the three interval test cases in pure For- tran 77 environments ............................ 60 61 63 67 70 71 104 125 126 126 129 132 148 171 188 1.3 6.2 Average execution times for the three interval test cases in COSY In- finity language environments. ...................... 6.3 Execution times in minutes and seconds for a typical mix of DA oper- ations on different platforms with different interfaces to COSY Infinity. 210 Al Orbital Elements of the major planets (angles in degree, distances in astronomical units) at the epoch J 2000 .................. 213 A2 Daily rates of change of the orbital elements {2, z', and w of the nine major planets . .............................. 214 A3 Daily rates of change of the orbital elements a, a, and M of the nine major planets ................................ 214 0.1 Defined combinations of COSY objects and standard data types for the addition ................................. 219 C2 Defined combinations of COSY objects and standard data types for the subtraction ............................... 219 C3 Defined combinations of COSY objects and standard data types for the multiplication. ............................ 220 0.4 Defined combinations of COSY objects and standard data types for the division ................................. 220 C5 Defined combinations of COSY objects and standard data types for the power operation. ........................... 220 0.6 Defined combinations of COSY objects and standard data types for the comparison .LT ............................. 220 C7 Defined combinations of COSY objects and standard data types for the comparison .GT ............................. 221 C8 Defined combinations of COSY objects and standard data types for the comparison .EQ ............................. 221 C9 Defined combinations of COSY objects and standard data types for 221 the comparison .NE ............................. C.10 Defined combinations of COSY objects and standard data types for the COSY concatenation .UN .. ..................... 221 H In. List of Figures 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 Commuting diagram for the elementary operations on "D” ....... 7 Commuting diagram for the elementary operations on intervals. . . . 14 Illustration of the wrapping effect caused by a rotation of 45 degrees in the plane ................................. 17 Illustration of how the wrapping effect can lead to arbitrary large over- estimations. ................................ 18 Taylor models of order one and five for the sine function over the domain interval [—1.5, 1.5] ......................... 25 Left: Enclosing a function by a Taylor model of order eight; right: Interval bounding of the same function with 50 sub—intervals ...... 25 Commuting diagram for elementary operations on Taylor models. . . 27 Commuting diagram for sine function on Taylor models. ....... 28 Commuting diagram for antiderivation of Taylor models. ....... 29 Interval bounding of a two-dimensional function by Taylor models of orders seven, eight, nine, and ten. .................... 31 Percentage of random functions that can be shown to be invertible as a function of dimensionality ........................ 47 Percentage of random functions that can be shown to be invertible as a function of domain size. ........................ 49 Percentage of functions that can be shown to be invertible as a function of the non-linearity 5 ............................ 50 Percentage of maximal domain size over which invertibility can be proven as a function of complexity and non-linearity in the partial derivatives. ................................ 52 Percentage of random functions that can be shown to be invertible as a function of functional complexity. ................... 53 Illustration of the domains in the definitions of inverse Taylor models. 59 Difference between the arcsine function and the reference polynomial of the 19—th order Taylor model shown in Table 3.2 ........... 62 Transformed Fixed-Point problem: the fixed point an of Eqn. (3.35) is the intersection point of the two graphs. ................ 67 xii 3.9 Taylor polynomial approximation of order 25 of the sine function over the interval [0, 12]. IIIIIIIIIIIIIIIIIIIIIIIIIIII 3.10 Charts for the cartesian coordinates x and y of the two-parameter 3.11 3.12 4.1 4.2 4.3 4.4 4.5 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 description Eqns. (3.44) for hyperbolas. Remainder bounds are in the order of 10‘s. ............................... Illustration of verified inversion with intervals. The sharpness of the 75 enclosures scales linearly with the magnitude of the domain sub-intervals. 77 Double-logarithmic plot of the accuracy of the first Taylor model N ew- ton step as a function of initial domain size; results for 6—th order Taylor models. .............................. Illustration of an electric circuit described by the index-1 DAE (4.69). Ub is the operating voltage and U8 is the driving voltage of the circuit. Illustration of a double pendulum with masses m1 and m2 connected by massless rods of lengths l1 and I2, respectively. ........... Coordinates 2:1 and y1 of the double pendulum for 901(0) = 902(0) = 5°,30°,90°, 273(0) = 235(0) = 311(0) = 315(0) = 0 and 0 S t _<_ 100 (integrated in order 6, h = 0.1, g = 1, l1 =12 =1,m1= m2 =1). . . . .732- and yz-phase space projections for the double pendulum with for 901(0) = ¢2(0) = 5°,30°,45°, 373(0) = $610) = 34(0) = 215(0) = 0 and 0 S t g 100 (integrated in order 6, h = 0.1, g = 1, 11 = l2 2 1, m1 = m2 = 1) ................................ Coordinates 2:1 and y1 of the double pendulum for 901(0) = 0, (p2(0) = 45°, cc’1(0) = 333(0) 2 y’1(0) = y§(0) = 0 and 0 S t S 100 (integrated in order 6, h = 0.1, g = 1,l1=l2 =1,m1=100, m2 =1). ........ Orbits of the inner planets Mercury, Venus, Earth, and Mars around the Sun. .................................. Orbits of the outer planets Jupiter, Saturn, Uranus, Neptune, and Pluto around the Sun. .......................... Illustration of the eccentric anomaly E ................. Residuals of Jupiter; two Very Long Baseline Interferometry observa- tions from the Galileo spacecraft; courtesy of [130]. .......... Cartesian positions (in AU) and velocities (in AU /TU) of 1997 XF11 during the ten year integration interval .................. Logarithmic plot of the diameters of enclosures for the positions (in AU) and velocities (in AU/TU) of 1997 XF11 during the ten year integration interval ............................. Evolution of the step size and one-step integration errors during the ten year integration of 1997 XF11. ................... Differences in the computed positions and velocities of the asteroid 1997 XF11 between the ten year integrations with and without rela- tivistic corrections. ............................ xiii 78 111 114 116 117 117 128 129 134 139 146 147 5.9 Illustration of the shrink wrapping method used in the verified inte- gration of asteroid orbits .......................... 5.10 One-step and total shrink factors obtained during the ten year integra- tion of the asteroid 1997 XF11. ..................... 5.11 Logarithmic plot of the diameters of the enclosures for positions ob— tained by AWA (upper) and the Taylor model based integrator. 5.12 Logarithmic plot of the diameters of the enclosures for positions ob- tained by AWA for an initial box reduced by a factor of 3.275 in each direction. ................................. 5.13 Logarithmic plot of the diameters of the enclosures for positions ob— tained by AWA for initial condition enclosures of 10‘15 ......... 5.14 ’Iiacking picture of the cubic two—dimensional symplectic map, and the box of guaranteed invertibility of the generator associated with S = 0. 5.15 Illustration of the Fermi-Pasta-Ulam system with three masses and four springs ................................. 5.16 (q1,p1) and (q1,q2) tracking pictures of the half period map of the Fermi-Pasta-Ulam system for a particle launched along the ql axis. The box of guaranteed invertibility of the generator associated with S = 0 extends to at least [—1, 1} in every direction. .......... 5.17 Illustration of the accelerator cell studied in Sec. 5.2.4 ......... 5.18 (x, a) and (y, b) tracking pictures of the one-turn map of the accelerator cell for particles launched along the spatial :r and y axes, respectively. The box of guaranteed invertibility of the generator associated with S = 0 extends to at least [—0.1, 0.1] in the spatial variables. ..... xiv 152 153 155 156 157 166 167 168 170 List of Algorithms 6.1 Interval addition with hardware support to set the rounding mode. . 180 6.2 Directed rounding of the two interval endpoints of [$1,132]. ...... 184 6.3 Portable determination of the machine constant 5 at run time. . . . . 184 6.4 Portable determination of the machine constant p at run time ..... 184 6.5 Reliable interval extension of the sine function .............. 186 XV H. Chapter 1 Introduction Over the last four decades, the use of computer programs has gained widespread acceptance in almost all scientific and engineering disciplines. An implicit trust in the computed results is based on the assumptions that computer programs are imple- mented correctly and perform as intended. However, the validity of these assumptions is questionable, to say the least. Considering the fallibility of humans, it is presumptuous to believe that software could be produced free of errors. Moreover, since verifying the correctness of imple- mentations is equivalent to the halting problem; it is generally impossible to prove that a given computer program does not contain bugs [117]. While elaborate testing procedures can reduce the number of errors, even the most important and extensively tested software systems can fail. As an example, consider the explosion of an Ariane 5 rocket in 1996 which has been traced to a very simple programming error [77]. However, even assuming that software could be written free of bugs, other sources of errors still exist. Most of these errors are caused by the transition from the perfect world of numerical analysis to the limited world of finite state machines. While numerical analysis often assumes the existence of infinitely many real numbers and infinitely precise computations, computers have only finite memory and floating point numbers are stored with only finitely many digits, resulting in computations that are only approximations to the mathematically correct results. Most modern computers implement double precision floating point numbers with approximately 16 significant digits. While this is sufficiently accurate in most ap- plications, the resulting rounding and truncation errors can accumulate and lead to arbitrary large deviations from the mathematically correct results. Moreover, while humans prefer decimal numbers, computers usually store information in binary repre- sentations [48]. The conversion from decimal numbers to binary numbers, combined with the limited accuracy of the number representations, can be a significant source of truncation errors. To illustrate this, consider the infinite binary representation of 0.1: 0.110 = 0.001 10012. (1.1) On a computer with 24 bit floating point numbers, the error caused by the truncation amounts to an absolute error of approximately 9.5 - 10'8. While this seems insignifi- cant, the high failure rates of the Patriot missile system during the Gulf War in 1991 have been traced to exactly this truncation error [141]. It is important to note that the problem of roundoff and truncation errors is a fundamental characteristic of finite state computers. While intimate knowledge about the underlying architecture and significant efforts in the implementation of computer programs could reduce the effects of these errors, such errors can never be completely eliminated. Moreover, the extra effort necessary would require a software development model contradicting the idea of using portable high-level programming languages that help the programmer focus on the correctness of the implementation. Thus, on the one hand we desire the portable programming of general purpose, limited precision binary architectures, on the other hand we know from theory and s»;- .1]. A 'I .g" .| .z‘ 0 . w. ; . *u \ l. . A“ V. . ‘3 . .1 h 's . practice that this combination can lead to significant computational errors. Nev- ertheless, computers are here to stay and their use is likely to increase. Moreover, computational results are increasingly used as the basis of important decisions, rang- ing from the operation of nuclear power plants to the world financial markets. In other words, potentially erroneous results are used to control systems of global importance. Somewhat contrary to naive intuition, the main problem with these numerical inaccuracies is that most computed results are actually sufficiently good approxima— tions of the mathematically correct results and only a very small number of cases exhibit significant errors that require further investigation. It is the goal of validated methods to solve this dilemma by providing the users with bounds that are guaran- teed to contain the mathematically correct result. While validated methods do not magically increase the accuracy of computations, they provide a self-validating ap- proach that computes both an approximation to the mathematically correct result and a rigorous upper bound on the error of the approximation. In most situations these are relatively tight bounds and only in a few cases will they be large, indicating that the computed results may be rather bad approximations of the correct results. While the fundamental ideas behind validated methods have already been intro— duced in the dissertation of R. Moore in 1962 [91], early methods have often computed overly large bounds. Only recently, starting with work by R. Lohner in 1987 [78], have validated methods been able to produce sharper and more usable bounds. In this dis- sertation we present new results from the field of interval methods [92], which lie at the core of the validated approach. The results are based on the theory and the ap- plication of the recently developed Taylor models [83, 82], which combine high order Taylor polynomials with intervals for validation. Taylor models can in some situations increase the sharpness of the computed bounds by several orders of magnitude over conventional methods, thereby improving the applicability of validated methods. Chapter 2 Background Information In this chapter we introduce the notational conventions that will be used throughout this dissertation. We also summarize the mathematical and computational theories that form the backbone of the material developed in this dissertation: the differential algebra flD,“ interval analysis, and Taylor models. Throughout this dissertations N denotes the set of positive integers, Q is the set of rational numbers, IR stands for the set of real numbers, and C is the set of complex numbers. For notational convenience, we will always assume that n, v, w E N. Following the usual conventions, IR” denotes the set of v-dimensional vectors with real entries; elements of R” are written as :1: = ($1,...,:1:.,), (2.1) with real numbers 2:,- 6 IR for z' = 1,. . .,v. If necessary, we denote vectors by if to distinguish them from real numbers 2:. However, in most cases, such a distinction will not be necessary. Finally, unless stated otherwise, all functions are assumed to be at least (n + 1)-times continuously differentiable over their domain. “1'31.“ I a. ‘-.r A“ v ‘: u . r- -‘t . ‘L '.l 1.1) l '9 .ll .1“ gr. I.‘ ix‘I ‘L‘~ .‘l 2.1 The Differential Algebra n.0,, Many of the results presented in the following chapters rest on the availability of accurate descriptions of high-order Taylor polynomials of sufficiently smooth functions on computers. To avoid the need for storing an infinite number of Taylor coefficients, it is convenient to consider the equivalence classes of functions with the same n—th order Taylor polynomials. Since these are finite objects, the resulting structures can be implemented on computers, and they can be used for the efficient and accurate modeling and representation of complicated functions. Let U C R” be an open set containing the origin and consider the space C "+1 (U, R”) of (n +1)-times continuously differentiable functions that map U into IR”. We define the relation of equality up to order n as follows. Definition 2.1. For f,g E C"+1(U,R"’) we say that f equals 9 up to order n if f(0) = 9(0), and all partial derivatives of orders up to n agree at the origin. If f equals 9 up to order n, we denote that by f =n g. It is easy to see that equality up to order n establishes an equivalence relation on the space C"+1(U,R“’) [17]. The resulting equivalence classes are called DA vectors, and the class containing the function f E C"+1(U,R"’) is denoted by [f],,. The collection of these equivalence classes is called ”B”. More details on this structure are given in [12, 17]. PropOSItion 2.1. For f E C"+I(U,R"’), the n-th order Taylor polynomial Tn( f) of f is contained in [f]n. This assertion follows easily from the basic definition of the equivalence classes. However, the fact that the n—th order Taylor polynomial of f can be used as a rep- 9 U» .“OH. 4 .J. C’ ‘I ll ~'l. .2 . .2an 1.1.1 4 uni l 1‘ K. 1.1 \ s a..." I‘M. resentative for the class [f],, opens the door for a computer implementation of the structure "D, by storing and manipulating the coefficients Of Taylor polynomials. 2.1.1 Elementary Operations Elementary operations like “+” and “x” can be lifted from C”+1(U, R”) in the usual way, and extend to the corresponding Operations “EB” and “(8” on "D, [12, 13]. Definition 2.2. Let f, g E C"+1(U,R‘”) be two functions. Then the sum of the DA vectors [flu and [g]n is given by lfln 69 [9]" = if + gln- (2.2) The product of the two DA vectors is defined by lfln ® [gin = if X Slim (23) Together with the scalar multiplication r - [f],, = [r - f],,, this definition of the elementary operations makes “D, an algebra [17]. Moreover, the diagram in Fig. 2.1 is closed and commuting. In other words, the extension of the elementary operations to "D, is transparent to the equivalent classes of DA vectors; i.e., knowledge of the values and derivatives of f and g at the origin is sufficient to obtain Taylor polynomials of their sums and products. Moreover, in [12, 13] the available Operations on "D, have been extended to include subtraction and, for a limited class of DA vectors, even the multiplicative inversion. From now on we omit the distinction between the Operations on C"+1(U,1R"’) and "D, and will always use the same symbols “+” and “X” for operations between numbers, functions, and DA vectors. It has been shown that the derivative Operation (derivation) can be extended from C"+1(U,R“’) to the algebra "D, in such a way that ”D, becomes a difierential algebra [12, 17]. While we will not use this intrinsic structure Of "D”, we will make frequent use of the antiderivation of DA vectors to be presented later. 11‘) fig 4 lfln’lgl" (+v ) (63, 6) ii i f(+. -)g a [flue ©>lgln Figure 2.1: Commuting diagram for the elementary Operations on "D”. 2.1.2 Contracting Operators and Fixed Points In the previous section we have shown how elementary operations on the function space C"+1 (U, R‘” ) can be extended to the differential algebra ”Du. Here we take a look at the more general concept of operators on “D, and summarize an important fixed point theorem. The availability of this powerful fixed point theorem for operators on nDu allows the use of the differential algebra ”D, in a large class of numerical applications, ranging from the analysis of dynamical systems [15, 23, 27, 18] to global optimization [68, 88]. Definition 2.3. For [f],, 6 "D”, the depth A ([f],,) is defined to be the order of the first non-vanishing derivative of f if [f],, 76 0, and n + 1 otherwise. By definition of the equivalence classes, this definition is independent of the choice of f E [f]. We note that any a E “D, with A(a) 2 1 satisfies the condition a"+1 = 0 and is therefore called nilpotent. Using the straightforward definition of the depth, contracting operatOrs on ”D, are defined as follows. Definition 2.4. Let 0 be an operator defined on M C "D”. O is contracting on M, ’31] 'l ‘v ,I.‘l '1. h I".. ‘. 5 | ,I‘IP I . iffor any two [f],,, [g]n E M, ’\(O(ifln) '- 0([gln)) 2 ’\(lfln — iglnl (2'4) with equality if and only if f 2,, g. This definition has a striking similarity to the corresponding definitions on stan- dard function spaces. Even more so, a theorem that resembles the Banach Fixed Point Theorem can be established on "D”. However, unlike in the case of the Banach Fixed Point Theorem, in “D, the sequence of iterates is guaranteed to converge in at most n + 1 steps [17]. Theorem 2.1 (DA Fixed Point Theorem). Let 0 be a contracting operator and self-map defined on M C "D”. Then 0 has a unique fixed point a E M. Moreover, for any a0 6 M it is 00‘“) (a0) = a. A proof and further discussion of the DA Fixed Point Theorem can be found in [17]. Here we just mention that, since )1(a + b) 2 min()\(a), A(b)), it follows easily that the sum and composition of two contracting Operators 01 and 02 defined on M is also a contracting Operator. 2.1.3 Functions and Antiderivation To fully utilize the differential algebra "Du, especially in numerical analysis and com- puter environments, it is necessary to not only define the elementary operations on an but also the standard mathematical functions commonly available on computers: square root, exponential, logarithm, trigonometric and hyperbolic functions. At a fundamental level, the basic functions on ”D" are simply defined in terms of the corresponding operations on the function space [17]. Definition 2.5. For f,g E C"+1(U,R"’) we define 9(lfln) = [9(fllrv Y . . r. .l. 1 X. L . l. Tru \ . v .l\t n\J “I 5 ”MI. I v m. 1 I 2.: ri , . u . 1 on; lulu ..w. 36 “.11 it '1. 'l n. .2; 7? .31 an.“ . u. o TA . .. .. .1. u . V .n . . . p. .ru. 1 i 1".“ N While this definition is straightforward and of great theoretical value, the actual computation of functions on "D” is based on addition theorems and the Taylor series of these functions. Although the full details of this procedure are beyond the scope of this summary, the following example illustrates the general approach. Let [f ]n 6 "D, be a DA vector and write [f],, = a0 + b with a0 = f (0) Then it is en) (lfln) = exp(ao + b) = exp(ao) exp(b)- (2-5) Using the definition of the exponential function as a power series and the fact that b is nilpotent, in fact mm b" exp ([fln) = exp(ao)' Z 75'" (2-6) k=0 ’ Further details on the implementation of functions of DA vectors can be found in [12, 13, 17]. The approach outlined above allows us for any function f, for which we have a code list or algorithm consisting of finitely many intrinsic functions and elementary Operations, to obtain the n—th order Taylor polynomial of f around the origin. By starting the evaluation with the identity DA vectors [id],,, we obtain Tn( f ) by evalu- ating the code list of f with the argument [id],,. This gives a convenient and powerful method of computing derivatives Of functions described by computer programs. Un- like conventional automatic differentiation [111, 52], which is often limited to first and second order, the Taylor series approach does not pose any arbitrary limits on the maximum order of derivatives that can be computed. We conclude this section on functions on "D, with an example of an operator that is unusual but, considering the structure of the differential algebra "D”, actually quite natural. The antiderivation; i.e., the integration with respect to any Of the v variables, turns out to be a contracting Operation on all;- Papa: v' ' ”mint“ “fill CL ... L '3 . .1 .‘q r:- , 4] 1L] '. éilir‘ . .5“.]' :1. . . 3L] 1 Proposition 2.2 (Antiderivation is Contracting). For k e {1, . . .,v}, the an- tiderivation 6,:1 : nD.,—>,,D,, is a contracting operator on "D”. The proof of this important result is based on the fact that if a, b E "D, agree up to order l, the first non-vanishing derivative of 6,:1(a - b) is Of order l + 1 [17, 59, 60]. It is important to realize that in the DA framework of "Du there is no fundamental difference between any of the standard mathematical functions like the sine and exponential functions on "D” and antiderivation. In fact, fully embracing antiderivation as a normal operation on DA vectors will enable us to develop a new and powerful method for the verified integration of ordinary differential equations (ODEs) and differential algebraic equations (DAES) in Chap. 4. 2.1.4 Summary For functions that are given by finitely many intrinsic functions and elementary Op- erations, the DA approach is equivalent to evaluating the code list with an n-th order automatic differentiation tool. However, the ease with which the DA approach com- putes Taylor polynomials, and therefore derivatives, of any given algorithm up to machine precision and virtually arbitrary order makes the differential algebra “D” an important computational tool for numerical analysis. Together with the powerful DA Fixed Point Theorem that guarantees convergence in at most n+ 1 steps, the method lies at the core of the map approach that has been used successfully in the analysis and the design of particle accelerators [10, 12, 14, 15, 17]. The differential algebra ”D, has been implemented in the arbitrary order code COSY Infinity [13, 23, 84]. And by storing only non-vanishing coefficients, the implementation can handle even high-dimensional problems to very high orders. By combining the symbolic nature Of operations on truncated polynomials with the 10 fl~ a“; fill Lin-‘- " ,. . ugh cl ¢ ‘1 in- "I D'"! w‘ t. t numerical operations on floating point coefficients, differential algebra based methods offer a unique combination of exact symbolic Operations with the efficiency of floating point computations. Moreover, the DA approach avoids the pitfalls of conventional computer algebra, since it does not suffer from the dramatic increase in storage that often plagues symbolic computer algebra tools in high-order applications. At the same time, propagating high order descriptions of computer code avoids many of the problems associated with traditional numerical methods and opens the door for rigorous high-order sensitivity analysis of dynamical systems [104]. 2.2 Interval Arithmetic The use of interval methods in computational sciences was started by R. Moore [91], who observed that, “Computers carry only a limited number of significant digits. Repeating a computer calculation with more significant digits does not necessarily increase the number of significant digits in the result.” To illustrate this problem, consider the sequence of numbers defined recursively by 1130 = 1 — 1.0—21 and 513".“ = 1:712- (27) If we compute this sequence on a computer with 10, 16, or even 20 digit accuracy, it is $o=$1=---=$75=1, (2.8) while in fact .1375 < 10““). Since the question in numerical computations is often how the result compares to some fixed number, this is actually an important problem and computational errors like the above are relatively frequent. These computational errors have many sources: floating point errors, rounding errors, truncation, and 11 *' *3]le i'“ . .,pl'" ,6 Lu. initial errors. And since they are hard to identify at the time of actual computation, they are often undetected and can have dire consequences. The previous example illustrates the main problem of naively using computers to implement results of numerical analysis. Due to the fact the computers can only represent a small subset of the rational numbers Q accurately, the results obtained by computer operations are generally only approximations of the mathematically correct results. While traditional methods of numerical analysis are very powerful and provide rigorous results in an ideal world, their transformations to computer systems suffer from the intrinsic limitations of finite state machines. However, modern interval techniques use knowledge of these limitations to com- pute results that are both accurate and reliable. In this section we summarize the fundamentals of conventional interval methods. The development of powerful new interval techniques is the main focus Of this dissertation. Interval Notations Before presenting details of interval analysis, we introduce notation and conventions that will be used in connection with intervals throughout the text. For two numbers a, b 6 R, the interval [a, b] is defined by [a,b]={a:€lR]a_<_:ch}. (2.9) Unless explicitly stated otherwise, intervals will always be assumed to be closed and bounded and therefore compact. Whenever we talk about intervals in general, we use boldface to distinguish them from regular numbers: X is a point, while X denotes a compact interval. To express containment of points and intervals within intervals, we use the stan- 12 a" U" 4‘4” l1; ‘1: .,.. _ J... _-x=_. 11‘ M4,. “4" ,1 dard notation xE[a,b]¢)anSb (2.10) [a,b]C[c,d]®cgagbgd (2.11) The union and intersection of two intervals are defined by [a, b] U [c, d] = {as E IRlz E [a, b] or :1: E [c,d]} (2.12) [a, b] (I [c, d] = {:12 E Rlx E [a, b] and a: E [c,d]} (2.13) The union of two intervals will generally not be an interval. The intersection of two intervals on the other hand is either an interval or empty. Lastly, we denote the midpoint and the width of an interval by m ([a, b]) and w ([a, b]), respectively. In the case Of IR” with v > 1, we define interval boxes to be vectors of intervals and all previously mentioned operations and relations are defined componentwise. The width w of an interval box is called the magnitude and is defined as the maximum of the widths of its components. 2.2.1 Elementary Operations and Functions We define elementary Operations on intervals set-wise such that they satisfy the fun- damental inclusion requirement: for * E {+, —, x, /} we demand that [a,b]*[c,dl={x*ylas:v_<_b,c5ysd} (2.14) In other words, the sum, difference, product, or quotient Of two intervals is the set of sums, differences, products, or quotients of pairs of real numbers, one from each of the two intervals. All elementary operations can be defined in such a way that the 13 any ‘ =- I,,,,,Iy (+,><, ) (as...) ii 1) $(+a X, )y (- 4 Ix(®1®t )Iy Figure 2.2: Commuting diagram for the elementary operations on intervals. results of Operations between intervals are again intervals: Ia: G Iy :: [$1 + 91,552 + y2] (2'15) Ian 9 I1, = [171— 312,332 — 91] (2-16) Iz®Iy = [min{$lyly$1y27$2y17$2y2}ymax‘f-lela-lebeI/la3729/2“ (2-17) If the interval 1,, does not contain 0, the multiplicative inverse is given by If = [1/y2,1/yll (2.18) and the quotient Of two intervals is defined in terms of the multiplicative inverse and multiplication: hegzhegi (mm With these definitions, the diagram in Fig. 2.2 closes and commutes under the inclu- sion relationship. From now on, we will denote the elementary Operations on intervals and real numbers by the same symbols. One problem of interval arithmetic is the dependency problem. It is best illustrated with the interval I = [-1, 1] and the computation I - I = [—2,2]. (2.20) 14 lu£»l‘ w- I Art. ‘4), ‘. I . .1 ‘1 :1 a“ l I '] [:pr. "‘~ ‘Jl - l \‘ ., L" l" 5"» , . .2, Thus, the width of I — I is twice as large as the width of the original interval I. This problem, known as the cancellation effect, is a manifestation of the more general problem that conventional interval arithmetic has no provision to distinguish between independent and dependent variables. As far as interval analysis is concerned, the two intervals on the left hand side of Eqn. (2.20) are two different entities with no relation between them, and the fact that they have the same endpoints is seen as a mere coincidence. The dependency problem especially affects large computations, since not all occurrences of the dependency problem are as easy tO spot as in this example. We note that work is under way to include techniques for resolving the dependency problem within Fortran and C++ compilers [136, 137], where the Opti- mizers already perform a dependency analysis that is capable of also reducing the dependency problem. We mention that even conventional floating point arithmetic is beset by cancellation effects: on a computer with less than twenty digits of accuracy (1020 +1) — 1020 = 0. (2.21) This exemplifies a general problem of floating point arithmetic: computing the dif- ference of two similar numbers Often results in large relative errors. And although interval computations always enclose the correct results, the cancellation effect can lead to significant overestimations. Similar to the way the elementary operations are defined on intervals to satisfy Eqn. (2.14), the standard mathematical functions for intervals are defined such that the basic requirement f (la, bl) = {f(rr) la S cc S. b} (2.22) is always satisfied. In many cases the computation of interval extensions of mathe matical functions is straightforward and details will be discussed in Sec. 6.1. 15 1')? a " flu. .. r ‘ , ‘ A‘ .... .. ift‘rh‘ _ll'. ill v 1 2‘5 "it; ' - . l)": 1,“,- l n. .I 1,, l 2.2.2 Interval Arithmetic and Set Theory Since interval arithmetic has its roots in the desire to compute rigorous results with limited precision floating point hardware, it has originally been used to model small intervals. However, by realizing that the interval [a, b] is actually the set of real num- bers between a and b, where a and b can be rational numbers with exact floating point representations, computer programs become able to accurately handle real numbers. Thus, by using intervals to represent sets of real numbers, intervals bring set theory and the ability to represent any real number to computers. Utilizing the elementary operations and functions defined in the previous section, we can Obtain the Fundamental Theorem of Interval Analysis [94]. Theorem 2.2 (Fundamental Theorem of Interval Analysis). Let I = [a, b] be a compact real interval. If f : I —> IR is a continuous function and f is its inclusion monotone interval extension, then 17 6 [a,b] => f(x) 6 f(I)- (2.23) The Fundamental Theorem of Interval Analysis allows the use of computers, de- spite their limitation, to rigorously answer the question: Given a function f described by a computer program, consisting of finitely many elementary Operations and intrinsic functions, and an interval I, find a set of real numbers that is guaranteed to contain the image of I under the function f. By allowing the computation of intervals that are guaranteed to contain the range 1' (I) of f over I, the use of interval analysis allows the use Of computers with finite precision to obtain rigorous and trustworthy results. 16 I [‘1 :14 Fist 1 j.;,.‘] ,l‘ldl \t') “.1 7‘, \‘su Some of the most powerful aspects of interval computations are tied to the Brouwer Fixed Point Theorem, which is the finite-dimensional version of the Schauder Fixed Point Theorem: Theorem 2.3 (Brouwer Fixed Point Theorem). Let D be a non—empty compact convex subset of IR” and suppose f is a continuous mapping such that the range 1‘ (D) C D. Then f has a fixed point in D, i.e. there is :1: E D such that f(at) : 1:. Since interval boxes satisfy the requirements on D, the Brouwer Fixed Point Theorem can be combined with the Fundamental Theorem of Interval Analysis to allow computer based proofs of the existence of solutions to linear and non-linear systems. One particular important application of this, the basic interval Newton method, will be discussed in Sec. 2.2.3. The Wrapping Effect Much like the dependency problem, the wrapping effect can lead to substantial over- estimations in interval computations. It is caused by the need to wrap the results of interval computations in interval boxes that are parallel to the axes. The effect is demonstrated in Fig. 2.3: a rotation by 45 degrees leads to an overestimation Of the magnitude of the result by a factor of \/2. Figure 2.3: Illustration of the wrapping effect caused by a rotation of 45 degrees in the plane. While efforts have been made to fight the wrapping effect by allowing rotated and even skewed interval boxes [78, 79], the fundamental problem remains and Fig. 2.4 17 p. .. 1', 5'1 ",. t 1.1“ ‘. c ‘ 1 ~ «.2 I J g?" crlr . it. \1 l 0. My, .q . illustrates how the fact that the edges of interval boxes are always given by straight lines can lead to arbitrary large overestimations of the image sets. Figure 2.4: Illustration of how the wrapping effect can lead to arbitrary large overes- timations. 2.2.3 Applications of Interval Arithmetic By utilizing numerical analysis tools that have been developed specifically for interval analysis, interval methods have been used successfully in a variety of applications over the last decades. In this section we take a closer look at some of the most important of these algorithms and applications: interval Newton methods, optimization, and integration of ODEs. Interval Newton Method Interval Newton methods are among the most fundamental and important interval algorithms. They are, in one form or another, part of almost any other interval algorithm. Since they utilize most of the concepts Of other interval methods, e.g. fixed point formulations, iteration, and intersections, they are outstanding examples of the more general interval analysis tools. Let f : [a, b] —-> IR be a C1 function such that either f’ > 0 or f’ < 0 on [a, b] and assume that 51:27 E [a, b] such that f (a) = 0. The goal of an interval Newton method is to compute a small interval enclosure Of :11. According to the Mean Value Theorem, for any y E [a, b] there is some 5 E [a, b] such that 0 = f(x) = f (y) + f’(€)(:r - y)- (2-24) 18 v_ .l l ,[.]o After choosing y = m([a, b]) and solving for 2:, it is M, , _muuz — (1 Jill f’(§) . (2.25) with some unknown 6 E [a, b]. However, if we use intervals to evaluate the interval extension of the derivative f’, we can guarantee that f’ (6) E f’ ([a, 5]). Assuming that 0 E f’ ([a, b]), it follows that f (m ([a, b1» "'f'(1a.b'1) ' (2'26) While f (m([a,b])) and m([a,b]) in the last equation are points, 1'" ([a,b]) is an 2: E m([a,b]) —— interval and the right hand side is therefore itself an interval. This last equation leads to the following definition Of the Newton operator N for an interval I C [a, b]: N (I) = m (I) — M (2.27) f’ (I) With the Newton operator, the foundation Of the standard interval Newton method is given by Thm. 2.4; more details of interval Newton methods can be found in [92, 94, 3, 57]. Theorem 2.4. For X 0 = I, define a sequence of intervals by Xn-l-l : N (Xn) n Xn- (2.28) If I contains a zero of f, X n will contain the zero for any n. Moreover, if for some n the intersection N (X n) n X n is empty, I did not contain a zero of f. It is important to note that intersecting the result of the Newton operator with the previous enclosure ensures that the elements of the sequence never increase in size. Thus, unlike in the regular Newton method where the sequence of iterates might diverge away from a zero, the sequence of magnitudes w (X n) is monotonically decreasing. 19 :3 '{3‘ 1166?. 9125151 I Y ‘zf "‘1- ‘.zl nun . .1 A; ‘-l‘-.,‘ ,IlM ' . a“. Optimization Interval algorithms for constrained and unconstrained optimization are based on an exhaustive search of the search initial region. They usually follow the general scheme for Lipschitz optimization [106], but use interval evaluations Of the Objective function instead of estimations based on Lipschitz constants. The general problem of constrained and unconstrained optimization is given by minimize: ¢(2:) (2.29a) such that: C(23) = 0 (2.29b) 9(2) 3 0. (2.29c) ct : D C IR" —+ IR is the objective function; c : D —+ IR’c and g : D —> IRl are constraint conditions. If l 2 2, condition (2.29c) is understood to apply to each of the components of g. If k = l = 0, the problem is considered to be unconstrained. The basic algorithm for the exhaustive search starts by dividing the initial search region D into sub-domains and maintaining a list of boxes where the minimum might be. The boxes in the list are repeatedly tested and split further if necessary. While this by itself leads to an exponential growth Of the number of boxes, the main idea of the algorithm is to speed up the search by rejecting boxes from the list; i.e., to rigorously determine that the minimum cannot be inside a particular box. Interval methods play an important part in this algorithm. They are used in storing information on possible regions in the list and, more importantly, they are used to compute verified range enclosures Of the objective function to possibly reject boxes from the list of candidate boxes or to verify that the box contains the global minimum. More details of such algorithms can be found in [114, 58, 56]- Sophisticated extensions of the basic branch and bound algorithms exist that use 20 u 1,” . .‘\ El. J. I‘l‘,“ 1 H ‘H. a *‘T 1,. .~-. ‘1", -.1- ‘| interval methods to improve the decision process Of rejecting boxes [115, 34]. However, while global optimization with intervals has found successful applications [37], it Often struggles due to the exponential growth in the number of candidate boxes. This is especially true for complicated multidimensional objective functions [21]. More details on constrained and unconstrained optimization can be found in [140, 29, 41, 55, 63, 95, 38]. We would like to highlight the GlObSol package [67, 68, 46] as one particular sophisticated implementation of global Optimization algorithms with intervals. ODE Initial Value Problems Consider the general explicit first order ODE initial value problem 22' :2 f(t,:r), 23(t0) = 2:0 (2.30) with a sufficiently smooth function f defined on a suitable subset of IR x R”. Inter- val methods for the integration Of ODE initial value problems provide enclosures for the final positions, such that at every time step t,- the actual solution to the initial value problem is contained in the solution interval. To achieve that goal, the integra- tion methods take extended initial conditions, mathematical truncation, and roundoff errors into account. However, mostly due to the wrapping effect discussed in Sec. 2.2.2, methods for the integration of ODE initial value problems are among the most demanding problems for designers of interval algorithms [98]. While the first attempts at interval based integration tools [92] have struggled with the prOpagation of large initial regions even over small time intervals, newer methods like AWA use changes of variables, ellipsoid interval boxes and other techniques [79, 101] to fight the wrapping effect with varying degrees of success. But only recently have interval-based techniques been developed that alleviate the problem of the wrapping effect in the integration of ODE initial 21 12.319 1 41.21.11 I [1 , .. l ‘ v- £915... «.11 at; .7 ."M r . 3.2.4 1" '11] 11.1 e I re "C u‘ ‘.;,_ \. ‘ ] gm) ~.‘ ad. 5"»... "Mid value problems [22, 82]. These methods will be summarized in Sec. 2.3.3 and their application to solar system dynamics will be presented in Sec. 5.1. By successfully avoiding the wrapping effect to high orders, the new approach allows the long term propagation of large initial condition sets without substantial overestimations. 2.2.4 Summary Interval arithmetic is arithmetic defined on sets of numbers instead of single numbers and its results have been connected to modern computing machinery and numerical analysis in 1962 by R. Moore [91]. The use of interval arithmetic allows rigorous computations Of verified results on computers by accounting for the effects of finite precision and roundoff errors, mathematical truncation, and extended sets Of starting values. As such, interval computations overcome the general limitations of finite state machines and bring the strict mathematical rigor of numerical analysis and set theory to the computational sciences. Good introductions and overviews to interval analysis can be found in [92, 94, 3]. Interval methods have been used successfully in a wide range Of applications: chemical and electrical engineering [102, 5], computer graphics [96], dynamical sys- tems and chaos [50, 124], computer assisted proofs [74, 75], and expert systems [70] to name only a few. However, the wrapping effect, the dependency problem, and the dimensionality curse prevent interval methods from widespread acceptance in the computational sciences. 2.3 Taylor Models Intervals are well suited to model enclosures Of real numbers. However, since they only propagate information about function values and neglect any further knowledge 22 w- ‘1‘."- ._.'. c. r‘ ‘\ “trig“A ”a. ., . a“ about derivatives, they are not well suited for the rigorous enclosure of functions. In fact, modeling functions with intervals usually requires splitting the domain D into a collection {D,-} of sufficiently small subdomains and evaluating the functions on each of the D,s. Then the “evaluation” of f consists of finding the appropriate interval enclosure 1' (D;) of the range of f over the subdomain D,-. While splitting a one-dimensional interval is usually acceptable, splitting higher dimensional interval boxes results in an exponential increase of the number of subdo- mains. This phenomenon is known as the dimensionality curse. It frequently plagues conventional interval techniques and limits their applicability to large, both in terms of dimensionality and domain size, problems. Taylor models Offer a remedy to this problem by combining high-order Taylor polynomials with real number coefficients and intervals for verification. Taylor mod— els allow the rigorous modeling of complicated multidimensional functions over do- mains that are usually several orders of magnitude larger than the ones over which conventional interval methods can work with the same sharpness. Definition 2.6 (Taylor Model). Let D C IR” be a box with 230 E D. Let P : D —+ IR‘” be a polynomial of order n and R C R” be an non-empty convex compact set. Then (P, 2:0,D,R) is called a Taylor model of order n with expansion point 2:0 over D. Following these notations, P is called the reference polynomial, 2:0 is the expansion or reference point, and R is called the remainder bound of the Taylor model. For the rigorous modeling of functions on computers, we usually view Taylor models as subsets of function spaces by virtue of the following definition. Definition 2.7 (Taylor Models as Sets Of Functions). Let T = (P, 2:0, D, R) be an order n Taylor model. Then, identify T with the set of functions f E C"+1(D, R‘”) 23 O—w 1'5 3,, If” r' - fl .1 § . 3"". r. , a: l "1» It". Ar- such that f (2:) —P(£L‘) E R for all 2: E D, and the n-th order Taylor series of f around 2:0 equals P. Furthermore, if a C"+1 function f is contained in a Taylor model T, we call T a Taylor model for f. To illustrate the concept of Taylor models, consider two Taylor models T1 and T5 of orders one and five given by T1 = (2:, 0, [—1.5,1.5] , (—1.122182, 1.122182» (2.31) T5 (2: — 23/31+ 25/51, 0, [—1.5, 1.5] , (~—0.015781,0.015781)) (2.32) The two Taylor models T1 and T5 have been computed as Taylor models for the sine function with COSY Infinity [23, 26] and are shown together with their remainder bounds in Fig. 2.5. While the first order Taylor model has remainder bounds that are similar in size to the domain interval; i.e., the first order reference polynomial is a rather bad approximation Of the sine function over the domain, the fifth order Taylor model encloses the sine function with a sharpness that is of the order of one percent of the domain size. We mention that the requirements on the Taylor expansion of the functions contained in a Taylor model are sometimes omitted, and Taylor models are then simply viewed as sets Of 6"“ functions that are R-close to the reference polynomial. However, these situations are rare and will always be mentioned explicitly. Similar to the way intervals allow the rigorous implementation of set theory and numerical analysis on computers, Taylor models allow the verified representation of arbitrary smooth functions within computer programs. As such, Taylor models com- bine the mathematically rigorous results of functional analysis with the approximative floating point representations on computers to obtain methods to rigorously compute analytical results. For purposes of illustration, Figure 2.6 compares the enclosure of a function by 24 Figure 2.5: Taylor models of order one and five for the sine function over the domain interval [—1.5,1.5]. Taylor models and regular intervals. In this case, the virtue Of the method lies in the fact that a single Taylor model over a relatively large domain box guarantees a sharp- ness that interval bounding cannot achieve, even with multiple smaller domains. This helps significantly in fighting the dimensionality curse inherent in interval bounding where it is of the utmost importance to avoid subdivisions Of the domain boxes, due to the otherwise exponential growth in the number of subdomains. Figure 2.6: Left: Enclosing a function by a Taylor model Of order eight; right: Interval bounding of the same function with 50 sub—intervals. We note that the actual implementation of Taylor models and related Operations in COSY Infinity is mostly due to K. Makino [82]. Since the definition of Taylor models is so intimately related to Taylor polynomials, the implementation of Taylor models rests on the foundations laid by the differential algebra package and the interval .25 1 «- [I 1'rl'. ind" u . 1., ‘ 1”,] ! riff], I“ ‘1', . Ll'gu'c l I . l I'M , V4 implementation to be discussed in Sec. 6.1. 2.3.1 Arithmetic and Operations Methods have been developed for arithmetic Operations on Taylor models that pre- serve the defining inclusion relationships. Hence these methods allow the computation of Taylor models for any sufficiently smooth computer function. In the following, let P and D be as above and denote by B(P, D) a guaranteed enclosure of the range of the polynomial function P over the box D. Moreover, for k E N, P0,) and P(k+) are the parts of P of orders up to k and greater than k, respectively. Using this notation, the next theorem summarizes results presented in [82] and shows how elementary operations can be defined on Taylor models. Theorem 2.5. Let T1 = (P1,2:0,D,R1) and T2 = (P2,2:O,D,R2) be two Taylor models of order n and define RP = R1 ' R2 + R1 ' B(P2, D) + B(Pla D) ' R2 + B((P1° P2)(n+)v D) (233) Obtain new Taylor models T5 and Tp by TS : (Pl + P2: $0,131 R1 + R2) 9 (2.34) TP ((PI ' P2)(n),$0, D1 RP) ° (2.35) Then, T; and Tp are Taylor models for the sum T5 and product Tp of T1 and T 2. In particular, for two functions f1 E T1 and f2 E T2, it is (f1 + f2) 6 Ts and (f1 - f2) 6 Tp. (2.36) Proof. If we define C"‘H functions 61 = f1 - P1 and 52 = f2 — P2, then 51(3) 6 R1 and 52(2)) E R2 for any 2: E D. Thus, for a given a E D ((1. + f2) — (P1 + 122)) (2:) = 3, (2:) + 62(2) 6 R1 + R2 .-_- Rs. (2.37) 26 in I Eb . nub . 7‘ how ".' N] a...“ Since the n-th order Taylor expansion of the sum f1+ f2 equals the sum of the Taylor series, T5 is indeed a Taylor model for the sum T1 + T2. Also, over the domain D ((f1'f2)—(P1'P2)(n)) = ((P1+61)(P2+62))—((Pl-P2)—(P1'P2)(n+)) = P1'52+51'P2+51'52+(P1'P2)(n+)- (2-38) Moreover, since the n—th order Taylor polynomial of f1 - f2 equals the polynomial product (P1 - P2)(n), T,» is a Taylor model for the product T1 -T2. E] The general theory of arithmetic on Taylor models has been developed in [83, 22, 82] and all elementary operations have been defined such that the basic inclusion properties are maintained throughout the arithmetic. In particular, all operations are inclusion monotone and lead to closed commuting diagrams as shown in Fig. 2.7. f) g (- T va T9 [*3 X] [39, ‘3] l, v fl+1']g r —’ Tfler ®ng Figure 2.7: Commuting diagram for elementary operations on Taylor models. As with the elementary operations, it is possible to extend the standard mathemat- ical functions to Taylor models such that the fundamental inclusions are maintained; the effect of this is illustrated by the closed commuting diagram for the sine func- tion shown in Fig. 2.8. The actual implementation of the intrinsic functions uses an approach resembling the way functions are defined on "D” to obtain Taylor models for the standard mathematical functions available on computers. More details on the 27 definition and implementation of intrinsic functions on Taylor models can be found in [83, 22, 82]. f 9 4' T1 [Sin( ' )1 [Sin( - )l l l sin(f) ‘ fi- sin(Tf) Figure 2.8: Commuting diagram for sine function on Taylor models. Last but not least, we define the antiderivation of Taylor models, which is a particularly powerful operation with important applications in verified Taylor model computations. Definition 2.8. For an n-th order Taylor model T = (P, .730, D, R) and k = 1, . . . , 1), let an. Qk 2/ P(n-1) (551,”-a$k—l,€k,$k+1a---,$v)d€k- (2.39) o The antiderivative 6,:1 of T is defined by 61:1(T) : (QM-730,1): (B(P(n) _ Phi-1), D) + R) ' B(xk: D» ' (2°40) Since (2;, is of order n, the definition assures that for a n-th order Taylor model T, the antiderivative 6;1(T) is again a n-th order Taylor model. Moreover, since all terms of P of exact order n are bound into the new remainder, this definition of the antiderivation is inclusion monotone and lets the diagram shown in Fig. 2.9 commute. As in the case of the differential algebra "D”, it is noteworthy that the antideriva- tion does not fundamentally differ from other intrinsic functions on Taylor mod- 28 al‘ H .3 ll} ’l’sgart. the .A I 'F‘wll Lil; L1“ ‘l T to“ :31 l l o“ f(€)d€k ‘ : 3"1(Tf) Figure 2.9: Commuting diagram for antiderivation of Taylor models. els. Moreover, since the corresponding DA operation is contracting and smoothness- preserving, it has desirable properties for computational applications. Finally, it should also be noted that the antiderivation of Taylor models is compatible with the corresponding Operation on the differential algebra ”DU. 2.3.2 Accuracy of Taylor Models As indicated earlier, one of the main advantages of Taylor models over conventional interval methods lies in the fact the propagating high-order information on the func- tions allows the accuracy of the enclosures to increase. In fact, the accuracy of the enclosures obtained by Taylor models scales with a high order of the domain size. This makes Taylor models particularly well suited for high-dimensional problems over large domains, since it reduces the number of necessary domain splittings. The following theorem states that in most cases the width of the remainder bounds scales in fact with the (n + 1)-st order of the domain size. Theorem 2.6. Let f be a function that is represented by a finite number of elemen- tary operations and intrinsic functions, and assume that B is an inclusion monotone ' palynomial bounder that scales at least linear with the magnitude of the domain. If 29 LJ‘ .1». T = (P, 2:0,D, R) is a Taylor model of order n obtained by evaluating the code list representing f, then the remainder bound R scales with the (n + 1)-st order of the magnitude of D. Proof. The proof follows by induction over the elementary operations and intrinsic functions that make up the code list for f. Firstly, the assertion is correct for the constant Taylor models and for the identity function, since the respective reference polynomials can be bound with an arbitrary precision. Since for two Taylor models A and B, R4“; = RA + By, the assertion also holds for the sum of Taylor models. As shown earlier, the remainder bound of the product A - B is given by RA-B = RA ° R3 + RA ° B (PB) + B (PA) R3 + B ((PA ° PB)(n+)) . (2.41) Since each of the terms scales with at least order (n + 1), so does their sum. Thus, the assertion holds for all elementary operations. (Note that division is defined and implemented in terms of the multiplicative inverse, which is conceptually an intrinsic function.) Finally, the computation of remainder bounds of intrinsic functions is exemplified by the exponential function. According to [82], RexpM) = (3 (PA — PA(0)) + RA)"+1 ' exp ([0,1] ' (3 ((PA * PA(0))) + RA))- (2-42) By inclusion monotonicity, the second term never increases with a decreasing domain size and thus, the product scales with the (n + 1)-st order of the domain size. Using the complete definitions [83, 82], similar arguments can be made for all the intrinsic functions of Taylor models, including the computation of the multiplicative inverse. Thus, the remainder bound of any finite Taylor model computation does indeed scale with the (n + 1)-st order of the domain. Cl 30 Figure 2.10 illustrates how the fact that the sharpness of Taylor models scales with the (n + 1)-st order of the domain allows us to obtain sharp bounds quickly, even in the multidimensional case. While the Taylor model of order seven has remainder bounds that are of the same order as the domain, the size of the remainder bounds quickly decreases with increasing order of the Taylor models. Finally, a Taylor model of order ten encloses the function with a remainder bound that is only a fraction of the underlying domain. Figure 2.10: Interval bounding of a two-dimensional function by Taylor models of orders seven, eight, nine, and ten. 31 2.3.3 3 'Fdilt if i'u'.’ la m m I w "i‘ , M ,1. tan},- 2.3.3 Applications of Taylor Models Since their first development in 1996 [83], Taylor model methods have been used for a variety of applications. While recent results will be discussed in later chapters, we summarize some of the fundamental Taylor model algorithms in this section to illustrate the general applicability of the approach. Bounding Ranges of Functions One of the most challenging problems in global optimization is the problem of deter- mining bounds on the ranges of functions over given domain boxes. While interval arithmetic can answer that question rigorously [113], the cancellation and dependency problems often lead to range bounds that are too pessimistic to be of practical value. Taylor models reduce the problem of bounding the range of arbitrary functions to the problem of bounding the range of polynomials. While this is still a hard problem, it is already substantially easier to solve than the more general problem [129]. These Taylor model based techniques have been used in several situations, including verified bounding of highly complex functions [21, 16] and global optimization [88]. Rigorous Quadrature Quadrature, the numerical computation of integrals, is an important area of numeri- cal analysis. Extensive theories provide formulas for computing approximative values of the integral and rigorous upper bounds for the integration errors. While these methods often converge with high-orders, they usually also require the computation of bounds of higher order derivatives. Although conventional interval methods and automatic differentiation can in principle be used to compute these bounds, the com- putational overhead for accurate results is often tremendous [92, 131] and the com- 32 n,‘ 1 ‘. 3“} Pl putations suffer from the well known problems of interval arithmetic and automatic differentiation: exponential inflation of code and dependency problems. Using the antiderivation of Taylor models, the computation of integrals becomes a mere application of the antiderivation. Simply modeling the function by a Taylor model over the region of interest and using the antiderivation yields, according to Def. 2.8, a rigorous enclosure of the primitive. Hence, subsequently computing a rigorous bound of the range of the Taylor model for the primitive gives the desired result. Details on using this approach for high-dimensional verified quadrature are given in [25]. Integration of Ordinary Differential Equations While conventional interval methods have long been used for the integration of ODEs [92, 66, 35], the wrapping effect and the dependency problem have prevented these methods from successfully integrating complicated systems over large time in- tervals. On the other hand, one of the most important applications of Taylor model methods lies in the computation of solution of ordinary differential equations (ODES) under avoidance of the wrapping effect for practical purposes [24]. To illustrate the approach, we consider the initial value problem 22' = f(t,:r), m(to) = :60. (2.43) It is a well established fact that the solution to Eqn. (2.43) can be obtained as the fixed point of the Picard operator ’P given by ’P(x) = 2:0 +/t f (T,:L')dr. (2.44) Using the antiderivation for Taylor models, the operator 17 can be extended to Taylor models and yields a method that allows the computation of verified enclosures of flows of Eqn. (2.43). This approach has been discussed in [24, 82] and has recently 33 been used successfully in a variety of applications ranging from solar system dynamics (of. Sec. 5.1 and [27]) to beam physics (c.f. Sec 5.2 and [81, 82]). Unlike methods that use regular interval computations to enclose the final conditions, the Taylor model approach avoids the wrapping effect to very high order and is therefore capable of propagating extended initial regions over large integration intervals. 2.3.4 Summary At the time of writing, Taylor models are a relatively new concept combining high or— der Taylor polynomials with floating point coefficients with intervals for verification. They allow verified computations while avoiding some of the difficulties inherent in normal interval arithmetic. The Taylor model approach guarantees inclusion of func- tional dependencies with an accuracy that scales with the (n + 1)-st order of the domain over which the functions are evaluated. In particular, as shown in [85], this method can often substantially alleviate the following problems inherent in conven- tional interval arithmetic: o Sharpness of the Result 0 Dependency Problem 0 Dimensionality Curse The remainder of this dissertation focuses on the development of new algorithms and computational methods utilizing the Taylor model approach to obtain guaranteed results. 34 Ch Chapter 3 Verified High-Order Inversion In this chapter we derive Taylor model methods to prove the existence of inverse functions and to compute Taylor models for them if they exist. The fundamental problem can be paraphrased as follows. Given a function f : D C R” —> R” (known only up to some accuracy) defined over a box D, is the function invertible over its range f (D)? And if so, find a representation of the inverse as accurately as possible. In real analysis, as opposed to interval analysis, it is usually sufficient to determine invertibility in the neighborhood of a point, and the Inverse Function Theorem es- tablishes a sufficient condition for this: Theorem 3.1 (Inverse Function Theorem). Let U C R” be open and assume that f : U ——> R” is of class C"+1 with n 2 0. If mflmyw4w an is a linear isomorphism, then there is a neighborhood V C U of $0 such that f : V -—> R” has a 6"“ inverse g : f (V) —> V such that g o f = idv and f 09 == idf(V)- (3.2) 35 While the Inverse Function Theorem is one of the most powerful theorems in conventional point analysis, its main limitations in practical situations are that it does not give any means of computing the neighborhood V over which the function f has an inverse, nor does it provide a method to actually determine that inverse function. This is usually not a problem in point analysis where the Inverse Function Theorem is used with great success. However, in numerical analysis and interval analysis we are often interested in both computing the inverse function and deciding the question of invertibility over extended regions. In this chapter we take a fresh look at the problem of computing rigorous enclo- sures of inverse functions. We first derive a new Taylor model based criterion for invertibility that requires only knowledge about first derivatives. While the method has important applications in studying dynamical systems [133], we will use it as a first step in computing Taylor models for inverse functions. In a second step, we show in Sec 3.2 how Taylor models for inverse functions can be computed from the origi- nal Taylor models. By studying several examples, we demonstrate that the method outperforms conventional interval approaches both in terms of accuracy and usable domain size. 3.1 Rigorous Invertibility Analysis The first step in the process of computing inverse functions is the determination of invertibility over the domain of interest. While the Inverse Function Theorem guarantees invertibility over neighborhoods of points, it gives no estimates on the sizes of these regions of invertibility. In this section we derive Taylor model based methods to rigorously prove invertibility of functions over given domains and demonstrate the applicability of the approach by comparing its performance with several conventional 36 interval methods. 3.1.1 Invertibility from First Derivatives There are a variety of ways to decide invertibility for a given function f, some of which rely on second and higher order derivatives. These methods are often computationally expensive since they require bounding complicated functions like the operator norm of the second derivative map over the domain of interest [2]. In this section we present a method that requires only bounds on first partial derivatives to decide the question of invertibility. Furthermore, and perhaps more importantly, the method exhibits some important structure regarding the points at which the derivatives actually have to be evaluated. We will later capitalize on this structure to significantly reduce cancellation problems in the necessary verified linear algebra. It should be noted that this method is not local in nature, but can indeed guarantee global invertibility everywhere in the given domain. The following theorem enables us to decide whether a given function is invertible over a domain by just evaluating its first derivatives. Hence it greatly reduces the computational overhead necessary. It seems to originate in works by K. Kovalevsky from the beginning of the twentieth century but has been “rediscovered” by many others (e.g. [103, E 5.3-4]). Theorem 3.2 (Invertibility from First Derivatives). Let D C IR” be a boa: and f I D —) IR” a C1 function. Assume that the matrix 3%(X1) 35:00) M: 3 a (3.3) %,,L‘;(x.) 35x.) is invertible for every choice of X1. . - - aXv 6 D ° Then f has a Cl-inverse defined on f (D); where f (D) denotes the range of f over D. 37 Proof. It suffices to show that f is injective over D. So assume that 3 y E f (D) and 331,2.) E D such that ff-Tl) = f(xz) = 31- (34) For t E [0,1] and i = 1, . . . , v, define auxiliary functions h.(t) = f.- (331 ' (1 - t) + $2 - t)- (3-5) Since h,(0) = h,(1) for all i = 1,. . . , v, the Mean Value Theorem asserts that for each i there exists a t,- E [0, 1] such that dh, EU»:(Vf)($1'(l-ti)+x2'ti)'(x2—$l)=0 (3.6) Then for i = 1,. . .,v, define Xi = $1 ' (1 — t,) + $2 ' ti- (3.7) By the convexity of D, x,- E D for all i = 1, . . . , v and moreover, (Vf,)-(:r2—$1)=0foralli=1,...,v. (3.8) However, the last equation is equivalent to M - (11:2 — $1) = 0. But by assumption M is regular and the kernel of the linear map associated with M is therefore just {0}. We conclude that $1 = 3:2 and therefore that f is one-to-one over the domain D. As such it has an inverse, which, by the Inverse Function Theorem, is also of class C1. [:1 Regularity of the Jacobian at each point in the domain is a necessary but insuffi- cient condition for invertibility. Thus, while Thm. 3.2 resembles the Inverse Function Theorem, it is in fact stronger in the requirements on f. To see that it is impossible to deduce invertibility from the fact that the Jacobian is non—singular at each point within the domain, consider the function x5 - 1011:3112 + 5503/4 f ( Z ) = ( y5 _ 10:32:13 +5134?! ) ’ (3.9) 38 rich 1 which is related to the complex function f : C —> C, f(2) = 25 via the identification z = x + iy. The derivative of f at any point (2:, y) E IR2 is given by "4 _ 2 2 4 3 _ 3 53: 301: 31 +53; 20:1: y 20:03; ) (310) M (x’ y) = ( 202:3/3 — 20x33; 511:4 — somzy2 + 53,4 The determinant det (M) : R2 —> R of M is a polynomial defined on the whole space R2 and det (M (:r, y)) = (5:1:4 — 305102;;2 + 5y4)2. (3.11) Hence, over the domain D = [01,1] x [01,1], the derivative of f is regular at every point. However, for (p1 = 1r/15 and 902 = 77r/15 f<:::::)=f(:::::) And since cos (p1 z 0.978, sin (pl :3 0.208, cos (pg z 0.105, and sin (pg 7:: 0.995, f satisfies the requirements of the Inverse Function Theorem at each point in the domain D but fails to be an isomorphism. This shows that a naive application of the Inverse Function Theorem is not sufficient to establish invertibility over extended domains. However, Thm. 3.2 shows that if we view the Jacobian of f as a function of v2 variables, instead of just v variables, its regularity at each point in the extended domain D” gives a sufficient criterion to guarantee invertibility of f. An immediate corollary of Thm. 3.2 is the following interval formulation of suffi- cient conditions for the existence of the inverse function, which in this version is very well known [100, 119]. Corollary 3.1. Let f and D be as in Thm. 3.2. For i,j = 1,...,v let pm C R be compact intervals such that 0f.- 6231' (:13) E pm- Va: E D. (3.13) If the interval matrix P = (Pm) is regular, then f has a Cl-inverse defined on f (D). 39 rise W “1‘ , _ . flint“. l. 3:;» a ,. u‘. "t '1‘ am: L 112V» 2%? .31. A practical algorithm based on Thm. 3.2 and Cor. 3.1 requires efficient and accu- rate methods to model partial derivatives and depends on the availability of efficient methods for the determination of regularity of an interval matrix with potentially large entries. Both these problems can become challenging for naive interval methods which have difficulties modeling complicated functions: they often give significant overestimations in bounding the derivatives of the functions, and suffer from the well known problems that come with an increasing number of variables and an increase in domain sizes. Moreover, interval methods are not particularly well suited for practical applica- tions of Thm. 3.2'via Cor. 3.1 because of difficulties in handling the fact that the individual rows of the Jacobian can each be evaluated at one point. A new Taylor model based method that circumvents this problem will be presented in the next sec- tion. But nevertheless, once the interval matrix P of Cor. 3.1 has been determined as accurately as possible, the next question is how to establish regularity of it. While in lower dimensions it may often be sufficient to compute the interval determinant of the matrix, more sophisticated methods are needed for higher dimensional problems. It has been shown that regularity of a matrix is equivalent to the componentwise distance to the next singular matrix being greater than one [122]. In the following paragraphs we present two advanced interval based methods to prove the regularity of the interval Jacobian P. In Sec. 3.1.3 we will compare these interval methods for the determination of invertibility with a new approach based on Taylor models. The first method of verifying the regularity of an interval matrix A = 431+ [—A, A] is based on the singular value decompositions of both the midpoint matrix ll and the A-matrix. If 0(M, i) denotes the i-th singular value of a matrix M (sorted in non- increasing order), then the following assertion holds [120, 121]: 40 ner“"’] [Lit A f L»: If!) mgr: -) , ”'1'": v_ k! l 'l; ,_ m, U. o;_ .5] ‘ Theorem 3.3. Given A = A + [-—A,A]. Then 0(A,1) < 0(A,n) implies that A is non-singular. The success of this method depends mostly on the sharpness of the models of the partial derivatives, and as such this method has difficulties scaling with dimensionality and complexity of the functions of interest. The next theorem has been presented in [100], and it is one of the most powerful approaches to proving the regularity of a given interval matrix A. Theorem 3.4. Given A = A + [—A,A]. If A is regular, let S be an approximate inverse of A. If the spectral radius p(I — S - A) is less than 1, then any matrix in A is non-singular. This method uses a preconditioning of the matrix A, which allows the method to work very well for medium sized problems, but as we will see below, in higher dimensions this method still suffers from the fact that the entries of the product matrix S - A are computed from addition and subtraction of intervals and as such are subject to significant overestimations due to cancellation and other dependency problems. We note that the most practical method to establish the contraction of an interval matrix A is to start out with an arbitrary non-empty interval vector 130 containing the origin and iterate n+1 = B - w,c with B = I — S - A. If for any k E N we can show that n+1 C (ck, we have proved that B is indeed contracting and therefore p(B) < 1 as required by Thm. 3.4. 3.1.2 Verifying Invertibility with Taylor Models We now come back to the question whether a given function f, defined over some interval box D, is invertible over its image. We Will combine the results 0f Thm. 3'2 41 with the regularity criterion T hm. 3.4 and Taylor model based techniques to derive a new algorithm to rigorously answer this question. Compared to the corresponding interval version of Cor. 3.1, we will be able to avoid an overly pessimistic behavior. The proof of Thm. 3.2 is based on proving regularity of the matrix gfi-(xl) 35m) 3,93%) 3530a) with points X1, . . . , Xv E D. Thus, any combination of entries from the same row of M can be evaluated using the same set of domain variables. While interval arithmetic cannot capitalize on this intrinsic structure, as we shall see in the following, Taylor models are particularly well suited for this task. Instead of modeling each partial derivative of the functions f,- by intervals, we now model each of the gradients by one single Taylor model, leaving us with the task of proving regularity of a matrix-valued Taylor model, which we will be doing using a method derived from Thm. 3.4. A naive application of Thm. 3.2 leads to matrix valued Taylor models depending on the v2 variables in D”. However, since the size of n-th order Taylor models in v variables is approximately (n+v)____(_m_)_!, (3.15) v v! n! a dependence on v2 variables would restrict the method to low orders even for rela- tively small problems. In the remainder of this section we present an algorithm that uses Thms. 3.2 and 3.4 to prove invertibility of functions modeled by Taylor models over extended domains and requires only Taylor models depending on v variables. Thus, the method avoids an inflation of dimensionality beyond the v variables that are required anyway to model the functions of interest. 42 Algorithm (1) For i, j = 1,. . . , v let Tm- be Taylor models for the partial derivatives of f; i.e., aft ‘ 5;; E Tm. (3.16) (2) Let T be an interval enclosure of the range of the Taylor models TM; i.e., g E Tm- => g(:v) E TM for all a: E D. (3) Let T be the midpoint matrix of T and similarly for the transpose T” of TT. (4) Let S be an approximate floating point inverse of TT. (5) Compute a new Taylor model N by Ni,j = 6i,j —' Z Si,kTIZ:j = 61d — ZSUCEJC, (3.17) k=1 k=1 where 6“ denotes the Kronecker Delta. (6) For i = 1, . . . , v compute Taylor models r,- with the Taylor models e,,+j that are defined over new domains, independent of those of the N,,J-’s: Ti = 2N3]: ' €u+j = ZN“ ° ev+j- (3.18) j=l 3'21 (7) Find interval enclosures I i for the Taylor models r,-. (8) If Ii C [—1, 1] for all i = 1, . . . , v, the original function f is invertible. Remarks In practice, the Taylor models Tm» needed in step (1) can be obtained by evaluating the first order automatic differentiation [111, 112] code belonging to f in Taylor model arithmetic. Additionally, extensions of the Taylor model approach are being 43 . . m w.“ developed that allow the determination of the TM directly from an extended Taylor model of f. Finally, if f is computed as the flow of a differential equation; i.e., the function that maps initial conditions to final conditions, Taylor models can be obtained as solutions of an augmented differential equation (c.f. Sec. 5.2). As outlined earlier, all Taylor models with the same index i; i.e., belonging to the same row of the matrix, can be evaluated using the same set of domain variables. The next step now is to compute a suitable preconditioning matrix for the transpose of T. However, instead of showing that the matrix T is regular, we will show the equivalent statement that the transpose of T is regular. By working with the transpose T7 in steps (3—5), the resulting entries of the matrix N are all computed using linear combinations of Taylor models that can be evaluated over the same set of domain variables. Hence each entry of N is again a Taylor model in v variables, and the columns of N are now modeled over the same set of domain variables. As a consequence, the sum 22:1 8.»chch can be evaluated in Taylor model arithmetic, and in particular any manifestations of the dependency problem are suppressed, as in other Taylor model computations. The final stage of the method in steps (6—8) is to show that the spectral radius of the resulting matrix-valued Taylor model N is less than unity. To establish this, the image of the unit ball (in the maximum norm) under the map N is considered. If we succeed in showing that the image is properly contained in the unit ball itself, invertibility is verified. Similar to the above, it proves advantageous not to work with the matrix N, but rather with its transpose NT. In this way, mixing of Taylor models from different columns of N is again avoided, which allows suppression of blow up in the necessary linear algebra. This is justified since the transposed matrix has the same spectral 44 9131. El radius as the original one. As a final application of Taylor models, we model the unit ball by Taylor models as follows. For any small 6 > 0, the identity function over the interval [—1, 1] is contained in (117,0, [-1, 1], (—€,e)). Thus, if we denote the i—th such Taylor model by e,-, a bound on the range of e,- contains the closed unit interval in z,. It should be noted that due to the special nature of the Taylor models r,-, the bounding is actually quite simple: each of the v terms of the sum that constitutes the Taylor models r,- has a different set of variables. Thus, there is no cancellation in the sum, and hence the sharpest possible bound of the sum equals the sum of the bounds of the individual terms. Moreover, the latter are just Taylor models multiplied by the interval [—1,1], and hence their ranges are just the ranges of the Taylor model multiplied by that interval. Summary The presented method utilizes Taylor models at two crucial points. Firstly Taylor models are used to model the partial derivatives of the function f. As will be seen later, this helps tremendously in fighting the issues of complexity and cancellation in modeling the derivatives. Moreover, the Taylor models are used to compute en- closures for the derivative matrices under consideration, and they are used for all arithmetic operations on these matrices which minimizes the overestimation due to cancellation, since the majority of the functional dependence is propagated in the reference polynomial and therefore not subject to the dependency problem [85]. Finally, we have been able to modify the computation in such a way that all arithmetic Operations on Taylor models can be performed using only v variables, which allows for a favorable computational complexity and allows the method to scale well to high dimensional problems. 45 3.1.3 )5 {rill} ten 0\' lill'illll [:lll Elli :33 we infra] 37% 1h :3 ~05 .U Ml 3.1.3 Comparison of Invertibility Criteria As outlined earlier, we are mainly interested in proving invertibility for a given func- tion over domains as large as possible. But not only the sheer sizes of the domains of invertibility matter, but an equally important question is how good we are able to fill out the maximum region of invertibility from within; i.e., how close to a critical point can we still prove invertibility. Moreover, in many applications we are likely to be confronted with high dimensional problems. The following examples will demonstrate how the new Taylor model based method performs in each of these areas as compared to standard interval methods. In each of the following examples we have tested a number of functions for invert- ibility over various domains. All tests have been performed with each of the methods listed below, and the results are presented in graphs and discussed. Note that the determinant test has only been used for problems of up to six variables. (1) Interval tests based on Cor. 3.1 and an interval determinant (2) Interval tests based on Cor. 3.1 and Thm. 3.3 (3) Interval tests based on Cor. 3.1 and Thm. 3.4 (4) Taylor model based tests as presented in Sec. 3.1.2 For the interpretation of the results discussed in this section it is important to note that in all cases the necessary bounds on range enclosures have been computed using regular interval arithmetic: no dedicated range bounders utilizing domain decompo- sition have been used either in the interval or Taylor model approaches to keep the computational overhead within reasonable limits. 46 Invertibility as a Function of Dimensionality The first example illustrates how the performance of the presented methods behaves as a function of dimensionality. Invertibility of one thousand 15-th order polynomials with uniformly distributed random coefficients in [-1, 1] has been tested with the methods (1) to (4). The domain of the random polynomials was [-0.005,0.005]” for v = 1,. . . , 10. Fig. 3.1 shows the percentage of these random polynomials that could be verified to be invertible as a function of dimensionality. To simulate realistic conditions, the Taylor models for these order 15 polynomials have been of order ten and the resulting remainder bounds are in the order of 10‘”. 100 fi— r 9 (r3 l r j— § u .. Z 0 g 80 _ a: o ‘] Q) I; 60 - x D l t I 2 5.: x ”5 a, 40 - ‘ l 0) £9 . x c u 8 + Interval Determmant 213 20 l x Singular Value Decomposition ' « ‘1 x Interval Spectral Radius . x a Taylor Model Based Spectral Radius i O ‘ ‘ 1 4V— l 4 E 0 2 4 6 8 10 Dimensionality v Figure 3.1: Percentage of random functions that can be shown to be invertible as a function of dimensionality. As expected, for increasing dimensionality the number of successfully established invertible polynomials is decreasing. But it is important to note that the Taylor model based method performs much better in establishing invertibility and suffers much less from an increase in the number of variables than the interval based methods do. 47 Unfortunately there is no good way to assess what fraction of the original functions are truly invertible; but it is to be expected that this fraction decreases with the dimensionality, accounting for the drop of predicted invertibility in all approaches. Invertibility as a Function of Domain Size For the next example we have investigated how the presented methods of establishing invertibility behave as functions of the domain size. We have restricted ourselves to the medium sized problem of 10—th order polynomials in six variables. For each polynomial we have applied the methods (1) to (4) over domains of increasing size. All domain boxes have been placed symmetrically around the origin and the abscissa in Fig. 3.2 shows the magnitude of the domains. For each of the different domain sizes we have plotted the percentage of polynomials for which invertibility could be proven using the four different schemes. As in the previous example, the computation was performed with 1000 random polynomials and the Taylor models were extended with artificial remainder bounds in the order of 10’”. The results of this experiment are illustrated in Fig. 3.2. In Thm. 2.6 we have shown that the accuracy of Taylor models scales with the (n + 1)-st order of the domain size, while interval arithmetic scales at most quadrat- ically with the domain sizes. That explains why the Taylor model based method can assert invertibility over domains that are much larger than the ones interval based methods can handle. Moreover, it should also be noted that the number of invertible polynomials naturally decreases with an increasing domain size, since the probability of singular points being contained in the domain D increases with the volume of D. 48 100 G F l l l— 1 I— *1 T I g a + Interval Determinant . a a x Singular Value Decomposition a, x ‘ E, n Interval Spectral Radius _5 80 - D u D Taylor Model Based Spectral Radius. Q ~ “’0 a U D Q X T D g 60 ' x l' D 0 ~ 11 D 2 . x " 0 .s x . D a "6 "I U 40 - x u ‘ 8, X I G D “'2' X I 0 q, + x n O I 5 20 " x ‘ .. o. + x x + ' . § x x ‘ X 0 4 4* +‘: .%£ a—t4: fiL—‘r :1:+T .Lg‘sé" 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 Domain size Figure 3.2: Percentage of random functions that can be shown to be invertible as a function of domain size. Invertibility as a Function of Non-Linearity This example demonstrates how the performance of the discussed methods changes with the non-linearity of the functions of interest. To that end we have considered functions f :2 (f1,...,f6) : [~1, 1]6 C R6 —+ R6 given by 6 _ a,- (I) 2;“ i" "2 (3.19) (I. — Zk=1bikak) + E2 with 5 ER and 6 x 6 matrices A 2 (age) and B = (bi,k) With coefficients in [“1:1l- fi($1‘l ' ' ' 12:6) : For small values of e and appropriate coefficients a”, and biJca this may give ill-defined functions because of vanishing denominators. However, in these cases we count the functions as not invertible. We have generated 500 functions by choosing the coefficients of the matrices A and B randomly in the interval [—1,1], and the presented methods have been used to prove invertibility over [—1,1]6. Fig. 3.3 shows the percentage of functions that can be shown to be invertible depending on the non-linearity 5. All Taylor model 49 fijrfilplll teal. .. LIEU. . 'l g fI-a‘! computations have been performed in order 6 and the value of 5 has been increased linearly from 0 to 250. 100 (D .5. so- ‘5 C .2 G) 3.5 60 ' t 9 .E ‘5 a, 40 - O) 5! C 8 x 5 20 - f + Interval Determinant ] o. _f + x Singular Value Decomposition x f u Interval Spectral Radius ,3 ,, +3». 0 Taylor Model Based Spectral Radius 0 ""13.” 4 n 1 1 1 O 50 100 150 200 250 Non-linearity s Figure 3.3: Percentage of functions that can be shown to be invertible as a function of the non-linearity e. The non-linearity of these functions is mostly determined by the quantity 5, such that for 5 >> 1, the functions are almost linear and their invertibility depends mostly on the invertibility of A. Since almost all computer generated random matrices are invertible, it is to be expected that all methods succeed in proving invertibility for sufficiently large 5. Fig. 3.3 indicates that this is indeed the case. For 5 ~ 1 on the other hand, the resulting functions show a large non-linearity and all methods fail in proving invertibility. It is likely that due to the size of the domain box and the high degree of non-linearity, these functions are truly not invertible. Between these two extremes, with increasing 8, the different methods become more successful in establishing invertibility. However, the Taylor model based method starts succeeding in proving invertibility very suddenly for e t: 20, while the success of the 50 Invert .,_, other methods sets in only for larger 5, and increases at a much slower rate. Invertibility in the Vicinity of a Critical Point As another example of how the presented methods scale to larger domain sizes and to illustrate how they behave in the neighborhood of a singular point, consider the function f = (fl, . . . , f6) : R6 ——> R6 with components f,- given by fi($1, - - - ,176) = (171' — 370)2 + ACOSHQR — $0)'($1r(i) — 550)) —A(:c,,(,-) — x0) sin(:r,- -— mg) (3.20) whererr(i)=i+1fori=1,...,5 and 1r(6)=1. At the point (2:0, . . . ,xo) the Jacobian of f vanishes and hence for symmetric domain boxes centered at the origin, 2:50 is an upper bound for the magnitude of the domain of invertibility. Fig. 3.4 shows the percentage of this maximal domain diameter 22:0 (for $0 = 0.1) over which the different methods can prove invertibility as a function of the parameter A, which is gradually increased from 0 to 2. The invertibility tests have been performed for the domains [-p - 0.001, p - 0.001]6 with p = 1,. . . , 100, and the figure shows the maximal value of p for which the various methods can determine invertibility as a function of A. The results indicate that the Taylor model based test can guarantee invertibility over a much larger region and is less sensitive to the perturbation A. It is important to note that all interval based methods perform equally poorly for increasing A. This indicates that the real problem of establishing invertibility in this example is with the use of conventional intervals for the modeling of the derivatives and not so much with the different methods to establish regularity of the interval matrix of derivatives. This example illustrates how Taylor models can enclose even complicated functions extremely accurate: the remainder bounds of the Taylor models 51 la: 100 We 0 m I F r— f r *1 w l D D a g 1‘ I! D a D o '5 80 " "‘ G U o D T .— I D (U n D D S ' .. ° 0 a 'U T n D D D D F: 60 - T n D D _] .g Ii " I! D U D n D 6 ‘ " l 0 D 0 El E " I . fi ‘ B 40 - ! " x ,. - 0 n T " l I O) . If n T " ll #3 + Interval Determmant 3 x Singular Value Decomposition 8 20 - n Interval Spectral Radius - o. D Taylor Model Based Spectral Radius 0 __1 1_ l l 1 L A __l_ 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Complexity and non-linearity A Figure 3.4: Percentage of maximal domain size over which invertibility can be proven as a function of complexity and non-linearity in the partial derivatives. for the partial derivatives are all in the order of 10‘”. Invertibility as a Function of Functional Complexity In this last example we study how the presented methods behave as a function of computational complexity of the original function f. To that end, we have emulated functional complexity using the following method. For the previously introduced random polynomials of order ten in six variables with coefficients in [—1, 1] we model computational complexity c by 1 c f(rv) — figs-(x) (3.21) where each of the f,- is a random polynomial as before and the scaling factor has been introduced to maintain the standard deviation of the coefficients of the resulting DOIynomials. All tests were performed over the [-—0.005, 0.005]6 domain box and the results are shown in Fig. 3.5. 52 It is important to note that the functions have been evaluated by adding the results of the individual polynomial evaluations. This is a realistic model of practical applications where the supplied functions often come as black boxes that do not permit any further simplifications to control cancellations. 100 I I 1— fi— I f I f I a a o E] a D 80 ["3 D a D “ n 60-} ~ .o. l Interval Determinant Singular Value Decomposition Interval Spectral Radius ~ Taylor Model Based Spectral Radius Percentage of invertible functions I DIX-O- 20 IX * . i i o s . . . _ o 5 1o 15 20 25 so 35 4o 45 so Complexity c .I. .I. I. .I. .I. Figure 3.5: Percentage of random functions that can be shown to be invertible as a function of functional complexity. The first observation worth noting here is that the Taylor model based method does not suffer from an increase in complexity, but rather even seems to improve with it. This is caused by the way we simulate computational complexity: while the coefficients in the original polynomials are uniformly distributed in [—1, 1], the coeffi- cients of the resulting “complex” one are not uniform anymore, and this distributional change leads to a simpler behavior with improved invertibility. As outlined in Sec. 2.3, Taylor models are particularly well suited to deal with the cancellation problems that plague conventional interval arithmetic, since the bigger Part of that cancellation takes place in the polynomial coefficient real number arith- 53 n‘ :- ‘1 use i metic and only a small fraction of it contributes directly to the remainder bounds. (This is in contrast to normal interval arithmetic where all the functional depen- dence is propagated in the remainder bound.) As such, Taylor models are expectedly much better in modeling the derivatives themselves and hence the Taylor model based methods can succeed in proving invertibility for computationally complex functions that cannot be properly modeled by intervals. 3.2 Guaranteed Enclosures of Inverse Functions Once the existence of an inverse function has been established for a function f con- tained in a Taylor model, the question of actually computing a Taylor model for the inverse function f ‘1 arises naturally. In the following we develop methods for com- puting inclusions of inverses of general functions in Taylor models. We show how Taylor model techniques can be combined with the methods for the verification of invertibility presented earlier to derive a new method that allows the computation of Taylor models containing the inverses of given invertible functions enclosed in Taylor models. 3.2.1 Polynomial Inverses If we start out with a Taylor model T for an invertible function f, the first step in computing a Taylor model S for f ‘1 is the computation of the reference polynomial of S. According to Def. 2.6, this requires the computation of the n—th order Taylor Polynomial of f ‘1. In this section we utilize the DA Fixed Point Theorem to calculate this n-th order Taylor polynomial of the inverse. Assume that f E C"+1(]R”,1R”) is an origin-preserving map; i.e., f (0) = 0. If the derivative M of f at 0 is a linear isomorphism, then the Inverse Function Theorem 54 guarantees that there is a neighborhood V of 0, such that f ‘1 exists and is also of class C“1 on f (V) Once the linearization of f has been shown to be invertible, within the DA framework it is possible to compute a representative of the equivalence class of [f 4],, from a single representative M of the equivalence class [ f1". Split the map Mas M=M+NM (3.22) where NM denotes the purely non-linear part of M. Composing from the right with the locally existing inverse results in I = MOM_1=M0M—1+NMOM_1 (3.23) => M-1 = M—1 o (I — NM 0 M“) = 0(M-1) (3.24) where I is the identity map. The last equation is actually a fixed point relation for the map M'l. Viewing this equation as a relation on equivalence classes, it turns out that 0 is contracting because (0 ([Aln) — 0 (Wall = M“ 0 (WMln ° [Bln - WMln ° [Alnl- (325) Thus, with [A],, and [8],, agreeing up to order k, the nilpotency of N implies that [NMln 0 [B]n and [NM],, 0 [4],, agree up to order k + 1. According to the DA Fixed Point Theorem, this assures the existence of a unique fixed point of the operator 0, which can be reached in at most n + 1 iterations. Moreover, this method generates all derivatives of up to order n of f"1 at the origin in finitely many steps through mere iteration. Because it only requires the iteration of a rather simple operator, this approach is particularly useful for computational applications. More details on this method can be found in [17]. While the presented method by itself allows only the computation of inverse poly- nomials for origin-preserving maps, we will show at the end of the next section how 55 m (””1“ . I‘leLt'J fir 3an [ll ilil‘l it can nevertheless be used for the computation of Taylor models for general inverse functions. 3.2.2 Inverse Taylor Models The computation of verified inverses starts by modeling the functional dependence of interest by a Taylor model and the a priori determination of invertibility. With the knowledge that the function contained in the Taylor model is actually invertible, the computation of an enclosure of the inverse follows the steps outlined in this section. First, it should be clarified what exactly constitutes an inverse Taylor model, since compatibility requirements of the underlying domains force us to consider two kinds of inverse Taylor models as given by the next definition. We denote by B (T) an inclusion of the range of all functions in the Taylor model T and by I a Taylor model (as, 0, D, (-—e,e)) of the identity function over the domain D (with some arbitrary small 8 > 0). Finally, for a polynomial P and a Taylor model T, P(T) denotes the result of the formal evaluation of the polynomial P with the argument T. This is well defined since it requires only additions and multiplications of Taylor models with Taylor models and real numbers. Definition 3.1 (Inverse Taylor Model). Let T = (Pm 1130, D, R) and S = (Gm yo, A,Q) be two Taylor models. S is called a left-inverse Taylor model for T if 1. GnOPnan, 2'Pn($0):y0; 3. fET=>f(D)CA, 4. B(G. (T) — I) g o. 56 5' _‘_l . .lkil ll ii'flldl Similarly, S is called a right-inverse Taylor model for T if T is a left-inverse Taylor model for S. It turns out that the rather general definition of left-inverse Taylor models en- sures the verified enclosure of left-inverses for all invertible functions contained in T. Moreover, the computation of left-inverse Taylor models is usually sufficient for prac- tical applications, their computation does not require a modification of the original domains, and their domains are easily computed. Theorem 3.5. Let T = (szro, D, R) and S = (Ga, yo, A, Q) be given Taylor models such that S is a left-inverse Taylor model for T. Assume that f E T is invertible over D. Then there is g E S such that g = f"1 on f (D). Proof. Define 6, : D —+ R” by (if (as) :2 Gn( f (cc)) — :12. Note that the Taylor expansion of 6; around 2:0 E D vanishes up to order n since the expansion of f around 2:0 agrees with R, up to order n and by assumption Gn(P,,(:r)) =n 3:. Since f is invertible over D, Vy E f (D), there is a unique 23,, E D such that f(IIJy) = y. Then, for y E f“(D) define the function {3(9) = (My) — 5r($y) = 07.01) — (Hf—101))- (3.26) Then, for any a: E D we have s (f (19)) = Gn (f (33)) - 5r (f‘1 (f (It‘ll) (3.27a) = a: + a, (x) — 6, (r‘ (f (z))) = e. (3271)) Thus 57 is the unique inverse to f on f (D)- Since lily) - Gn(y)l = l5f(f"l(y))| E 9, it follows that g is contained in the Taylor model 5; := (Gmyo, f (D) ,Q). As the continuous image of a compact set, 1’ (D) is itself compact and hence the function Q can be extended (n + 1)-times continuously differentiable to a function g 57 defined on the whole set A. Moreover, since (2 is a proper superset of g (f (D)), it is possible to satisfy (g — G") (A) C Q. Since the n-th order Taylor expansion of g around yo = f (.720) E D equals the n—th order expansion of g, which in turn equals G", it follows that g E S. Finally, for :r E D, g(f(:r)) == g(f(:r)) = 2:. E] The first step in translating the previous theorem into a practicable algorithm lies in the efficient computation of the inverse polynomial C". To this end, let T 2 (Pa, 2:0, D, R) be a given Taylor model and assume that P" is locally invertible around 2:0. Then consider the polynomial P(x) = P(a: + 2:0) — P(rro). (3.28) Apparently P(O) = 0, and hence, according to Sec. 3.2.1, [P‘l] can be calculated in at most n + 1 steps. With yo = P(ro) we obtain 012(9) as (My) = are + 13‘1 (y - 310) (3.29) Once the reference polynomial of the inverse Taylor model has been computed, the next step is the determination of an appropriate domain A such that either f E T => f (D) C A (left-inverse Taylor model) or g E S => g(A) C D (right-inverse Taylor model). Fig. 3.6 illustrates the various domains involved in the computation of inverse Taylor models and their relationships X=Uf(D) andY=nf(D)a (3-30) feT feT and hence X and Y are the smallest and largest sets that could be used as domains for left- and right-inverse Taylor models, respectively. Since a general computation and representation of X and Y is impossible, a practical application strives to find a small overestimation W for X and a large underestimation Z for Y such that 58 Ii‘vv] .1. AAA‘ , II. . r} 11' g) . mm.“ «31.4. [f 'i \[i‘] ‘_——---- 'o I I \ . \ ‘s 0’. 'I. ‘0' .-.-I-’ Figure 3.6: Illustration of the domains in the definitions of inverse Taylor models. Z C Y C X C W. Our implementation of this method chooses A = T(D) as the domain for the left-inverse Taylor model. Computation of the latter consists of bounding the range of the reference polynomial Pu over D and adding the remainder bound R. While the bound on the polynomial does not necessarily have to be accurate to obtain a left-inverse Taylor model, the practical application of left-inverse Taylor models benefits from overestimations as small as possible. 3.3 Examples of Taylor Model Inversion In the following we demonstrate the performance of the inversion methods presented in the previous two sections in examples. All computations have been performed using the Taylor model objects in the arbitrary order code COSY Infinity [23, 26, 82]. 59 ..-. 3.3.1 One-Dimensional Function As a first example, consider the one—dimensional sine function over the domain D = [—0.5, 0.5]. The functional dependence is modeled by a 19-th order Taylor model over that domain, and the result is presented in Table 3.1. The 19—th order Taylor model encloses the sine function with an accuracy that is close to the machine epsilon of approximately 10‘”. RDA VARIABLE: NO8 19, NV= 1 I COEFFICIENT ORDER EXPONENTS 1 1.000000000000000 1 1 2 -.1666666666666667 3 3 3 0.8333333333333333E-02 5 5 4 -.1984126984126984E-03 7 7 5 0.2755731922398589E-05 9 9 6 -.2505210838544172E-07 11 11 7 0.1605904383682162E-09 13 13 8 -.7647163731819817E-12 15 15 9 0.2811457254345521E-14 17 17 VAR REFERENCE POINT DOMAIN INTERVAL 1 0.000000000000000 [-.5000000000000000 ,0.5000000000000000] REHAINDER BOUND INTERVAL R [-.1085432243394823E-014,0.1085432243394823E-014] *t**#*****************##tt*********tttttttttttttttttttttt Table 3.1: Taylor model of one-dimensional sine function computed with COSY Infin- ity. Shown are Taylor coefficients, reference point, domain information, and remainder bound. Invertibility of the sine function is easily verified using the method presented in Sec. 3.1.2. The next step is the computation of a polynomial approximation of the inverse of the reference function. Not surprisingly, the methods of Sec. 3.2.2 obtain the 19—th order Taylor expansion of the arcsine function up to machine precision. Proper computation of the remainder bounds and the domain gives the left-inverse Tiii/let model shown in Tab. 3.2. Fig. 3.7 shows the difference between the arcsine function and the reference poly- 60 '2 l “I RDA VARIABLE: NO= 19, NV= 1 I COEFFICIENT ORDER EXPONENTS 1 1.000000000000000 1 1 2 0.1666666666666667 3 3 3 0.7500000000000000E-01 5 5 4 0.4464285714285714E-01 7 7 5 0.3038194444444444E-01 9 9 6 0.2237215909090909E-01 11 11 7 0.1735276442307693E-01 13 13 8 0.1396484375000000E-01 15 15 9 0.1155180089613970E-01 17 17 10 0.9761609529194068E-02 19 19 VAR REFERENCE POINT DOMAIN INTERVAL 1 0.000000000000000 ['.5210953054937487 ,0.5210953054937487] REMAINDER BOUND INTERVAL R [-.7707363654262549E-008,0.7707363654262549E-008] ******$*****¥*******************************¥*****¥**#*#* Table 3.2: Left-inverse Taylor model for the Taylor model shown in Table 3.1. Shown are Taylor coefficients, reference point, domain information, and remainder bound. nomial over the interval sin (D) = [-0.479425...,0.479425...]. As expected, the error stays within the remainder bounds (illustrated by the horizontal lines), and the remainder bounds represent an overestimation of the true error by less than a factor of four. It should be noted, however, that the arcsine function as such is not contained in the resulting left-inverse Taylor model. But, as outlined earlier, the part of the arcsine function necessary to act as a left-inverse of the sine function over the domain sin (D) is contained and can be extended such that the resulting function is defined over the whole domain of the left-inverse Taylor model. 3.3.2 Six-Dimensional Function The following example involves a six-dimensional exponential function that has a dense 8—th order Taylor polynomial with 3003 non-vanishing coefficients in each of the six reference polynomials. 61 le ' r. ‘r. . limit DU with 10-08 1 fl f I fl“ -5909 . 3 Error of Arc Sine Approximation Remainder Bounds ------- _1 908 1 Mr 1 1 1 1 1 1 1 -0.4 03 -0.2 -0.1 0 0.1 0.2 0.3 0.4 Figure 3.7: Difference between the arcsine function and the reference polynomial of the 19-th order Taylor model shown in Table 3.2. Let A = (a,,-) be an invertible 6 x6 matrix and consider f = (fl, . . . , f6) : R6 —-> R6 defined by 6 \ f,(:r1, . --,$6) = exp (201347)“ I — 1. (3.31) i=1 / If B = (b,,-) is the inverse matrix of A, the inverse function g = (g1, . . . , g6) of f is given by 6 91(91,- - - 1316) = 2% log (”l/j + 1) - i=1 Using the matrix 1 -—1 -1 —1 1 1 p—tb—iHi—‘HH p—Lh—lF—‘HI—‘H HHHHHH 1 1 1 -1 —1 1 1 ) —1 1 -1 —1 _1 ) (3.32) . (3.33) the 8-th order Taylor model method for verified inversion has been applied to the 62 firth: ll'f w Ffldlfl Still 1 L' ~ “All”? .7? in Arl 1"“ function f over the domain box D :2 [—-0.01,0.01]6. The resulting left-inverse Taylor model is defined over the domain [——0.061686, 0.061686]6, and it has been verified that it contains the true inverse function g. The remainder bounds of the individual components of the left-inverse Taylor model are listed in Tab. 3.3. Since the range of the exponential function and its inverse are in the order of 10“, the relative accuracy of the remainder bounds is about 10‘”). In Sec. 3.6.1 we will present a comparison of these results with conventional interval techniques. Component Remainder Bound 1 [-0.4190638646976846E—011, 0.4184087823912867E—011] [-0.2791908825275360E—011, 0.2791908821988238E—011] [-0.2791908824574486E—011, 0.2791908821987869E—011] [~0.1396454411975411E—011, 0.1396454410994258E—011] [—0.1396454411909750E—011, 0.1396454410994186E—011] [-0.1396454411225902E—011, 0.1396454410994267E—011] GWAOJM Table 3.3: Remainder Bounds of the left-inverse Taylor model for six-dimensional exponential function given by Eqns. (3.31) and (3.33). 3.4 A Superconvergent Newton Method While the conventional interval Newton method presented in Sec. 2.2.3 is a very powerful and versatile tool for the rigorous enclosing of zeros, its relatively slow quadratic convergence often results in the need for domain subdivisions. We will now apply the methods derived in the previous section to a higher order version of the interval Newton method. For practical root finding problems the new method offers a larger domain of convergence, and is more robust for complicated functions because of the decrease of the dependency and cancellation problems. Since the method converges with an order higher than two, it is also called superconvergent. 63 C7 5’1 While the conventional interval Newton methods are essentially based on first order Taylor approximations and inversions, the availability of high order Taylor models and the corresponding inversion tools enables us to overcome the limitation of slow convergence of the otherwise powerful interval Newton method. Based on the idea of approximating the function f by a high order Taylor model and computing a left inverse Taylor model for it, the Taylor model based computation of zeros is elegant and straightforward, and can easily be implemented using the COSY Infinity [23, 84] language environment. The algorithm takes the following form: (1) Set k = 0 and let D”) be a box that contains a zero of f and assume that f is invertible over 0(0). (2) Model the functional dependence of f over DU“) by an n-th order Taylor model Tl") with the reference point m (D02). (3) Compute a left-inverse Taylor model S (k) for T (k). (4) 001111311“? a new domain via DOC“) = D“) F) S (k) (0). If the enclosure DOC“) of the zero is not accurate enough, replace 1: by k + 1 and repeat from (2). Obviously, for n = 1 this method is equivalent to the traditional interval Newton method with centered forms (c.f. Sec. 2.2 and [92, 94, 3]). However, the high order version of the interval Newton algorithm allows for a fast and accurate determination 0f guaranteed enclosures of the zeros of a given function. Since the accuracy of n-th order Taylor models scales with the (n + 1)-st order of. the domain size, this method converges very rapidly once the domain size has become sufficiently small. This fact will be illustrated further in Sec. 3.6.2. Using the notation of the previous algorithmic description, the following is an important consequence: 64 Theorem 3.6 (Interval Newton for Taylor Models). Given a Taylor model T for f as before and assume that f is invertible over the domain D. If f has a zero 33" E D then, for 13(0) = D, it is 33" E DU“) for each k E No. On the other hand, if D“) 0 S0“) (0) = 0 for any k E No, then f has no zeros in D. Proof. We prove the first part by induction. For k = 0 the assertion is true by assumption, since :5“ E D = 13(0). Now we assume that the assertion is true for some k E No. Then, since .v“ E DU‘) and f E T”), 0 is contained in the domain of Sl"). Moreover, since 506) is a left-inverse Taylor model for TV“), El g E S (k) such that g( f (:r*)) = ac”. However, since f (22") = 0, that implies that x“ E SW (0). Consequently, the first part of the theorem holds for all k E No. The second part of the assertion follows from the first part: if there was a zero in DU“) then, according to the first part, it would also be in DU“) 1') S“) (0). However, if the intersection is empty, there was no zero to begin with. E] It is important to note that, just like for the conventional interval Newton method, it is generally not possible to give rigorous estimates on the speed of convergence of this method. However, once the requirements of Thm. 3.6 have been satisfied, it will usually converge with the (n + l)-st order of the size of the inclusion D; see also Sec. 3.6.2 for a further discussion. It should also be pointed out that unlike the ex- tended Interval Newton method, the presented algorithm cannot handle local extrema of f within the domain. However, this can easily be circumvented by performing ex- tended Interval Newton preconditioning steps that split the original domain until the Taylor models on the smaller boxes satisfy the requirements of Thm. 3.6. 65 3.4.1 vir’tm ‘ {infulf 0190 w - .2b 933' me .\e' a or ¢39L~ 4.5.le , l. .. with: 3.4.1 Examples The following examples demonstrate the Taylor model based high-order extension of the Interval Newton method. In particular, they show that the method converges extremely fast over relatively large domains. One-Dimensional Fixed Points This example compares the presented DA and Taylor model high-order extensions of the Newton method and the floating point and interval versions of that algorithm in a one-dimensional setting. It shows how the convergence rate of the high-order methods can indeed outperform the traditional algorithms by a wide margin. The floating depth h of a cylindrical trunk of density p = 0.66 and radius r in water is given by h = r (1 — cos (do/2)), (3.34) where do = 3.655403079564624 is the unique fixed point of a = sin(a) + 21rp. (3.35) This fixed-point relation for do can readily be transformed into a root finding problem and the resulting function is shown in Fig. 3.8. We used several Newton methods to determine the fixed point by solving the corresponding root finding problem: a traditional point Newton algorithm, an interval Newton method, a Taylor model based Newton method as presented in this section, and a non-verified point-polynomial version of the latter. The interval [3.3, 4.3] has been the initial enclosure of the zero, and the midpoint 3.8 has been used as a starting value for the non-verified methods. The computations have been performed in order 19 with an accuracy goal of 10‘”. The final approximations and enclosures of do are 66 4 . .— r r 3 _ sin(x)+1.32n-x 2 _ . 1 _ - 0 -1 . -2 . .. -3 . I . . 1 2 3 4 5 6 Figure 3.8: Transformed Fixed-Point problem: the fixed point do of Eqn. (3.35) is the intersection point of the two graphs. shown in Tab. 3.4 together with the number of iterations and CPU time required by the different methods. Newton Method Result Steps Time Floating Point 3.655403079564624 4 000001663 Polynomial 3.655403079564624 0.00162493 1 Interval [3.655403079564622, 3.655403079564625] 4 000003908 Taylor model [3.655403079564619, 3.655403079564628] 1 0.00484088 Table 3.4: Comparison of different Newton methods for the determination of the fixed point do of Eqn. (3.35). It should be noted that the high-order methods reach the desired accuracy after just a single step. The traditional algorithms on the other hand need four iterations to achieve the desired precision for this rather simple function. However, since the computational complexity of Taylor models exceeds the one of intervals, and since all methods work without domain splitting, the traditional algorithms are still more 67 efficient. Polynomial Functions The following one-dimensional example uses the verified Newton methods to enclose the value of 7r to almost machine epsilon — in fact the desired accuracy goal has been set to 10‘”. To that end, the sine function has been approximated by its 25—th order Taylor polynomial P, and the two Newton methods have been used to determine the unique zero of P in the interval [1.8, 4]. The order of the approximation has been chosen in such a way that the value of that unique zero and 7r agree up to machine epsilon. $2k+1 12 P(x) = §(-1)km (3.36) Due to the alternating sign of the terms, the naive implementation and evalua- tion of P exhibits strong cancellations. While it is clear that more efficient methods for the approximation of the sine function exist, this particular setup is typical for applications where the functional dependence is often given by complicated and re- dundant expressions that exhibit strong cancellation. While interval methods are known to have problems with this, the Taylor model approach can successfully avoid these problems. Fig. 3.9 shows the sine function and its 25-th order approximation between 0 and 12. In the region around 3 that is most interesting for this problem, the two functions agree up to machine precision. This example illustrates how Taylor models can successfully control cancellation and how the high-order methods converge much faster than the regular first-order Newton methods: over the same initial domain, it takes the regular interval Newton method eight iterations to enclose the zero with the desired accuracy. On the other hand, Tab. 3.5 shows that the Taylor model method achieves a comparable accuracy in just a single step. Moreover, the interval Newton method had to resort to extended 68 25"“ order TaWar-Approximation of Sine Function Sine Function --------- I— l l -1.5 ' 4 A o 2 4 6 3 1o 12 Figure 3.9: Taylor polynomial approximation of order 25 of the sine function over the interval [0, 12]. interval divisions and subsequent domain splitting in the first four of the eight steps. A Multidimensional Function This example illustrates the ability of the Taylor model based Newton method to read- ily extend to higher dimensional problems. The example shows the determination of the zeros of the six-dimensional exponential function defined in Sec. 3.3.2 with initial enclosures of the zeros given by [——0.25,0.25]6 for Taylor models and [—0.02,0.02]6 for the conventional interval Newton method. The computation has been performed with 8-th order Taylor models and the desired accuracy has been set to 10"”. The enclosures returned after each step of the iteration are listed in Tab. 3.6. Unlike the conventional interval method, where the computation and inversion of the interval matrix for the derivative poses a significant challenge, the Taylor model based Newton method generalizes easily from the one-dimensional case to 69 Table I v] , M~L t 1,. \ “pun!- '. n, . LIV] Method Step Enclosures of 7r Interval Newton 0 [1.800000000000000, 4.000000000000000] [2 .91907 767271 1 509, 4.000000000000001] [2.919077672711509, 3.427717355896160] [2.919077672711508, 3.165511666111006] [3.085798724266786, 3.16551 1666111007] [3.136629538119056, 3.154844930481172] [3.141102949424404, 3.141988434766514] [3.141592414707353, 3.141592894936869] [3.141592653589790, 3.141592653589801] [1 .800000000000000, 4.000000000000000] [3.141592653589765, 3.141592653589826] Taylor models i—aoooacfio‘fiww" Table 3.5: Enclosures of 7r, obtained by the interval Newton and Taylor model meth- ods. higher dimensional problems. Furthermore, the rapid convergence of the method over a relatively large domain demonstrates the power of the high-order approach. It converges in only two steps to an extremely accurate enclosure of the zero. But even more striking, the second step improves the sharpness by about 11 orders of magnitude. Once again, this is a vivid example of how the accuracy of Taylor models scales with the (n + 1)-st order of the domain size. The interval Newton method on the other hand, reached the desired accuracy in eight steps; since the computational complexity of the Taylor models is about three orders higher than that of intervals, this improvement does not yet justify the use of Taylor models. However, the interval Newton method required a much smaller initial enclosure of the zero to converge. Combining this with the dimensionality of the problem and assuming that the domain has to be split uniformly in each coordinate direction, the Taylor model method turns out to be two orders of magnitude more efficient than the interval Newton. We will revisit this particular aspect of the method in Sec. 3.6.2. 70 E3“ Method Step Magnitude of Enclosures of Zero Interval Newton 0 [-0.20000000000000E—001, 0.20000000000000E-001] [0.200000000000003001, 0.20000000000000E-001] [-0.20000000000000E—001, 0.200000000000003001] [~0.20000000000000E—001, 0.200000000000003001} [-0.87410175378563E-002, 0.87410175378563E—002] [-0.96852367808498E—003, 0.96852367808498E—003] [-0.54933148362165E-005, 0.54933148362165E-005] [-0.82475281503554E-010, 0.82475281503554E-010] [-0.94355649697386E—014, 0.94355649697386E-Ol4] [0.250000000000003000, 0.250000000000003000] [-0.47478831445046E—003, 0.47478831445046E—003] [-0.60171167482408E—014, 0.60171167482408E-014] Taylor models [\DD—‘OCN‘QCDCfii-bww'" Table 3.6: Enclosures of the zero for the six-dimensional exponential function, com- puted with an interval Newton method and the high order Taylor model approach. 3.5 Implicit Equations In this section we develop a method to compute Taylor model enclosures of functions that are described by implicit relations. This method has interesting applications in prescribed path control [58] and can even be used to obtain solution enclosures for implicit ordinary differential equations [59]. While the new approach uses the methods for computing inverse Taylor models discussed earlier, the method is based on the Implicit Function Theorem. Theorem 3.7 (Implicit Function Theorem). Given two open sets U C IR” and V C R‘" and afunction f E C"+1(U X V,R“’). Write f(iE, y) With 33 E U and y E V and suppose that (2:0, yo) 6 U x V and that f (270, yo) = 0- If D(fl({a:o}lew)nW) ($0ay0) 3 Rm “i Rm (3.37) is a linear isomorphism, then there is an open set U’ C U C IR” containing 1120, an Open set V’ C V C R” containing yo, and a unique 0"“ function g : U’ —+ V’ such 71 that for By tr inc Elflff iZELUH sigh] that for all :1: E U’ 9(130) = 310 and f (x, g (3)) = 0- (3-38) By combining this fundamental result of calculus with the presented methods for enclosing inverse functions, we are able to compute Taylor models for implicitly described functions 9 as above. Suppose that f satisfies the conditions of the Implicit Function Theorem. Then we define a function (I) : U’ x V’ :-—> R‘” by <1> ( 3 ) = ( “If,” ). (3.39) Since f satisfies the Implicit Function Theorem, it follows that (I) satisfies the con- ditions of the Inverse Function Theorem around (1:0, yo) 6 U’ X V’. Thus, there are neighborhoods U” of 2:0 and V” of go such that the inverse <1)“1 : U” x V” ——> R'“ exists and is well defined. Then we use "l to compute an explicit expression for the implicitly described function g: (;)=¢-‘(3)=(,E‘.))- (3'40) Thus, in the semi-algebraic framework of Taylor models, we can use an inver- sion algorithm to obtain rigorous enclosures of implicitly prescribed functions. The only drawback of this method is that augmenting the system requires additional vari- ables that might not be available due to memory constraints. However, in many important applications like prescribed path control [58], the variable 3: represents one-dimensional time. Thus, the augmented system is often enlarged by only a single variable, which is usually within the workable range of the Taylor model methods. We also mention that the technique of augmenting implicit systems has been applied extensively to first order in bifurcation theory [72] and to higher order for symplectic integration [11]. Other combinations of the Implicit Function Theorem 72 and fir are dis ‘ _‘ tr ‘4 . (11.1, i . ll [l EUHj 0 Eli and first order automatic differentiation [111, 52, 112] in the field of control theory are discussed by Evtushenko [42]. 3.5.1 Curves, Surfaces, Manifolds Hyperbolas are conic sections with an eccentricity e > 1. Denoting the foci of a hyperbola by f1 and f2. Hyperbola with a semi-major axis of a are described by where the two different signs denote the two branches of the hyperbola [61]. If we denote the distance between the origin and the foci by e; i.e., e = [ f1], the eccentricity e and the semi-major axis a are related by: a - e = e. (3.42) Given these relations and notation, the semi-minor axis b and the parameter p are connected via b2 = 32 -— a2 and pa 2: b2. (3.43) We denote the angle between the zit-axis and the vector connecting the focal point f1 with the point (x, y) on the positive branch of the hyperbola by oi. Then in the special case of a = e = 1, the cartesian coordinates x and y and the parameters 5 and o are connected by h, (5,45,33,31) = x2 — 52": 1 — 1 = 0 (3.443) hz (e, d), 3:, y) = (1 — ecos(¢)) - (flit — e)2 + y2 — (52 — 1) =: 0 (3.44b) In the remainder of this section, we use the approach outlined in Eqns. (3.39) and (3.40) to compute Taylor models containing charts for the cartesian coordinates 73 l0 [ll Si] : W‘iw maul The 1 mail 1". I: "y.. JON :1: and y as functions of the parameters 5 and d. The results show that the method can indeed be used to solve implicit equations and to compute parameterizations. Besides its applicability to static problems as discussed here, this method has an im- portant application in prescribed path control, where the implicit relation connecting coordinates and control functions is obtained as the result of an ODE integration step [58]. To compute charts for :1: and y from Eqns. (3.44), we start with initial conditions of so = 2 and (to = 7r/ 2 and use a standard Newton method to determine consistent initial conditions 2:0 and go such that hl (50195011301310) = h2 (50, $0,150, yo) = 0- (3-45) Once the consistent initial conditions 360 = 2 and yo = 3 have been determined to machine precision, we model the functions hl and I12 by 19-th order Taylor models T1 and T2 with reference point (so, (to, $0,110) and domain D = [199,201] x [1.56079633,1.58079633] x [1.95, 2.05] x [288,312]. (3.46) The resulting Taylor models have remainder bounds with a width in the order of the machine precision 10'16 and their reference polynomials have 59 and 4284 coefficients, respectively. After augmenting the system and defining the function (1) according to Eqn. (3.39), we use Taylor model inversion and Eqn. (3.40) to obtain Taylor models for a: and y as functions of s and (b. The resulting Taylor models are defined over the domain D’ = [199,201] x [1.56079633,1.58079633], (3.47) have reference polynomials with 210 and 205 non-vanishing coefficients, and remainder 74 6111115 fig. 3.? 25“ Eli - 138' bounds R1 = [—0.635775672959 x 10-8,0.635776786515 x 104], (3.483) R2 = [—0.537150447777 x 10—3, 0.537151554117 x 1078]. (3.48b) Fig. 3.10 shows the resulting surfaces 3: (e, (1)) and y (e, 45) over the domain box D’. Implicit x(e.¢) Implicit y(c.0) 1.56 ’51-99 1'56 15 0052.01 ' . 2.005 ' . _ . .01 Angle 0 Ewentricity e An9'9 4’ Eccentricity a Figure 3.10: Charts for the cartesian coordinates x and y of the two-parameter de- scription Eqns. (3.44) for hyperbolas. Remainder bounds are in the order of 10‘s. This example demonstrates how the Taylor model inversion tools can indeed be used for the computations of implicit functions and charts of smooth manifolds. While the Implicit Fimction Theorem can guarantee the existence of local parameteriza- tions and is therefore of great theoretical value, the presented method allows rigorous computations of these parameterizations. As illustrated by the example, the result- ing Taylor models offer a highly accurate description of the charts of the implicitly described manifolds. The ability to rigorously solve implicit equations with small overestimation is of great interest to fields as diverse as prescribed path control [58] and differential algebraic equations [59]. 3.6 Analysis of Inversion Algorithms Before we close this chapter on Taylor model based inversion methods, we use the previously presented examples to study some of the complexity and performance 75 312036 niques ahead 11150 3.6.] The C if Elli. the s: E characteristics that distinguish conventional interval methods and Taylor model tech— niques. While the theoretical limits that govern Taylor model computations have already been discussed in Sec. 2.3, here we focus on practical results to underline and illustrate these aspects. 3.6.1 Space and Time Requirements The computational complexity of individual Taylor model operations exceeds the cost of the corresponding standard interval arithmetic by a wide margin. For example, the storage size S of an n-th order Taylor model in 1) variables scales approximately as S=(n:v)=(n+v)!' (3.49) v! n! Moreover, while the cost of simple operations like addition and assignment scales linear with S, multiplication and intrinsic functions scale with higher orders of S [83, 22, 82]. However, since Taylor models can often avoid domain splitting, their real strength lies in high-dimensional problems with large domains. As an example, con- sider enclosing the inverse of a one-dimensional function as illustrated in Fig. 3.11. Using only intervals, this requires domains of the size of the desired accuracy, since the number of sub-intervals has to scale with the inverse of the desired accuracy. As another example, consider the six-dimensional exponential function presented in Sec. 3.3.2. The computational cost of the 8—th order Taylor models is about 1500 times higher than that of conventional intervals. However, achieving a uniform in- clusion accuracy of 10‘11 with intervals, as is often needed for root finding problems, would require approximately 109 sub-intervals in each dimension: resulting in 1054 separate interval evaluations, which would make it impossible to use conventional in- terval methods for that task. Thus, the advantages of Taylor models often outweigh their additional costs over naive interval techniques if enclosures of high precision are 76 l J I L I l T l I _L f l fl 1 T J_ I l I— 7 l Figure 3.11: Illustration of verified inversion with intervals. The sharpness of the enclosures scales linearly with the magnitude of the domain sub-intervals. desired. As a general rule of thumb, if conventional interval methods and Taylor models compete over the same domain and box-splitting is not required to achieve the neces- sary sharpness, the simple interval methods will almost always be more efficient. But if using intervals requires splitting of the domains, as is frequently the case, the use of Taylor models will generally become favorable. 3.6.2 Order Dependence This final example of the chapter illustrates how the accuracy of Taylor model ap- proaches indeed scales with the (n + 1)-st order of the domain size as shown in Thm. 2.6. For the six-dimensional exponential function defined in Sec. 3.3.2, Fig. 3.12 shows the magnitude of DU) after just one iteration of the 6-th order Taylor model Newton method as a function of 0“” = [—0.25, 025l- The data illustrate how that accuracy of the method scales with the 7-th order of 77 the initial domain 13(0), For purposes of illustration, a line of slope seven has been fitted to the data points. Deviations from that theoretical line for small domains are due to limits imposed on the verified computations by the machine epsilon. For large initial domains on the other hand, D“) is bigger than the theoretical line would imply, caused by a decreasing accuracy of the 6-th order Taylor polynomial approximations over the initial domain. Hence, for large domains the Taylor models operate in a region where the rules of (n + 1)-st order convergence fail to apply. 0 I T I I + + X ‘2 r + :1,” l ,x. + t” .4 r- t 4"”, _] 13;, ,v —6 ,. yr l -8 h ’3’”, .. ,4- I, -10 - ,«J « ’2' (I -12 )- A» + .4!" +’+’+ '14 7 +++++++++jxv - ’ + Accuracy vs. Domainsize -16 . ----- Linear Approxnmation wuth slope 7 - _18 " 1 I l 4— -2.5 -2 -1.5 -1 -0.5 0 Figure 3.12: Double-logarithmic plot of the accuracy of the first Taylor model Newton step as a function of initial domain size; results for 6-th order Taylor models. 3.7 Summary In this chapter we have presented new Taylor model based methods to rigorously answer the question of whether a function is invertible over a certain domain of interest. This method uses only first derivatives and as such it can be implemented with little computational expense. 78 ls intent it] an 00:15 01rd) [161' 9120 I fill] lite? As a next step, we developed a new method for computing rigorous enclosures of inverse functions of arbitrary complicated multidimensional functions. Its applicabil- ity and performance has been demonstrated in the computation of verified enclosures of the inverse function for given invertible functions. The method combines Tay- lor models with the method for proving invertibility and efficient differential algebra algorithms that allow the determination of exact n-th order polynomial inverses in finitely many steps. Using Taylor models, we have been able to overcome some limitations of conven- tional interval based methods to determine invertibility. Namely we have successfully shown that Taylor model based methods can model computationally complex func- tions much more accurately than interval methods. Since the accuracy of Taylor models scales with the (n -l- 1)-st order of the domain size, we have shown that the new methods work over significantly larger domains and scale better to high dimen- sional problems than interval based methods. Moreover, the newly developed method can handle complicated functions by using Taylor models which have been shown to control the dependency and cancellation problems inherent in regular interval arith- metic [85]. This is especially important in the establishment of invertibility prior to the computation of the inverse Taylor model itself. As an immediate application of the new method, an extension of the conventional interval Newton method has been presented. It has been demonstrated to converge much faster than the regular interval-based version of that algorithm. Moreover, combining the inversion tools with a fresh look at the Implicit Function Theorem resulted in a new method for the rigorous computation of implicitly described func— tions. Important applications of the new methods are in the computation of solutions of differential algebraic equations (c.f., Chap. 4 and [4, 107, 59]), the determination of existence of generating functions (c.f., Sec. 5.2), and the stability analysis of asteroid 79 (villi: orbits near the Lagrangian points [133]. 80 Chapter 4 Differential Algebraic Equations Under certain conditions, the solutions of ordinary differential equations (ODEs) and differential algebraic equations (DAEs) can be expanded in Taylor series in the independent variable and the initial conditions. In these cases, we can obtain good approximations of the solutions by computing the respective Taylor series [31, 32]. Moreover, the Taylor model approach often allows us to compute rigorous enclosures of the solutions of initial value problems [24, 82]. By using structural analysis [105, 108] and differentiation, it is often possible to transform a given DAE into an equivalent system of implicit ODEs. If the derived system is described by a Taylor model, representing each derivative by an independent variable, verified inversion methods discussed in the previous chapter can be utilized to solve for the highest derivatives as functions of lower order ones. The resulting Taylor model forms an enclosure of the right hand side of an explicit ODE initial value problem that is equivalent to the original DAE. While this explicit system is suitable for integration with Taylor model solvers [24], the intermediate inversion often requires a substantial increase in the dimensionality of the problem, limiting the approach to relatively small systems. An application of this inversion-based integration of differential equations has been discussed in [59]. 81 In this chapter we derive a method for the verified integration of general implicit ODEs that is based on the observation that solutions can be obtained as fixed points of a certain operator containing the antiderivation. We show that this Operator is particularly well suited for practical applications in Taylor model settings since its restriction to DA vectors is guaranteed to converge to the exact solution in at most 11 + 1 steps, where n E N is the order of the Taylor models. While sophisticated methods have been developed for the numerical integration of DAEs [51, 30, 53, 4], they are usually based on multistep methods that generally do not provide any means of verification and validation. However, since the new method also allows the computation of the differentiation index Vd of DAEs [4], it can be utilized for a verified index analysis. Combining this with a more traditional structural analysis of DAEs [105, 108], we get a scheme for transforming DAEs into implicit ODEs, which can then be solved with the new method. Since this combination uses verified Taylor models at every stage of the computation, it can be used to compute Taylor model enclosures of the solutions of DAEs. Additionally, by utilizing high order Taylor models (n > 20 is not uncommon), the scheme can even be applied to high-index problems that pose serious challenges to existing non-verified DAE integration methods. 4.1 Background and Motivation The standard first order initial value problem is usually associated with an explicit ordinary differential equation of the form x’ =: f(t,:L‘). (4.1) However, a more general first order ODE problem is given by the implicit form F(t, 516,:6’) = 0, (4.2) 82 when ‘5 as {0 ill . [1 HIGH 5'10. Slilll‘ l0 1 where the Jacobian 6F (t, u, v) 01) (4.3) is assumed to be non-singular for all arguments in an appropriate domain. According to the Implicit Function Theorem, it is then possible to turn the problem into an explicit first order system that is as smooth as the original problem. Thus, if F was sufficiently smooth to begin with, the implicit system is guaranteed to have a smooth solution. However, the solutions of implicit ODEs are not guaranteed to be unique. To illustrate this, consider the following initial value problem for an implicit ODE (:6’)2 = 3:2, 12(0) = 1. (4.4) ‘t are both valid solutions to this problem, Obviously 23+(t) = e‘ and 2:-(t) = e illustrating that smooth implicit systems may indeed have multiple smooth solutions. We will later address the implications of this phenomenon by introducing the concept of consistent points, which will enable us to guarantee the local uniqueness of solutions of implicit ODEs and DAEs. 4.1.1 First Order Differential Algebraic Equations ODEs with constraints are another extension of the regular eXplicit first order ODE problem 23’ = f(t, x, z) (4.5a) 0 = g(t, 27, 2)- (4.5b) Unlike in the case of standard ODEs, m(t) is not only a solution to the differential Dart (4.5a), but is also forced to satisfy the algebraic constraint condition (4.5b). While the system (4.5) is given in the so—called semi-explicit form, it can easily be 83 rewritten in an implicit form by introducing the variable 5 == (2:, z)T. I _ x'_f(t1$iz) _. 46 F(t7€a€)_( g(t,$,Z))-O. (') However, unlike in the case of implicit ODEs, the resulting Jacobian matrix 6F(t,u,v) __ I 0 Ti. .) is no longer regular. Generalizing the previous example and the definition of the general first order implicit ODE (4.2), we write the general first order differential algebraic equation (DAE) as F(t, :12, :r’) = 0, (4.8) where F is a sufficiently smooth function defined in an appropriate domain. In the case of DAEs we do not impose any additional restrictions on the regularity of the Jacobian (4.3). In fact, the rank and the structure of the Jacobian may even depend on the solution m(t) However, for simplicity we will always assume that the rank and the structure of the Jacobian along a single solution m(t) are independent of t. The class of DAEs contains all ODEs and all algebraic equations, as well as a large number of problems between these two extremes. As such, DAEs can generally not be solved by simple algebraic means. Moreover, unlike the case of ODEs, even smooth DAEs do not necessarily have solutions. It is interesting that the first re- sults on the existence and uniqueness of solutions of general DAEs have only recently appeared [109, 110]. While these results are based on theoretical arguments of dif- ferential geometry, in Sec. 4.3 we will present results that guarantee the existence of solutions to a large class of DAEs based on direct computations. The most common approach of determining the solution of a given DAE is to differentiate the system until we can pick 1) equations from the enlarged system such 84 that the Jacobian (4.3) of this new system is regular. The resulting implicit ODE system can often be solved for a solution of the original DAE problem. This basic idea is formalized by the definition of the index [4]: Definition 4.1 (Differentiation Index). For the general DAE system (4.8), the differentiation index 11,, along a solution m(t) is the minimum number of differentia- tions of the system which would be required to solve for 33’ uniquely in terms of the dependent variable a: and the independent variable t. Thus, the index Va) is defined in terms of the overdetermined system F(t,a:,:r’) = 0 F’(t, :6, 27’, :6”) = 0 (4.9) F(p)(t,:1:, ,m’, . . . ,x(p+1)) : 0 to be the smallest integer p so that (4.9) can be solved for the highest derivatives in terms of lower order derivatives and t. In the following, the index will indicate the differentiation index unless we compare to other concepts of indices in which case we will use the fully qualified name. It might be somewhat of a surprise that the definition of the index explicitly refers to a solution of (4.8). However, while we will always assume that the index does not change along a single solution; i.e., is independent of the independent variable t, a given system can have different indices along different solutions. To illustrate this, consider the system 0 = 13’ — 2 (4.103) 0:30-9) (4mm 0 =xy+z(1—y)—t (4.10c) 85 Assuming continuous solutions, the second equation permits exactly two solutions: ya(t) = 0 and ye (t) = 1. Depending on which of the solutions is compatible with the initial conditions, we get two different scenarios: a: The system reduces to at’ = z and 0 = z —- t; it has index 1 with the solution (180 + t2/2, 0, t)T. 6: The system reduces to 56’ = z and 0 = :1: — t; it has index 2 with the solution (t, 1,1)7'. While this example shows that the index does indeed depend on the solution of the DAE, it also illustrates how the solution manifold of a DAE can vary with the index. While any sufficiently smooth first order ODE system of size v has exactly v degrees of freedom, the solution manifold of the DAE system (4.8) can have any dimensionality between 0 and v. By permitting only discrete values for y, the previous example highlights how the constrained nature of DAEs even extends to the initial conditions. While the initial conditions for ODEs are usually free to vary over large ranges, the initial conditions of DAE problems have to be consistent with the constraints. However, frequently the consistent initialization of a DAE problem is not obvious: problems with an index larger than one have hidden constraints which are only revealed by the intermediate derivatives in (4.9). Thus, while the highest derivatives are needed to transform the system into a solvable ODE, all lower order derivatives are describing the constraint manifold and the consistent initial conditions. Standard DAE Problems Short of testing all possible combinations of derivatives, Def. 4.1 provides us with no efficient algorithm for computing the index of the DAB (4.8). While we will revisit 86 this issue in Sec. 4.3, many practical DAEs can be written in standard Hessenberg forms with known index, separating the ODE part and the algebraic constraints. Definition 4.2 (Hessenberg Index-1). The difl'erential algebraic equation IL" = f(t,a:,z) (4.113) 0 2 g(t, 2:, z), (4.11b) with the Jacobian matrix function Bg/Bz assumed to be regular for all t, is in Hes- senberg Index-1 form. DAEs in Hessenberg Index-1 form are also called semi-explicit index-1 systems. Since the constraints can in principle be solved for the algebraic variable z, these systems are closely related to implicit ODE problems. Another important class of DAEs is given by the Hessenberg Index—2 problems. Definition 4.3 (Hessenberg Index-2). The differential algebraic equation 33’ = f(t,:L‘, z) (4.12a) 0 = g(t, 2:), (4.12b) with the product of the Jacobian matrices (fig/6:6) (6 f / (92) assumed to be regular for all t, is in Hessenberg Index-2 form. We note that in many index-2 problems the algebraic quantity 2 plays the role of a Lagrange multiplier. Thus, many problems of classical mechanics can easily be written in the form of index-2 DAEs. An example of this will be discussed in Sec. 4.4.2. 4.1.2 General Differential Equations Up to this point, we have exclusively dealt with systems of first order differential equations. This is commonly justified by the observation that any n—th order differ- 87 ential equation in v variables can be transformed into an equivalent first order systems with up to n x v variables. This approach is known as order reduction and has many advantages for the theoretical analysis of general ODE problems. However, the in- creased size of the problems often challenges verified numerical integration methods. In the remainder of this chapter we will therefore focus on methods for the direct numerical integration of general high-order differential equations. For a definition of the most general DAE problem, we consider the independent scalar variable t, the v dependent variables 2:, = xj(t), and sufficiently smooth func- tions f,- for i = 1, . . . , v, and form the system f1(t,x1,...,a:(1£”),...,a:v,...,a:,(,£l”)) = 0 (4.13) fv (t,:c1,. ..,:1:[£"‘),...,2:,,,...,xffm) = 0. For a given j, the i-th equation of this system does not necessarily depend on 1:], the derivatives 33;”) for n < (,5, or even any of its derivatives. However, if there is a dependence on at least one of the derivatives (including the O-th derivative 3]” = 317,-), we denote the order of that highest derivative by f,g. 4.2 Verified Integration of Implicit ODEs While sophisticated general-purpose methods for the verified integration of explicit first order ODEs have been developed [79, 78, 82, 24, 97, 99], none of these can be readily used for the verified integration of implicit ODEs, let alone differential algebraic equations. As a first step toward the verified integration of general DAEs, we consider the problem of integrating the implicit ODE initial value problem (4.2). Unlike in standard ODE integration methods, we will not limit ourselves to first order 88 problems, but will eventually be able to integrate the arbitrary order problem F (t,:r,:r’,...,x(”)) , (4.14) where the derivative of F with respect to the highest derivatives in each of the compo- nents is assumed to be non-singular for all argument values in an appropriate domain. The new integration method is based on an extended use of antiderivation. To motivate the combination of high order methods and antiderivation for the verified integration of general differential equations, consider the explicit second order ODE initial value problem :6” = f (11:,33’, t), m(te) = 1:0, .r’(t0) = 3:2,. (4.15) While the conventional approach to solving this system is based on order reduction to a two-dimensional first order problem, in the framework of Taylor models we can use the intrinsic antiderivation and substitute f = :r’; i.e., 270) = $0 + [toté(r)d'r- (4.16) After inserting the expanded expression for m(t) into (4.15), we obtain an explicit first order ODE for 5. t E’ = F(€,t) = f (2:6 + avast). (4.17) to This system can readily be integrated with the existing Taylor model based integration scheme discussed in Sec. 2.3.3, and the solution m(t) of Eqn. (4.15) can be computed from the Taylor model for f (t) by a final application of antiderivation. c :1: = x0 + €(r)dr. (4.18) to While this approach seems obvious from a theoretical point of view, numerical analysis has traditionally avoided explicit references to antiderivation in its algo- rithms since antiderivation is usually not available in the classical computational 89 frameworks. Only recently, with an increased use of symbolic and semi-symbolic tools like Mathematica and COSY Infinity, has the explicit use of antiderivation in numerical computations become feasible. From a practical point of view, this ap- proach has the advantage of reducing the computational complexity of the right hand side of the ODE (4.15), often allowing for favorable computations. 4.2.1 High Order Taylor Model Solutions In this section we present a Taylor model based algorithm for the verified integration of the general first order ODE initial value problem F(t, IL‘, 113’) = O, (B(to) 2 $0. . (4.19) While we will later extend the algorithm to higher order ODEs, for now we assume that the problem is stated as an implicit first order system with an (n + 2)-times continuously differentiable v-dimensional function F and regular Jacobian matrix 6F(t, u, v) 6v (4.20) in appropriate domains. The method discussed in this section uses techniques and tools based on the dif. ferential algebra "D, and the Taylor model approach presented in Chap. 2. Utilizing DA vectors and Taylor model methods for the verified integration of initial value problems allows the propagation of initial conditions by not only expanding the so- lution in time, but also in the transverse variables [24]. By representing the initial conditions as additional DA variables, their dependence can be propagated through the integration process, allowing Taylor model based integration schemes to reduce the wrapping effect to high orders [86]. Moreover, in the context of the new algorithm, expanding the consistent initial derivative in the transverse variables further reduces 90 the wrapping effect and allows the system to be rewritten in a derivative-free, origin preserving form suitable for verified integration with Taylor models. Algorithm: Taylor Model Integration of Implicit ODEs For notational convenience, we write the initial condition of the implicit ODE prob- lem (4.19) as are, and denote functional dependence on the initial condition 2:0 in a suitable domain by y. With these conventions, a single n-th order integration step of the implicit first order ODE (4.19) consists of the following sub-steps: 1. Using a suitable numerical method like the Newton method discussed in Sec. 3.4, solve the implicit system F(th y’xl) = 0 (421) for a consistent initial condition :c’(t0, y) = 2:6(y). 2. Utilizing antiderivation, rewrite the original problem in a derivative-free form t (t.y,€) = F (t,v+/ €(T.y)dr,€) = 0, (4.22) to where 5 = g (t, y) = :r’ (t, y) has been substituted for the derivative of 2:. 3. Translate the problem into an origin-preserving form by substituting C (t, y) = E (t, y) — :1:[,(y) to obtain a new function ‘1’(t. y. C) = <1) (t. y. (0. y) + $6(y))- (4.23) 4. Since \Il(to, are, 0) = 0, within the differential algebraic framework it is possible to write the first order truncation of ‘1! without the constant part as ‘1’“, y, C) :1 LC(<) + LR“) g), (4.24) where L( and L R denote the linear parts in C and (t, y), respectively. 91 . If LC is regular, transform the previous expression into an equivalent fixed point formulation for C. ((t. y) = MC) = -L;1(4(t,y,g)— 14(4)). (4.25) . Using the operator ”H, define a sequence (av) of DA vectors in "DH, by a0 = 0 and au+1 = 71(au). (4.26) Then define the polynomial P(t, y) = on“. . Construct a Taylor model T with the reference polynomial P over an appropriate domain T x D C R x R” containing the reference point (to, 2:0) such that ’H(T) c T. (4.27) . Compute a Taylor model X from T by using the relation may) = y + f (<14. 4) ”some. (4.28) to As we will see in the remainder of this section, the final Taylor model X is in fact a Taylor model for the flow x(t,y) of the original ODE problem (4.19) for the particular choice of the initial derivative. 4.2.2 Mathematical Background The algorithm presented in the previous section rests on several non-trivial assertions that will be proven here. We provide its mathematical foundation and establish the basis for the discussion in the next subsection. The proofs can be split into two groups: the differential algebraic part of the algorithm and the Taylor model results. Implementation and user interface issues of the algorithm will be discussed in the next section. 92 Differential Algebraic Results In this section we discuss the differential algebraic results needed to justify the pre- sented algorithm for the verified integration of implicit ODEs. First, we show that the operator ’H introduced above is well defined and DA-contracting. We then show that up to an additive constant, its unique fixed point lies in the same equivalence class as the derivative of the flow of the original ODE problem. Lemma 4.1. The operator ”H given by Eqn. 4.25 is a contracting operator and self map on the set M = {a 6 "DH, [ Ma) > 0}. Proof. To first order regularity of L4 is equivalent to the assumed regularity of the Jacobian 8F(t, u, v) /0v. g(teme = g(tmyofl‘l‘ $6) =1 aa—iaoflmo (429) Thus by assumption, the linear map LC is regular in a neighborhood of the initial conditions and the operator ’H is therefore well defined on all of "DH”. To show that ’H is contracting on the subset M of nilpotent DA vectors, let a,b E M be given and assume that a and b agree up to order k. Since L21 is invertible, it suffices to show that A ((‘I’ftoa Elma) — Lc(a))—(‘1’(to,yo,bl— Lc(b))) > k- (4.30) Since ‘I’(te, 51:0, 0) = 0, the map \Il(te, 2:0, C) is origin-preserving and can be written as ‘I’(t0, 5130, C) = L((C) + LR(t03$O) +N(t0, 53010, (4.31) where N is a purely non-linear function. Thus, it suffices to show that N is contract- ing. However, if a and b agree up to order 16, their images N(to, are, a) and N(to, :60, b) trivially agree up to order k + 1. Finally, since \II is origin-preserving, ’H is indeed a self map of M, and therefore a contracting operator on M. E] 93 While this lemma guarantees the existence of a unique fixed point of the operator ’H, the next theorem summarizes the main result of the DA part of the presented algorithm. Theorem 4.1. If we denote the flow of the implicit first order ODE initial value problem (4.19) by m(t, y), then the fixed point of ’H is a representative for [540431) - xtfylln- (43?) around the expansion point (to, 9:0). Proof. In principle, this assertion follows from the construction of the operator 71. However, in the following we summarize observations on existence, uniqueness, and smoothness of the solutions to the IVP (4.19) that are required for a full justification. Since the original function F is of class 6"” and its Jacobian matrix is assumed to be non-singular over a suitable region containing (to, 2:0, 363), at least one solution to the initial value problem exists. Moreover, once a consistent initial derivative 276 has been fixed, the Inverse Function Theorem guarantees the existence of a unique solution m(t, y) in a neighborhood of (to, 1:0, 3:3). The solution is the unique solution to an explicit first order system, which can in principle be derived from the original implicit problem. Since the Inverse Function Theorem guarantees that the explicit system is also of class 6"” in a neighborhood of the consistent initial conditions, Thm. 5.2 ensures that the flow m(t, y) is a C"+2 function of its variables. Thus, the equivalence class of its derivative is well defined in C’“r 1, and since the solution C is unique, the fixed point of ’H is indeed a representative for that class. Cl Taylor Model Results In this section we prove the main Taylor model result needed for the presented al- gorithm: the self-inclusion of the Taylor model in (4.27) is a sufficient condition to 94 guarantee the enclosure of a fixed point of Eqn. (4.25). Definition 4.4 (L-Taylor Model). Let T = (P, 11:0, D, R) be an n-th order Taylor model and let L > 0 be a Lipschitz constant for P over D. Then the L-Taylor model TL is the set of all functions f E C0 (D,R“’) such that 1. f(x) — P(x) 6 R for all a: E D, 2. f(xo) = P(IBO), 3- [f(rv) - f(y)l S le - ill for all 23.9 E D- First, it should be noted that this gives a well defined and non—empty set of functions, since P E T 1,. And while there is an obvious connection to normal Taylor models, L—Taylor models contain a different set of functions: 1. Taylor models are “more restrictive” than the corresponding L-Taylor models, since the latter may contain non-differentiable functions. 2. Taylor models are also “less restrictive” than the corresponding L—Taylor mod- els, since they do not pose any limits on the Lipschitz constant of its members. Finally, unlike standard Taylor models, L-Taylor models are subsets of a Banach space, namely C0 (D, R”)- Lemma 4.2. The L-Taylor model TL defined as above is a convex subset of the Banach space C0 (D, R”) [82]- Proof. Given two functions f0 and f1 in TL and t 6 [0,1], define ft = t- f1 + (1 —t) - f0. Since the sum of continuous functions is continuous, ft 6 C0 (D, R”) for any t E [0, 1]. Moreover, by convexity of R Mas) — P(x) = (t - 11(4) + (1 - t) - 3(4)) - P(x) 6 R (4.33) 95 for any t E [0, 1]. Since UM?) " ft(y)l = lt' f1($) + (1 ‘tl 'f2(17) —t°f1(f/) — (1 “ t) ' 12(9)] St: |f1($) -f1(y)|+(1—t) ' Ifo(1‘) -fo(y)| (4-34) gt-L+(1—t)-L=L, TL is indeed a convex subset of the Banach space C0 (D, R”). E] Lemma 4.3. The L-Taylor model TL defined as above is a compact subset of the Banach space C0 (D, lR‘”) [82]. Proof. It suffices to show that every sequence (fu) in TL has at least one limit point in TL. According to the Arzela—Ascoli Theorem it suffices to show that (fu) is uniformly bounded and equicontinuous. We will first show that the sequence is uniformly bounded. Since P is a polynomial and therefore continuous on D, [P] assumes its finite maximum M1 over D at some point 5:: M1: [P(zit)[= max{|P(:r)| [z E D}. (4.35) Since R is bounded, there is a constant M2 > 0 such that y E R => |y| 3 M2. (4.36) After defining 5,, = f,, — P, 6,,(x) E R for all a: E D. Thus for any a: E D and any :2 e N |f,,(:r)| _<_ [P(cc)| + |6,,(:r)[ < M1 + M2 < oo. (4.37) Hence, the sequence is indeed uniformly bounded. Moreover, since by definition the sequence is also Lipschitzian with uniform Lipschitz constant L that is independent is u, it follows that it also equicontinuous. Thus, the L-Taylor model TL is a compact subset of the Banach space C0 (D, R'”). C] 96 The main result of this section is that the self-inclusion given in Eqn. (4.27) in indeed sufficient to guarantee that the fixed point of ’H is contained in the Taylor model T. The proof of this assertion is based on the next theorem, which is an almost immediate consequence of the previous two lemmas. Theorem 4.2. Let T and D be interval domains containing the points to and 51:0 respectively. Let L > 0 be a Lipschitz constant for P and define the L-Taylor model TL = (Pa “0.580) ,T X D4Rl- (438) If 7-l(TL) C T L, then the L-Taylor model TL contains a fixed point of ’H. Proof. Since ”H is a continuous operator on C0 (T x D, R”), the assertion follows from the last two lemmas and the Schauder Fixed Point Theorem. El According to this theorem, a fixed point of H is contained in the L-Taylor model TL. However, as indicated earlier, one of the fixed points is in fact equal to $’(t, y) — r[,(y) and is therefore of class 6"“. Moreover, according to Thm. 4.1, the n-th order Taylor expansion of the fixed point around (t0,:1:0) equals P. Thus, if [H] < 1, the fixed point is also contained in the regular Taylor model T = (P, (40,420) ,T x D, R). (4.39) We also note that since the reference polynomial P is already a very good ap- proximation of the mathematically correct fixed point. In practice the self-inclusion of the image can almost always be achieved by choosing sufficiently small domains T and D around the reference points to and 2:0. 4.2.3 Discussion of the Algorithm After having established the mathematical foundations of the algorithm in the previ- ous subsection, we will now comment on the individual steps of the basic method. In 97 the following discussion we focus on how each stage can be performed automatically in a computer environment, with only limited need for manual user interventions. 1. In the integration of explicit ODEs, the initial derivative 11:6 2 23’ (to) is uniquely determined by the systems initial conditions. In the context of implicit ODE on the other hand, the solution to an initial value problem is not guaranteed to be unique. Thus, the consistent initial condition 2:6 has to be obtained during a pre—analysis step. And since the consistent initial condition may not be unique, verified methods have to be used for an exhaustive global search. To simplify this, the user should be able to supply initial search regions for 326. In addition to the example given in Eqn. (4.4), the following example illustrates how the solutions of implicit ODE initial value problems are not necessarily unique. (45(4))2 + (4.2111(4))2 =1 and 4(0) = 0 (4.40) Obviously, s:._(t) = — sin(t) and 23+ (t) = +sin(t) are two distinct smooth solu- tions of the problem with different initial derivatives -1 and +1. Since the following stages of the algorithm require an expansion of the consistent initial velocity in the transverse variables, the use of high-order Newton methods or the inversion schemes discussed in Sec. 3.5 are warranted at this point. This will compute both a Taylor model and a suitable DA vector for the following stages of the algorithm. As part of the first step, the regularity of the J acobian matrix (4.20) should also be verified, and rigorous bounds on the sizes of the domains of regularity have to be established for t, 2:, and z’. Naturally, all subsequent steps then have to ensure, by possibly shrinking the domains for the independent variable t and the initial conditions, that all computed quantities stay within these upper bounds. 98 2. From a user’s perspective, it is important that by combining a suitable user interface with a dynamically typed runtime environment like COSY Infinity, the substitution of the variables with antiderivatives can be performed auto- matically, and there is no need for the user to rewrite any of the equations by hand. 3. By shifting to coordinates that are relative to the consistent initial condition 2:3, the solution space is restricted to the set M = {a E "D, : /\(a) 2 1} of nilpotent DA vectors. In step 6, this allows the definition of a DA—contracting operator, and the application of the DA Fixed Point Theorem. Again, this coordinate shift can be performed automatically within the semi-algebraic DA framework of COSY Infinity. 4. As in the previous two steps, the semi-symbolic nature of the DA framework allows the linear part L; to be extracted accurately and automatically, since the linearizations are computed automatically and are readily available in DA and Taylor model based methods. It should also be noted that the existence of the inverse over the domains of interest has already been established in the first step, since the regularity of LC follows from the regularity of the Jacobian (4.20). 5. Frequently the solution of a system can be specified as the fixed point of a sufficiently smooth operator. Here we have the derivative of the solution of an implicit ODE initial value problem given as the fixed point of such an Operator. Moreover, according to Thm. 4.1, the fixed point does indeed represent the desired solution. 6. According to Lem. 4.1, ’H is well defined and DA-contracting on the set M of nilpotent DA vectors. Hence the DA Fixed Point Theorem guarantees that the sequence (a,,) of iterates converges in at most n + 1 steps to the n-th order 99 solution polynomial. an+1 = P(t, v) =4 C(t, v). (4.41) It should however be noted that within a computer environment the implemen- tation and iteration of ’H finds only a floating point polynomial which is a fixed point of the floating point operator ’H. While the coefficients of the polynomial might differ from the mathematically exact n—th order expansion by the ma- chine epsilon, it is in fact sufficient to find a fixed point of ’H only to machine precision, since deviations from the exact result will be accounted for in the remainder bound of the Taylor models. . It has been shown that for explicit ODEs and the Picard operator ’P defined in Eqn. 2.44, inclusion is guaranteed if the solution Taylor model T satisfies P(T) C T. Although ’H differs from ’P, based on the proof sketched in [82], we have been able to establish a similar result in Thm. 4.2. While all previous steps are guaranteed to succeed whenever at least one consis- tent 2:[, can be found for which the linear part is regular, practical applications of this step of the algorithm can fail if no suitable Taylor model can be con- structed, although decreasing the size of the domains will generally lead to a successful inclusion. Details on suitable strategies for the construction of the Taylor model T can be found in [82, 24]. Lastly, it should be noted that this step of the algorithm requires a guaranteed computation of ’H, using Taylor model enclosures for 2:3 and L4- . This final step computes an enclosure of the solution to the original problem from the computed Taylor model containing the derivative of the actual solu- tion, and it relies on antiderivation being inclusion-preserving. To maintain full verification, the Taylor model enclosure of 2:3 has to be added to the result of 100 the integration. Similar to the example presented in (4.15), the new method also allows the direct integration of higher order implicit ODE initial value problems without explicit order reduction. While this approach does not reduce the dimensionality of the Taylor models if dependence on initial conditions is desired, it can however improve the sharpness of the final remainder bounds. Since the magnitude of the time domain T is usually smaller than one, Def. 2.8 shows that the remainder bounds of the actual solution will often be smaller than the ones of the computed highest derivative. To illustrate how the algorithm can be adapted to higher order ODEs, consider the general second order implicit ODE initial value problem 0(1), 2:, 29,42") = 0, 23(to) = 2:0, 2:’(t0) = 2:2,. (4.42) If we assume that the Jacobian matrix 8G (t, u, v, w) / 0w is non-singular in a suitable domain, this can be written as t r t (I)(t,§) = G (we +/to (2:6 + to {(o)do) dr,2:[, + to £(r)dr,§,t) = 0, (4.43) and the algorithm works with only minor adjustments. Similar arguments can be made for more general higher order ODEs. The performance of this approach will be illustrated in the next section with the direct integration of a second order system. 4.2.4 A Second Order Example To illustrate how the new algorithm can be used for the direct integration of higher order problems and to demonstrate how the method works in practice, consider the 101 implicit second order ODE initial value problem ex" + 2:” + :r = 0, (4.44a) 2:’(0) = 2:3 2 0. (4.44c) While the demonstration in this section uses explicit algebraic transformations for illustrative purposes, it is important to note that the actual implementation uses the DA framework and does not rely on such explicit manipulations. We also mention that for purposes of illustration, this example does not expand the solution in the transverse variables. 1. Compute a consistent initial value for 233’ = ”(0) such that e30 + x3 ’+ 2:0— — 0. A simple interval Newton method, with a starting value of 0, finds an enclosure of the unique solution 23’ = -1.278464542761074 in just a few steps. 2. Rewrite the original ODE in a derivative-free form by substituting C = 2:”. t r @(C, t) = 35“) + C(t) + (2:0 +f (2:3 +/ C(o)do) dr) 2 O. (4.45) 0 o 3. Define the new dependent variable C as the relative distance of C to its consistent initial value and substitute C = C — 2:3 in <1> to obtain the new function \II given by n t r ‘IJ(C, t=) C+ 2:3 ’+ eerC+1+ 1:2th +/ / C(o)dodr = 0. (4.46) o o 4. The linear part L((C) of ‘II is 1+ 633 where the 1 is the constant coefficient and e33 results from the linear part of the exponential function eC. 5. With LC from the previous step, the solution C is a fixed point of the DA- contracting operator ’H defined by 1 1 + 6%, (e353 (( — (3C) — 2:3’— —$23 -—t2-— [0 /0( C o() dodr). (4.47) 102 71(4): . Starting with an initial value of C (0) = 0, the n—th order expansion P of C is obtained in exactly n steps. 0"“) = a (0“) . (4.43) . The result is verified by constructing a Taylor model T with the computed reference polynomial P such that ’H(T) C T. With the Taylor model T = (P, 0, [0, 0.5], [—10-14,10—14]), (4.49) ’H(T) = (P, 0, [0, 0.5], [—0.659807722 . 10‘”, 0.659857319 - 10-14]) . (4.50) Since P is a fixed point of H, the inclusion ’H(T) C T can be checked by simply comparing the remainder bounds of T and ’H(T); the inclusion requirement is satisfied for the constructed T. . A Taylor model for 2: is obtained by using the antiderivation of Taylor models according to t r 2: (t) = 2:0 +/ (233 +/ (163’ + C (0)) do) dr. (4.51) o o The listing in Tab. 4.1 shows the resulting Taylor model of order 25 computed by COSY Infinity This example demonstrates how the new method can be used for the verified integration of implicit ODE initial value problems to high accuracy. In this context it is important to note that the width of the final enclosure of the solution is in the order of 10‘14 for a relatively large time step of h = 0.5. 4.3 Verified Integration of DAEs In this section we revisit the question of how a given DAE initial value problem can be transformed into an equivalent system of implicit ODEs. Since no general purpose 103 RDA VARIABLE: N03 25, NV= 1 I COEFFICIENT ORDER 1.000000000000000 -.6392322713805370 .4166666666666668E-01 .1993921404777223E-02 .6314945441169959E-04 .2635524930464548E-05 1 .4411105791086625E-06 12 .1533094467519992E-07 14 0.8104707776528831E-08 16 10 -.3384116382961162E-09 18 11 -.1389729003787960E-09 20 12 0.1981078695604361E-10 22 13 0.1549987273495670E-11 24 IOOIO OGGPN’O 1 2 3 4 5 6 7 8 9 VAR REFERENCE POINT DOMAIN INTERVAL 1 0.000000000000000 [0.00000, 0.50000] REMAINDER BOUND INTERVAL R [-.2500253775762034E-014,0.2500000000000003E-014] *ttttttastestastiestattests*ttttastes*ttttttttttttttttttt Table 4.1: Taylor model for the solution of the implicit second order ODE initial value problem given by Eqn. (4.44). methods for the direct integration of DAEs exist, this transformation is in general the only available method of computing solutions to the initial value problems of differential algebraic equations. We first summarize the main result of the E-method by J. D. Pryce [107, 108], which uses a structural analysis of DAEs to establish the existence of solutions. Ad- ditionally, the method finds an upper bound on the problem’s differentiation index Vd- Moreover, and more importantly, it also determines an explicit scheme for the transformation of the DAE system into a solvable system of implicit ODEs. However, within the previously presented algorithm for the integration of ODEs, the regularity of L4 already offers a sufficient criterion for the solvability of the derived ODEs. While the linear map LC will generally be singular, by repeatedly differenti- 104 ating the individual equations of the DAB, we will eventually obtain a regular linear map LC- Additionally, once the consistent initial derivative is computed, the minimum number of differentiations gives the differentiation index 114 of the DAE. 4.3.1 Structural Analysis Here we summarize the main results of the structural analysis of DAEs deveIOped by J. D. Pryce [105, 107, 108]. The signature method, or E-method, is used to decide whether a given DAE possesses a unique solution and to transform it into an equivalent system of implicit ODEs. The method considers the general DAE problem f1(t,2:1,...,2r[E“),...,:r,,,...,rff‘”) = 0 (4.52) f.) (t,2:1, . . . ,xfif'”), . . . ,2:,,,...,2:,(f””)) = 0. We recall that for a given j, the i-th equation of this system does not necessarily depend on 2:], the derivatives 2)") for 77 < 5,], or even any of its derivatives. However, if there is a dependence on at least one of the derivatives, including the zeroth order derivative 56].") = 22,-, we denote the order of that highest derivative by Cij. We formally define the v x v matrix 2 = (0.3:) by —oo if the j-th variable doesn’t occur in f: (4.53) 01'2“ — 5,.)- otherwise. The matrix 2 is called the signature matrix of the problem (4.52) and it forms the foundation of the structural analysis. We denote the set of permutations of length v by P” and consider the assignment Problem [9] of finding a maximal transversal T E P, of the signature matrix 2. Maximize ”T” = 23:10:2(1) With 0‘1',T(i) Z 0- (4-54) 105 Generally, an assignment problem is the task of matching v workers with v jobs. The entries of the assignment matrix 2 measure the effectiveness of the person i for the job j. An effectiveness of —00 indicates the inability of performing a particular task. By requiring the entries of the transversal T to be finite, we effectively limit the problem to the sparsity pattern S={@flbe>—my (4%) It should be noted that this analysis of the sparsity structure can be ill-posed and the assignment problem does not have any feasible solutions. More importantly, it has been shown by Pantelides that the assignment problem might be ill-posed even for solvable DAE problems. While this is one of the largest drawbacks of the structural analysis, the method works for many important DAE problems, including DAEs of index 0 and DAEs in Hessenberg form [108]. If a maximal transversal exists, we can consider the linear programming problem in the variables (a), (dj) 6 Z” defined by Minimize Z = Ed,- — Z c,-, (4.56) j=l i=1 subject to the restrictions dj-q20ijV(i,j)€S, c,20‘v’i=1,...,v. (4.57) While this problem does not have a unique solution, there is a uniquely determined smallest solution [108], and the smallest (Ci), (dj) 6 Z” are called the oflsets of the problem. Once the offsets have been determined, we form the v x v system Jacobian matrix J = (J,,-) with entries as 458 Jr = '————.'— ( ' ) J 106 where J,, = 0 if the derivative is not present or if d,- < cg. We will later see that in the framework of structural analysis the system Jacobian plays the role of the Jacobian matrix (4.20) in conventional index analysis of DAEs. Next we consider the collection of equations obtained by taking derivatives of the f, with respect to the independent variable t: fl“) = fl” = =- ff“) = o. (4.59) fl”) = fl” = ... = fl“) = 0. By definition of 0,,- and Eqn. (4.57), the derivatives of the 2:,- that appear in the equations (4.59) are all in the set {25(10),2:[1),.. ”2:90,. ..,2:S,0),2:§,1),. ..,2:(1d”)} . (4.60) If we write the set (4.60) as a vector X and collect all the equations (4.59) in the function F, we can write the derived system as F (t,X) = 0, (4.61) where F and X have M = v + Zc, and N = v + Z d,- components, respectively. Within the framework of structural analysis of DAEs, a consistent point of (4.52) is defined as a scalar to and a set of scalars C), indexed over a finite set I of indices (3', I) such that there is a unique solution of (4.52) near t = to such that (cg-(Rte) = 35,, However, the set I is not required to be a minimal index set, and it has been shown that it generally will not be minimal [116]. The concept of consistent points allows a simple success check of the structural analysis. If we view (4.61) as M algebraic equations in N variables and assume that the system has a solution (t‘, X *), then if the system Jacobian (4.58) is regular at the point (t‘, X ‘), it is a consistent point and 107 the system (4.52) has a unique solution near t = t“. Moreover, in a neighborhood of (t*, X ‘) the system has D degrees of freedom, where D is the common optimal value of (4.54) and (4.57), U D = urn = 2 = Zd, — 26,. (4.62) j=1 i=1 The structural analysis also provides an upper bound for the differentiation index Vd of the system. 0 ifalldj>0 1 otherwise. (4'63) udgmaxq+{ While it has been suggested that the index 14, might equal the given expression, examples in [116] show that (4.63) provides only be an upper bound. To utilize structural analysis for the numerical integration of the DAE (4.52) we use the subset of equations belonging to the highest derivatives in (4.59) and solve it as an implicit ODE initial value problem. Thus, we solve the v-dimensional ODE problem 1“" = --- = fl“) = 0 (4.64) to find the unique solution of the original DAE problem. Moreover, all the interme- diate equations (0) _ = f1(°1_1)=0 1 (4.65) is”) = = fir-1E 0 express the obvious and hidden constraints. These equations can in principle be used to verify initial conditions and make a posteriori adjustments to the computed final coordinates. Finally, it should be noted that the structural analysis is not infallible, since it relies on the presence of a substantial sparsity pattern. To illustrate this, consider 108 some DAE system F = 0 for which the structural analysis succeeds and the c,- are not all zero. If we then premultiply the system by a random non-singular matrix, the resulting DAE is equivalent to the original one. However, the signature matrix of the new system will generally be constant in each column, since the system is unlikely to exhibit any significant sparsity patterns. Thus, all the offsets will be c,- = 0 and the system Jacobian is identically singular. 4.3.2 Verified Index Analysis The applicability of the structural analysis summarized in the previous section is of- ten limited by its inherent limitation; namely, the E-method fails if the system has no obvious sparsity pattern. Moreover, it has been shown that index one problems can have arbitrary large structural and Taylor indices, rendering the method unpractical in these cases [116]. More formally, the fact that the index set I is not necessarily minimal, may lead to arbitrarily large indices and unnecessary computational over- head. Within the DA framework on the other hand, the regularity of LC in the neigh- borhood of a consistent point is a sufficient criterion for the existence and uniqueness of solutions. Moreover, since the regularity of LC is equivalent to the regularity of the Jacobian (4.20), this approach allows a rigorous determination of the differen- tiation index 11.1. In practical terms, this approach is simplified by the availability of the differentiation in the DA framework of COSY Infinity. After determination of a consistent initial velocity, a conventional DA analysis can easily determine the differentiation index of the system, provided it is in fact smaller than the order of the DA vectors. High order Taylor model methods can then be used to rigorously guarantee regularity of the resulting system Jacobian and determine the existence of a consistent point for the initial conditions. 109 Clearly this approach benefits from the availability of high order derivatives in the DA framework. However, to fully automate this method and avoid the need for hand- coded derivatives, it requires the additional availability of the derivative operation on Taylor models. Its utilization can completely eliminate the need for manual user intervention from the whole algorithm, which has repeatedly been recognized as one of the most important aspects of successful user interface design for verified integration schemes [135, 34, 65]. 4.4 Examples and Applications To further illustrate the concept of DAEs and to show the wide range of applications in which they arise naturally we discuss two important representatives of general DAE problems in this section. Moreover, as an interesting application of DAEs, we will demonstrate how prescribed path control problems can be formulated as DAE problems. In light of the Taylor model approach the method has distinct advantages over other existing methods for verified feedforward control problems; namely, a re- duction in computational complexity and a broader range of applications, since it is not limited to explicit first order systems. 4.4.1 Kirchhoff ’3 Law and Electric Circuits Differential algebraic equations are a natural way of modeling systems that obey intrinsic conservation laws. Frequently these systems can be written in either Hes- senberg Index-1 or Index-2 form. While we will encounter an index two problem in the next section, here we illustrate the concept of DAEs with an index one example from of electrical engineering. Consider the electric circuit shown in Fig. 4.1, which Combines resistors, capacitors, and transistors. 110 Figure 4.1: Illustration of an electric circuit described by the index-1 DAE (4.69). Ub is the operating voltage and U, is the driving voltage of the circuit. Ohm’s Law describes the simple linear relationship between the voltage U and the current I at a resistor of resistance R by U = R1. (4.66) At a transistor on the other hand, the voltage and the currents at base, emitter, and collector are connected by nonlinear relationships of the form IE = f(U) = 5 (BU/UP — 1) (4.67a) IC = —GIE (4.67b) I3 = (0l — 1) IE, (4.67c) where U p, 02, and H are characteristic quantities of the transistor. Across a capacitor C the connection between voltage and current is given by a simple first order differential equation I = CU’. (4.68) If we apply Kirchhoff ’s Law, which formalizes the conservation of current, to the 111 five points marked in Fig. 4.1, we get the following set of equations: U1-—Ue 0 = — (a — 1) f (U1 —- U3) (4.6931) R1 C(U; — Uj) = UbgzUz + “AU? — af(U1— U3) (4.6%) 0 =.- U3 ’ U0 —- f (Ul — U3) — f (U, — U3) (4.69C) R3 0 (U1 — 04) = (a — 1) f (U. — U3) — U4 1:12 (4.69d) 0 2 U4 — Uh + Of (U4 — U3). (4.698) R5 Thus, the response of the system to the driving voltage Ue(t) is determined by a differential algebraic equation with index one. The system consists of four algebraic constraints. Since the two differential equations are linearly dependent, they amount to one differential and one algebraic equation. It is easy to see that the problem has only one degree of freedom, since only the difference U2(0) — U4(0) can be chosen freely. This circuit is a typical example for a large class of engineering problems that have conservation laws at their core. Most, if not all of these dynamical systems are naturally described by differential algebraic equations, providing an almost infinite number of interesting and important DAE problems, many of which can and should be solved with verified methods. 112 4.4.2 Constrained Mechanical Systems Constrained mechanical systems can often be written as differential algebraic equa- tions and are typically of the form M (2:) -2:” := f (t,2:,x’) + G (2:) /\ (4.70a) 5(4) = 0 (47%) where 2: E R” is the state vector, A 6 R‘” is the Lagrange multiplier, and ’ = BQ/aq = G”. We will generally assume that M is positive-definite symmetric; so its diagonal is non-zero. Moreover we will also assume that ’ has full row rank, implying that v 2 w. Extensive theories have been developed that give cookbook recipes for solving these problems by using the constraint conditions (4.70b) to introduce new variables that reduce the problem’s dimensionality [49, 61]. However, these schemes usually rely, to a certain degree, on the user’s intuition and often require substantial arithmetic to reorganize and simplify the resulting ODEs into explicit first order systems. Moreover, while the use of the constraints (4.70b) does lead to simplified ODEs, it generally does not reveal the hidden constraints of the system. Within the context of differential algebraic equations, the regular and the hidden constraints are easily accessible. Moreover, the availability of automated integration schemes for DAE initial value problems removes the need for intuition and physical insight into the problems and allows the automated integration of these important Systems without user intervention. Example: Double Pendulum As a prototypical example for constrained mechanical systems, consider a planar pair of connected pendulums in a frictionless environment. Assume that the pendulums 113 l, Figure 4.2: Illustration of a double pendulum with masses m1 and me connected by massless rods of lengths l1 and [2, respectively. are massless and inextensible, with point masses on the ends, as illustrated in Fig. 4.2. If we denote the tensions in the strings by A1 and A2, they take the roles of the Lagrange multipliers in this problem and the equations of motion, expressed in the Cartesian coordinates 2:1, y1,2:2, ye, are given by 771133,; + MIDI/l1 — 42(332 " 370/12 = 0 7712in + 4191/51 - 42(92 — yll/l2 — mlg = 0 7112.72"; + /\2(.’132 — 330/12 2 0 (4.71) mzvé' + 42(92 - y1)/lz - 0229 = 0 2:? + y? — l? = 0 (2:2 — 11:1)2 + (92 — yl)2 — 13 = 0. While it is common practice in classical mechanics to reduce this problem to a two dimensional ODE in the variables (t, and 452, here the focus is on treating this problem as a differential algebraic initial value problem in the six variables 2:1, y1,2:2,y2 and A1, Ag. The 2 matrix of the system (4.71), with the entries forming a maximal transversal 114 in bold face, is given by 2—10—10 0( —12—1 0 0 0 0—1 2—1—10 —-1 0 -1 2 —1 0 ' (4'72) —1 -1) O O —1 -I -—1 Further analysis gives the offsets c = (0, 0, 0, 0,2, 2) and d = (2,2,2,2,0,0) and (0000—1 reveals that the system has four degrees of freedom, a structural index of three and a differentiation index Vd of two. Based on this analysis, the implicit ODE problem that can be derived from the DAE (4.71) is given by mir’,’ + Alan/l1 — A2(2:2 — 2:1)/12 = 0 ”123/’1’ + Mill/l1 — 42(112 — fill/l2 — 77119 = 0 mgx’z’ + /\2(2:2 —- 2:1)/lg = 0 N (4.73) 7712312 + 42(92 — yll/lz — ng = 0 2 (2’12 + 2:12:3’ + y’12 + yiyi’) = 0 2 ((4; -- xi)’ + (y; — yl)’ + (42 — 2:1)(23’2’~— 43') + (42 — y1)(y3’-—- 3,9) = 0, Additionally, the complete set of obvious and hidden constraint conditions for this system is given by 2:? + y? — if = 0 ($2 - 51302 + (312 — 3102 - lg = 0 (4.74) 2($13314- 3131i) = 0 2 ((422 — 2,) (2’2 — 2’1)+(yz — 91) (vi _ 91)): 0 The first two constraints express that the lengths l1 and l2 are constant, and the 115 -o.4 ' Y I . l -o.2~ -. ” l .i I: :2 5‘. I~ - ..., -‘ .. .. . -. .- :. .3 .. .. .. .. .. _. "_ ‘. .. ,. .. .. .. ;. ,4 1;; (~; ..1l .. >1 .. ,n .- , - -, .. 1- a .. ., .. _. .. _. . ,, ., .1.l ., .. .. .- ,:'-‘.>‘ .... 41.. :. .. _. 1; '...._‘ .. .1- .. . .. .. . '..:.,,V ., ., .~ .. .. .. _. .. '.. . ,. . .. .r .u - W! :3.:.-:... .. 4.. ~. a 5‘ h .. A. '.~ ‘. ~04 l 004 -v-1-‘.- 4. .. .. :..~. .4 ...,;:,: .. .:.. ,,: fifj.v'-11 l-?. 'l'. :'. .v'. .. v ; ’ '. '- ., .‘ . ~ a. ....:.,»:. _,‘. .--. ,.‘- .,.~ .... ... _, _,‘. . . .. .,.- 111.... .-- U. ..._ _.:. .,.- . ..,..~-.--. ,'. ‘vv.. _... .,.v ...-.~-.--. .'. -.8.. 4H -.' ..,.,.-....-.<, -, ... '....'.. .,,,,..‘..,-.-, ~, ..,4,.._-.. ..,...-....;,-. .....--,.... ~. '~ - .:.«.. ,-. 1 .. - : ' - _A__P_A_ ;_L_L _A_—l g L4 0102030405060708090100 0102030405060708090100 Figure 4.3: Coordinates 2:1 and y1 of the double pendulum for 221(0) = rp2(0) = 5°,30°, 90°, 2’, (0) = 23(0) = y[(0) = y3(0) = 0 and 0 S t _<_ 100 (integrated in order 6,h=0.1,g=1, l1=12=1, m1=m2=1). hidden constraints reveal that the velocities of the masses are tangent to the circles on which the masses move. To illustrate the motion of the double pendulum, Figs. 4.3 to 4.5 show results obtained from various integrations of the double pendulum with different initial val- ues. All plots highlight the oscillatory motion of the system with an energy transfer between the two masses. Fig. 4.3 shows the coordinates of the upper mass ml for various initial energies. While the system follows an almost periodic path for small energies, it shows signs of chaotic motion for high energies. The phase space projections of the second mass m2 shown in Fig. 4.4 indicate that the system moves on almost periodic orbits on its oscillating orbit. Finally, Fig. 4.5 shows the time evolution of the coordinates 2:1 and y1 for the case where the upper mass is 100 times heavier than the lower one. In this case, we see slow oscillations superimposed with the small and fast oscillations of the lower mass. Moreover, the evolution of y1 indicates slow but noticeable energy transfers between the two masses. 116 X2092 Yz’Y'z 2 ~ 1 ’ l 1 l 0.5 r .......... .. .“_...:-.;;‘7$;”#"”' ii“, .3 52),”: ° ‘ or (91/ "l§?‘<;\«;;\ ' unit“: -1 h '0.5 L ‘‘‘‘‘‘‘ “3. 2:1...” .2 . l .1 _ -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 Figure 4.4: my and yg-phase space projections for the double pendulum with for 901(0) = 092(0) = 5°,30°,45°, 221(0) = 103(0) = y{(0) = y§(0) = 0 and 0 5 t g 100 (integrated in order 6, h = 0.1, g 2 1, l1 =12 =1, m1 = m2 = 1). X1 Y1 0.015 f f . 0.9999 . . . . . . . n 0.01 lv 1 0.99992 U U U n n n U n °-°°5 0.99994 1 ° ' 0.99996 l , -0. 005 l 0.99998 + . -0.01 l I U U U U w 1 0.015 g a 4_- e 4 0102030405060708090100 0102030405060708090100 Figure 4.5: Coordinates $1 and yl 0f the double pendulum for 901(0) = 0’ “22(0) = 450’ 51(0) = 332(0) —_— y'1(0) _—_ 315(0) = O and 0 _<_ t g 100 (integrated in order 6, h = 0.1, 9:1,llzl221,m1=100,m2=1) 117 4.4.3 Optimal Prescribed Path Control Prescribed path control problems for dynamical systems are traditionally solved by using feedback control methods. Feedback control is based on the idea that we can let the system operate with a pre—determined control function A for a small time interval h and adjust the control function at the end of the time interval, depending on the deviation of the system from its prescribed path. While feedback control is usually easy to implement, its applicability is often limited to non—stiff problems that vary only slowly over time. However, it is also possible to handle prescribed path problems in feedforward mode, where the appropriate control function is determined according to the desired future behavior of the system. Feedforward control tends to be more accurate in following the prescribed path with larger time steps h. Thus, feedforward control can often result in more accurate modeling with less computational cost. If we denote the independent variable by t E [t,-, t,+1], the state of the system by m(t) and the control input by A(t), the general prescribed path problem can be written as F(t,x,x', A) = 0 Vt E [t,-,t,~+1]. (4.75) Using the Taylor model approach for the integration of implicit ODES, problems of this form can be solved for the system’s state as a function of both the independent variable t and the control function A. Thus, we can obtain a Taylor model with an explicit dependence on the control function A, m(t, A). (4.76) Then, given a. prescribed path Mt), we define a new function, A = A(t, A), as the difference of the solution m(t, A) and the prescribed path 112(15). Mt, A) = m(t, A) — ¢(t)- (4.77) 118 The function A can be viewed as an implicit description of the control function A(t) such that the system follows the prescribed path. We may assume that m(t,) = 2/2(t,-) and that we have an initial control A(t,~) such that (éég—Al) (t,, A(t,-)). (4.78) is regular. Then, the Implicit Function Theorem guarantees the existence of a function A A(t) such that A(t, 5\) = 0 (4.79) in a neighborhood of ti. Since the partial inversion discussed in Sec. 3.5 allows the de- termination of a Taylor model for the optimal control function A from this expression, at this point the problem of verified prescribed path control can be considered solved. However, the potentially large number of additional variables required in solving the implicit equation (4.79) often renders this method unfeasible. With the availability of general DAE integration tools on the other hand, we can view (4.75) as a DAE problem and use methods for the verified integration of high order DAE initial value problems to compute appropriate control functions for the system (4.75). In comparison to traditional methods of prescribed path control, the virtue of this method lies in the fact it requires only a single ODE integration for a rather long time step. Existing methods for prescribed path control optimize the control function through a series of repeated integrations or use extremely short time steps with the assumption that the control function stays constant over a single integration interval. The DAE approach on the other hand, allows the determination of a smooth and explicit expression for the control function to very high order. Moreover, by utilizing Taylor models at each step of the algorithm, the new method even allows the computation of rigorous error bounds on the computed control functions, providing the operators of the system with guaranteed and trustworthy results. 119 Example: Two-Link Robot Motion To illustrate the formulation of prescribed path control problems as DAE initial value problems, the following system describes a slightly simplified version of the equations for the prescribed path control of a two-link robot arm [107]. With sufficiently smooth and highly non-linear functions f1, f2, and f3, the system has a differentiation index of 5 and is given by x’l' = fl (t, 2:2, :33, 3’1, 223,111,112) (4.80a) 233' = f2(t,:c2,:r3,:r'1,zg,u1, U2) (4.80b) 17;: = f3(t, $2, $53,333, $3,111,112) (4.80c) 0 = cos(a:1) + cos(:l:1 + 11:3) — cos(e‘ — 1) — cos(t — 1) (4.80d) 0 = sin(xl) + sin(scl + x3) -— sin(l — e‘) — sin(l - t). (4.80e) The prescribed path is expressed in the constraint conditions (4.80d) and (4.80c) and requires valid solutions 3:1 and $3 to be of the form 271(1) = 1 — et ' (4.81a) x3(t) = e‘ —t (4.81b) 4.5 Summary We have developed a method for the verified integration of general ODE initial value problems. It is based on a combination of high order Taylor models and the an- tiderivation, and it allows the automatic computation of tight solution enclosures for a large class of solvable initial value problems. The new approach avoids the wrap- ping effect to very high orders by expanding the solution not only in the independent variable t, but also in the transverse variables. Thus, the method computes accurate 120 enclosures of the flow of ODE systems. Additionally, by using the method, or alter- natively a more traditional structural analysis, for the index-analysis of differential algebraic equations, it can be extended to the integration of the large and important class of DAE initial value problems. Finally, we point out that although we deliberately restricted the discussion to DAE problems with constant indices, many interesting practical problems have in fact non-constant index: valves in chemical processes may open and close, surfaces in mechanical devices may touch. However, a change in index is usually also associated with a non—smooth change in the equations, which is something the Taylor model approach cannot handle well. But the method can still be used in these cases by dynamically changing the equations between time steps and using the Taylor model describing the final coordinates at the end of the previous step as the initial condition for the next step. 121 Chapter 5 Applications of High Order Taylor Model Methods In this chapter we present applications in which the availability of Taylor model meth- ods allows the use of computers to obtain rigorous results. These applications require the high order Taylor model approach for two reasons: standard floating point meth- ods fail to provide the required mathematical rigor, and conventional interval tech- niques fail to adequately address the problems of wrapping and dependency. Since the systems discussed in this chapter are not necessarily available for direct experiments, the Taylor model approach is the only method to study them rigorously. First we study the verified integration of near-earth asteroid orbits. This problem has long been seen as a significant challenge for interval techniques, since it involves large initial sets, long integration intervals, and highly non-linear equations of motion, resulting in significant wrapping and dependency problems. However, we will show that the Taylor model approach can overcome all of these obstacles and return prac- tical results. Unlike in the case of conventional floating point methods, the computed results are guaranteed enclosures of future positions of the asteroids, accounting for the effects of uncertain initial conditions, measured data, and mathematical trunca- tion. 122 In a second application we use high order Taylor models to prove the existence of generating functions for Hamiltonian systems over extended domains. The approach combines the mathematical rigor of the invertibility tests discussed in Sec. 3.1 with Taylor model based integration of initial value problems. This novel combination pro- vides the rigorous justification for the method of symplectic tracking of Hamiltonian systems, which has important applications in the study and the design of particle accelerators. 5.1 Solar System Dynamics and Asteroid Orbits Spurred by recent media attention, the possibility and implications of asteroid im- pacts on Earth have recently come to the public’s attention. Indeed, the potentially devastating consequences of such events justify research in preventing them, including possibly by altering the orbits of near-earth asteroids on collision courses. But pro- tective measures cannot be taken unless near-earth asteroids are detected and their orbits are predicted. Moreover, for an increased response time it is essential to detect hazardous asteroids as early as possible and perform long-term integration of orbits from known observational data. Since the results of these predictions will potentially justify substantial preventive efforts, it is of utmost importance to actually get them right. Thus, the integration of asteroid orbits from initial measurements is an obvious application of verified ordinary differential equation (ODE) integration methods. Due to the computationally challenging models and the large uncertainties, however, so far all such computations have been performed using conventional high order floating point integration schemes with subsequent non-verified error analysis [73]. The integration of the motion of asteroids in the solar system poses several chal- lenges to verified integration methods, mostly due to the large uncertainties in the 123 initial conditions, measured physical constants (e.g. the gravitational constant is one the least accurately known constants [62]), and strong non-linearity in the mathemat- ical models. Therefore, verified tools have not found widespread application for these problems. However, due to the potential seriousness of the outcome of the integration, they are important applications and test cases of verified methods [90, 135]. The major problem of integrating solar system dynamics lies in the relatively large uncertainties in the positions and velocities of the objects involved. This problem particularly hurts conventional verified integration schemes that are prone to the wrapping effect discussed in Sec. 2.2.2. It leads to a substantial overestimation on the range of final positions and velocities, and this overestimation often results in situations that do not permit reasonable conclusions, since the bounds are overly pessimistic. The other main obstacle to verified integration of solar system dynamics lies in the strong non-linearity of the equations of motion (5.3). Together with the uncertainties in the positions and velocities, this lets conventional verified integration suffer from cancellation problems and makes it almost impossible to obtain usable long-term integration results. We will demonstrate, however, that the use of Taylor model based ODE integration methods [83, 85] can successfully alleviate these problems. In fact, the Taylor model approach can indeed be used to perform verified long-term integration of asteroid orbits in the solar system. We will demonstrate the performance of the method with several examples obtained from orbit integrations of the asteroid 1997 XF11, which is predicted to have several nearby passages of Earth between the years 2000 and 2040. 5.1.1 Physical Background Before we discuss any details of rigorous asteroid integrations, in this section we summarize the physical principles that govern the motion of planets and asteroids in 124 the solar system. We give an overview of the solar system and present an accurate formulation of the equations of motion for a point mass in it. Some of the notation introduced deviate from the standard conventions used in the International System of Units (SI) and will be used throughout the text. Physical Units The basic physical quantities relevant for the mathematical description of planetary motion are length, time, and mass. They are usually expressed in terms of the SI units Meter, Second, and Kilogram. Quantity SI Unit Symbol Length Meter In Time Second 5 Mass Kilogram kg Table 5.1: SI units for length, time, and mass. For practical calculations it is sometimes convenient to use different units than the ones listed in Tab. 5.1, mostly to avoid excessively large or small values of phys- ical quantities. The following list gives a short description of the commonly used astronomical units of time, mass, and length: 0 The astronomical unit of time is the time interval of one Earth-day; an interval of 36525 days is called a Julian century. 0 The astronomical unit of mass is the mass of the Sun. 0 The astronomical unit of length (AU) is the radius of a circular orbit in which a body of negligible mass, and free of perturbations, would revolve around the Sun in 365.25692 days, and its value is almost equal to the mean distance between the Earth and the Sun [125]. 125 The rules required for a conversion of these units into the SI system are given in Tab. 5.2. Several other physical quantities, namely the velocity, the acceleration, and the force), are expressed in derived units; Tab. 5.3 lists them in terms of the fundamental quantities length, mass, and time. Quantity Unit Symbols Conversion Time Day d 1 d = 86400 3 Julian Year J 1 J = 365.25 (1 Mass Sun Mass M 1 M = 1.9891 x 1030 kg Length Astronomical Unit AU 1 AU 2 1.496 x 1011 m Table 5.2: Astronomical units and their conversions to SI units. Quantity Unit Velocity Length/ Time Acceleration Length/Time2 Force Length x Mass/Time2 Table 5.3: Units of the derived quantities velocity, acceleration, and force. Astronomical calculations also make frequent use of angles, which are commonly measured in degrees (360 deg is equivalent to a full circle) or radians (27r rad E 360 deg). Finally, practical computations often require the use of ad hoc units to satisfy requirements of the algorithms used, and these units will be introduced as needed. Gravitation Gravitation is an attractive force between masses, and for the purpose of gravitational interactions in the solar system, all masses can be treated as point masses; i.e., as if their mass was concentrated in the center of mass. Classically, the force inflicted by a point mass mg on a mass m1 acts along the line connecting the two points and is 126 given by (5.1) where G denotes the gravitational constant and F12 is the vector connecting ml with m2. Eqn. (5.1) is due to Newton, who derived it from observational data of planetary motion and Kepler’s Laws. Since E12 2 -—F‘21, gravitation satisfies Newton’s Third Law: actio = reactio. Predicting the motion of two point masses m1 and m2 under their mutual New- tonian gravitational influence, neglecting all other forces and interactions, is known as the two-body problem and has been studied extensively. It can be shown that the barycenter of the two masses moves with constant velocity and that the two masses move on one of the conic sections circles, ellipses, parabolas, or hyperbolas with the common barycenter in the focus [49, 61]. Consequently, in a barycentric coordinate system, the two masses move on planar conic sections around the origin. Details of the motion on conic sections will be discussed in Sec. 5.1.2. Within the framework of classical celestial mechanics, the superposition of the gravitational force is linear. Hence, the total force ET on a mass M induced by N other masses m,- is merely the vector sum of the individual forces 13,-z N _‘i M 3 -0 FT = 21:14}: 2:1|:G|TTUL—3T. (5.2) A more general and more complete model for the force on a point mass under the influence of other masses will be discussed in Sec. 5.1.1. The Solar System The solar system is the dynamical system consisting of the objects that are confined to orbits around the Sun. The largest members of the solar system are the Sun, 127 the nine major planets, some of which are orbited by moons, and a large number of asteroids and comets. Fig. 5.1 shows the orbits of the four innermost planets: Mercury, Venus, Earth, and Mars. As shown in Tab. A.1, the orbital ellipse of Mercury has a noticeable ec- centricity. Fig. 5.2 shows the orbits of all nine major planets; the m/y-plane coincides Y a , ye Figure 5.1: Orbits of the inner planets Mercury, Venus, Earth, and Mars around the Sun. with the ecliptic; i.e., the plane of Earth’s orbit, of the year 2000. The up-to-scale inclusion of the orbits of the inner planets illustrates the size of the solar system: Pluto’s perihelion is almost 50 astronomical units away from the Sun, and it takes the light from the Sun almost seven hours to reach Pluto. The Sun is by far the heaviest and largest object in the solar system, as its mass of approximately 1.9891 x 1030 kg exceeds the masses of the other planets by at least three orders of magnitude. As an important consequence of this dominance, the maximal distance between the solar system barycenter and the Sun’s center of mass is only 1.5 x 109 m, or about twice the radius of the Sun. For most practical purposes it is therefore sufficient to assume that all planets orbit the fixed Sun. Tab. 5.4 lists the masses of the nine major planets with their moons as fractions of the Sun mass, their mean distance to the Sun, and the average times needed for one revolution around 128 \t ...) Figure 5.2: Orbits of the outer planets Jupiter, Saturn, Uranus, Neptune, and Pluto around the Sun. the Sun. It should be noted that, according to Kepler’s Third Law, the period T and the mean distances a are related via T = am. Planet Mass (m) Dist. a a3/2 Time T Mercury 1 / 6023600 0.39 0.2436 0.240 Venus 1/408523.5 0.72 0.6109 0.615 Earth 1/328900.5 1.00 1.0000 0.999 Mars 1 / 3098710 1.52 1.8739 1.881 Jupiter 1/ 1047.355 5.20 11.858 11.856 Saturn 1/ 3498.5 9.54 29.466 29.424 Uranus 1 / 22869 19.19 84.064 83.747 Neptune 1 / 19314 30.07 164.89 163.723 Pluto 1/ 3000000 39.48 248.07 248.021 Table 5.4: Masses, mean distances, and periods of the major planets. Masses are given as fractions of the Sun mass M, mean distances in astronomical units, and times in Julian years [125]. In addition to the Sun, the nine major planets and their moons, the solar system contains numerous other small objects, mostly asteroids and comets, the largest of which are commonly referred to as minor planets [134]. The sixteen largest asteroids in the solar-system have diameters between 240 km and 1000 km and are part of 129 the main asteroid belt, which consists of several thousand objects of non-negligible mass and size and is located between the planets Earth and Mars. During past space craft missions, the asteroids Ceres, Pallas, Vesta, Iris, and Bamberga, all of which are orbiting the Sun in the main asteroid belt, have been found to have a noticeable gravitational effect on the motion of objects in the Earth-Mars range [125]. Hence, an accurate formulation of the equations of motion of objects on near-earth orbits has to include their gravitational effects. We will revisit this when we discuss an accurate mathematical model for the motion of near-earth asteroids in Sec. 5.1.1. Finally, it should be noted that the term asteroid belt is misleading, since it is by no means a dense object. In fact, several spacecrafts have passed it undamaged and the mean distance between non-negligible objects in the main asteroid belt is several million kilometers. The Equations of Motion In the previous sections we have given an overview of the solar system and the classical description of gravitation. Now we present an accurate mathematical model [139, 125, 44, 76] describing the motion of point masses in a solar system barycentric coordinate system by a second order differential equation. The model includes Newtonian effects of the Sun and the major and minor planets and relativistic corrections to the classical Newtonian description of gravitation of the Sun and the major planets. r}, 73,-, and if,- denote the solar-system barycentric position, velocity, and accel- eration of the j-th major planet, where j = 0 denotes the Sun. 7 and B are the post-Newtonian parameters measuring space curvature produced by the unit rest mass and nonlinearity in the superposition of gravity; in accordance with general relativity, both are assumed to be one. The masses have been combined with the gravitational constant into the coeflicients u,- = G'mj, and the position of the point 130 mass is given by 7'". The distances between point mass and major planets are expressed by r,» = Hr —- 773]]2, the respective velocities are denoted by v, = H153”2 and v =2 “F”, N is the number of minor planets that are included in the computation and c is the velocity of light. 9 N 3 4 u .1:— + + 7 affixing? (5.3) This differential equation shows that the force on a point mass is essentially given by the Newton force inflicted by all other major objects in the solar system and some correction terms in the order of 1/c2. Although Eqn. (5.3) gives only a partial description of the acceleration exerted on a point mass in the solar system, it provides a highly accurate description of the motion, and it has been used for the preparation of the most accurate ephemeris system so far [125]. 5.1.2 Geometric Description of Kepler Orbits As mentioned earlier, the Kepler problem of predicting the motion of two point masses under their mutual Newtonian gravitation permits only conic sections as orbits, with the barycenter being the focus. The orbits lie in planes that are perpendicular to the SYstem’s angular momentum and pass through the barycenter. In general, the shape 131 of conic sections is determined by the polar equation, which relates distance to the barycenter r and the angle in the plane of motion «p via r (1 + scos((,0)) = a (1 -— 52). (5.4) Depending on the value of the eccentricity e, the orbit is given by one of the following geometric figures: Eccentricity Orbit s = 0 Circle 0 < s < 1 Ellipse e = 1 Parabola s > 1 Hyperbola Table 5.5: Geometric shape of conic sections, depending on the eccentricity 5. In this section we summarize how Kepler ellipses can be described in terms of orbital elements and how this description can be used to compute the positions and velocities of planets in a barycentric cartesian coordinate system whose :r/y-plane coincides with the ecliptic of the Earth. Orbital Elements The shape and orientation of Kepler ellipses in the solar system can be described un— ambiguously by six parameters, the so-called orbital elements. While several equiva- lent choices for these parameters exist, the following list shows one particular selection that significantly simplifies the subsequent computations: O Longitude of ascending node (2: The angle, along the ecliptic, from the vernal point to the ascending node, which is the intersection between the orbit and the ecliptic, where the planet moves from south of to north of the ecliptic; i.e., from negative to positive latitudes. 132 o Inclination i: The tilt of the orbit plane relative to the ecliptic. Argument of perihelion w: The angle from the ascending node to the perihelion, along the orbit. Semi-major axis a: The mean distance to the Sun. Eccentricity e: The parameter describing the shape of the orbit (c.f. Tab. 5.5). Mean anomaly M: The product of the mean motion of the orbiting body and the time since the body passed the pericenter at some specific point in time. The mean anomaly is best described if one imagines a small circle within the orbit, and centered on the focus f. Then M is the angle between the perihelion P, the focus f, and the position the planet would occupy if it was traveling at a constant angular speed. While the last three elements describe the orbit ellipse themselves, (2, i, and w describe the orientation of this ellipse in space. Other equivalent choices of the orbital elements include the mean longitude L = M + w + Q, the time at perihelion, and the mean daily motion 17.. For reference, Appendix A lists the orbital elements and daily rates of change of the major planets. Cartesian Coordinates Since the equations of motion (5.3) describing the acceleration of a point mass in the solar system are only valid in a cartesian system of inertia, for practical computations it is usually necessary to convert orbital elements into cartesian coordinates. This transformation is based on the eccentric anomaly E, which is due to Kepler and connects the cartesian coordinates of a planet with its orbital elements. Consider a large circle of radius equal to the semi-major axis of the orbit. Let the center of this circle be at c. Now, on the diameter of the circle which passes through 133 the long axis of the ellipse, raise a perpendicular which meets the planet at r on the ellipse. Continue this perpendicular until it crosses the large circle and call the point where the perpendicular crosses the large circle Q. Then the angle E is called the eccentric anomaly (c.f. Fig. 5.3). The mean and eccentric anomalies are related by Figure 5.3: Illustration of the eccentric anomaly E the implicit equation MzE—esmm. am If the eccentricity 5 is small, a Newton method can easily find E for given values of M and 6. Once E has been determined, the cartesian positions F = (51:, y, z) and velocities i" = (is, g, z") are computed from the following equations [118, 93]: G(M+m) 1 n = a (l — scos(E)) (5.6a) a = a (cos(E) — 5) (5.6b) fl = a\/1 — 52 sin(E) (5.6c) fy = —nsin(E) (5.6d) 6 = nx/l — 52 cos(E) (5-59) 134 p. = cos(w) cosm) —- sin(a) sin(w) cos(i) p, = cos(w) sin(Q) + cos(fl) sin(w) cos(i) p, = sin(w) sin(i) q. = — sin(w) cosm) — sinm) WM 0080‘) q, = — sin(w) sin(a) + cos(fl) cos(w) cos(i) Q, = cos(w) sin(i) x=a~px+59x y=a-py+6-qy z=a-pz+fi-qz i=7-pz+6-qx UZV‘Py'l’a'Qy 2:7'pz'l'6‘Qz (5.7a) (5.71.) (5.7c) (5.7a) (5.7e) (5.70 (5.8a) (5.8b) (5.8c) (5.8a) (5.8e) (5.8f) The origin of this coordinate system coincides with the barycenter of the combined two-body system of the Sun and planet, but the coordinates can easily be transformed into a. solar system barycentric coordinate system by translation. Planetary Motion and Perturbations For the two-body problem, the six orbital elements are constants and are determined uniquely by the six initial conditions of the position i" and the velocity i". Due to the overwhelming mass of the Sun and the large distances between the planets, even in the presence of other planets, the orbits of the planets are very close to being ellipses. Minor corrections to the two-body Kepler ellipses of the planets are due to: o influences from planets other than the Sun, 135 o deviations of the force from the simple inverse square law (relativistic correc- tions) . It is convenient to treat the orbits as instantaneous ellipses whose parameters are defined by the instantaneous values of the position and velocity vectors, since for small perturbations the orbits are approximately ellipses. However, these perturbations cause the six formerly constant parameters to vary slowly, and the instantaneous perturbed orbits are called osculating ellipses: the osculating ellipses are the elliptical orbits that would be assumed by the bodies if all the perturbing forces and corrections were suddenly turned off. This perturbation theoretical approach of describing planet orbits relative to slowly varying elliptical reference orbits has proven itself extremely successful in parti- cle optical systems [17, 87], and based on the arbitrary order approach of COSY Infin- ity [23, 26], sufficiently accurate smooth description of planet orbits can be computed even for long-term integrations. Moreover, in light of the Taylor model approach, this strategy has the particular advantage that the orbital elements of planets are almost linear functions of time, which eases the description by Taylor models and increases the accuracy of the enclosures. 5.1.3 Verified Orbit Integration In the previous sections we summarized the concepts that govern the motion of masses in the solar system and presented the principle of describing nearly closed orbits by perturbations to osculating ellipses. Here we show how these concepts can be used for the verified integration of asteroid orbits. Since this involves highly non-linear equa- tions, large sets of initial conditions, and substantial uncertainties in the mathematical description, we have found conventional verified integration schemes [79, 99] unable 136 to handle the problem of creating a guaranteed ephemeris for asteroids in the solar system. As indicated earlier however, Taylor model based integration schemes [83, 85] can offer accurate long-term verification even under these challenging conditions. To illustrate this, we have based the integration on Taylor model methods, Eqn. (5.3), and the ephemeris DE405 [130]. The basic idea of the integration is that the presence of a relatively small asteroid will not change the orbits described by the DE405 ephemeris. While this is clearly only an approximation, due to the mass ratio between the planets and asteroids and due to the relatively large distances separating them, the deviations from this approx- imation are in fact extremely small and will be accounted for through general error bounds on the planets’ positions and velocities. Thus, instead of integrating the full (60 + 3 x N )-dimensional system involving the Sun, the nine major planets, and the N minor planets as Eqn. (5.3) would suggest, we use osculating ellipses together with the corresponding errors of the ephemeris to obtain Taylor models for the positions and velocities of the planets. These Taylor models describing positions and velocities of the planets and interval enclosures for effects not described by Eqn. (5.3) are then used for a Taylor model based integration of the six-dimensional equations of motions for the asteroid alone. This particular approach allows a substantial reduction of the problem’s dimen- sionality and makes its verified integration feasible. Utilizing high-order Taylor mod- els limits the wrapping effect and enables us to obtain long-term predictions, even for large sets of initial conditions. Moreover, the approach does not lead to a dramatic inflation in the size of enclosures of the final coordinates. In the remainder of this section we present a detailed discussion of how the verified integration of Eqn. (5.3) can be achieved and results will be presented in Sec. 5.1.4. 137 Accuracy of the Ephemeris The mathematical models and concepts presented so far allow the computation of accurate planet orbits. For the computation of the ephemeris DE405 [125, 130], which has been used to compute the motion of asteroids, rockets, and satellites in the solar system [142], the following effects have also been included: 0 Newtonian gravitation of the Sun and the major planets o relativistic effects of orders up to 1 / 62 o Newtonian effects of the 300 largest minor planets Since planets are not rigid spheres with a uniform mass density, an accurate inte- gration includes the deviations from the simple point-mass concept. The following factors have been included in the preparation of the DE405/LE405 ephemeris: o the difference between geometric centers and barycenters (figure effects) 0 Earth tide effects 0 physical librations of moons Finally, the results of the numerical integration have been adjusted to past measure- ments. Measurements are, among other things, based on the following techniques: 0 meridian transits (i.e., optical triangulation), 0 satellite astronomy, 0 radar ranging to planet surfaces (including surface models), 0 radar ranging to spacecrafts, 138 0 laser ranging to lunar reflectors. Residuals of the ephemeris DE405/LE405 have been presented in [130] and it has been shown that the ephemeris DE405 is indeed suitable to accurately predict the positions and velocities of massive objects in the solar system [142]. Since Jupiter has the second largest mass of all objects in the solar system, its positional errors translate into potentially large errors in the integration results of asteroids and its residuals with respect to DE405 are therefore of great importance. Fig. 5.4 shows these residuals with respect to DE405, obtained via Very Long Baseline Interferometry observations from the Galileo spacecraft. VLBI residuals of Jupiter w.r.t. DE405 0.05 r *— fi— *fi 1— f v— r“ 1 mm f fi—fi— fi— fi— fi fir —' .. J, l '3 ] [ ] l l g o Mfigfi f_‘L1 M] l r; rd j1_ hf Ir— L[‘ll fl M4 f LL L i l T l [ U1 11 l ’ . 1 8+ _L l M M l M l M 1* l ‘? 96 1996.5 1997 1997.5 1998 Figure 5.4: Residuals of Jupiter; two Very Long Baseline Interferometry observations from the Galileo spacecraft; courtesy of [130]. The residuals shown in Fig. 5.4 translate into positional uncertainties of only about 100 km. Combining this with even higher accuracies for the inner planets (the ephemeris system DE405 is believed to be accurate to 0.001 arc seconds [130]), it is possible to predict the positions and velocities of the Sun and the major and minor planets quite accurately for the next 50 to 100 years. As a result, the data of DE405 can be used to obtain time-dependent descriptions of the orbital elements. These can then be used to compute Taylor model descriptions of the barycentric cartesian coordinates of the planets that are sufficiently accurate for verified long- term integrations. Consequently, Taylor model integration schemes can be used to 139 obtain verified long-term results for the motion of small objects in the solar system. In addition to the effects included in Eqn. (5.3), our integration also considers the acceleration due to the rotation of the solar system barycenter around the galaxy; i.e., the deviation of the solar system from an inertial frame of reference. While this contribution amounts only to about 5-10"8 [7], it is in the same order as the relativistic corrections and can therefore not be neglected. By determining the direction of this acceleration, we include it in the accurate model of determining the acceleration on a point mass in the solar system given by Eqn. (5.3). The center of our galaxy is believed to coincide with Sgr A * and its ecliptic latitude fl and longitude A at the year 2000 are given by [127]: s = 266.85 deg, (5.9) A = -5.61 deg. (5.10) Translating these into cartesian coordinates allows the determination of a vector de- scribing the acceleration of objects in the solar system due to the rotation of the solar system barycenter around the center of the galaxy. Since this vector varies over time, and since its direction and length are based on measurements, it is necessary to in- clude error bounds on it. However, since these errors are several orders of magnitude smaller than other errors, they do not challenge the rigorous integration. Additional Effects Traditional numerical integration schemes, as used in the preparation of the ephemeris DE405, can only include quantifiable eflects. Moreover, in the framework of tradi- tional floating point integration schemes it is useless to account for quantities that have smaller effects on the result than the known uncertainties in other more in- fluential parameters. However, some of these omitted effects are actually indirectly 140 accounted for, since the results of the ephemeris have been fitted to observational data by modifying the underlying physics [125]. The most important of the omitted forces in the preparation of the ephemeris DE405 and Eqn. (5.3) are: relativistic effects of orders higher than 1/62, relativistic corrections for the minor planets, non-gravitational forces (e. g. friction), static background forces. Verified integration methods, however, cannot ignore these contributions, regardless of their actual influence. The easiest, and usually sufficiently accurate, way of including them in the mathematical models is to add error intervals to the right hand side of the equations of motion (5.3). Taylor Models for Orbits To obtain an enclosure for the acceleration of a point mass in the presence of the Sun and the planets, Eqn. (5.3) requires verified descriptions of the positions, velocities, and accelerations of the Sun and planets. As outlined before, we obtain these Taylor models from the time dependent descriptions of the osculating ellipses. The orbital elements of the planets are slowly time-dependent functions with a dominating linear part. As such they can easily be modeled within the differential algebraic framework presented in Sec. 2.1. Within the differential algebraic framework of COSY Infinity, these polynomials can be carried through Equations (5.6 — 5.8) to obtain time de- pendent polynomials describing the positions and velocities of the individual planets and the Sun. Taylor models for the cartesian coordinates are then constructed by 141 utilizing bounds obtained from the residual analysis of the DE405 ephemeris [130] and an accuracy analysis of the polynomial descriptions. It is important to underscore that we currently do not have long-term data for the orbital elements for which the described procedure yields polynomial descriptions that are equivalent in accuracy to DE405. However, the system provides the framework necessary to include time dependent descriptions of the orbital elements that allow the computation of such orbits. Currently, the implementation assumes a linear time- dependence of the orbital elements as shown in Appendix A. Thus, while the results presented in Sec. 5.1.4 are not verified with respect to the actual solar system, they are consistent and accurate within the given description of the orbits. Moreover, our results indicate that a more accurate description of the orbital elements will neither degrade the performance nor the sharpness of our integration method. Initial Conditions and Observational Data A major challenge for verified integration of asteroid orbits is the problem of obtain- ing guaranteed sets of initial conditions for the asteroid at some time. Since actual measurements are the main source for initial conditions, they do not have rigorous error bounds, but come with confidence intervals instead. This is based on the as- sumption that, in the absence of systematic errors, the measurements have a Gaussian distribution around the correct value. Hence, the quantity that describes the width of the confidence interval most naturally is a, the parameter describing the width of the Gaussian curve. While this description leaves the possibility of undetected sys- tematic errors, measurements are frequently confirmed using independent methods. Thus, it is usually appropriate to assume that initial conditions obtained from a 30 neighborhood are good enclosures of the actual initial conditions. Due to the numerous measurements of the inner planets, their positions and ve- 142 locities are known to a very high degree of accuracy. For asteroids on the other hand, often only a few measurements are available. However, if the asteroid ever comes close to Earth, high accuracy measurements using radar Doppler effects and radar delays will yield accurate results with small uncertainties for the initial conditions [142]. In the computation of the results shown in Sec. 5.1.4, we assumed that the ini- tial positions and velocities of the asteroid of interest can be enclosed in cartesian boxes of approximately i,75 km and 21:2 . 10‘5 km/s. Since current sigma position uncertainties for near Earth asteroids tend to be between 10 km and 103 km [142], obtaining enclosures this sharp poses a significant challenge for astronomers. How- ever, the actual results shown in Fig. 5.6 indicate that the Taylor model methods can successfully handle initial condition sets as wide as 10"5 AU, or approximately 1500 km. Thus, for the performance of the actual integration with algorithms based on Taylor models, the width of the enclosure boxes is not a significant factor. 5.1.4 Results and Discussion After combining the theory of Taylor model based ODE integration with the detailed model of the underlying physical concepts presented earlier, we have performed several verified integrations of asteroid orbits using Taylor model methods. Here, we show results from these integrations to illustrate the performance and applicability of these integration tools. All computations presented in this section are based on special units for the dis- tance and time that aid the integration process: 0 Time is measured in time units (TU), where one time unit equals 365.25/27r z 58.131 days. 0 Distance is measured in astronomical units (AU). 143 Since the Earth moves on an almost circular orbit, it travels approximately 2n astro- nomical units in 271’ time units. Hence, in these units the velocity of the Earth on its orbit is approximately v = (21r AU) / (21r TU) 2: 1 AU / TU. Since in these units both the distance of the Earth to the coordinate center and its orbital velocity are close to unity, all components of the right hand side of the differential equation describ- ing the asteroid’s orbit are of similar magnitude. This simplifies the error analysis of the integration tools, since it removes the necessity of weighting the errors of the individual components against each other. Moreover, this aids the step size control of the integrators, since it turns the potentially stiff system, with different step sizes for different components of the system, into a non-stiff one. Integration Results for 1997 XF11 In this section we present results obtained from integrations of the near-Earth aster- oid 1997 XF11, which has a high eccentricity of e = 0.48, low-inclination orbit with i = 4.10, and a diameter of several kilometers. It orbits the Sun with a period of approximately 1.7 years. Since one of the two intersection points of its orbit with the ecliptic is in the vicinity of the Earth’s orbit ellipse, there is a non-vanishing probabil- ity for close encounters between the asteroid and Earth. Tab. 5.6 lists predictions for the close encounters between the years 2000 and 2040: the closest approach distance is predicted to be 0.006 astronomical units in the year 2028. While this is about one quarter of the distance between the Earth and the Moon, it is close enough to warrant further research into the future orbit of the asteroid 1997 XF11. In fact, 1997 XF11 is an important test case for verified solar system integrations, since a successful integration would be able to verify the non-impact in the year 2028. The integrations discussed in this section are based on initial conditions obtained from the HORIZONS system [128]. As indicated in Sec. 5.1.3, we assumed the initial 144 conditions to be in boxes of 10“7 and 10“6 for the positions and velocities. The exact position r?) and velocity r?) at the initial time to are given by: to = JD2450465.5 z January 17,1997 (5.11) as e (—1.772691...,+0.14s722...,—0.0792s4...)4.0.5-(lo-7,104,104) r0 6 (+0.237203...,—0.612525...,+0.045832...)$0.5-(104,104,104). The cartesian positions and velocities of the asteroid 1997 XF11 during the integration interval starting at to are shown in Fig. 5.5 and were obtained from the relativistic integration discussed in Sec. 5.1.4. Since this problem is dominated by the two-body interaction between the Sun and the asteroid, the evolution of the coordinates exhibits the periodic behavior of elliptical orbits with high eccentricity. In the remainder of this section we present results obtained from two different integrations of 1997 XFll’s orbit. The first integration uses Eqn. (5.3), while the mathematical model of the second integration does not include the relativistic cor- rections to the classical Newton gravitation. As such, the two integrations allow us to measure the actual influence of the correction terms on the asteroid’s final coordi- nates, which gives an estimate on the computational cost and final benefit associated with these corrections. Integration with Relativistic Corrections In this section we present results from a ten year integration of the asteroid 1997 XF11, using Eqn. (5.3) with relativistic corrections and initial conditions boxes given by Eqn. (5.11). The asteroid’s positions obtained from this integration are shown in Fig. 5.5. As indicated earlier, the results demonstrate that the asteroid is moving on an almost elliptical orbit around the Sun. More interestingly though, Fig. 5.6 gives the size of the enclosures for the coordinate a: and the velocity i2. It is important 145 Position x E Position Y 2 _ -------------- Fashion 2 "I, “\,\ :/'\\ l, “A: "I, 1 . .' ' ' l f ’- "f . 1‘. .‘ l‘ [ o .1 .................... .7 ......... .......... . . ’ :1 , -1 l " I 1 997 1 999 2001 2003 2005 2007 Velocity v? .......... Velocity v 15 ~ -------------- Velocity v; MM -1 4 . . 1 997 1 999 2001 2003 2005 2007 Figure 5.5: Cartesian positions (in AU) and velocities (in AU/TU) of 1997 XF11 during the ten year integration interval. to note that the sizes of these enclosures do not just grow monotonically with time, but even decreases at certain points. Moreover, the enclosures have a very small average growth rate, which is the main reason for the ability of the Taylor model based integrator to propagate initial conditions over large time intervals. Tab. 5.6 lists the closest approach distance between 1997 XF11 and the Earth in the year 2028 as 0.006 AU. By extrapolating the enclosures that can be obtained from the Taylor model integration on the other hand, we observe that the enclosures at that point will likely be smaller than 0.006 AU. Since this is sufficiently accurate to guarantee a non-impact, it shows that the Taylor model approach can indeed be used 146 0'001 Enclosure 01x 19-04 : ‘1 16-05 r 1 1 l l 1e-06 ~ l 1 16-07 4 t 4 1997 1999 2001 2003 2005 2007 0-001 3' Enclosure?! v, f r u l l l 16-04 g l l l l 1 l 1 16'05 .- '] E i r l l l 1606 + + g a * 1997 1999 2001 2003 2005 2007 Figure 5.6: Logarithmic plot of the diameters of enclosures for the positions (in AU) and velocities (in AU/TU) of 1997 XF11 during the ten year integration interval. to obtain meaningful results from verified long-term integrations of real world systems. In Sec. 5.1.5 we will compare this with similar results obtained from conventional verified integration tools. To illustrate the technical aspects of the integration method, Fig. 5.7 shows how the automatic error and step size control of the Taylor model integration method work in practice. The maximally allowed local error was set to 10‘9 and Fig. 5.7 shows that the integrator has been able to meet this goal by reducing the step sizes accordingly. A comparison between Figs. 5.7 and 5.5 shows a strong correlation 147 Date Min. Dist. (AU) October 31, 2002 0.064 June 10, 2016 0.180 October 26, 2028 0.006 June 8, 2035 0.174 November 4, 2040 0.150 Table 5.6: Dates and minimal distances of the predicted closest Earth approach distances for the near-Earth asteroid 1997 XF11 between the years 2000 and 2040 [73]. between the local error and the current location of the asteroid on its orbit: the local error tends to increase, with a corresponding reduction in the step size, whenever the asteroid is closest to the Sun. In these regions its acceleration is at a maximum and the asteroid moves with its largest orbital velocity. It should be noted that a significant contribution to the local errors stems from the uncertainties of measured quantities and the planets’ orbits. Thus, the main source for local error can therefore not be subject to verified error control. Consequently, the desired local error is close to the optimum for this particular problem. Integration without Relativistic Corrections Here we present results of a ten year integration similar to the one discussed in the previous section. However, this integration is based on the classical Newton gravitation without any of the relativistic corrections included in Eqn (5.3). This allows us to quantify the actual effects of the relativistic corrections on the computed orbit of 1997 XF11. Hence, the analysis determines the importance of a complete model including the relativistic corrections given by Eqn. 5.1 and puts a measure on the computational overhead of these corrections. For this integration, we used the same initial conditions as in the previous example. Moreover, the automatic step size control of the Taylor model integrator chooses 148 r h —— Stepsize mgr time o 1 _1 1_ 1 1 997 1 999 2001 2003 2005 2007 1,2 . . a fi 6‘09 -—-— One-step integration error 16-09 [- 86-10 F i 6e-10 - J 4e-10 l 26-10 [ M 1997 1999 2001 2003 2005 2007 Figure 5.7: Evolution of the step size and one-step integration errors during the ten year integration of 1997 XF11. the same step sizes as in the previous computation, indicating that the relativistic corrections have only a small influence on the final result. Firstly, it should be noted that the integration without the relativistic corrections requires only about one tenth of the CPU time necessary for the accurate integration discussed before. Secondly, the enclosures of the computed final positions and veloc- ities are virtually indistinguishable from the ones shown in Fig. 5.6. The differences between the positions and velocities computed by the corrected and the uncorrected integrations are shown in Fig. 5.8. Since the differences are much smaller than the sizes of the enclosure boxes, the performance of verified solar system integration could 149 39-05 28-05 7 19-05 - o . -1 e-05 1 -29-05 > ——-— Difference—Ax ......... -. Difference Ay .............. Difference A2 1* - .' - -2 J . w::::.‘:_". I. .2- - L -3e-05 1997 39-05 2.56—05 26-05 1 .56-05 1 e-OS 56-06 0 . -56-06 -16-05 1 997 Figure 5.8: Differences in the computed positions and velocities of the asteroid 1997 XF11 between the ten year integrations with and without relativistic correc- tions. be improved significantly by finding guaranteed enclosures of the relativistic correc- tions in Eqn. (5.3) that are easier to compute than the corrections themselves and are not overly pessimistic. The net effect of this approach would be a relatively small tradeoff in accuracy for a substantial gain in integration speed. However, such rigor- ous estimates on the size of the relativistic corrections are not known and it is not clear if they could be rigorously computed with an effort that is significantly smaller 1999 2001 2003 2005 2007 r Difference sz -—— 'DiffererTce Av, ........... Difference Avy r l 1 J 1999 2001 2003 2005 than the one needed to compute the corrections proper. 150 Shrink Wrapping An important aspect of the Taylor model based integration method is the ability to dynamically shrink the size of the remainder bounds during the integration. This is extremely useful since it delays the exponential inflation of errors by several orders in time. This dynamic reduction of errors is called shrink wrapping [28]. The main idea behind shrink wrapping is illustrated in Fig. 5.9: at certain points during the integration, the linear part of the Taylor model M propagating the initial conditions is enlarged, resulting in a new Taylor model M5. The polynomial part of this new Taylor model maps the set of initial conditions to a set that properly contains the original enclosure of the final coordinates plus the remainder bounds. At this point, the remainder bounds of MS can be shrunk to zero width, since the poly- nomial map of M 5 already ensures containment of the image of the initial conditions in the final enclosure. However, this parts with the Taylor model approach viewing the reference polynomial as the Taylor expansion of the flow of the differential equa- tion. On the other hand, shrink wrapping enables Taylor model based integration methods to obtain long-term results without an overly pessimistic enclosure of the final positions and velocities. Since the errors of the right hand side of Eqn. (5.3) are relatively large, shrink wrapping has been performed after each integration step in the computation of the results presented in Sec. 5.1.4. The shrink factor is the factor by which the linear parts of the Taylor models’ reference polynomials are enlarged. We distinguish two kinds of shrink factors: 0 The one-step shrink factor determines the amount of scaling in a single shrink wrap Operation. 0 The accumulated shrink factor gives the amount by which the linear parts have 151 Figure 5.9: Illustration of the shrink wrapping method used in the verified integration of asteroid orbits. been scaled over the whole integration at this point. Fig. 5.10 showsvthe values of both shrink factors, computed during the ten year integration with relativistic corrections presented in Sec. 5.1.4. The figures illustrate two important and somewhat surprising aspects of the shrink wrapping method: on average, the one-step shrink factor decreases exponentially over time and the total shrink factor grows only linearly and not exponentially. This behavior is particularly helpful for the computation of long term results, since it prevents an exponential inflation of the enclosures of the final positions and velocities over the integration interval. Comparisons between integrations with and without shrink wrapping show that the availability of the method is a key ingredient of an overall strategy geared towards verified long-term integration results [86]. 152 1.0010 1.0009 ~ 1.0008 ~ 1.0007 - 1.0006 - 1.0005 - 1.0004 1.0003 1.0002 1 .0001 ‘ ‘ ‘ 1997 1999 2001 2003 2005 2007 One-step shrink fafztor 3.5 . Total shrink factor 3.0 - 2.5 r 2.0 t 1.5 * _L J _1 1 .0 ‘ L ‘ 1 997 1 999 2001 2003 2005 2007 Figure 5.10: One-step and total shrink factors obtained during the ten year integration of the asteroid 1997 XF11. 5.1.5 Comparison with AWA During the last decade, AWA [79, 80] has become a standard benchmark for verified ODE integration methods. While it tends to be overly pessimistic for complicated problems, the method works reasonably well for converging linear problems. In this section we compare the performance of AWA with the Taylor model based integra- tion schemes used before, and we will demonstrate that the Taylor model approach performs favorably when compared with the conventional interval techniques used in 153 AWA. Since AWA lacks advanced control structures which would aid the implementation of the right hand side as given by Eqn. (5.3), and cannot handle uncertainties in the right hand side of the differential equations very well, we have only used it for the integration of the two-body problem of the Sun and the asteroid 1997 XF11. 2°. ’7 .. r = ———-—-r. (5.12) HT”3 '7 =2 0.9986 is the value of the product of gravitational constant C and the Sun mass M in the units that have been used for all previous integrations. It should be noted that Eqn. (5.12) represents a much simpler problem than the one previously discussed, since it contains only one planet orbiting the Sun, no relativistic corrections, no corrections accounting for the solar system being a rotating coordinate system, and no uncertainties for the positions and velocities of the Sun. As in the previous cases, the initial conditions are given by Eqn. (5.11), and AWA has been set to integrate the system in order 18 with an initial step size of 0.0001. For comparison, the same integration has been performed with Taylor models of order ten and an initial step size of 0.1. In both cases the maximum local errors were set to 10‘“. It should be noted that the model discussed in this section is an extremely simplified version of the real problem described by Eqn. (5.3) and has therefore no relevance for the verified integration of solar system dynamics. Also, while we would have liked to perform the integrations over time intervals spanning more than one year, AWA was unable to integrate the system further, since the wrapping effect led to an exponential inflation of the enclosures. (The time evolution of the position and velocity enclosures is shown in Fig. 5.11). The results given in Fig. 5.11 show that the Taylor model approach allows the com- putation of much sharper bounds on the final coordinates than conventional interval 154 0.1 , Y— T fi 1 '— ' + Enclosure of x , . . x Enclosure 01y ,3 . 0.01 r - Enclosure of z #3.. , . , x . 1 [ +/M xx I 0.001 1' +3. x .- * If“, I l ’ fi 16-04 E + ,th / , + + x I ++' x ' 1 6'05 1' W. x .. l #yl" X I 1 ‘4'” Mix... l r I [ 1e-06 - £39593" , . i , ' l 1 6’07 :— 9 J; 1* g 1997 1997.2 1997.4 1997.6 1997.8 1998 1998.2 0'1 l . ' Enclosfire ofx‘ # ( .- : 52333329: 3 0.001 [ l 1e-04 [- [ 16'05 I Minamzflll i 1e-06 “‘1“. i 1e-07 [ ‘I ] 1e-08 ’ + A 4 1997 1997.2 1997.4 1997.6 1997.8 1998 1998.2 Figure 5.11: Logarithmic plot of the diameters of the enclosures for positions obtained by AWA (upper) and the Taylor model based integrator. techniques, represented by AWA. At the point where AWA terminates the integration, the Taylor model enclosures are almost four orders of magnitude sharper than the corresponding results obtained by AWA. Moreover, at the final point of integration, the enclosures obtained by AWA exhibit a large exponential growth while the Taylor model results show a much smaller growth rate and even decrease at certain times. This smaller rate of increase gives the Taylor model methods the ability to integrate over much larger time intervals than conventional interval based tools. 155 0.1 ‘1 I F— + Enclosuré of x J x Enclosure of y , .1! J 0-01 F - Enclosure of z :" ‘ 0.001 L Wt: _' 3 . . . 18-04 7 3' i . xi ’ . ++X 1e-05 :- .. [ . ”W” x 'l 1 19-06 :- .+*w W 1‘"- 1 I 5.1%}: l 1e-07 :- _‘ 1 . ' l 1e-08 ‘ + L E 1997 1997.2 1997.4 1997.6 1997.8 1998 1998.2 Figure 5.12: Logarithmic plot of the diameters of the enclosures for positions obtained by AWA for an initial box reduced by a factor of 3.275 in each direction. Since the Taylor model approach is much more computationally complex than regular interval arithmetic, it is natural to ask if we could split the original problem for AWA in such a way that the total CPU time used is the same for both AWA and the Taylor model integration. Using AWA, the computation of the results shown in Fig. 5.11 took 4.5 seconds of CPU time. On the same machine, the Taylor model based integration took more than 92 minutes, or about 1234 times longer. Since (3.275)6 z 1234, we then split the boxes of initial conditions into boxes that are 3.275 times smaller in each coordinate direction and ran AWA with initial conditions given by these smaller boxes. Fig. 5.12 shows the resulting enclosures of the final positions for the same integration period as before. It is striking that even if AWA uses the same amount of CPU time as the Taylor model based integrator, there is hardly any change in the resulting size of the enclosures. Moreover, the sharpness of the results given by AWA is still not comparable with the one computed from Taylor models. As a final test, we tried to determine the size of initial condition boxes for which 156 AWA can compute a final enclosure comparable in sharpness to the Taylor model results. However, even the smallest possible initial sets; i.e., boxes with a magnitude of the machine epsilon of 10‘15 in each of the six coordinate directions, do not allow the computation of final enclosures comparable to the Taylor model results. Details of this integration are shown Fig. 5.13 and it is important to note that even with these tiny initial boxes, the final enclosures are virtually indistinguishable from the ones obtained with larger initial sets. This indicates that the main challenge for the verified integration of this system lies in the wrapping effect, which is known to limit the applicability of conventional interval techniques. 0.01 w— r fl r + Enc[osure o: x z , s x no osureo y ,.;:. . 19'04 ' - Enclosure ofz ,W..-‘ J l 9.13“?" 1 16-06 F “,5“? 5 ' W;: :1 II. [ 9"” .- 16-08 ' WW II" ~ l K . 1e-10 [ ,. l . Ix ] 1e-12 l :i' l . . [ 19-14 l t‘ ( [ f 1e—16 " . 1997 1997.2 1997.4 1997.6 1997.8 1998 1998.2 Figure 5.13: Logarithmic plot of the diameters of the enclosures for positions obtained by AWA for initial condition enclosures of 10‘”. While normal interval techniques will almost always outperform Taylor model methods for problems that do not require any domain splitting, it has been shown in Sec. 3.6.1 that if domain splitting becomes necessary, the additional expense of using Taylor models is more than compensated for by the increased accuracy. Here we even have a problem for which the normal interval methods, represented by AWA, cannot compete with Taylor models in terms of sharpness regardless of the amount of domain splitting. Lastly we mention that the integration of a two—body system could 157 have been implemented more favorably for AWA, since the problem can be reduced to four dimensions. However, the main purpose of these computations was to gauge the applicability of conventional interval techniques to the six-dimensional problem of solar system dynamics and not to fabricate a simple test case. 5.1.6 Summary The results presented in Sec. 5.1.4 and Sec. 5.1.5 show that the Taylor model approach can indeed be used for the verified integration of solar system dynamics. And it has been shown that the Taylor model approach is capable of propagating large sets of initial conditions of large time spans without falling prey to the wrapping effect that often leads to significant overestimations of the final enclosures in verified integrations of ODE initial value problems. Moreover, by comparing the results from the Taylor model based integration method with similar results obtained from tools utilizing conventional interval tech- niques, it has to be concluded that the Taylor model approach is the only method that allows the verified long-term integration of asteroid orbits and other real world problems. Compared to conventional interval techniques, the Taylor model approach gives several orders of magnitude sharper results and can propagate initial conditions over orders of magnitude larger time spans. 5.2 Existence of Generating Functions In this section we discuss a particularly interesting application of Taylor model meth- ods to the theory of generating functions of canonical transformations. We show that using high-order Taylor models enables us to rigorously establish lower bounds on the size of the domains of definition of any type of generating function. Moreover, 158 we will see that the computed domains often enclose the dynamic apertures of the examples studied. The results of this analysis have far reaching implications in the study of particle accelerators and symplectic tracking. 5.2. 1 Introduction It is well known that coordinate transformations of Hamiltonian systems that keep the Hamiltonian structure intact are called canonical transformations, or, in more modern terminology, symplectic maps [17, 49, 123]; the time evolutions of Hamilto- nian systems are also symplectic. Any symplectic map can be represented in terms of a single scalar function, called the generating function. Until recently, only the four conventional Goldstein—type generators were well known [49]. However, following the introduction of extended generating functions, it has been shown that to each sym- plectic map, infinitely many generating function types can be constructed [40]. In certain applications, such as long term simulation of accelerators and other Hamilto- nian systems, it is important to maintain symplecticity during tracking [47, 1, 40]; one available method to achieve this is the generating function based symplectic track- ing[14,40,138] In principle, any of the valid generators could be used for tracking of a given Hamiltonian system. Although it can be shown that for any globally defined sym- plectic map, global generators can always be found [40], their construction is however usually not straightforward. Indeed, as shown in [17, 40], the representations of sym- plectic maps through generators commonly available in practice are often only locally valid. Frequently, the purpose of tracking simulations is to estimate the region of space where stable particle orbits can be found, the so—called non-wandering set, or dynamic, aperture in the beam dynamics terminology. Hence, a sizable phase space region must be covered by tracking, and if the generating function method is used, 159 the generating function must be defined at least in that region. However, so far nothing has been known about the size of the domains of definition of generating functions, and the necessity to study this interesting problem has been recognized in [14, 138, 33, 40]. 5.2.2 Theory of Extended Generating Functions Following the exposition of [40], we regard every map as a column vector. Let a = ( Z: ) (5.13) be a diffeomorphism from the Open set U C IR” onto its image; a; and (12 are first and second 2w components of 0:, respectively. In other words, for i = 1,2, a,:U-—>V,~C1R2“’. L613 A B Jac(a) = ( C D ) (5.14) be the 420 x 4w Jacobian of a, split into 2w x 220 blocks and define ” _ J2w 0211) A,” ), (.15 where \ 0w 1,, = ( -1. 0w ). (5.... In the last equation, Iw denotes the unit matrix of appropriate dimension. We call a map a conformal symplectic if its Jacobian satisfies (Jac (0))T ' J4w ' Jac (a) = p ' j4w, (5.17) with some non-zero real constant p [6]. Finally, a map M is called symplectic if its J acobian M satisfies the symplectic condition [39] MT-J-Mz J. (5.18) 160 Henceforth, we will consider only symplectic maps that are origin preserving; i.e., maps around fixed points. We call a map a gradient map if it has a symmetric J acobian N, and it is well known that, at least over simply connected subsets, gradient maps can be written as the gradient of a function, which is called the potential of the map [17]. With these preliminaries, the main result of the extended generating function theory is best expressed by the following theorem [40]. Theorem 5.1 (Existence of Extended Generating Functions). Let M be a symplectic map and denote the identity map by I. Then for every point z in the domain of M, there is a neighborhood of z such that M can be represented by functions F via the relation W..-(..(A;».(..(g»a where a = (0:1, a2)T is any conformal symplectic map such that det (C (M (z) ,z) - M2 + D (M (2) , 2)) 7t 0. (5.20) Conversely, for any (32 function F the map .M defined by M = (NC — A)“1 (B — ND) (5.21) is symplectic, where the matrices A, B, C, D, M, and N are defined as above. The function F is called the generating function of type a of M, and denoted by Fam. Theorem 5.1 says that, once the generator type a is fixed, locally there is a one-to-one correspondence between symplectic maps and scalar functions, which are unique up to an additive constant. The constant can be normalized to zero Without loss of generality. Due to the fact that there exist uncountably many maps of the form (5.17), to each symplectic map one can construct infinitely many generating 161 function types. But in any case, according to Eqn. (5.19), for a given M the domain of definition of the generator of type a equals the domain of invertibility of (mm)- Thus, finding the domain of definition of a generator is equivalent to finding the domain of invertibility of (5.22). A large class of practically used generator types are the generators associated with linear maps oz. These can be organized into equivalence classes [S], associated with —J L"1 J ) a = _ , (5.23) (§(I+JS)Ll §(I—JS) and represented by arbitrary symmetric matrices S [40] where the matrix L denotes the linear part of M. Finally, we note that the size of the domains of definition of the extended generat- ing functions is directly related to the problem of optimal symplectic approximation of Hamiltonian flows [40]. A very general theoretical argumentation based on Hofer’s metric states that in general the generating function associated with S = 0 gives the Optimal approximation. In the examples shown in Sec. 5.2.4 we will therefore concentrate on this particular choice of S. 5.2.3 Enclosing Derivatives of Flows In order to prove the existence of a generating function over extended domains with Taylor model methods, we need to obtain Taylor models for the derivative of Eqn. (5.22). In the case of linear maps oz, this reduces to the computation of Taylor models for the J acobian of the map M. While simple systems may allow the explicit determination of the symplectic map and its derivative, in the more common case of practical problems, the exactly known information are the Hamiltonian functions, or 162 equivalently, the corresponding ordinary differential equations. Since closed form so- lutions of these are generally not available, numerical integration schemes are utilized to compute approximations to the mathematically correct solutions. If we use the Taylor model integration method discusses in Sec. 2.3.3, we can enclose the mathe- matically correct solution by a Taylor model. However, since the derivative operation does not extend to Taylor models, we cannot use this Taylor model to compute Tay- lor models for the derivates. But the following result from the theory of ordinary differential equations gives us the means to compute Taylor models for the J acobians of flows [54]. Theorem 5.2 (Smooth Dependence on Initial Conditions). Suppose that f (t, 11:) is a 6"“ function defined on an open set 52 C R x IR". Let p(t, u, y) be the unique solution to the initial value problem 17' = f(t,x), m(u) = y. (5.24) Then the function ¢(t, u, y) is a 6"“ function of the variables (t, u, y) and the partial derivative of the solution 239% with respect to the spatial initial conditions satisfies the matrix initial value problem I” = (3;: (t,¢(t,u,y))) -Y, Y(u) = id. (5.25) Thus, by using Taylor model based integration techniques to rigorously enclose the solution of the enlarged ODE initial value problem 23' = f(t, 1:), (5.2%) Y’ = (g (m) ~Y, (5.26b) m(to) = x0, (526(2) mo) = I, (5.26d) 163 we can compute Taylor models containing not only the map M(tf, to, .710) at the final time tf, but also the Jacobian g—fifflbtmxo) of the flow. In light of Eqn. (5.26), the symplectic condition (5.18) loses its differential algebraic nature and becomes a purely algebraic constraint. Thus, similar to energy conservation, the symplectic condition expresses in fact a conservation law, which can open the door for symplectic integration in the framework of as differential algebraic equations. It should be noted that the augmented ODE system (5.26) is less smooth than the original problem, resulting in solutions that are also less smooth than before. However, since the systems of interest are commonly analytic over their domains, this is usually not a problem. Finally, it is important to note that while the augmented system is v(v + 1) dimensional, there is no need to propagate any dependency on the additional initial conditions in the augmented system. Thus, the integration of the augmented system can still be done with Taylor model in v variables, resulting in a system that is as manageable as the original problem. 5.2.4 Examples To illustrate how the existence of generating functions can indeed be proved with full mathematical rigor, we will now discuss several examples of symplectic systems. We show that the Taylor model based invertibility tests can often guarantee invertibility of the maps over regions that properly enclose the dynamic aperture. An Exactly Symplectic Map As a first example we consider a two dimensional cubic polynomial, which is exactly symplectic and given by M = N 0 £, (5.27) 164 where __ cos§ sin§ Lfl(—sin§ cosg)’ (5'28) and —3((1+I))3 N(q)=(q ). 5.29 p P+3(Q+P)3 ( ) The main advantage of this relatively simple example over the more complicated examples discussed below is it allows an analytical treatment. According to Thm. 3.2, extended generating function exist over extended regions D C IR2 if the following polynomial expression for the determinant in the four variables ql, p1, qg and p2 has no zeros in D x D C R4: [(1--\/§)ql+(1+~/®101]2 1 + 9 8(1+\/§)2 [(2+\/3)(312+1)+322] (5.30) - 1 x/— 1 + 9l(1 @fl:%; 3)p]2 [(2+\/-)(312-1)+ (7+4\/§)311], where S = ( S“ 312 ) (5.31) 812 822 is any symmetric matrix with real entries. It is easy to see that Thm. 3.2 guarantees the existence of globally defined inverses if the entries of S satisfy: ’322< (2 + f) (312 + 1) (5.32a) 2+\/§ —311$ 7+ 4f (312 — 1) (5.325) Thus, we can analytically prove the existence of a large class of. globally defined generators for this symplectic map. However, while global domains of definition are of great theoretical value, in practical applications we are interested mainly in proving invertibility in finite regions that enclose the dynamic aperture. For these cases, the Taylor model based methods are well suited. For example, employing the Taylor 165 Figure 5.14: Tracking picture of the cubic two-dimensional symplectic map, and the box of guaranteed invertibility of the generator associated with S = 0. method for S = 0, the tracking picture shown in Fig. 5.14 was obtained. The box surrounding the tracking is the region for which invertibility of the corresponding generator can be proved using the Taylor model approach. Clearly, it extends well beyond the dynamic aperture, which is approximately 0.2. The Fermi-Pasta-Ulam System Since the map of the previous system was given by an exactly symplectic polynomial, we were able to derive an analytic expression for the derivative of the map and compute a Taylor model for the Jacobian. This Taylor model allowed us to apply the invertibility test to prove the existence of the generating function over extended domains. Here and in the next example we study problems in which the maps are obtained by integrating an ODE initial value problem, requiring the application of Thm. 5.2 to obtain Taylor models for the Jacobians of the flows. The Fermi-Pasta—Ulam system [45] is an extensively studied Hamiltonian system which is known to exhibit many interesting features and has stimulated important developments in nonlinear dynamics. The Fermi-Pasta-Ulam system consists of a 166 finite series of Springs; regular springs obeying a linear force law alternating with stiff non-linear springs. The general Hamiltonian function of the system is given by d 1 n k n n H ((1,15) = E 2 (13327—1 + 1033') + 5 Z (921' " Q2i—1)2 + 2012:“ "' (1204- (5-33) i=1 i=1 ':0 For this example we considered the problem with n = 2, k = 5000, and q4 = q5 :- p4 = 0, resulting in the system shown in Fig. 5.15. V \3 ql q2 q3 Otttouattttauo Figure 5.15: Illustration of the Fermi-Pasta—Ulam system with three masses and four springs. By using the Taylor model based integration scheme, we obtained a Taylor model for the half-period map of the flow of (5.33), and then tracked several particles by iterating this map. Typical phase space plots are shown in Fig. 5.15. The domain of the generator associated with S = 0 extends to at least [—1, 1] in every direction of the phase space. Although the concept of dynamic aperture does not apply in this example, since the system is bounded, the size of the resulting domains are comfortably large. A Particle Accelerator Cell We conclude the section with an example of practical interest in accelerator physics. The Hamiltonian representing an accelerator magnet with the arc length 3 as inde- pendent variable and on-energy particles is given by [17] H(:r,y,a,b)=—(1+%)m—%(1+%)As(x,y). (5,34) 167 Figure 5.16: ((11, p1) and (q1, Q2) tracking pictures of the half period map of the Fermi— Pasta—Ulam system for a particle launched along the q] axis. The box of guaranteed invertibility of the generator associated with S = 0 extends to at least [—1, 1] in every direction. Here p is the curvature radius of the magnet, A, is the 3 component of the vector potential, and ((m, a) ; (y, b)) are two pairs of canonically conjugate variables. For the sake of computational simplicity, we assume that the magnetic fields are piecewise independent of the arc length, in particular we use the sharp cutoff approximation for fringe fields, as is commonly done in beam physics. To find the Hamiltonian functions for a specific magnetic element or a field free regions of space, all that is necessary are the specific values of p and A,. In field free 168 space is obviously A, = 0 and p 2 00. For a homogeneous dipole magnet, the vector potential and curvature radius are connected via B A. = «,9 (a: + p), (5.35) where qBo/Po = 1/ p. For a quadrupole magnet on the other hand p z co and the vector potential is given by A, = -—-;E (2:2 — yz) , (5.36) where k is the quadruple strength. The case It > 0 denotes a magnet that is focusing in the :r. direction, and correspondingly k < 0 indicates a magnet that is defocusing in the a: direction and focusing in the y direction. Lastly, the Hamiltonian function of a sextupole magnet is determined by p z co and A8 = —g (a:3 —- myz) , (5.37) where h is the sextupole strength. Combining these elements, we set up a cell consisting of the following sequence of elements: drift, defocusing quadrupole superimposed with a sextupole, drift, bend- ing dipole, drift, defocusing quadrupole superimposed with a sextupole, and drift. The defocusing quadrupoles have the same strength k = —0.0085 and the sextupole strength is h = 0.06. The lengths are 1 meter for the drifts and 0.5 meter for the magnets; the dipole’s curvature radius is p = 2.5 meter and the reference particle is a proton with an energy of 1 MeV. The resulting arrangement is illustrated in Fig. 5.17. Integration of the system yielded a Taylor model containing the true solution of the corresponding differential equations. For illustration, Tab. 5.7 shows the first ten terms of the Taylor model containing the x-component of the flow, which is the function mapping the initial conditions :ro, a0, yo, and b0 to the at position at the end of the accelerator cell. The full Taylor model is of order 17 and its reference 169 E:g 0. With the exception of the number 0, all numbers are normalized such that the leading digit 2:1 is non-zero. While the IEEE standard 754 requires B to be two and allows T = 24 and T = 53, the newer standard IEEE 854 also permits B = 10. However, in the remainder of this section we will always assume that all floating point numbers are represented in a fixed binary format; i.e., we assume that B = 2 and that the integers 176 T, Emin, and Emax are fixed. In the following small letters at, y, and 2 denote real numbers, and their corre- sponding representations on the computer are denoted by 2:, y and z. The elementary operations +, —, x, and / on the set of floating point numbers are denoted by EB, 69, ®, and Q, respectively. The most fundamental and most limiting property of floating point numbers is that only a finite subset of rational numbers has an exact floating point representation. We denote this set of floating point numbers by T f={x€Qi(Z$12-i) 2EaEminSESEmaXazlzl} (6'2) i=1 It should be noted that within any given floating point system .73, some important real numbers have exact floating point representations: —2, —1, 0, 1, 2 E .7. The next proposition highlights another two important floating point numbers which will be used throughout this section. Proposition 6.1. For a given binary representation of floating point numbers, define the numbers 6 and ,u by e = 21‘T (6.3a) p = 23min“. (6.3b) Then 5, u E .7: and the two numbers can be characterized by e = min{x|$6}'and1€Bm>1}, (6.4a) p = min{x|z€fand:c>0}. (6.4b) Proof. To prove that the two numbers have exact floating point representations, we write 5 = 21‘T _-: 21 ~2‘T. (6.5) 177 Thus, according to Eqn. (6.2), it is indeed e = 6:. Similarly, it is ’1’ : 28min+1 : 21 . 23min (6.6) Hence, the number it also has an exact floating point representation. Next we show that e is the smallest floating point satisfying the condition 1 ED :1: > 1: 1 69 e (21 - 2-1) 99 (21 . 2") = 21 ® (2-1 + 2-T) (6.7a) II = 21-(2—‘+2-T)=1+e>1=1. (6.7b) On the other hand if we assume that there is y E F such that y < e and y 69 1 > 1, then 3; is of the form T y = 21‘T - 23 . 2.2.2“ (6.8) i=1 with Emin + T — 1 g E S 0 and 2:1 2 1. Consequently it would be T 1+3; = (21.2-1)ea(21434-2334) (6.9a) i=1 T = 21- (2‘1 EB 2244“?) (6.95) i=1 = 21-2*1=1= 1, (6.9c) in contradiction to the assumption that y ED 1 > 1. Thus, the floating point number 5 is indeed characterized by (6.4a). The last part of the assertion follows from the definition of p and the basic require- ments of floating point representations: clearly ,u > 0, and any positive floating point number smaller than p would have either E < Emin or $1 = 0, both violating the definition of the set .7: of floating point numbers. Thus, it can indeed be characterized by (6.4b), making it the smallest positive floating point number. (:1 178 6.1.2 Directed Rounding In Sec. 2.2 we showed that the sum of two intervals is given by [131,332] + [91,312] = [331 + 132,311+ 312] (6-10) However, if we naively implement this interval addition on a computer with floating point numbers satisfying T = 53, the result of the following addition would violate the basic inclusion property of intervals: [—1,1]+ [264,264]. (6.11) While all interval endpoints have exact floating point representations, the floating point sums needed for the interval addition would lead to truncation: —169264 = 264, (6.12a) 169264 = 264. (6.125) Clearly, this situation is not acceptable in an interval library used for rigorous com- putations and the necessity of finding a solution to this problem has already been recognized in [92]. The easiest way of alleviating this problem is by using directed rounding for the result of interval computations: at the end of each elementary inter- val computation, the lower bound of the resulting interval is replaced by a number that is, by an appropriate amount, smaller than the originally computed floating point number. Similarly, the computed result for the upper bound is shifted by an appro- priate amount towards +00. Since floating point numbers are not uniformly spaced, the results cannot be computed by simply adding fixed numbers to the endpoints. Instead, directed rounding towards —oo and +00 is used. The most accurate and most efficient way of implementing directed rounding in an interval library would be to rely on intrinsic routines built into IEEE 754 compliant 179 floating point hardware: the standard requires the hardware to provide at least the following three rounding modes: 1. rounding towards —00, 2. rounding towards +00, 3. rounding to the nearest floating point number. With the ability to set the rounding mode, Alg. 6.1 could be used for the addition of two intervals [11:1, 12:2]. Algorithm 6.1 Interval addition with hardware support to set the rounding mode. save current rounding mode set rounding towards —0o 21 4: 101 EB 91 set rounding towards +00 22 ¢ $2 + 312 restore saved rounding mode However, this strategy defies the goal of writing a portable interval library, since it would require hardware-specific code to change the rounding modes. Moreover, even on the same platform, different compilers use different conventions of embedding hard- ware instructions, leading to an unacceptable platform dependence in the resulting code. Finally, it should also be noted that on many current hardware platforms a change in the rounding mode results in a severe performance penalty, since it leads to pipeline flushes that render this approach unacceptably slow [65, 135]. Only recently have interval researchers been able to convince hardware manufacturers of the need for more efficient rounding mode switches, and the resulting interval implementations have shown substantial increases in performance [136, 137]. In order to implement directed rounding in software, it is necessary to know by how many Units in the Last Place (ULPs) the computed result of any of the floating 180 point Operations differs from the mathematically correct result. We will denote this number by N and note that most Of the current computer manufacturers guarantee N to be at most one. Based on the number N, the following theorem allows the use of directed rounding for verified interval results: Theorem 6.1. Leta: E .77 be given and assume that it is sufficiently large ( c. f. Cor. 6.1 for a definition of how large a: has to be). Provided that none of the following opera- tions leads to a floating point overflow, the number (1 EBsgn (m)®N®e) ®x (6.13) is by at least N ULPs bigger than :6. Similarly the number (1 esgn($)®N®6) 82: (6-14) is by at least N ULPs smaller than 2:. Proof. The assertions will be proved for N = 1 and for positive numbers 2:. Similar arguments hold for higher values of N and negative numbers at. It is (1 + e)a:— — 21 (2‘ + 2”) 2E-(: av.- 2") (6-159) = 21+E (2 +2”) (22.2 “) (6.15b) T T = E. (Z x12"‘+Z:r,-2‘(‘+T‘l)) (5-150) T- 1 =]ZE ((iafl” + (xr+:r1) )2 T+Z$ '2 (”T 1)))J (6115(1) = E. (: +($T + $1)2_T) (6.15e) 181 Since $1 =- 1, the result is indeed by at least one ULP bigger than the original number (a. A similar computation shows that the second multiplication does indeed round towards —00: (1 - e)m = a: 9 ea: (6.16a) T T = 23. (211.2“) — 21*T - 253- (21,24) (6.165) i=1 i=1 T T = 23 - (211,2-i — Z ass-(”F”) (6.16c) i=1 i=1 T—l T 2 ]2E ' (2 “322—1. + (331" " $1)2_T + Z$i2_(i+Ti3))] (616d) i=1 i=2 Since 131 = 1, the result is indeed by at least one ULP smaller than the original number m, which finishes the proof. [:1 A Portable Algorithm for Directed Rounding Based on the results derived in the previous section, we can now present an algorithm that, under mild assumptions on the accuracy of the hardware and system software, allows the computations of verified interval Operations with directed rounding. The next two corollaries are immediate consequences of Thm. 6.1 and they provide the mathematical foundation for software implemented directed rounding: Corollary 6.1. The number 6 = p/ (1 — Ne) is the smallest positive floating point number that can be rounded towards zero without remaining fixed or resulting in an underflow result of 0. Corollary 6.2. Let x,y E R be such that min(|x|,ly|) 2 6 and assume that the following multiplications do not result in a floating point overflows Then the following 182 inclusion holds: [23, y] C [(1 — sgn(:c)Ns) - m, (1 + sgn(a:)Ne) - m]. (6.17) Corollary 6.2 excludes two special cases that require exception handling in an im- plementation of an interval library that uses software implemented directed rounding: o The multiplications might result in overflow errors and since this cannot be handled in a general and portable way, the exception handling is left tO the runtime environment. 0 Rounding Of floating point numbers a: with 0 _<_ [2:] < 6 need special attention because a multiplication may result in an underflow error or a fixed point. From now on we will assume that N = 1, and we define the following two numbers: M1, 2 1+5, (6.18a) Md 2 1—8. (6.18b) With these numbers and the previously defined 6, Alg. 6.2 describes portable software simulated directed rounding: if [$1,232] is a floating point interval, then Alg. 6.2 performs the necessary rounding of that interval. On average, Alg. 6.2 is too pessimistic since it assumes that all floating point Operations have systematic errors, even the ones with exact floating point Operands and results. Moreover, both endpoints of the intervals are always rounded outwards, although on average half the computed endpoints are already properly rounded. How- ever, considering the result Of Thm. 6.1, the presented algorithm is optimal for soft- ware simulated directed rounding. In order tO implement Alg. 6.2, it is necessary to know the machine dependent numbers )1 and e, which can be calculated with the algorithms 6.3 and 6.4. 183 Algorithm 6.2 Directed rounding of the two interval endpoints Of [2:1, 272]. if 2:1 2 6 then $1 4: $1 ® Md else if (121 5 —6 then $1 <= $1 ® Mu else if x1 > 0 then $1 4: 0 else $1 <= —0 end if if $2 2 0 then $2 ¢ $2 ® Mu else if x2 3 —6 then $2 4: $2 ® Md else if x2 _>_ 0 then $2 4: (5 else $2 <= 0 end if Algorithm 6.3 Portable determination Of the machine constant e at run time. 8 <= 1 while 1 + 6/2 > 1 do 8 <= 5 / 2 end while A_lgorithm 6.4 Portable determination Of the machine constant p at run time. u 4: 1 while u/2 > 0 do H <= M2 _end while 184 Thus, by combining Alg. 6.2 with the two algorithms 6.3 and 6.4 we can indeed implement a truly portable interval library with fully rigorous directed rounding. The performance Of the combined algorithms will be discussed in Sec. 6.1.4. 6.1.3 Implementation of Interval Operations In Sec. 2.2 we noted that the intrinsic functions can be extended to intervals in such a way that the basic inclusion property is always maintained: f([a1bll = {f(ib‘) I a S a: S b}- (6-19) Here we discuss for some Of the intrinsic functions how the extension to intervals can be defined such that inclusions are indeed propagated. The interval extensions of the standard mathematical functions are not based on their respective Taylor expansions, but rely on domain decompositions such that the functions are monotonic over the subdomains. The following two proposition show how the monotonicity Of certain functions can be used to compute simple and accurate interval extensions for monotone functions: Proposition 6.2. For any interval [a, b] C R, the interval extension of the exponen- tial function is defined as exp ([0131]) = [exp(a), exp(b)l (5-20) and satisfies the basic requirement (6.19). Proposition 6.3. If we assume that a > 0, then the interval extension of the logo- rithm is defined as 10909.0) = [10800110800] (6.21) and satisfies the basic requirement (6.19). 185 Alg. 6.5 shows an algorithm for extending the sine function to interval arguments [a1,a2]. Special attention is given to arguments that are too large for an accurate determination of the sine function. However, since the interval [—1, 1] is always a correct, albeit overly pessimistic enclosure for the correct result, even this case can be handled rigorously. Algorithm 6.5 Reliable interval extension of the sine function. n2 <= 2n if max(|a1|, [a2|) < 106 then perform outward rounding on the interval [a1, a2] and store the result in [b1, b2] r1 4: nearest integer to ((b1 ® W2) 6 1/ 4) 89 7T 2 if ()1 > T1 then 7‘1 <= 71 EB 7T2 end if r2 4: nearest integer to ((bg ® fig) 03 1/4) (8 W2 if b2 < r2 then 7‘2 4': T2 8 7T2 end if if b1 3 r2 and b2 2 r1 then C1 4: —1 C2 <= +1 else if b1 5 r2 then c1 4: min (sin(bl), sin(b2)) C2 <= +1 round cl towards —00 else if b2 2 r1 then C1 <= —1 02 4: max (sin(bl), sin(b2)) round c2 towards +00 else c1 <= min (sin(bl), sin(bzl) c2 42 max (sin(bl),sin(b2)) perform outward rounding on the interval [61. C2] end if else C1 ¢ —1 C2 4: +1 end if return the interval [c1, C2] 186 6.1.4 Benchmarks and Results Here we present several benchmark results Of the interval implementation in the COSY Infinity [26] language environment, which is based on the previously presented algorithms. The results show that the implementation is very efficient and introduces only a small overhead in the interval Operations Of COSY Infinity. We chose the multiplication Of intervals as the benchmark Operation, since its computational cost is similar to the one Of the presented software-based rounding algorithm. For the actual implementation, we have measured the execution times Of 9 - 106 interval multiplications that test all nine cases of the interval multiplication. Thus, the results allow a very good comparison between the execution times and the imposed overheads. For additional comparison, we also compared the new code tO the previous implementation Of the interval routines in COSY Infinity, which did not provide any rounding and. were therefore not fully verified. Altogether, we tested the following three cases: A: The code has been executed with an unmodified version Of the original interval routines Of COSY Infinity. B: The code has been executed with the new interval routines, but rounding was disabled. C: The code has been executed with the new interval routines and outward round- ing was enabled. The results Of A and B allow a comparison between the new and the Old interval libraries in COSY Infinity. A comparison between B and C on the other hand, puts a measure on the performance overhead imposed by the software rounding. 187 First we tested the new routines in a pure Fortran 77 program, to eliminate the overhead Of the COSY Infinity runtime environment. The program had tO be compiled without any optimization, thus the execution times should be viewed relative to each other, rather than in absolute numbers. Platform, Compiler Test A Test B Test C A /B C/B VMS/Alpha, Digital F77 9.16 s 4.63 8 13.39 s 1.98 2.89 Linux/Intel, G77 9.75 S 8.44 8 14.77 S 1.16 1.75 Solaris/Sparc, Forte F77 8.29 s 6.22 s 11.03 s 1.33 1.77 Digital Unix/Alpha, Digital F77 2.80 s 1.40 s 3.10 s 2.00 2.21 Table 6.1: Average execution times for the three interval test cases in pure Fortran 77 environments. Next we computed the interval products in the COSY Infinity language environ- ment and the corresponding execution times are given in the next table. This time we used the highest Optimization levels available. Thus the numbers should not be compared to the previous benchmarks. Platform, Compiler Test A Test B Test C A /B C/B VMS/Alpha, Digital F77 28.97 8 25.55 S 58.28 S 1.13 2.28 Linux/Intel, G77 24.71 3 26.24 8 31.72 s 0.94 1.21 Solaris/Sparc, Forte F77 17.27 8 18.32 8 19.81 s 0.94 1.08 Digital Unix/Alpha, Digital F77 7.24 s 7.63 s 8.77 s 0.95 1.15 Table 6.2: Average execution times for the three interval test cases in COSY Infinity language environments. It should be noted that the new interval routines have a different calling scheme in the COSY Infinity language environment than the Old interval routines: compared to the case A, B requires one additional subroutine call, while C requires two additional calls. The plain Fortran programs on the other hand have only one additional call in the case Of C and the same number Of subroutine calls in the cases of A and B. Since the new call graph leads to an improved maintainability Of the code, the additional 188 overhead is however acceptable. It should also be noted that the actually computed interval results Of A and B are always equal, but unverified. The case C on the other hand produces correctly rounded enclosures for the mathematically correct results in all situations. 6.1.5 Summary In this section we have presented the basic ideas behind rigorous interval arithmetic implementations. In particular, we introduced the concept Of directed rounding, which can be used to eliminate truncation and underflow errors from conventional floating point Operations. We also outlined how hardware support can in principle be utilized for outward rounding. However, the main focus of the discussion has been the presentation of a portable algorithm for the rigorous implementation Of directed rounding. To illustrate how the use Of directed rounding can avoid the problem illustrated in (6.11), consider the following similar example with directed rounding: ([264, 264] + [_1,1]) _ [264, 264] z [264 J" 264 T] _ [264, 264] (622a) : [264 _ 2111 264 + 211] _ [264,264] (622b) = [—211 — 2'42, 211 + 2‘42] . (6.22c) Although the last interval is a rather pessimistic enclosure Of the mathematically correct result of [—1,1], it is nevertheless a guaranteed enclosure. And it is easy to see that the result is in fact the best enclosure one can Obtain for this problem on a computer with T = 53. Generall, overly large final interval enclosures Often indicate that the program, although mathematically correct, is stated in a numerically unstable form and should be reformulated. 189 6.2 Language Independent Programming Scientific software is commonly deve10ped in high-level programming languages like C [69], C++ [132], and Fortran 90/ 95 [89]. More recently the use of Java has been discussed for the use in scientific computing. Despite its performance problems, the main advantage of Java lies in its combination of the object oriented programming paradigm with the promise of portability across a wide range of platforms; while almost all other languages have undergone standardization efforts, none has ever fully achieved true portability between systems. Although all compiled programming languages with support for floating point numbers are more or less equivalent, for the purposes of computational mathematics, the portability of software is an issue of great concern and often guides important design decisions: compilers are not necessarily available all platforms and legacy applications need to be extended. Thus, it is natural to ask for a language independent programming model. The obvious approach to this problem is to write software in a meta-language and use code generators to derive native language programs from this high-level description. While we are not aware of any large scale system developed in this model, we have been able to ad0pt a similar approach for the development of the C++ and Fortran 90/95 interfaces to COSY Infinity. The basic ideas, results, and experiences with this approach will be discussed in this section. 6.2.1 The Least Common Denominator Approach The COSY Infinity language environment offers an object oriented approach to ad- vanced numerical data types. However, access to these data types has traditionally required using the COSY Infinity programming and run time environment. This restriction has often made it difficult to interface COSY Infinity’s data types and 190 algorithms with existing software packages, which are likely to be written in compiled languages like C++ and Fortran 90/95. The C++ and Fortran 90/ 95 interfaces to COSY Infinity offer a solution to this problem: the flexibility of a modern object- oriented language combined with the power of the high performance data types and algorithms of COSY Infinity. While it can be argued that Fortran 77 this is not the most appropriate language for the deveIOpment of a large-scale software package — version 8.1 of COSY Infinity consists of more than 35,000 lines of highly optimized source code — the use of For- tran 77 also has undeniable advantages over the deployment of more recent languages like C++, Java, or Fortran 90/ 95: Availability Together with C, Fortran 77 is one of the most widely available pro- gramming languages. Virtually all computer platforms have a Fortran 77 com- piler. Performance Fortran 77 is a very simple language with low runtime overhead, re- sulting in fast execution times and outstanding resource utilizations. Legacy Many existing and well tested software packages are written in Fortran 77, providing programmers with a large pool of reusable software components. Fortran 77 is in fact the least common denominator of most modern programming languages, since it has no support for operator overloading, dynamic memory allo- cation, and low level hardware operations like bit manipulation and pointers. Thus, any program written in Fortran 77 can, at least in principle, be mapped one-to-one to languages like C, C++, and Fortran 90/ 95. Since the Fortran family of programming languages is backwards compatible, this mapping is trivial for the Fortran 90 / 95 en- vironments. Moreover, with the availability of the F2C converter [43], we also found 191 the appropriate tool for a conversion from Fortran 77 to C and C++ source code. COSY Infinity is a programming environment that uses a byte-compiler approach similar to Java: source code is translated into an internal byte code, which is then exe- cuted by an interpreter. Thus, the COSY Infinity system consists of three main parts: the compiler, the interpreter, and the support libraries implementing all the available language features. Not much would be gained by simply converting Fortran 77 code to C++ source or renaming the Fortran 77 source files into Fortran 90/ 95 source files. However, both C++ and Fortran 90/ 95 offer object oriented language features, which allow the definition of new data types, and have been used for the implementations of the interfaces to COSY Infinity. This approach offers: 0 native language support via classes (C++) and Modules (Fortran 90/ 95), a thin and lightweight object oriented wrappers with low overhead, 0 embedding of the COSY interpreter into Fortran 90/ 95 and C++ programs, 0 operator overloading for transparent transition from built-in numerical data types to COSY objects. The interfaces disregard the byte-compiler part of COSY Infinity and embed the interpreter to gain access to the support libraries of the system. Since all advanced data types of COSY Infinity are implemented as part of the support libraries, they are immediately available to C++ and Fortran 90/ 95 programmers. We also mention that the system has been implemented in a way that completely removes the need for any manual code manipulations in any language other than Fortran 77. Thus, at the core of the least common denominator (LCD) approach lies a set of Fortran 77 source files that implement all algorithms and data structures and all code for the interfaces is generated from this uniform code base in an automated 192 way. Therefore, the interfaces immediately benefit from any enhancements in the COSY Infinity system. In the remainder of this section we discuss the C++ and the Fortran 90/ 95 in further details and present their application programming interfaces (APIs). 6.2.2 The C++ Interface to COSY Infinity The C++ interface is implemented through the Cosy class, which offers access from within C++ to the core of COSY Infinity. This interfacing is achieved by embedding the COSY Infinity execution engine into a C++ class. Since the glue that holds the two systems together is a very thin wrapper of C++ code, the performance of the resulting class is comparable with the performance of COSY Infinity itself and exceeds that of other approaches (c.f. Sec. 6.2.4 for further details). The COSY Infinity programming language [20] uses an object-oriented approach to programming which centers around the idea of dynamic typing: all Cosy objects have an internal type, which may be real, string, logical, etc., and the exact meaning of operations on Cosy objects is determined at runtime and not at compile time. The Cosy class attempts to be compatible with the C++ double precision data type. In most cases, it should be possible to convert an existing numerical program to a Cosy-based one by simply replacing the string “double” with the string “Cosy” in the source. However, using this approach would underutilize the Cosy class, which shows its real strengths if the advanced data types like intervals, DA vectors, or Taylor models are used. For example, replacing the double precision numbers in an existing program with Cosy objects that are initialized to DA vectors would allow high-order sensitivity analysis of the original program. Other benefits lie in the automatic veri- fication of existing programs by using intervals or Taylor models. 193 Memory Management The Cosy class manages its own internal memory and does not use dynamic allocation of memory by either malloc or new. To a large extent, this is the reason for the performance advantage that COSY Infinity has over languages like C and C++. As a consequence of this, every Cosy object requires a small portion of space in some non-dynamic memory region. While this is never an issue with global and local variables, this becomes an issue when Cosy objects are created dynamically by using new or new []. Consequently, dynamic allocation of Cosy objects should be avoided whenever possible. If Cosy objects really have to be created dynamically, care should be taken to delete the objects as soon as possible, or the COSY system will exhaust its internal memory. Constructors To allow an easy conversion of existing code from the double data type to the Cosy data type, several constructors have been defined that should accommodate this through a variety of implicit constructions. Together with the built-in type con- versions of C++, this mechanism should be able to handle almost any situation correctly. Cosy( ); The default constructor creates a Cosy object with enough internal space to store one number or character. The object’s type is initialized to RE (real) and its value is set to zero. Cosy( const double val, int len = 1 ); Create a Cosy object with enough internal space to hold len numbers or characters. The parameter 1en is Optional and defaults to 1. The object’s type is initialized to 194 RE and its value is set to val. Cosy( const int val, int len = 1 ); Create a Cosy object with enough internal space to store len numbers or characters. The parameter 1en is optional and defaults to 1. The type of the object is initialized to RE, since COSY Infinity does not have a dedicated data type for integers; its value is set to val. Cosy( const bool f ); Create a Cosy object with enough internal space to store one number or character. The object’s type is initialized to L0 (logical or boolean) and its value is set to the boolean value f. Cosy( const char *str ); Create a Cosy object from a C string str. The object’s type is set to ST (string) and enough internal memory locations are allocated to hold the string (without the terminating NULL character, which is not needed in COSY). The object is initialized with the string str. Cosy( const Cosy& src ); Create a new Cosy object from an existing one. The new object is initialized with a deep copy of src. Cosy( integer 1en, const int 11, const int dimI] ); This special constructor creates a Cosy object that represents a Cosy array of dimen- sionality n. The length of each of the dimensions is given in the array dim. And each entry of the array has internal space for len numbers and is initialized to zero with tYpe RE. For further details on Cosy arrays, refer to Sec. 6.2.2. 195 Assignment Operators The Cosy class supports all assignment operations available in C++. Moreover, all the assignment operations that are commonly used with floating point numbers are implemented in a way compatible with the standard C++ definitions for floating point data types. Cosy& operator =(const Cosy& rhs) Assign a deep c0py of rhs to the object and return a reference to it. Cosy& operator+=(const Cosy& rhs) Add rhs to the object and return a reference to it; equivalent to :1: = a: + rhs. Cosy& operator-=(const Cosy& rhs) Subtract rhs from the object and return a reference to it; equivalent to a: = :c — rhs. Cosy& operator*=(const Cosyk rhs) Multiply the object with rhs and return a reference to it; equivalent to a: = a: a: rhs. Cosy& operator/=(const Cosy& rhs) Divide the object by rhs and return a reference to it; equivalent to a: = 112/rhs. Cosy& operator&=(const Cosy& rhs) Unite the object with rhs and return a reference to it. For numerical Cosy objects, the result of a union is usually a vector. It should be noted that this implementation of this operator is not compatible with the default behavior of this operator in C++. Unary Mathematical Operators The Cosy class supports all unary operators available in C++. The operators are compatible with the default implementations for floating point variables. 196 Cosy operator+() Return the positive of the object. This is in fact an identity operation and is included only for completeness. Cosy operator-() Return the negative of the object without modifying it. Cosy operator++() Add one to the object and return the result. Cosy Operator--() Subtract one from the object and return the result. Cosy operator++(int) Add one to the object and return a c0py of the object before the operation. Cosy operator--(int) Subtract one from the object and return a copy of the object before the operation. Array Access Cosy get(const int coeffI], const int n) Obtain a copy of an array element. The element is described by the n-dimensional array coeff. More details on Cosy arrays are provided in Sec. 6.2.2. void set(const Cosyh arg, const int coeff [] ,const int n) COpy the Cosy object arg into an array. The target element is described by the n-dimensional array coefl'. More details on Cosy arrays are provided in Sec. 6.2.2. 197 Printing, IO, and Streams As indicated earlier, the code for the Cosy class is automatically derived from For- tran 77 code by using the F 2C [43] converter. Consequently, the IO handling of the underlying C code is conceptually closer to the printf-type ideas of C than it is to the streams of C++. However, by using temporary files, the Cosy class has partial support for the stream based IO of C++. This mechanism uses the file COSY.TMP in the current working directory as a translation buffer. This allows the Cosy class to be compatible with output streams. friend ostreamh operator<<(ostream& 3, const Cosy& src) Print a representation of the object src onto the ostream s. The printing uses the formats specified in [20]. Type Conversion While the implicit type conversion mechanisms of C++ allow a transparent transition from the default C++ data types to Cosy objects. The conversion of Cosy objects into standard C++ data types on the other hand requires use of the dedicated conversion functions listed below. friend double toDouble(const Cosy& arg) Return a double precision variable that represents the result of calling the function CONS on the Cosy object arg. friend bool toBool (const Cosy& arg) Return a boolean variable that contains the boolean value of the Cosy object arg. If arg is not of type LO, the return value is undefined. 198 friend string toString(const Cosytz arg) Return a C++ string object that contains the string contained in the Cosy object arg. If arg is not of type ST, the result is undefined. Elementary Operations and Functions The COSY Infinity environment has a large number of operators and functions built into its language. The C++ interface to COSY Infinity aims to give transparent access to these functions by trying to be compatible with both the notations of C++ and of COSY Infinity. To that end, the Operators are compatible with the C++ no- tations, and the elementary functions are compatible with the standard C++ naming conventions, and almost all functions defined in math.h for double precision floating point numbers are supported for Cosy objects. As a general rule, all functions in C++ are named with the lower case version of their corresponding COSY Infinity identifier. However, whenever COSY Infinity uses a name for a function that does not exist in C++, e.g.,the absolute value function is called obs in COSY Infinity, while it should be called fabs in C++, both names are made available. Whenever the name of a COSY function clashes with reserved words of C++, the first letter of that function’s name is capitalized: e.g., the COSY Infinity function REAL is called Real in the C++ interface. For the operators defined in COSY Infinity, the following deviations from these general rules exist: 0 While the exponentiation is an operation in COSY Infinity, C++ uses the func- tion pow(. . .) for this. 0 The operator # of COSY Infinity is not defined in C++ and has been replaced with the standard C++ operator !=. 199 o The operator & does not follow the standard C++ conventions and computes the union Of two Cosy objects. However, since the Cosy class is meant to be used for the development of new programs, or as a replacement for double variables, overloading this operator is unlikely to cause any problems. The signatures Of all Operators and functions derived from the COSY operators are: Cosy Operator+ Cosy Operator- Cosy operator* Cosy operator/ Cosy pow bool operator< bool Operator> bool operator== bool Operator!= Cosy operatort bool Operator<= bool operator): The standard functions defined for the Cosy class are listed in Appendix 8.1. These functions are also referred to as intrinsic functions for Cosy Objects. COSY Procedures The COSY Infinity language environment has several procedures built into its lan- guage. These procedures range from diagnostic tOOls (e.g., MEMFRE) over file han- dling tO complex tasks (e.g., POLVAL). For a complete interface from C++ to COSY Infinity it was necessary to make these procedures available as void functions. The C++ interfaces to the procedures all have a standardized signature: void (...); All procedures take at least one argument, and all arguments are either Of type Cosy at or const Cosy &. The list in Appendix B.2 shows the names of the COSY Infinity procedures in the first column, and the declaration Of the corresponding C++ 200 functions. Further details on the corresponding COSY Infinity procedures are given in [20]. Cosy Arrays vs. Arrays Of Cosy Objects In the COSY Infinity language environment, arrays are collections of Objects that may or may not have the same internal type. Thus, within COSY Infinity, it is conceivable to have an array with entries representing strings, intervals, and real numbers. In that sense, the notion of arrays in COSY Infinity is quite similar to the notion of arrays Of Cosy Objects in C++. However, there is a fundamental difference between the two concepts: a C++ array Of Cosy Objects is not a Cosy Object. Due to this difference, the C++ interface does not use C++ arrays Of Cosy objects (although the user obviously has the freedom to declare and use them). As a consequence, the interface provides two different (and slightly incompatible) notions of arrays. Arrays of Cosy Objects are C++ arrays and they can be used wherever C++ permits the use of arrays. Cosy Arrays, on the other hand, are individual Cosy Objects which themselves contain Cosy Objects. Since several important procedures Of COSY Infinity assume their arguments to be Cosy arrays, Cosy arrays are quite important in the context of COSY Infinity and its C++ interface. Since the C++ interface to Cosy does not use the [] Operator for the access to elements, users should use the utility functions Cosy get(const int coeffI], const int n) and void set(const Cosy& arg, const int coeffI], const int n) described in Sec. 6.2.2 to access the elements Of a Cosy array. To simplify the access 201 to individual array elements, we suggest that users use inheritance or external utility functions for convenient access to the elements Of Cosy arrays. For two-dimensional arrays, such a function could be imlemented as: Cosy get(Cosy ta, int 1, int j) { int c[2] = {1+1. j+1}; return a.get(c, 2); Since Cosy arrays start at one, (as Opposed to C++ arrays that start at 0), these utility functions could also be used to mask this implementational detail from the user. However, since the user’s requirements on the dimensionality Of Cosy arrays vary widely, the distribution Of the C++ interface does not provide any of these convenience functions. Finally, we point out that the two different concepts of arrays lead to the possibility Of having C++ arrays Of Cosy arrays — although it might be challenging to maintain a clear distinction between the various indices needed to access the individual elements. 6.2.3 The Fortran 90/ 95 Interface to COSY Infinity The Fortran 90/ 95 interface to COSY Infinity gives Fortran 90/ 95 programmers easy access tO the sophisticated data types Of COSY Infinity. The interface has been implemented in the form of a Fortran 90/ 95 module, which has recently been used as part of the GlObSOl [46] system for constrained global Optimization. Special Utility Routines The Fortran 90 / 95 interface tO COSY Infinity uses a small number Of utility routines for low-level access to the internals. In this section we describe these routines in detail. 202 SUBROUTINE COSYJNIT [] [] [] Initialize the COSY system. This subroutine has to be called before any COSY objects are used. NTEMP sets the size Of the pool Of temporary Objects and defaults to 20. This pool of variables is used for the allocation of temporary COSY Objects. Since Fortran 90/ 95 does not support automatic destruction of Objects, it is necessary to allocate all temporary Objects beforehand and never deallocate them during the execution of the program. The pool is organized as a circular list; and in the absence of automatic destruction Of Objects, if the number of actually used temporary variables ever exceeds NTEMP, memory corruption will occur. It is the responsibility Of the user to set the size appropriately. NSCR defaults to 5000 and sets the size Of the variables in the pool. Additionally, the subroutine SCRLEN is called to set the size of COSY’s internal temp variables. MEMDBG may be either 0 (no debug output) or 1 (print debug information on memory usage). It should never be necessary for users Of the Fortran 90/95 module to set MEMDBG. Neither the size Of the pool, nor the size Of the variables in the pool can be changed after this call. More details on the pool Of temporary Objects are provided in Sec. 6.2.3. SUBROUTINE COSY-CREATE [] [] [] [] Create a variable in the cosy core. All COSY Objects have to be created before they can be used! This routine allocates space for the variable and registers it with the COSY system. SELF is the COSY variable to be created. LEN is the desired size Of the variable SELF, which determines how many DOU- BLE PRECISION values can be stored in SELF, and defaults to 1. If VAL is given, 203 the variable is initialized to it; otherwise VAL defaults to zero. Independent of the parameters LEN and VAL, the type of the variable is set to RE. This routine can also be used for the creation Of COSY arrays, which are discussed in Sec. 6.2.3. If NDIMS and DIMS are specified, the variable SELF is initialized to be an NDIMS-dimensional COSY array with length DIMS(I) in the i-th direction. Each entry of the array has length LEN and is initialized to VAL with type RE. SUBROUTINE COSYDESTROY (SELF) Destruct the COSY Object SELF and free the associated memory. If SELF hasn’t been initialized with COSY_CREATE, the results are undefined. SUBROUTINE COSYJRRAYGET Return a copy of an element of the array SELF. NDIMS specifies the dimensionality Of the array and IDXS is an array containing the index Of the desired element. More details on COSY arrays are provided in Sec. 6.2.3. SUBROUTINE CUSYJRRAYSET (SELF) Copy the COSY Object ARC into an element of the NDIMS-dimensional array SELF. The target is specified by the NDIMS—dimensional array IDXS which contains the index of the target. More details on COSY arrays are provided in Sec. 6.2.3. SUBROUTINE COSY_GE'I'I'EMP Return the address of the next available temporary Object from the circular buffer of such Objects. While the value Of the returned variable is undefined, the type is guaranteed to be RE. Refer to Sec. 6.2.3 for more details. SUBROUTINE COSYDDUBLE (SELF) Extracts the DOUBLE PRECISION value from the variable SELF by calling the function COSY function CONS. 204 SUBROUTINE COSYLLDGICAL (SELF) Extracts the logical value from the variable SELF. If the type of SELF is not LO, the result is undefined. SUBRDUTINE COSYLURITE [] Writes the COSY variable SELF to the unit IUNI'I‘, which defaults tO 6. This func- tion uses the same algorithms employed by the COSY procedure WRITE, which is discussed in [20]. SUBRUUTINE CUSYQTMP Return a temporary COSY Object initialized with the value ARG, which may be either Of type DOUBLE PRECISION or INTEGER. The main purpose Of this function is for the temporary conversion Of parameters to COSY procedures. As an example, consider the following two equivalent code fragments. They illustrate that the use of the function COSY.'I'MP leads to simpler and less error prone code: TYPE(CDSY) :: A,B,X CALL COSY_CREATE(A) CALL CDSY_CREATE(B) CALL COSY_CREATE(X,2) A 8 2 B = 5 CALL INTERV(A,B,X) CALL COSY_DESTROY(A) CALL COSY_DESTROY(B) TYPE(CDSY) :: x CALL COSY_CREATE(X.2) CALL INTERV(CDSY_THP(2),COSY_THP(5),X) Operators The Fortran 90/ 95 interface to COSY Infinity Offers all Operators that the standard COSY system offers. For the convenience of of the user, additional support functions are provided that allow mixed Operations between built-in data types and the COSY 205 Objects and all Operations behave as expected. A complete list of all the defined operations between COSY Objects and built-in types is given in Appendix C. It should be noted that all Operations involving COSY Objects return COSY Objects. Assignment The Fortran 90 / 95 interface to COSY Infinity provides several assignment Operations that allow an easy transition between built-in data types and COSY Objects. This section lists all the defined assignment Operators involving COSY Objects. COSY LHS = COSY RHS Copies the COSY object RHS to LHS. If LHS hasn’t been created yet, it will be created automatically. DOUBLE PRECISION LHS = COSY RHS Converts the COSY Object RHS to the DOUBLE PRECISION number LHS by calling the function COSY_DOUBLE. LOGICAL LHS = COSY RHS Converts the COSY Object RHS to the LOGICAL variable LHS by calling the function COSYLLOGICAL. COSY LHS = DOUBLE PRECISION RHS Copies the DOUBLE PRECISION variable RHS to the COSY Object LHS. If LHS hasn’t been created yet, it will be created automatically. The type Of LHS will be set to RE. COSY LHS = LOGICAL RHS Copies the LOGICAL variable RHS to the COSY Object LHS. If LHS hasn’t been created yet, it will be created automatically. The type Of LHS will be set to L0. COSY LHS = INTEGER RHS 206 Copies the INTEGER variable RHS to the COSY Object LHS. If LHS hasn’t been created yet, it will be created automatically. The type of LHS will be set to RE. Functions The Fortran 90/95 interface to COSY Infinity supports most of the functions sup- ported by the COSY language environment. The following list shows all the functions that are supported for COSY objects by the Fortran 90 / 95 interface to COSY Infinity. The left column shows the Fortran 90/ 95 name of the function and the right column shows the name Of the function in the COSY Infinity environment. All functions have the exact same names as the corresponding COSY infinity functions discussed in [20]. Subroutines All the standard procedures of the COSY Infinity language environment are available as subroutines from the Fortran 90 / 95 interface to COSY. The names and parameter lists of the subroutines match the names and parameter lists Of the normal COSY Infinity procedures. Automatic argument conversion is not available. That means that all arguments have to be either previously created COSY Objects or temporary COSY Objects Ob- tained from calls to COSYJI'MP. Memory Management The COSY Fortran 90/ 95 module is based on the standard core functions and algo- rithms of COSY Infinity. As such, it uses the fixed size memory buffers Of COSY Infinity for storage of COSY Objects. While this fact is mostly hidden from the user, understanding this concept helps in writing efficient code. When a COSY Object is created by using the routine COSY_CREATE, memory is 207 allocate in the internal COSY memory. This memory is not freed until the routine COSY_DESTROY is called for this object. Moreover, since COSY’s internal memory is stack based (and not garbage collected), memory occupied by one Object will not be freed until all Objects that have been created at a later time have also been destroyed. Since Fortran 90/95 does not have automatic constructors and destructors, all objects have to be deleted manually. While this is generally acceptable for normal Objects, this is impossible to guarantee for temporary Objects. TO allow temporary Objects in the COSY module, a circular buffer Of temp. objects is created when the COSY system is initialized with COSY-INIT. As an example on how the pool of temporary Objects should be used, consider the following fragment Of code that implements a convenience interface to the COSY procedure RERAN. Internally, the function CRAN obtains one Object from the pool for its return value. This avoids the obvious memory leak that would result if it was creating a new COSY Object. FUNCTION CRAN() USE COSY_MODULE IMPLICIT NONE TYPE(COSY) :: CRAN CALL COSY_GETTEHP(CRAN) CALL RERAN(CRAN) END FUNCTION GRAN However, it has to be stressed that the fixed size Of the pool Of temporaries bears a potential problem: there is no check in place for possible exhaustion Of the pool. In other words, the pool has to be sized large enough to accommodate the maximum number Of temp. objects at any given time during the execution of the program. Since this number is easily underestimated, especially for deeply nested expressions, the buffer should be sized rather generously. 208 COSY Arrays vs. Arrays of COSY Objects In the COSY Infinity language environment, arrays are collections Of objects that may or may not have the same internal type. Thus, within COSY Infinity, it is conceivable to have an array with entries representing strings, intervals, and real numbers. In that sense, the notion Of arrays in COSY Infinity is quite similar to the notion Of arrays of COSY Objects in Fortran 90/ 95. However, there is a fundamental difference between the two concepts: a For- tran 90/ 95 array of COSY objects is not again a COSY Object. Due to this difference, the Fortran 90/95 module does not use Fortran arrays Of COSY Objects (although the user Obviously has the freedom to declare and use them). As a consequence, the interface provides two different (and slightly incompatible) notions Of arrays. Arrays of COS Y Objects are Fortran 90/ 95 arrays and they can be used wherever Fortran permits the use Of arrays. COS Y Arrays on the other hand, are individual COSY Objects which themselves contain COSY Objects. Since several important procedures of COSY Infinity assume their arguments to be COSY arrays, COSY arrays are quite important in the context Of COSY Infinity and its Fortran 90/95 interface modules. TO access the elements Of COSY arrays, users should use the utility routines SUBROUTINE COSYJRRAYGET (SELF) and SUBROUTINE COSYJRRAYSET Finally, we point out that the two different concepts of arrays lead to the possi- bility Of having Fortran 90/95 arrays of COSY arrays — although it would be quite challenging tO maintain a clear distinction between the various indices needed to access the individual elements. 209 6.2.4 Performance Analysis To demonstrate the success Of the LCD approach in the creation of C++ and For- tran 90 / 95 interfaces to COSY Infinity, Tab. 6.3 lists the execution times Of a mod- erately complicated DA algorithm in six variables and order ten on various hard- ware/compiler combinations. We have timed the exact same algorithm with the different interfaces and it is important to note that: a) the Fortran 90/ 95 interface is not significantly slower than the standard COSY Infinity system b) the C++ interface is only about two to three times slower than the correspond- ing Fortran code While the C++ interface seems slow, most Of the computational overhead can be attributed to the implementation of the Cosy class that is Optimized for maintain- ability and not performance. Moreover, to the best of our knowledge all other DA package implemented in C++ are at least ten times slower than COSY Infinity. Thus, the C++ interface to COSY Infinity is in fact the fastest general purpose DA pack- age available; and gives users from all areas of the computational sciences access to advanced Taylor model algorithms and high-order verification. COSY, C77 C++ COSY,Fort77 Fort90 Alpha/666, Linux 1:26 4:15 1:11 1:09 1:17 4:18 1:09 1:10 Intel PII/333, Linux 5:36 10:31 5:51 10:31 5:54 10:21 Alpha, True64 1:53 1:54 1:53 1:54 Table 6.3: Execution times in minutes and seconds for a typical mix Of DA Operations on different platforms with different interfaces to COSY Infinity. 210 6.2.5 Summary In this section we have presented the least common denominator (LCD) approach to language independent software develOpment. We have demonstrated its applica- bility by summarizing its utilization in the development Of C++ and Fortran 90/ 95 interfaces to COSY Infinity. While we refrain from recommending the LCD approach with Fortran 77 code in new software projects, our experiences show that the method can be used for the maintenance Of large existing legacy codes that cannot easily be rewritten in new programming languages. However, it has proven invaluable in the design and development Of the C++ and Fortran 90/ 95 interfaces to COSY Infinity, broadening the applicability and availability of advanced Taylor model methods. 211 APPENDIX 212 Appendix A Orbital Elements of Major Planets For the computation Of asteroid orbits presented in Sec. 5.1, the orbital elements of the nine major planets play an important role: they are the basis of the computation Of the planets’ positions and velocities. For reference purposes, the following tables list the orbital elements and their linear time dependences. It should however be noted that approximating the time dependence Of the orbital elements by linear relations is not suflicient to reproduce the accuracy of the ephemeris DE405 [130]. Tab. A.1 lists the orbital elements Of the nine major planets; Tab. A.2 and A.3 list the daily rates Of change Of the orbital elements. The numerical values are listed in [125] and are accurate at the epoch J2000. While the number Of significant digits is reduced in the following tables, the computations have been performed with the full accuracy available. Planet 0 i w a e M Mercury 48.33167 7.00487 29.12478 0.38710 0.20563 174.79439 Venus 76.68069 3.39471 54.85229 0.72333 0.00677 50.44675 Earth —— 11.26064 0.00005 114.20783 1.00000 0.01671 357.51716 Mars 49.57854 1.85061 286.46230 1.52366 0.09341 19.41248 Jupiter 100.55615 1.30530 —85.80230 5.20336 0.04839 19.65053 Saturn 113.71504 2.48446 —21.28310 9.53707 0.05415 317.51238 Uranus 74.22988 0.76986 96.73436 19.19126 0.04717 142.26794 Neptune 131 .72169 1.76917 —86.75034 30.06896 0.00859 259.90868 Pluto 110.30347 17.14175 113.76329 39.48169 0.24881 14.86205 Table A.1: Orbital Elements Of the major planets (angles in degree, distances in astronomical units) at the epoch J 2000. 213 Planet (2' i’ w’ Mercury —0.3394174461 - 10‘5 —0.1787968669- 10'6 0.7756255234- 10'5 Venus —0.7581489085 - 10‘5 —0.2175070345 - 10‘7 0.6754049719- 10‘5 Earth —0.1386284128° 10‘3 —0.3569853221 - 10’6 0.1477415013- 10‘3 Mars —0.7758688874 - 10‘5 —0.1937029431 - 10"6 0.1962864091 - 10‘4 Jupiter 0.9256749562 ' 10'5 —0.3156133568- 10‘7 —0.2868963411 - 10“5 Saturn —0.1210015971 - 10‘4 0.4646741210 - 10‘”7 —0.2721423688 - 10"5 Uranus 0.1278728420- 10‘4 —0.1589474479- 10‘7 —0.2805080227- 10’5 Neptune —0.1150277600 - 10’5 —0.2768271345 - 10‘7 —0.5271731676- 10‘5 Pluto —0.2838999240° 10‘6 0.8418891347- 10‘7 —0.7218799621 - 10"6 Table A.2: Daily rates Of change Of the orbital elements (2, i, and w Of the nine major planets . Planet a’ e’ M’ Mercury 0.1806982342- 10‘10 0.6918549067~ 10"9 4.092334433949377 Venus 0.2518818487- 10’10 —0.1351950719- 10’8 1.602131301695948 Earth —0.1368904989- 10‘ll —0.1041478438~ 10‘8 0.985599987451460 Mars -0.1977002118- 10‘08 03258590009 10‘8 0.524021165107646 Jupiter 0.1662888405- 10‘07 —0.3526351815- 10'8 0.083080374325026 Saturn —0.8255441486- 10‘07 —0.1006488706- 10"7 0.033485450148293 Uranus 0.4162217593- 10‘07 —0.5242984255o 10‘8 0.011721311354533 Neptune —0.3427679829 - 10‘07 0.6882956878- 10’9 0.005987479199909 Pluto —0.2105736030 - 10'07 0.1770020547 ~ 10’8 0.003976577306242 Table A.3: Daily rates of change of the orbital elements a, e, and M Of the nine major planets. 214 Appendix B The COSY C++ Interface In this section we describe the intrinsic functions and procedures of COSY Infinity that have been exported to the C++ interface. The definitive reference for these information is given by [20]. B.1 Available COSY Functions The following table lists all COSY Infinity functions that are available from the C++ interface. The first column lists the name of the COSY function as described in [20] and the second column shows the complete signature of the corresponding C++ function. TYPE Cosy type ( const Cosytz x ); LENGTH Cosy 1ength( const Cosytz x ); VARMEM Cosy varmem( const Cosytz x ); VARPOI Cosy varpoi( const Cosytz x ); NOT Cosy not ( const Cosytc x ); EXP Cosy exp ( const Cosy& x ); LOG Cosy log ( const Cosy& x ); SIN Cosy sin ( const Cosytt x ); COS Cosy cos ( const Cosyl’t x ); TAN Cosy tan ( const Cosyh x ); ASIN Cosy asin ( const Cosytt x ); ACOS Cosy acos ( const Cosytc x ); ATAN Cosy atan ( const Cosytr x ); SINH Cosy sinh ( const Cosytz x ); COSH Cosy cosh ( const Cosy& x ); TANH Cosy tanh ( const Cosytr x ); SQRT Cosy sqrt ( const Cosy& x ); ISRT Cosy isrt ( const Cosytc x ); SQR Cosy sqr ( const Cosy& x l; 215 ERF WERF ABS ABS NORM CONS IN WIDTH REAL IMAG INT NIN T NIN T DA CMPLX C ON J Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy Cosy erf werf fabs abs norm cons re in width Real Imag Int rint nint da cmplx conj const const const const const const const const const const const const const const const const const Cosyt Cosy& Cosy& Cosyt Cosy& Cosyk Cosy& Cosyh Cosyh Cosy& Cosyl Cosy& Cosy& Cosyk Cosy& Cosy& Cosy& NNfiNflNNNNNNNHNNNN VVVVVVVVVVVVVVVVV .0 .0 .0 .0 .0 .0 .0 .0 no no .0 90 no no .0 .0 B.2 Available COSY Procedures The following table lists all COSY Infinity procedures that are available from the C++ interface. The first column lists the name Of the COSY procedure as described in [20], while the second column shows the complete signature Of the coresponding C++ function (for readability, c denotes a const Cosy A argument, while v stands for Cosy & arguments). MEMALL void memall ( v ); MEMFRE void memfre ( v ); OPENF void openf ( c, c, c ); CLOSEF void closef ( c ); REWF void rewf ( c ); BACKF void backf ( c ); CPUSEC void cpusec ( v ); QUIT void quit ( c ); SCRLEN void scrlen ( c ); DAINI void daini ( c, c, c, v ). DANOT void danot ( c ); DAEPS void daeps ( c ); DAPEW void dapew ( c, C. c, C ) DAREA void darea ( c, v, c ); DAPRV void daprv (c, c, c, c, c ); DAREV void darev (v, c, c, C. C ); DAFLO void daflo (C C. V. C )3 (SIIFIA) IEUERUAIJ [MAFLAIJ IlAJDIIJ IlAJJEfli IIAJIQTT IIAJILIJ I1AJ?EHE l)AdPEU\ IlAJPEH? IlAJNCTVV CDF2 CHDPJF‘ CHDPIFTlA. CHDPJFTJS AdlfllEflfl IJTQ‘I IJDEXP MBLOCK EHJEMTITL SONSIEEI IiEKDS]? VELSET ’VTELC3EXF VUEZHERK) VUELBJACK VEFILL IhdLHNITT LEUILHE LJWALéHE IDVFEHIV' IIJSIENH) IPIL() IPJLHP I‘FVTELC) I‘fVOELHP IN/SEKF IDVCHEIT OIELEM CHILI (TVOELJELJ ()VODV INTPOL ILILANHAIL ‘EUD‘DAIL void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void void cdflo reran daran dadiu dader daint daplu dapee dapea dapep danow cdf2 cdnf cdnfda cdnfds mtree linv 1det mblock substr stcre recst velset velget vezero velmax vefill imunit ltrue 1false interv insrnd inlo inup ivvelo ivveup ivset ivget Oielem Oiin ovelem oviv intpol rdavar rdvar AAA/\AAAAAAA/‘AAAAAAAAAAAAAAAAAAAAAAAAAAAAA/\AAA c, c, v, c ); v ); v, c ); c, v, v ); c, v, v ); c, v, v ); v, c, c, v ); c, c, v ); c, c, c, v ); v, v, v, v ); c, c, v ); v, v, v, v v ); v, v, v, v v, v, v, v, v, v v, v, v, v, v, v v, v, v, v, v, v v, v, c, v, c, c, v ); c, c, c, v ); c, v, v, c, c ). c, c, c, v ); c, v ); c, c, v ); v, c, c ); c, c, v ); v, v, v ); c, v ); v, c, c, c, c ), v ); v ); v ); c, c, v ). c ); c, v ); c, v ); c, v ); c, v ); v, c, C )i c, c, V ); c, v ); c, c, v ); c, v ); c, c, v ); v, c )3 c, c, c, c, c, c, v ); c, c, c, v ); 217 RDANOT RDREA DAEXT RDRBN D RDAREF RDADOM RDITIG RDNPNT RDIN T RDPRIS CLEAR GRMOVE GRDRAW GRDOT GRCURV GRPROJ GRCHAR GRCOLR GRWDTH GRMIMA RKCO POLVAL POLSET FOXDPR MEMWRT VARMAX void void void void void void void void void void void void void void void void void void void void void void void void void void rdanot rdrea daext rdrbnd rdaref rdadom rditig rdnpnt rdint rdpris clear grmove grdraw grdot grcurv grproj grchar grcolr grwdth grmima rkco polval polset foxdpr memwrt varmax AAA/NAAA/\AAAAAAAAAAAAAAAAAA . . . . . VVVVV no 0.- no no .- ~<<<<