. I _ - . ‘ l l H . . , _~ _ - ~ - . . ' ~ - ' ' . ‘ ' 1 ~ 1 . . , V . . i A ' ., , , , A . . '- . ' ‘ .’ - ~ . . ' , _ ‘ :- ‘ . . ‘- . ‘ , - , V - , . ‘ ‘ y . - . .' _ ' 7 4 ‘ . . ‘ ‘ r , . ‘ .v‘ - , I : , ‘ ‘ , 4 0 - . ' . . . o ' ‘ A ' V . I . n A ' . . - I, . ‘ YA . .‘ > ‘y a .... n..-':;.,.:|ALAMIL'J- q“ .HHJ‘UH" ‘“ THESIS .‘J VERSITY LIBRARIE Illiilll’ilill‘nllu mm" Illlll 3 1293 01682 6301 This is to certify that the dissertation entitled RIGOROUS ANALYSIS OF NONLINEAR MOTION IN PARTICLE ACCELERATORS presented by Kyoko Makino has been accepted towards fulfillment of the requirements for Ph.D. degree in Ehygics Major professor Date 3!“:(01 9 MS U is an Affirmative Action/Equal Opportunity Institution 0- 12771 LIBRARY I Mlchlgan State University PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE M52 0 2002*? 1/98 animus-nu RIGOROUS ANALYSIS OF NONLINEAR MOTION IN PARTICLE ACCELERATORS By Kyoko Makino A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Physics and Astronomy 1998 ABSTRACT RIGOROUS ANALYSIS OF NONLINEAR MOTION IN PARTICLE ACCELERATORS By Kyoko Makino Methods to obtain rigorous descriptions for the motion of ensembles of particles relative to a reference curve are developed. To account for the most general refer- ence curves, a Lagrangian and Hamiltonian formulation of the relativistic motion of charged particles in electromagnetic fields in three dimensional curvilinear coordinates including torsion and with arclength as independent variable is derived. To allow for the use of the most general fields of particle optical devices, a wavelet-based method to include measured field data in the equations of motion in an appropriate manner is discussed. The method of transfer maps is a powerful tool for the study of weakly nonlinear dynamical systems, especially for the case of large phase space acceptances as in modern particle Spectrographs, and the intricate dynamical behaviour of repetitive systems as in circular accelerators and storage rings. The Differential Algebraic (DA) techniques have proven fruitful for various computational problems in beam physics, including the determination of high order Taylor transfer maps. A new approach, the Remainder—enhanced Differential Algebraic (RDA) method, is presented, which extends the method to allow the determination of remainder bounds for functional dependencies and solutions of ODEs. First, the basic theory of the method is developed and applied to problems of verified optimization and quadra- ture. Next, schemes are derived that allow the construction of numerical integrators of arbitrary order with rigorous verification of the error for both the integration of individual initial conditions as well as Taylor transfer maps. The methods are based on a differential algebraic fixed point problem which is studied using Schauder’s the- orem and other functional analysis tools. Employing various compactness arguments on suitable function spaces in combination with the RDA tools, in each integration step a proof of existence of a solution within a tight inclusion is performed. The resulting computational tools are implemented in the arbitrary order beam physics code COSY INFINITY, and their behavior and performance is studied. Using integration orders around ten and suitable step sizes, rigorous remainder bounds in the range of 10‘10 for transfer maps of orders around ten are obtained. Copyright by KYOKO MAKINO 1998 To the ultimate truth and those who love it ACKNOWLEDGMENTS Seeing my Ph.D. thesis work summarized in front of me, my thoughts float around in the memories of these last years. Probably merely a short period of one’s life, it has been a dramatical one for me, who had not imagined to work in the academic society again. The chance was opened by my academic advisor Prof. Martin Berz, who believes in the nobility of science, and who makes great efforts to support and encourage his students. While I had been exposed to scientific research earlier as a graduate student in high energy physics and afterward as a staff scientist working on computer codes, it was quite shocking to me when I saw the way Prof. Berz attacks a scientific question. In 1995 we were working on the problem of the relativistic equation of motion of spin. Prof. Berz questioned the meaning of the first component of the spin four vector, challenging us by asking whether he even could “put his grandmother” there. It then turned out that the axiomatic foundation of this time-like component is indeed rather shaky, and I had to realize that such a pure critical attitude is an important catalyst for real progress in science. The riddle of the spin four-vector continues to have several open subtleties, and unfortunately my dissertation does not discuss this interesting topic. I would like to use the occasion of the completion of my thesis work to thank, first of all, Prof. Berz who has guided me with pure scientific rigor through various vi interesting topics, as well as various practical opportunities to gain exposure to the wider scientific community at conferences and other occasions. My respect is not only tied to his scientific rigor, but also to his continuous support and warm encouragement of his students. In the period of the spin, we were having greatly exciting discussions in the group every day, and it was one of the wonderful moments in my life. The enjoyment of exciting scientific discussions is only possible with the right kind of people. I would like to express my appreciation to those people who share their time in the group of Prof. Berz with me: Dr. Ralf Degenhardt, Dr. Georg Hoffstiitter, Dr. Weishi Wan, Silke Rolles, Dr. Vladimir Balandin, Dr. Nina Golubeva and Meng Zhao, in the past; and at the moment, Khodr Shamseddine, Jens Hoefkens, Bela Erdelyi, Lars Diening, Jens von Bergmann and Michael Lindemann. My recent thanks go to a long-term warm friendship with Khodr Shamseddine, and mighty help from Jens Hoefkens. My life at Michigan State University is supported and encouraged also by those professors who kindly have served on my advisory committee; Prof. J erzy Borysowicz, Prof. Julius S. Kovacs, Prof. Jerry Nolen at Argonne, Prof. Jon Pumplin and Prof. Brad Sherrill. My thesis work has received substantial help from various people at the National Superconducting Cyclotron Laboratory at MSU. In particular, the group working on the S800 spectrograph has supplied valuable data and discussions, and I would like to thank Prof. Sherrill, Jac Caggiano, and Dr. Daniel Bazin, also for his last minute help supplying valuable pictures of the 8800. Several pictures in this dissertation were made using programs by David Johnson. Since my main work has been focused on computer simulations, I have received a lot of help from the computer department. For financial support of this research I owe thanks to the US Department of vii Energy, the Alfred P. Sloan Foundation, and the National Science Foundation, and I am grateful for their continuous support for basic science and their contribution to the search for the ultimate truth. Outside MSU, I have also encountered strong encouragement from various people. Especially I would like to thank Dr. Carol Johnstone at Fermilab, Dr. Christian Bischof at Argonne, the US Particle Accelerator School, Japanese visiting scientists I had the pleasure to meet in the USA, and several scientists at KEK who had shared their time with me back in Japan. My special thanks go to Prof. Seigi Iwata at KEK and Prof. Ryouichi Kajikawa, who were my academic advisors back at Nagoya in Japan, and have strongly supported and encouraged my continued work in the academic society. Before closing, I would like to thank my family who has always loved me. Un- fortunately I had to confront the loss of my father, who had believed in the power in science, in the last December. My mother, who has shared the grief of the loss of my father, has been a big mental support for me through her unconditional love since my childhood. I appreciate my relatives, my sister and my friends for their encouragement and support. Finally my wish is that my beloved son Masayuki and daughter Kazuko may enjoy reading this thesis some time in the future. viii Contents LIST OF TABLES xii LIST OF FIGURES xiv 1 Introduction 1 2 The Particle Optical Equations of Motion 8 2.1 Nonplanar Curvilinear Coordinates ................... 9 2.2 Differential Operators in Nonplanar Curvilinear Coordinates ..... 15 2.2.1 Gradient .............................. 17 2.2.2 Divergence ............................. 19 2.2.3 Curl ................................ 20 2.2.4 Laplacian ............................. 22 2.2.5 Velocity Vector in Nonplanar Curvilinear Coordinates ..... 23 2.3 Dynamics in Nonplanar Curvilinear Coordinates ............ 24 2.3.1 Electromagnetic Fields and Lorentz Force ........... 24 2.3.2 The Lagrangian and Lagrange’s Equations ........... 29 2.3.3 The Hamiltonian and Hamilton’s Equations .......... 33 2.3.4 Arclength as Independent Variable for the Hamiltonian . . . . 39 2.4 Planar Motion ............................... 43 3 Transfer Maps for Elements Characterized by Measured Fields 47 3.1 Wavelet Representations ......................... 49 3.2 Elements Characterized by Measured Fields .............. 53 3.3 Transfer Maps using Linear Field Data ................. 54 3.3.1 The Homogeneous Dipole Field ................. 55 3.3.2 Inhomogeneous Dipole Field ................... 58 3.4 Transfer Maps using Analytical Field Data ............... 61 3.5 Transfer Maps of the S800 Spectrograph ................ 68 ix 3.5.1 Measured Field Data of the S800 Dipole ............ 69 4 Remainder-enhanced Differential Algebra (RDA) Methods 76 4.1 Introduction ................................ 76 4.2 Differential Algebra and Interval Arithmetic .............. 77 4.3 Remainder-enhanced Differential Algebraic Operations ................................ 80 4.3.1 Addition and Multiplication ................... 83 4.3.2 Intrinsic Functions ........................ 85 4.3.3 Derivations and Antiderivations ................. 91 4.4 Examples ................................. 93 4.4.1 A Simple Function ........................ 93 4.4.2 Bound Enclosures of Functions ................. 98 5 Implementation of Remainder-enhanced Differential Algebraic Op- erations in COSY INFINITY 103 5.1 Supported Elements and Features .................... 104 5.2 Data Structure of Taylor Models ..................... 108 5.3 Operations on Taylor Models ...................... 110 5.3.1 Addition and Subtraction, etc ................... 111 5.3.2 Multiplication ........................... 112 5.3.3 Intrinsic Functions ........................ 114 5.3.4 Integral .............................. 115 5.4 Methods to Tighten Order Bound Intervals ............... 116 5.4.1 Implemented Methods ...................... 116 5.4.2 Examples of Performance .................... 118 5.4.3 Further Tightening Methods ................... 123 5.5 Examples of Computation ........................ 130 5.5.1 A Small Multidimensional Function ............... 131 5.5.2 Multidimensional Integrals .................... 133 5.5.3 Bounds of Normal Form Deviation Function .......... 136 6 Verified Integration of ODEs and Flows 139 6.1 Verified Integration with Taylor Models ................. 139 6.1.1 Schauder’s Fixed Point Theorem ................ 140 6.1.2 Strategy to Satisfy the Requirements .............. 141 6.1.3 Schauder Candidate Sets ..................... 142 6.1.4 Convexity, Compactness, and Invariance ............ 142 X 6.1.5 Satisfying the Inclusion Requirement with Differential Alge- braic Methods ........................... 145 6.1.6 Iterative Refinement of the Inclusion .............. 146 6.2 Example: Remainder Bounds for a Dipole of the 8800 Spectrograph . 148 xi List of Tables 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 4.1 4.2 4.3 4.4 5.1 5.2 5.3 5.4 ' Accuracy of the Gaussian wavelets representation for one dimensional functions. ................................. 51 Taylor transfer map of a homogeneous dipole of the radius 2.15m with a bending angle of 26.65° computed geometrically by the default pro- cedure in COSY INFINITY. ....................... 56 Difference between the Taylor transfer maps of a homogeneous dipole computed by the default procedure in COSY INFINITY and computed using the new element based on measured field data. ......... 57 Taylor transfer map of an inhomogeneous dipole of radius 2.15m and bending angle 26.65° with inhomogeneity 0.5 computed geometrically by the default procedure in COSY INFINITY .............. 59 Difference between the Taylor transfer maps of an inhomogeneous dipole computed by the default procedure in COSY INFINITY and computed using the measured field element. .................... 60 Transfer maps of a dipole with fringe field effects computed by default in COSY INFINITY, and using analytically generated field data. . . . 64 A part of the fifth order Taylor transfer map of a S800-like dipole with the fringe field effects computed by the default procedure in COSY INFINITY, and the error of transfer map using analytical field data. . 67 Design parameters of the S800 spectrograph ............... 69 Elementary properties of interval arithmetic; 11 = [(11, b1], 12 = [0.2, (22]. 79 The remainder bound interval I 1 /1=+m for various orders; x0 = 2, [a, b] = [1.9, 2.1] ................................... 95 The width of the bound interval of f (:13) = 1 / :1: +1: by various methods; 11:0 = 2, [a, b] = [19,21]. ......................... 96 The total number of F P operations required to bound a simple function like f(:i:') = 2,-(1/223- + 23,-). ........................ 98 Data types related to RDA in COSY INFINITY. ........... 104 Binary operations related to RDA in COSY INFINITY ......... 105 Intrinsic functions available for RDA in COSY INFINITY. ...... 106 Data structure of a Taylor model (RD data type) in COSY INFINITY. 109 xii 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 6.1 6.2 6.3 6.4 Number of applications of the lower dimensional quadratic function bounders to bound a higher dimensional quadratic function at the boundaries. ................................ 127 Remainder bound intervals of a small multidimensional function for various orders ................................ 132 Bounds estimates of a two dimensional integral with the interval method. 134 Bounds estimates of a two dimensional integral with the RDA method. 134 Bounds estimates of a four dimensional integral with the interval method. 135 Bounds estimates of a four dimensional integral with the RDA method. 135 Bounds estimates of a normal form deviation function with the interval method. .................................. 137 Bounds estimates of a normal form deviation function with the RDA method. .................................. 138 Error bounds of the Taylor transfer map of a S800-like dipole with initial condition within {—0.02, 0.02]4 and step size 1° .......... 149 Error bounds of the Taylor transfer map of a 8800-like dipole with initial condition within {—0.01, 0.01]4 and step size 1° .......... 149 Error bounds of the Taylor transfer map of a S800-like dipole with initial condition within [—0.01,0.01]4 and step size 0.5° ......... 150 Computational behavior of the iteration to find the first solution of the inclusion requirement. .......................... 151 xiii List of Figures 1.1 1.2 2.1 2.2 2.3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 Upper: Tracking pictures of a repetitive system. The left is in particle optical coordinates x and a = pm/po, and the right is in normal form coordinates. Lower: Deviation from a normal form invariant circle. . . Bound enclosures of a function. Upper: RDA method by 7th (left) and 8th(right) order. Lower: Interval method with 25 (left) and 200 (right) subintervals ............................. Reference curve and the locally attached Dreibeins. .......... Non-uniqueness of curvilinear coordinates. ............... Uniqueness of curvilinear coordinates within a tube ........... Measured field of a bending magnet of the S800 spectrograph. The picture shows the field at 65 x 74 points. ................ Gaussian wavelets representation for f (3:) = 1 (upper left), f (2:2) = x+1 (upper right), f (3:) = cosar+1 (lower left) and f (x) = exp(—:r ) (lower right). ................................... Derivatives of the function f (2:) = exp(—:r2) when represented by an ensemble of Gaussian wavelets. The interpolated function (upper left), and the first (upper right), the second (lower left) and the third (lower right) order derivatives of the interpolated function ........... Specification of measured field data for a particle optical element in COSY INFINITY. ............................ Analytical dipole field distribution with fringe field effects. Bend: 26.65°. Radius: 2.15m. Aperture: 10cm. In the upper picture, a beam enters from the center at the left rightward then goes clockwise. The lower picture is viewed from a lower angle. ............... Analytical S800-like dipole field distribution with fringe field effects. In the upper picture, a beam enters from the upper left rightward, then goes clockwise. The lower picture is viewed from a lower angle. . . . . Layout of the S800 spectrograph of the National Superconducting Cy- 10 11 12 48 50 52 54 62 66 clotron Laboratory at Michigan State University. Courtesy Daniel Bazin. 69 Two rectangular areas of the measured field data of the two S800 dipoles. Courtesy Daniel Bazin. ..................... xiv 72 3.9 Field distribution of the 8800 dipole measured data. Maximum: 0.965 Tesla. In the upper picture, a beam enters from the lower left rightward and goes counterclockwise. The lower picture is viewed from a lower angle. ................................... 3.10 Top part of the field distribution of the 8800 dipole measured data. The picture shows the field values ranging from 0.962Tesla to the maximum 0.965Tesla .................................. 3.11 Further detail of the 8800 dipole measured data in a small area 34.5cm x 34.5cm on the top region shows the step structure due to the lim- ited data digits (top). The data is smoothed (middle), and compared (bottom). ................................. 4.1 Function f (3:) = exp(—1/a:2) if :c 75 0 ; 0 else, and its Taylor polyno- mial, which vanishes identically. ..................... 4.2 One dimensional function and its bound enclosures. From top to bot- tom right: the function, the bounds by the interval method with the 25 and 50 divided domain intervals, the bounds by the 7th and 8th order Taylor models. ........................... 4.3 Bound enclosures of a two dimensional function by the interval method. From top to bottom right: The domain is divided to 10 x 10, 20 x 20, 40 x 40 and 80 x 80 subintervals. .................... 4.4 Bound enclosures of a two dimensional function by the RDA method. From top to bottom right: Computations in 7th, 8th, 9th and 10th order Taylor models. ........................... 5.1 The order bounds I 2 = [a_2,b_2] and I 3 = [a_3, b_3] around the linear line co + c_1:r 6 [co + a_1, co + b_1] ..................... 5.2 An upper bound of 153 is above Co + b_1 + a_2 + a_3. Find a point 23(1) to specify a new domain .......................... 5.3 Bounds estimates of a definite integral of a function f by the interval method (left) and the RDA method (right). .............. 6.1 Convexity; a set of numbers M (left) and a set of functions inside a Taylor model (right). ........................... 6.2 Functions on [t0, t1] are uniformly equicontinuous and uniformly bounded. (The Ascoli-Arzela Theorem) ...................... 6.3 Finding an inclusion set as a Taylor model satisfying Schauder require- ment ..................................... 6.4 Iterative refinement of the inclusions as Taylor models. ........ XV 74 75 78 100 101 102 128 129 143 144 146 147 Chapter 1 Introduction Following the general scientific approach, the study of the motion of particles is performed in two steps. The first step is to establish the equations of motion, the second step is to solve the equations and analyze the results. The work in this thesis covers both aspects for the motion of an ensemble of particles in beam physics. The standard way to treat such a system is to pick a typical particle, the reference particle, and to describe the motion of the other particles relative to the reference motion, after which the problem is cast into the framework of weakly nonlinear dynamical systems. This thesis work focuses on the rigorous approaches to analyze such systems. The first elaboration in this thesis is to derive the relativistic equation of motion of a charged particle in electromagnetic fields relative to a general three dimensional reference orbit. In chapter 2, Newton’s equation as well as a Lagrangian and a Hamiltonian describing the motion in nonplanar curvilinear coordinates are derived. The various forms were checked for consistency, and an interchange of the independent variables from time t to arclength s is performed. Besides providing the complete set of equations, the importance of this work lies in the fact that the relative motion with nonplanar reference orbit can still be recognized as Hamiltonian, which has key consequences in the study of this weakly nonlinear dynamical system. One of the by-products of the elaboration is the representation of the Laplacian operator A in nonplanar curvilinear coordinates, which enables the usage of the Differential Algebraic (DA) fixed point theorem [20] to solve the three dimensional potential problem. Another by-product is a set of equations of motion in the more conventional curvilinear coordinates in which the reference orbit is restricted to a plane. At this point, the electromagnetic fields are described by E and B, or the scalar potential (I) and the vector potential A in the equations of motions. When we study and design a beam optical system at a laboratory, the detailed knowledge of the fields sometimes becomes a crucial point to determine the quality of the device. Modern large acceptance spectrographs represent one such example. In chapter 3, a new method to take measured data of the fields into the equations of motion in an appro- priate manner is discussed. The method has to fit to the principal work horse of our study in beam physics, the Differential Algebraic (DA) technique [3] [4] [5] [8] [20], so the interpolation method of the field distribution function has to supply differentia— bility to any order. For this purpose, we used the Gaussian wavelet approach. The method is studied to check the numerical performance, then is applied to the magnetic dipoles of the S800 spectrograph at the National Superconducting Cyclotron Labo- ratory at Michigan State University [60], where the accuracy of the fields is required to about 10"4 in order to achieve its high energy resolution of 1 / 10, 000 [16]. The method of transfer maps has various advantages to study beam optical sys- tems compared to the conventional ray tracing method. It offers direct access to aberrations for single pass systems, and to chromaticities, amplitude tune shifts, and resonance strengths for repetitive systems such as circular accelerators and storage rings. Another topic where the map method is especially powerful is the analysis of the long-term behaviour of particles in repetitive systems. The DA techniques have offered a very elegant and accurate way to obtain high order Taylor transfer maps of the action on phase space. Typically derivatives of up to order ten in six variables are needed to study the long-term behaviour of particles in repetitive systems, so other methods are far from providing a robust way to study the weakly nonlinear behavior of beams. The problem to estimate the long-term stability of weakly nonlinear systems finds its origin in the detailed study of the solar system. Many perturbative methods for repetitive motion have been developed from this question. Recently, ideas of Lyapunov, Nekhoroshev and others triggered an analysis of stability in particle accel- erators based on approximate invariants[66] [14] [39], and the question of long-term stability can be re—cast into a highly complicated optimization problem [14] [15] [39]. The upper left picture in Figure 1.1 shows a tracking picture of particles launched from five different locations in a repetitive system described by a six dimensional Taylor transfer map in the sixth order in actual coordinates a: and a = pin/pg. The upper right picture shows a tracking picture of the same system in normal form coordinates after a nonlinear normal form transformation, where the motion is seen approximately on invariant circles. The deviations from invariant circles in the normal form coordinates directly relate to the number of stable turns, and thus to the time for particles to be lost. The sharpness to which these deviations can be bounded is the key to guaranteeing a large number of stable turns. The deviation functions are multidimensional polynomials up to roughly 500th order, and consist of about 105 floating-point operations. The functions have a very large number of local extrema, many large terms cancel each other, and very small fluctuations have to be estimated. To be useful, the maxima have to be sharp to about 10‘6, and for some applications to 10‘”. The lower picture in Figure 1.1 shows the Figure 1.1: Upper: Tiacking pictures of a repetitive system. The left is in particle px/po, and the right is in normal form coordinates. Lower: Deviation from a normal form invariant circle. optical coordinates x and a deviation from a normal form invariant circle as a function of two phase space angles. The problem of finding rigorous bounds on the extrema of functions has con- tributed to the development of many methods in numerical analysis. In our case, the irregularity of the deviation functions as well as the high dimensionality makes the question very troublesome for conventional optimization methods. Interval methods give a mathematically rigorous estimate of the bounds, but complicated functions like ours cannot avoid the severe blow-up problem, the control of which is the key to get a practical estimate. Usually the division of the domain of interest into many small subintervals helps to suppress the severe blow-up problem. The deviation functions require about 104 subintervals per dimension to be sharply bounded in a small phase space volume 0.026 as discussed in subsection 5.5.3, and the total number of subin- terval boxes goes up to 1024 to cover the whole six dimensional small domain, which limits the practicality. A new technique presented in the following chapters, the method of Remainder- enhanced Differential Algebras (RDA) [52] [9] [12] [54], combines the DA technique to express the model function by a Taylor polynomial, and interval computations to evaluate the bound of the Taylor remainder. The resulting error bounds are usually rather sharp, in particular at higher orders. Figure 1.2 shows bound enclosures of a one dimensional function using a seventh order (upper left) and eighth order (upper right) RDA method, in comparison with the interval method applied to the divided domain, using 25 (lower left) and 200 (lower right) subintervals. The RDA method can be used for rigorous global optimization of highly complex multi-dimensional objective functions. Now a practical answer to the bounds of the deviation functions can be given. Chapter 4 discusses the development of the method as well as some example cal- Figure 1.2: Bound enclosures of a function. Upper: RDA method by 7th (left) and 8th(right) order. Lower: Interval method with 25 (left) and 200 (right) subintervals. culations. The implementation of the method in the arbitrary order general purpose beam optics code COSY INFINITY [11] [15] [51] [10] is discussed in chapter 5. Exam- ple computations in the chapter include a new scheme to compute multi—dimensional integrals, as well as the normal form deviation function. To complete the question, the other important aspect is to have rigorous repre- sentations of transfer maps to describe weakly nonlinear systems. Chapter 6 develops the theory for a rigorous integration scheme for flows of ordinary differential equa- tions within the framework of RDA, which enables to obtain rigorous remainders of Taylor transfer maps[18]. The task requires to start from the foundations in func- tional analysis and Schauder’s fixed point theorem was utilized. The method has been implemented in COSY INFINITY, and some example calculations are shown. Chapter 2 The Particle Optical Equations of Motion In this chapter, we will derive the relativistic equations of motion of a particle in an electromagnetic field in the so—called curvilinear coordinates. These coordinates are measured in a moving right-handed coordinate system that has one of its axes attached and parallel to a given reference curve in space; furthermore, usually the time is replaced as the independent variable by the arclength along the given reference CUI‘VB. While this approach seems to complicate the description of the motion, it has several advantages. Firstly, if the chosen reference curve in space is itself a valid orbit, then the resulting transfer map will be origin preserving, because the origin just corresponds to the reference curve itself. This fact then opens the door to the use of perturbative techniques for the analysis of the motion in order to study how small deviations from the reference curve propagate. In particular, if the system of interest is repetitive and the reference curve is closed, then the origin will be a fixed point of the motion. If the arclength is used as the independent variable, then after one turn around the reference orbit, the system is repetitive, and perturbative techniques around fixed points can be employed to study the one-turn transfer map, which here corresponds to the Poincare map of the motion. Secondly, the method is also very practical in the sense that beams themselves are usually rather small in size, while they often cover large territory. So it is more convenient to describe them in a local coordinate system following a reference parti- cle instead of a Cartesian coordinate system attached to the laboratory. Expressing the motion in terms of value of the transfer map at a given arclength very directly corresponds to the measurements by detectors, which usually determine particle coor- dinates at a fixed plane around the system instead of at a given time. And expressing the motion in terms of the arclength as independent variable directly provides a nat- ural scale, since it is more natural to measure in meters along the system instead of nano— or microseconds. The following sections describe in detail the derivation of the motion in curvilinear coordinates. We will study the transformations of Maxwell’s equations and the re- sulting fields and their potentials to the new coordinates, and then derive the explicit forms of Newton’s equations as well as the Lagrangian and Hamiltonian with time as the independent variable. Finally, a special transformation on Hamiltonians is applied that replaces time as the independent variable by the arclength, while maintaining the Hamiltonian structure of the motion. 2.1 Nonplanar Curvilinear Coordinates Let {61,62,63} denote a Dreibein, a right-handed set of fixed orthonormal basis vec- tors, which defines the so—called Cartesian coordinate systems. For any point in space, let (11:1, 11:2, 11:3) denote its Cartesian coordinates. In order to introduce the curvilinear coordinates, let R(s) be an infinitely often differentiable curve parameterized in terms of its arc length 3, the so—called reference curve. For each value of 3, let the vector é}, 10 _) e y —’ ex R(s) S E’. x 3 Reference Curve X2 Figure 2.1: Reference curve and the locally attached Dreibeins. be parallel to the reference curve, i.e. _d1'i 68(3) — 213-. (2.1) We now choose the infinitely often differentiable vectors é’x(s) and 611(3) such that for any value of s, the three vectors {63, é}, 63,} form a Dreibein, a right-handed orthonormal system. For notational simplicity, in the following we also sometimes - - . —+ .. .. -C -C -C denote the curv111near baSlS vectors {83,835,831} by {e1 ,82 ,e3 }. Apparently, for a given curve 18(3) there are a variety of choices for é'$(s) and é'y(s) that result in valid Dreibeins since (2(3) and é'y(s) can be rotated around é}. A specific choice is often made such that additional requirements are satisfied; for example, if the curve [2(8) is never parallel to the vertical Cartesian coordinate 53, one may demand that é'x(s) always lie in the horizontal plane spanned by 51 and 62. The functions R(s), 61(3), and 631(3) describe the so—called curvilinear coordinate 11 Figure 2.2: Non-uniqueness of curvilinear coordinates. system, in which a position is described in terms of s, a: and y via Apparently the position 7" in Cartesian coordinates is uniquely determined for any choice of (s,:r,y). The converse, however, is not generally true: a point with given Cartesian coordinates F may lie in several different planes that are perpendicular to the curve R(s), as shown in Figure 2.2. The situation can be remedied if the curvature rc(s) of the reference curve R(s) never grows beyond a threshold, i.e. if 1'1 2 1/msax IK.(S)[ (2.2) is finite. As Figure 2.3 illustrates, if in this case we restrict ourselves to the inside of a tube of radius r1 around R(s), for any vector within the tube, there is always one and only one set of coordinates (s, 2:, y) describing the point 7". Let us now study the transformation matrix from the Cartesian basis (61,52, 63} to the local basis of the curvilinear system {63,613, é'y} = {éfi ég, €30}. The transfor- 12 (s) ‘ / ~. Figure 2.3: Uniqueness of curvilinear coordinates within a tube. mation between these basis vectors and the old ones is described by the matrix O(s) '51) '52) ) - (2-3) . 63) which has the form D—i 8 p—a It: AAA 051 m m 8 5,3131 ('01 vvv as 7"st ca1&1 l 8 A CM 551 m 3131 ('01 vvv Because the system {65, 63,, 63,} is orthonormal, so is O(s), and hence it satisfies O(s)-O(3)t=f and O(s)t-O(s)=f. (2.4) Since both the old and the new bases have the same handedness, we also have det(O(s)) = 1, (2.5) and hence altogether, O(s) belongs to the group 80(3). We remind ourselves that elements of 80(3) preserve cross products, i.e. for O 6 80(3) and any vectors ('1', 5, we have (2.6) One way to see this is to study the requirement of orthonormality on the matrix elements of O. The elements of the matrix O describe the coordinates of the new 13 parameter dependent basis vectors in terms of the original Cartesian basis; explicitly, we have [53h = 01:1, [Exlk = 0km [gylk = 0k3- (2.7) The demand of the right-handedness then reads the} X6 365m 3 22E:6 Ehnnena where 6,5,, is the common totally antisymmetric tensor of rank three defined as 1 for (2', j, k) = (1, 2, 3) and any cyclic permutation thereof éijk = —1 for other permutations of (1,2,3) 0 for two or more equal indices and boils down to a condition on the elements of the matrix O 3 3 Z €ijk0i20jm = Z fzmnokn- (2-8) t,j=1 n=l We remind ourselves that the symbol Cijk is very useful for the calculation of vector cross products; for vectors ('1', 5, we have 3 :ijzlé Etjka’ibj Using condition (2.8), we readily obtain (2.6). For the following discussion, it is useful to study how the transformation matrix O changes with s. Differentiating (2.4) with respect to the parameter s, we have z—Ot.O:_ Ot— 0 ( ) 0+ (13 d3 d ~ - dO’ dO_ Ot— dO +O’— dO d3 d8— 80, the matrix T = O‘ - dO/ds is antisymmetric; we describe it in terms of its three free elements via 0 A (10 A -73 7”2 0t ' — = T 2 7'3 0 —7'1 . (2.9) d5 —T2 T1 0 14 The three elements we group into the vector ?, which has the form ‘11 TI T2 7'3 We observe that for any vector (1', we then have the relation T . The components of the vector ”F, computed as dzfxi and hence the elements of the matrix T, can be T1 = 8 ' Z —81 ' ’ y d3 d3 _, dé‘, 4 d5, 72 = e, - — = —e - (13 y ds T3 2 ex- = —e, - —. (2.10) (13 ds These relationships give some practical meaning to the components of the vector 7": Apparently, 7'1 describes the current rate of rotation of the Dreibein around the reference curve 13(8); 72 describes the current amount curvature of 13(9) in the plane spanned by 5,, and 6,; and 73 similarly describes the curvature of R(s) in the plane spanned by 6E and 6,. In more mathematical terms, because of we have dé's d3 dé’x d3 dé'y d3 as successive multiplication with é}, —-o dex dcy 83; ° —d_9_ : , 6y ' d3 = 0, (2.11) = —T3€s + Tley = T263 — 7161, (2.12) é}, and 6,, and comparison with (2.10) reveals. 2.2 Differential Operators in Nonplanar Curvilin- ear Coordinates As the first step in the transformation of the Maxwell’s equations as well the equations of motion to the curvilinear coordinates, it is necessary to study the form of common 15 differential operators in the new coordinates. From (2.7), which has the form 3 3 r: :3 xkét = Z {R - a. + $01.2 + 901:3} a, k=l k=l we see that the Cartesian components of F are $1,, = if Bk +$Ok2 +y0k3, for [(2 21,2,3. Hence the partial derivatives of 5r,c with respect to s, :1: and y are % 68 9121: am where (2.1) and (2.7) have been used (3‘ dI-f(s) (13 0k?) : (8($19$Qa$3) 06,22,111) d012 /'—\ COO OOH ’° d3 y d3 6:1: and 55k "—"- Ok3, =Ok1+$ d0k2 + dOk3 (13 y d3 . Thus, the Jacobian matrix C is (81131 (9.732 (95133 \ 33 593 593 011 021 031 ) = :1 a: :3 = ( 012 022 032 31:1 3:2 2:3 013 023 033 \ By By 6‘31 / (1013 $61022 + 61023 $d032 + d033 y d3 ds 7’ ds as y ds 0 0 0 0 9 d0’ . 0 f y d0“t ~ ~ 0 S 0 0 0 3 y A . 1—ngr+7'2y —7'1y 711: A 0 -T’ -0‘= 0 1 0 -0‘. 0 0 0 1 16 It is convenient to denote the first part of the Jacobian matrix C by A, i.e. . 1— T311: + my —T1y 71:1: A = 0 1 0 0 0 1 then the Jacobian matrix can be written as A A 0:30“. The inverse matrix of A is found easily; we obtain I le — T113 A" _ 1—73x+72y 1—T3:r+72y 1—T3:r+7'2y _ 0 1 0 0 0 1 For later convenience, it is advantageous to introduce the abbreviation 05:1—7'3113‘l‘7'2’y. (2.14) We note that for :r and y sufficiently close to zero, a does not vanish and is positive. Hence besides the restriction for the motion to be inside a tube of radius r1 imposed by the need for uniqueness of the transformation to curvilinear coordinates in (2.2), there is another condition; defining I . I 1 7‘2 2 — m81n( , >, then if we restrict 11:, y to satisfy |:z:|, [y] < T2, the quantity a never vanishes. In this case, for the inverse matrix of the Jacobian matrix C we have — T151} .1. 23/. —l_ ‘-1_‘ a a C —O‘A —O’ 0 1 0 0 (2.15) (2.16) N ow all the necessary preparations are made for the calculation of partial differ- ential operators such as gradient, divergence, curl and Laplacian in the curvilinear coordinate system. 17 2.2.1 Gradient Let f be a scalar function, expressed either in the Cartesian coordinates (2:1, 3:2, 2:3), or the curvilinear coordinates (s, 2:, y). From the chain rule, we have _ _ — _ 8 833 5’5. 3’3. 5’3: 5" (T)— f: a: a: :r . 2317— f : COV f, 9: .5133; 3532:. 25313.3 62 \ 8y / \ 3y 3y 5y / \ 55; ) where V“ is the Cartesian differential operator vector. Multiplying with O”, we obtain the expression of this vector in terms of partial derivatives with respect to the particle optical coordinates as a y a {1 a a a \ { 5:613: { 36—3 \ a (5g +lea—x' — T111359) €707: 3;— f=C‘"- @— f=0' 5‘?— f, _._2 i 3 where (2.16) was used. We now define the vector differential operator V0 as ( .1. 2+ 2. a ) V10 Vs a as 713/82: — T1235; v67: vg? f = v, f = 3 f. (2.17) V? V. 0.9” \ 5 ) Then we have 22; = 0 - 60; and 60f = 0‘ - vetf, (2.18) or, for later use, in components, 3 3 V? f = 20w? f and VS f = Z ORV,“ f. (2.19) (:1 [=1 18 Let us consider a vector function A; we express it in both Cartesian and curvilinear coordinates: 21' = A161 + A252 + A363 = 21,3; + Ag; + Aye, We denote the component vectors in Cartesian and curvilinear coordinates with A“ and AC, respectively, and have A1 - .410 A, A“ = A2 , AC _ A20 = A, A3 A5,? A, Then, because of (2.3) A A = A353 + A252 + Aggy = As(0€1)+ Ax(O€2) + Ay(O€3) = 0' (A351 + ’4ng + Asia), and so we have A“ = O - AC as well as AC 2 O’ - Aft . (2.20) As a first step, we now want to determine the form of the gradient operator in curvilinear coordinates. In the Cartesian system, the gradient operation is f 21:. \ 6:1: gradctf zed; = 9— : it. +fia +39; €911? 8231 1 8.732 ,2 61133 3' lag) As in the situation with the more common coordinate systems, the gradient operator in the curvilinear system should determine the Cartesian gradient of a function, and then express it in terms of curvilinear coordinates; so we must have gradCf 2 Of - gradC‘f. 19 We find that V0 f defined in (2.18) satisfies this demand; so the gradient Operation in the curvilinear system is gradszVsz if that is 1 6 8 gradCf = [a (51+T1y5— — TliE— a ) f] 63 ‘f‘ iii—.51; + 91:6,).(221) 2.2.2 Divergence The divergence in the curvilinear system is calculated as follows. In the Cartesian system, . " 8A1 8A2 8143 d A = . 1v 8:131 + 62:2 + 8:133 Our goal is to express it in terms of the curvilinear system. For this purpose, we apply (2.19) to the three components Act :2 gm Ohm/IS, div/I = (66‘ - A“) = 2 Vi: ~ ACt = Z Z 0W? (Z OMAS.) k = Z: OkzokaCAC + Zi0;(( V, COkm)AC. kl,m l,m Since O = O(s), we have 1 dOkm- d3 v? 01...: 3,,v 0,.m— _ (2,— (2.22) Using this relationship, we obtain div/I = Z Z amokmvaC + Z Z 0,,1613—1 ad—Q—gf'" Am 20 = :[O OhmVFAS, + Em 6,,—[O‘—],mAS, l,m . 1 = Z (SImVlCAgl + Z aTsmA-C: Z VCAC+—T'1:3A + TQAy ) l,m m 1 a a a 0A, 0A,, 1 — a (5; + T1y5; — T111553?) A3 + 8513 + a + 8(‘7'3/41: + T2Ay)° Thus, the divergence expressed in the curvilinear system is obtained as . ~ 1 8 8 8 8 8 leA— ;{(OS+ lea—‘L‘—Tl$5—)As 3+0$ (OA$)+‘a—y(aAy)}. (2.23) 2.2.3 Curl The derivation of the curl in the curvilinear coordinates is a little more involved. In the Cartesian system, curlC‘A .—. vctxfict=wctx AC‘]1€1+[VC’XAC‘]2€2+[VC’XAC’]3€3 (6:19.212) 61132 (9133 _ fiflfi — 6133 01131 @352 \ am, 822 / The curl in the curvilinear coordinates, which we denote by curlCA and which has the components curlCA 3 curlCA = curlCA x [curlCA y has to satisfy the condition curlC‘A = O - curlCAo. (2.24) 21 First, let us express each component of curlC‘A in terms of the curvilinear system, again using the transformation rules for derivative components (2.19). We obtain [CHEVY],c —.= [6“ x A], = Z eijkat/igt = 2 ZeijkoiNfloijg) 1,] i,j l,m : Z Z EijkOilemvchg + Z Z 6ijkOil(VICOJ'TH)I4g l,m i,j l,m i,j : Z (Z élmnokn) VCAC + Z Z: 6ijk0i1615_ 1 adggm Am l,m l,m i,j = :Okn (IE: ElmnVICAC)1+—' a: (Zj CijkOis— —%:2) AC where (2.8) and (2.22) are used to obtain the third line. Making use of the fact that dO-m déC ZjGijkOisT —_-€[.3 X —dS—-]k as well as the relationships in (2.12) which entail -o dé’s deg, de e, x —— — —7'3é'y + T263, 6, x — = —T1€I, and e, x —y- = —7'1é'y, 3 d3 d3 we obtain .. .. .. 1 [curlC‘A],c = Z 0,.,,[VC x A0], + 3{(T30"3 + 720mm, — now/1,. — now/1,} _. —+ 1 = Zoknwc x AOL, + a{0k2(7-2As - 71/12:) + 01:3(73As - TlAy)}- So 6 X AC 0 _. - .. s 1 1 curlCtA = O- V0 x AC 2 + — ( 72A. — 71A: ) [60 x [1'0] 0 73A, — 71A,, 3 22 Now transforming the curl vector to curvilinear coordinates according to (2.24), we have 60 X EC 1 0 curlC/i = O‘ - curlazéio 2: {70 x KC 2 + —1- 72A, —- 71A,; 60 X EC 0 T3143 — TlAy 3 vay — VyA, 1 = VyAs — VsAy + 3(72As — 71A,) 1 V5141; ‘- V3143 + 3(7'3/43 — TlAy) So altogether, expressed in terms of partial derivatives with respect to the curvilinear coordinates, we have r (La—A— 1 8m 6y curlC/f 2: 66:3 — El!- (56; + 71.71;; — fizz-6%) A3, + éths — 71/12:) \ 31(5); + figs—m — mg?!) A, — % + 373A, — n.4,) ) 1 954.6; 1 z i £11.43) — 71A,, - (g + T1959; — 7,15%)A, . (2.25) ( —(%(aAs) + (g; + 11.11% — 7122593) An: - TlAy ) 2.2.4 Laplacian The Laplacian operator in the Cartesian system is Actf : (fiat ofict)f : (ficty '6ctf : __2_ + __2_ + ___E In terms of the curvilinear system, utilizing (2.19), we have ACf = (fiery . fictf = Z: OHVflOkagilf k l,m 23 = Z Z 01101....V,C(V,€,f) + Z Z 0.,(Vf okmeg‘, f) l,m k k l,m : ZélmVflVCf) f)+:: loksd W m :. Vs(st)+ (V ——) (V f)+ (V y +%)(Vyf), where (2.4) and (2.22) are used from the second to the third line. Thus, the Laplacian operator in the curvilinear system, expressed in partials of curvilinear coordinates, has the form c _ l 3 2._ 2 31’ 3f (if A f “ a(0.s+le0$187-1x8y_aa_s+ly8xfixay 18 aaf 8f 2.2.5 Velocity Vector in Nonplanar Curvilinear Coordinates The final differential quantity we want to express in terms of curvilinear coordinates is the velocity vector 17. It is expressed as ”(7 = @151 + U262 + U363 1‘ U353 + Uzgx + Uygy, and similar to before, we define and we have 17“ = O - 170. To determine the velocity expressed in curvilinear coordi- nates, we differentiate the position vector F with respect to time t; from (2.13), we have {R‘ - a, + 3301.2 + yOk3}é'k 24 : 2 {011.18 + 01,3211} + 05;?) + S k=l s . 0 s . 0 . d3 . d3 31 y y y A S A 0 A S (1—T3$+T2y) O- i: +sT- a: 20- it—s'rly , 3.] 31 31+ $7133 where (2.1) is used from the first line to the second line. Comparing with the previous d0k2$ + é61(01c3 .. (13 (is y 6" ll Q. equation, we see that the velocity expressed in terms of curvilinear coordinates is given 1), 5': ~ (1 — 733: + 723/) so {I _ —_ o . — . . v — v1 —— :1: — srly — :1: — srly , (2.27) vy (Hem: y'+.é7'1:r where a = 1 — 732: + 723; as (2.14). For future reference, we note that because of the by orthonormality of O, we also have the relationships ”2:27am. :1 t = 17C - 270 (2.28) V 3K ,l-j'ct . Act 210. (2.29) 2.3 Dynamics in Nonplanar Curvilinear Coordi- nates 2.3.1 Electromagnetic Fields and Lorentz Force In the Cartesian system, the electric field E and the magnetic field B are expressed in terms of the scalar potential (I) and the vector potential 4 as aid — _6ctq) _ 81‘?“ act _____ _ ct _ E grad (I) 8t at B“ = curlCtA = V“ x A“. 25 Our goal here is to express these fields and the Lorentz force in terms of the curvilinear coordinates. We need to find EC and BC and their relationships to the potentials such that E“ : O - EC, and g“ = O - EC, # E1 - E, _. B1 g B, Ect 2: E2 , EC 2 E, , and Bct = 32 , BC = B, . E3 E, 33 B, Using the differential operators from the last section, we have where 615d “ct .. .4 .. 0A ) = 42m _ 0t.— at -o EC ___ O‘tE’ct : 0‘1 (_§ctq) _ 8t ts. (——) ) EC = -33 _ 5A3 . (2.30) 331% 3‘31, \ "('53 _ 7’97 / The magnetic field Eccan be determined in a straightforward way from the transfor- mation rule for the curl (2.25), and we have VxAy — VyAx 1 BC 2 curlC/I : VyAs — VsAy + 3(72143 — Tle) 1 VsAx — VxAs + 3(T3A3 — TlAy) 26 ( (9A,, _ BA, 8m 8y \ 8A., 1 (‘3 6 6 1 1 8 8 6 6A, 1 \ a (a; “We: ”1’55le ‘ a. + EM ”A” 1 In the Cartesian system, the Lorentz force per unit charge, denoted as f, is given fct : Ect + ,I-J'ct X Bet, and we want to find f0 such that f.“ = O . ffC, where we write the components as 4 f1 f3 f“: 12 , f0: 1. . f3 fy Because of the orthonormality of O, we first observe A -o f0 : 0tf'ct : Ot(E'ct + ,D'ct X éct) : Ot(OE'C + (0170) X (080)) _. _, Es + vay — '0sz = EC+17CXBC= Ex-l—vyBs—vsBy , Ey + vsBx - ’Ust where the invariance of the cross product under 80(3) transformations (2.6) was used. Expressed in terms of curvilinear coordinates, we have explicitly for the components of f0. f3 = E, + vay — vyBx 6A, 6t I —v_,, {VyAs - VsAy + 8(7'2143 - 71/42)} 2 —V, — I + ”U; {Vs/ix - VxAs + 3(7'3As — 711411)} dA3+,6A3+,8A3+,%+ VA—v% dt 883 $851: gay '03; s I 38:1: = —V,<1> — 27 v1 8A +;(T3As — TlAy)_ ”y 83/ + vyvs Ay a —y.9(T2A _ Tle ) (1A, : — dt — VS(<§C=— 1 d d d _. +5 {Asa—N”— T-zy)+Ax dt —y(7'19)+A dt—( 7137)}‘esa (2-32) knowing the Lorentz force is the first step towards determining Newton’s equations of a charged particle in the curvilinear system. The momentum of the particle [9' is expressed as 15': 10151 + P262 + psé‘a = psé‘s + 1015; + pyé'y, and similar to before, we define [’1 -C 1’3 15‘“ = 192 , p = pm , 173 pg and have as before that 15” = O - 13C. In the Cartesian system, a particle with the charge 6 obeys Newton’s equation dfitt "ct Now we merely have to express the momentum derivatives in terms of curvilinear coordinates: dfi“_ da{;_Ade dOA€_A dp tdOfiC 71? * az<0pl—O"d7+fip ‘0 d7 + Odf A (115C A 6115C —O(d—t-+$'Tp): 0(dt+STXp), and then (2.33) can be written as O - (de/dt + S? x 13C) 2 e O - F}, or directly df , ~ 0% + 5'? x 15" = 6 f0. (2.34) 29 Explicitly we then have d p8 Tl p3 E3 ”3 BS 8; pm ' 72 X p2: = 6 ° Ex + ”I X Bit pg T3 py Ey v31 By P K — + T 6 — T 23—0— \ d A. a (7313/3331 a z e —a-t' ( Ax “‘ 3 ((D — [DC . EC) A, 3g - \ 3y l 1 d d d _, +3 {A 212”” — 722/) + A d —(T1y) + A, d—,(-m)} «2] . (2.35) 2.3.2 The Lagrangian and Lagrange’s Equations Now we are ready to develop Lagrangian and Hamiltonian methods in curvilinear coordinates. Following the transformation properties of Lagrangians, it is conceptu- ally directly possible, albeit practically somewhat involved, to obtain the Lagrangian in curvilinear coordinates. To this end, we merely have to take the Lagrangian of a charged particle in an electromagnetic field in the Cartesian system 62/ 2 I C I v a L(:r1,:r2,at3;a:1,:r2,:r3;t)=—m 21——2—e<1>+effC-A6t c and express all Cartesian quantities in terms of the curvilinear quantities. In this re- spect, it is very convenient that the scalar product of the velocity with itself and with A is just the same in the Cartesian and curvilinear systems, according to (2.28) and (2.29). So the Lagrangian in the curvilinear system is obtained completely straight- forwardly as 202 L(s,:c,y; 33m); t) = —mc2 1 — UC—,2 — 6(1) + 62) ~AC, (2.36) where 2702 2 of + 112+ v5, and 17C - AC 2 USA, + vax + vyAy. 30 Here (I) and AC are dependent on the position, i.e. {3,x,y} and the time t. The quantities O, T, and hence T1, 72, 73 used below are dependent on s. The derivatives of 213,253,223, with respect to s,:z:,y,é,:i:,g] are useful in order to determine the explicit form of Lagrange’s equations. (92), 8U, 3223 as ‘a’ 82': ‘0’ 0;) ‘0’ Box __ _T 6’01. _ 1 av, _ 0 598 _ I?!) 88.2% "‘ a Bay "" ) .1)! = _v_y_ = 12 = as “‘5’ 6i: 0’ 6g 1’ (2.37) 31,8 _ “—919.34. 172 ) 8218 _ , (9v, _, (93 _ d3 d3 ’ 8:1: _ ST?” 83/ _ 8T2, 82),, _ —,§ fit 61),, _ 0 Box _ . 03 — (18 J’ 030 — ’ 6y — 871’ Q‘Ly _ 5. 1:1,, % _ 3, fl _ as ‘ d3 ’ 0:1: _ 1’ ay " The Lagrange equation for a: is derived as follows. Using the derivatives of 113,113,223, in (2.37), we have 8122 8'08 (92),, Buy a; — ”*5 + 2%; + 2W5 — ”’1 8(27C . A‘C) 8U 82) 0v —— : AS—S 1413—1: ___y_ = 1: 02': 0:13 + as: “131 8i: A 829 82), 82):, av , , 63—23 = 221351:— + 2vI—8—x— + 2vya—ai’ = —2v33'r3 + 22131371 : —2.§ (73v, — fivy) = —2.§ [7" X 17cl2 So altogether we have L , 8 = ——T——vx + 6A1. 2 pm + 6A2, (2.38) 811°: /1 _ ”02/62 where p“ = mod/(n — 122/02 and correspondingly 15C 2 mic/M1 — v2/62 was used. 31 We also have 9% = —————.1__Tn;mé[?x6‘c]2—ea%(—17C-AC) — _s[7x,5012_e§5(q>_50 A”) Thus, the Lagrange equation for a: is %—+S[?Xp€]2=e {—12, — %(¢—17C-AC)]. (2.39) The Lagrange equation for y is derived in the same way, and it is %+S[?xp0]3=e[—%—6Ey(@—17C-AC)J. (2.40) It is a little more complicated to derive the Lagrange equation for 3. Using the derivatives of 2),, 22$, 2),, in (2.37), we obtain ,2 85:— : 222,0 — 2v$71y+2vyrlzr «c.1427 6(1) 65‘ ) = Asa — Aley'i'AyTlIE 8'02 dT3 dTQ drl dT1 “=29. __ _ _2x.___ ._____ as us< dsx+dsy) vsdsy+2vy3d8$ d" d" : —2s$ [lxfic] —2.éy [1x66], d3 2 ds 3 and so 8L m —-—,- 2 ~-—————————- . (vsa —— only + vyrlx) + e (Asa -- Axrly + A 713:) as /1__ 212/02 3’ = (1),. + 8A5) a — (Pa: + eAI) mg + (py + eAy) 7'12: (2.41) as well as 6L m d? _C d? _,C 3 _,C ~C — = ————- —— — — — — — (I>— A as (1-v2/c2 { SIB [dsxv L 3y [dsxv]3} €83( ) 32 Thus, the Lagrange equation for s is d( a 72+ Tx)+s'a: de'C +’ d?x'C _ s -’ x 1 _ 3 _ dtp P 1.1 le (18 P 2 31 d3 P 3 -' e —i(AC¥*-AT +AT:1:)—-£9—(—17C AC) (242) — dt 8 1: 131 y 1 ('38 . . The left hand side is modified as follows dps (1px dp . . . . 0 dt -T1’J—‘ dt +T117—df+Ps(- T3$+T27J)‘PxT1?J+PyT1£E + (— —sd—T3$+s 1T3)— .9@ + éfix+éx (fix '0 +' de'C Ps (13 (13 y Pm dsy Py ds d3 P 2 33/ d3 PC 3 dps dpx dp adt-Tiydt+T1$-d—:—P1TXP12— leXPCls ll d’ s , _, d 1: . _. d . _. oz (d1: + 8[T X pC]1>— 71y (32— + s[7' x fich) +713: (7?; + 8[7' x 150k) —"U3[7_'. X 15611 - 14:1:[77 X lack — 'Uy[7:' X fich d 3 , _, d d . _. a (—§£—+3[T xpC]1)—le (2+sh xfich) +T1x (%+S[T xpc]3), where (2.27) is used from the second step to the third step, and U3[?XI3C]1+Ur[?XfiC]2+1)y[?X]5C]3=2_)C°(7?X15C)=0 is used in the last step. So, the Lagrange equation for s simplifies to alps . a dpx d . .. 01(71t—+31TXPC11) —T1?J(dt + l Xlflz) +T1$ (%+SlT X15013) [ dA, dAx dAy = e —a d dt+nwa‘ “37+Asdt( Tax - 722/) d d a _,C C +Ad0w%h%a(nfl-5j@ A4 The equations for :r and y, (2.39) and (2.40), can be used to simplify the above equation. Doing this, we obtain dps .1 {3 _ dA, 3 a a 4C 40 a(dt +S[TXp]1) — e[ adt (83+T1y05r leay)(<1> v A) 33 d d d +Asa(T3x-sz) +A d —+A.d—,(— 71.22)] and with the above requirement that :1: and y are small enough such that a = 1 - T3113 + 72y > O, the equation can be written as (M, 1 8 0 (9 dt _ —_Co‘-‘C dt (163+ leazz: —T1:1:ay)(<1> v A) d —”- +S[T'Xp]1= 6(— 1 d d d +—{A d—,(m— sz>+A d(ny)+Aydt(-T1x)}]- Thus the set of three Lagrange equations can be summarized as below; it appar- ently agrees with Newton’s equations in curvilinear coordinates (2.35). d p3 . Tl p5 (7t pa: ' 72 X p2: pg 7'3 pg d 6A, 6 8 88 63/ _. : _a— 6A; —"'" a— (Q—IUCOAC) t eAy a Bar \ “BE / e d d d _, +3 {Asd—t(T3:c— Tg’y) + A xtd —(T1y)+ Ay ytd —(— —T1:1:)} 6,. (2.43) 2.3.3 The Hamiltonian and Hamilton’s Equations To obtain the Hamiltonian now is also conceptually standard fare, although practi- cally it gets somewhat involved. We adopt the curvilinear coordinates {s,a:,y} as generalized coordinates, and we denote the corresponding generalized momentum by 136 == (Pf, Pf, PyG) The generalized momentum is obtained via the partials of L with respect to the generalized velocities; using (2.38) and (2.41), we obtain 31. P? = a? w(ps+eA)a—(px+eA)ny+(py+eAy-')m: 34 g P? = . = pa + 6A.. 0:1: 8L Pf = a?) = p, + eAy. (2.44) It is worthwhile to express the mechanical momentum fiech, namely 13C, in terms of the generalized momentum 130 2 (PE, Pf, PyG). By combining the above expression (2.44), we have PE 2 (p, + 63/13) a - Pley + Plex, and so 103 + 6.43 = (P? + any — P5712“), QIH and altogether I — (P? + Pley — Plex) — 6A, a filech = [90 : Pf _ 8A1. ' (245) Pyc’ — eAy Squaring 15‘2" = 7772.270 2 172.270 / \/1 — (170)2/02 and re-organizing yields (27C)2 _ 62(pc)2 ‘ (:50)? + m2c2 and because 170 and pf are parallel we even have _C c 270 = p (2.46) \/(]5C)2 + m2c2. We also observe that 1 (15")2 m6 T : \/1_(IUC)2/02 : J1_ (if)? + m2c2 : \/(25C)2 + "L262. (247) The Hamiltonian in the curvilinear system H is defined from the Lagrangian L (2.36) and the generalized momentum 160 (2.44) via the Legendre transformation H = s'Pf+:i:Pf+g)Pf—L = sPSG+ziIPf+fo+mc2 1——+e—efchAC, c 35 and the subsequent expression in terms of only 3, x, y, PG, PG, PyG and t, if this 18 possible. Using (2.45), (2.46) and (2.47), we have from (2.27) that s 1 1 1 . 3 = 1 = —— {— (P? + 1926le — PyGTlx) — 6A8} , 02 r1170 01 1 1 i‘ = vx-l—Slez— [PG— eA + —yT1{—(PSG+Pley—Plex)-eAs}] my a a 1 1 0+ 2 2 2 G 2 G 2 y = ’Uy—STIIL‘ 1 1 : mgr—2{—T1$PSG_T12$fo+ (012+sz172) Pya+€T1$CYAs “802 Ag}: where we used the abbreviation 2y from (2.47), which is in terms of the generalized coordinates and the generalized momenta 1 C 1;)._ PG+PGT —PG7'.’E—C¥6A32 . \/< . . 22 02” 2 l + (P? — 6.41..)2 + (12.2 — eAy>2 + m2c2 (2.48) We also have 50.2? ::.JL7G A? m7 1 1 : __ [{_ (1380 + P107131 -— PfT1$)- 6143} As my a +(PxG - 6’41le + (Pg/G _ eAylAy]a and in particular it proved possible to invert the relationships between generalized velocities and generalized momenta. Hence the Hamiltonian H can be expressed in curvilinear coordinates, and it is given by 1 11 = -——[$ (Pdepffiy— Pan)PWG4. m7 T13} —PEPG a2 a y + {1 +(71y)2}(Pf)2-L1$P8GPG+{1+(712)2}(PG)2 36 2, , _27’1.I,‘y 71:1: P—Ge,A I 1 . PGPyG — an’eAs — 71—fo eA— PGe A +— 02 —PfeAy — :7 (P? + any — Pu 0713:) 8.43 + 62/13 —(P.f' — eA x)eAx — (PyG — eAy)eAy +1m2c2] + e +(7[—(:—2( (P? + any - PyGTISE)2 — 2— :3(PG + PGle— P0712“) 814 + 62/42 —eA)2 +22(PyG—eAy)2+mc]+e PIE PG 1 G G G 2 = —:y- +({a(Ps +Px7‘1y—Pynz)—6As} G m +(P — eAx)2 + (PyG —— eAy)2 + 771202] + 6(1) 1 '2 —-(mc7)2 + 6(1) 2 mc27 + ed). m7 Explicitly, the Hamiltonian in curvilinear coordinates is H \/(PSG + Pf'rly — PyGT1513 — cue/13V - c + (Pf -— 6A$)2 + (PyG — (3Ay)2 + mzc2 a2 +e, (2.49) where again a = 1 — 73:13 + my. Thus we derive Hamilton’s equations as follows. H 11 1 5‘ = a — -—{ (PSG+PGT1y—P3/GT1$)—€AM} BPSG— m'ya a . _ 8H _ 1 Try 1 G G G G :1: — 8P5 _Hv_l—{E(Ps +Px le—Py T1$)—8A3}+Px —eAx], . 1 G G G C y :: W:m[—?{E(Ps +Pley—PyT1$)—€A3}+Py —€Ay], PE = —%—:I — 7721; [ {a (PG+PIGTy PGT1$)—6A } {.<— _..._)(.. -3. +$ (Pf—(33y Pfflfix) 685:3} +6(PI 84$)ag? + €(PyG eAy)%£S€] — leg—i), (2.50) 37 - 0H 1 1 . Pf : ——0—; = ;7- [— {a (P86 + Pley — P3107133) — 6A3} . & ~{E(PSC’+Pley—Pf7-l$)—:1-PyG"ea } a 02 0:1: +e(Pf — magi + 6(1):; — 61409;] - 62%, (2.51) P316 2 _00—:/I = iii—r; [—- {i (P? + any — PyGT11E)- 6A3} .{_(_:25 (p80 + PIGTIy _ prlx) + ER? _ 863:3) +e(P$G — 63/1395?ij + e(PyG — eAy)%/:—y] — 6%, (2.52) where the abbreviation (2.48) is used. To verify the derivations, we check Hamilton’s equations to agree with previous results. It is shown easily that the first three equations agree with (2.27). The last three equations are shown to agree with Lagrange’s equations (2.39), (2.40) and (2.42). We have from (2.51) - 1 T G _ 3 G PG G G Pa: — ——’U3{—a (P3 + :Ble—Py T1$)—T1Py} 06 6A3 a4. 8Ay 8(1) -—-+e'v —+ev ——e— x33: 3’81: 6:1: ' .éTT 7’17 =-—4fn4fifi+m(i+0fi C! a 0CD 2) 63/43 v BA;r 8A,, 8— _ 213— — v — ‘ ('93: 3:1: 8:1: y 8:1: Expressing the equation in terms of the mechanical momentum 130 rather than the generalized momentum I30 according to (2.44) and using (2.37), we have , dA éT: Pm + e d: = ——a—3 (193 + eAs)oz — (pm + eA$)T1y + (pg + eAy)T1;z:} éflfiy W (1r)I + eAx) + 3T1 (“a—r +1) (pg + 64.) 38 (9 (9113 311$ 32) = —.é(Tp —Tp)—e—g—(@—fiC-§C) 3 8 1 y 011: , which is in agreement with the first Lagrange equation (2.39). The Hamilton equation for y is similarly modified from (2.52) . $72 . , 72y STl’rgsc p; z .42; + (_ -1)pg _ __pyc a a a (9y say 3333/ gay ' In terms of the mechanical momentum 150, we have , (1A , 8 .. py + e—dt—y = _S(Tlpx — T2193) - 68—31((1) - 170 - AC), and it again agrees with the second Lagrange equation (2.40). Similarly, the Hamilton equation for s is modified from (2.50) -. ' d d ' d (1 pg z 2(_ Ta +2)psc+{§_(__t2x+_sz)le_s-fiy}pf a 8 S S —:r oz (13 (13 y 3 dm (1T2 ,dTl G +{ a( dsx+d8y)T1x+sdsx}Py —-e @_U%—U%—U%— 8.9 583 x03 3’63 ' In terms of the mechanical momentum 13C, it takes the form d — {(1). + eAs)a — (p. + 614$)le + (py + eAy)Tlx} dt —§—gT-3 +fl {(+A)a—()+A)T +( +A)} "' a (13$ dsy ps 6 3 I11: 8 11y pg 8 317,117 ' d. d d + 3 __2. +—T—2y w—é—Ey (px+e4.) a (13 (19 ds 39 s (1T3! dTQ' .6171 +{ ;( (18 .1: + (13 y) T12: + 38—8-17} (py + eAy) 5E _ Us??? — “fa—s— ‘- v, as _e 8.4, 8A,. (9A,) and a little reorganization leads to d( + )+' dfx‘c — .sa— IT Ta: 3:1: —- d 3 _. = e {—“(E (Asa — Aley + AyTlfL‘) — 5;( — 17C - AC)] , which agrees with the third Lagrange equation (2.42), as it should. 2.3.4 Arclength as Independent Variable for the Hamiltonian As the last step, we perform a change of the independent variable from the time t to the space coordinate s. For such an interchange, there is a surprisingly simple procedure which merely requires viewing t as a new position variable, —H as the associated momentum, and —PsG as the new Hamiltonian, and expressing it in terms of the new variables, if this is possible. Then the equations are dzr a(—P,G) @_a(—pf) dt_8(—PG) _ 8 E‘ 639’ ds— an’ EE"a(—H)’ de: a(—PG) dPyG__8(—PG) d(—H) a(—PG) _ S 8 d3 8:1: ’ d3 — 8y ’ d3 — 8t To begin, let us try to express —PSG in terms of 13,23, y, —H P0 Pyc. From (2.49) we 7 x 3 obtain that 2 {5; (P80 + Pf'rl'y — PyGTlcr) — 6A5} + (Pf — eA$)2 + (10;; — eAy)2 + m2c2 1 = —(H — 6(1))2, C2 40 SO . 2 (P? + Pley — Plex - 08A,) 1 = a2 {—2(H — 6(1))2 — (Pf — 6.43)2 — (PyG — eAy)2 —— m2c2}. 0 Considering the case that [f = 0 and :z: and y are small, we demand p5 should be positive (and stay that way throughout); we also remind ourselves that a > 0, and hence the choice of sign is done such that P36 = —PIGT1y + PyGTlx + 06A, 1 +0 E(H — 6(1))2 — (Pf — eA$)2 —- (PyG — 63/11,)? — 777.202. Thus, —PSG and hence the new Hamiltonian H 3 is obtained as H3 = —PSG = Pley — PyGTlrr — ore/13 1 —a c_2(H — 8(1))2 — (Pf — 6A3)? — (PyG — eAy)2 — m2c2. Here, for later convenience, note 1 —2(H — 8(1))2 — (Pf — eAx)2 — (PyG — eAy)2 — m2c2 C 1 z 2; (PE + pgfiy _ p57...) _ 6A, = 1),. (253) Then, the equations of motion are 513 _ a<—P.G> _ T + aw? — eA.) d3 ’ 3P5 _ ‘y 1 G , §(H - 8(1))2 — (Pm — eAx)2 — (PI/G — (3143,)2 — m2c2 (2.54) __ G PG _ A g :: 8(0PI2 ) = —T1$+ 1 a( y e y) 1 s y ;(H — 6(1))2 — (Pf — eAx)2 — (PyG — eAy)2 — mzc2 dPIG 6(-—PSG) G 3 d8 — —- 8:1: —— Py 7'1 —- 673.4, + (16 8:1: —T3 37(H — 6(1))2 — (Pf - eAx)2 — (PyG — eAy)2 — mzc2 -15(H — «ea-83 — (Pf — (24)an — (Pf — wail-”- —08 C 1 6.73 01E 61E ’ (257) \[c—2(H — 8(1))2 — (Pf — €A$)2 — (PyG — eAy)2 — m2c2 dPG 8(—PG) . (9A 3! __ _ s : _ C, __3_ d3 — 6y Pm T1 + BTQAS + (16 By +T2\/—C-1—2-(H —— 6(1))2 — (Pf —— eAx)2 — (PyG - eAy)2 -— m202 i(H — 8406;? — (PG — 8A1)6A$ — (PG -— eA )a_A£ C2 0y “7 8y y y 8y —ae 1 , (2.58) 25(H — 6(1))2 - (PIG — (214$)2 — (PyG - eAy)2 — m262 d(—H) _ _a(_PsG) d3 _ 8t 1 6(1) 0A 8A ——H———PG—— A, x— G— —1 :ea 8548—0“ e)6t (x 6 )6t (Py eAy)0t t 1 \/c_2(H — 6(1))2 — (Pf — eA$)2 — (PyG - (3A,)? — m202 (2.59) For the sake of convenience and checking purposes, we replace Pf, P3? and H by [BC using (2.44) and (2.49) and with the help of (2.46) and (2.53). Then we have from 42 (2.54), (2.55) and (2.56) d3? pa, E — 713! + 01;): (2.60) d—y = —m: + 01& (2.61) (18 p, dt 1 (Pf)2 + m262 — = a— . d3 p, c And we have from (2.57) (1121” dA:1C 3143 d8 8 d8 : (py + €Ay)7'1 — €7'3As + area—z — T3ps ‘0‘“ c 55 We? ”ya—.5 e {\/(I3C)2 + m2c2 :39 8A1; BAy} p3 ’ and organizing the expression using (2.37), (2.27) and (2.46) we find dpl‘ _. ‘C _ dAI 1 a "C "C ds +[TXp L—e[— d3 —§5:;(—v -A ) . (2.62) In a similar way, we obtain from (2.58) d])y _. ‘C' dAy 1 a "C "'C — = ———‘ — —— (I) - . A . ds +[sz2 ]y e[ (13 say( '0 ) , (263) and from (2.59) This concludes the derivations of dynamics in curvilinear coordinates. In partic- ular, we have succeeded to derive the equations of motion of a particle moving in an electromagnetic field in curvilinear coordinates, with the arclength s as the indepen- dent variable. Moreover, we know that these equations of motion are Hamiltonian in nature, which has important consequences for theoretical studies. 43 2.4 Planar Motion As an application of the concepts just derived, let us consider a particularly important special case, namely the situation in which the reference curve stays in the $1412 plane. This so-called planar curvilinear system occurs frequently in practice, in particular if the reference curve is an actual orbit and the fields governing the motion have a symmetry around the horizontal plane. The basis vectors in this 2D curvilinear system can be expressed by the Cartesian basis vectors via —o —0 6y 2 63 63 = cos 061 — sin 962 e; = sin 661 + cos 662, where 6 depends on the arclength s; denoting its derivative by h, i.e. : d6(s) h = h(s) ds From (2.7), all the elements of the matrix 0 are determined as . cos 0 sin 6 0 O = — sin 6 cos 6 0 0 0 1 So, the antisymmetric matrix T of (2.9) has the form - .. d0 c030 —sin6 0 —sin0-h cosO-h 0 T = O‘-—= sin6 c050 0 - —cos€-h —sin0-h 0 d5 0 0 1 0 0 0 O h 0 = —h 0 0 , 0 0 0 44 thus the elements of T and hence 7" are given as finally, we have azl—T3x+'r2y=1+ha:. The partial differential operators in this 2D curvilinear system are, from (2.17), (2.21), (2.23), (2.25) and (2.26) { 1 0 Vs 1+}533g 60f = (V.)i= @— f v, at \ 'a‘y / 1 8f_, (9fa 8f_, C _ _ _ __ grad f — 1+hzr 8383+8xe$+0yey , _. _ 1 6A, 1 8 6A,, dWA _ 1+ha: 83 +1+hx6;{(1+hx)Ax}+W ( 25% ) 62: By curlC/I = 0A,_ 1 % 0.48?! 1+h$683 1 a, 1 )1+hx-6?—_1+h25§{(1+hx)A3}) 1 8 1 6f 1 0 6f 82f C _ _ _ _ _ __ A f _ 1+hcc63 <1+h$ 63)+1+hxax{(1+h$)3x}+8y2' The velocity expressed in this system is, from (2.27 ) (Us) ($(1+h$)) 17C: ’11,; = j: . vy 3) 45 The electromagnetic fields and the Lorentz force expressed in this system are, from (2.30), (2.31) and (2.32) 1 6(1) 8A,\ magnum ‘C \___/ H + P‘ (*3 C13 C13 98 / 24an ) Bs 8.1: By Bx 8A8?) 1+h:1(c9 63 'J 1 I 1 — — — —{<1 + mm.) ) \1+h$ as 1+h$823 _' f, E, + vay — 22sz E, + iBy — QB; f0 = fx 2 Ex+vyBs —v3By = Emmi-QB, —s(1+h:r)By fy E, + vsBx — va, Ey + 5(1+ h$)BI — 2'38, 1 8 \ d A3 1 + m5; A,d/dt(—hx) = ——(Ax)— 3 (Q—zTC-KC)+ ”if”: alt A 6 . g 0 K 83/ ) Thus, the equations of motion expressed in this system are, using the Newton’s equa- tion derived at (2.35) H r 1 2) d A, 1+151: as Asd/dt(—h:1:) = e ——(Ax)— _ ((15-170.A'C)+ 1+bhx dt A 01' y 2 0 ) By ) 46 Furthermore, the equations of motion after space-time interchange in this system are, from (2.60), (2.61), (2.62) and (2.63) dill p1: — = 1 I — . d8 ( + mt)ps, (2 64) (1?! pt _ = 1 } _J , V (13 ( + If)“, (2 60) d“ h — e —dA“ Jinn—50 2(0) (13 p3 _ ds 5* (9:1: 6 1 (1y _ gf“ _ e [EE‘E + EB, — (1 +Iz:1:)By] , (2.66) (11),, (1A,, 1 0 _,C ~C _ = __ _ __ (p _ . A (13 e [ d3 3 8y( 2) ) e 1 d2: _ gfy _ e [gEy + (1+ har:)Bat — 5B,] . (2.67) It is customary to have the equations of motion regarding the momentum slopes a = px/po and b = pb/po, instead of p, and py, where p0 is the initial momentum of the reference particle. Then the equations (2.64) and (2.65) can be expressed as d_2: d8 11% d3 where ps/po = \/(p2 - p2 - p§)/p3 = (1+ hm) ps/po’ (2'68) (1+ hzr) 198/190, (269) \/(p/po)2 — a2 -— b2 can be utilized. Because p0 is s-independent, the equations (2.66) and (2.67) can be expressed as do. I d TEx + -_yB3 — P3 6 — 2 h— — 1 h. B . d3 pa + 190 [3 ds ( + r) y] , (2 70) db 6 1 (1:1: _ = _ —E 1 I ‘ Bx _ —BS 0 0 d8 pg [3- ,+( + m7) d3 ] (2 71) Chapter 3 Transfer Maps for Elements Characterized by Measured Fields The use of transfer maps is widespread in the design and study of particle optical systems such as accelerators, spectrometers, beamlines, electron microscopes as well as glass optical system. Combined with the Differential Algebraic (DA) techniques [3] [4] [5] [8], computation of Taylor transfer maps to high order are now performed extensively. The higher order aberrations obtained in this approach offer a more efficient analysis of the particle optical systems, and hence often simplify the design of a system. One of the crucial particle optical systems which require the knowledge of high order aberrations are modern high resolution spectrographs for nuclear physics. Their large phase space acceptance sometimes turns out to demand a careful study of higher order aberrations up to seventh order [16]. To this end, the most critical question is that the simulation treats the field as precisely as possible. Since such spectrographs use large aperture magnets, a careful consideration of the fringe fields is indispens- able. Some efficient methods to take account of the fringe field were discussed in the references [36] [37] [39] [38] in detail. In this chapter, we will discuss a technique to compute Taylor transfer maps for the 47 48 fields from measured data in the framework of the DA methods. We limit ourselves to the common situation in which measurements only exist in the midplane of the device, based on the understanding that the out-of-plane fields can then be calculated. The analysis in this chapter is based on this paradigm, and is applicable whenever the out-of-plane expansion does indeed represent the true fields. The method has been implemented in the the DA based code COSY INFINITY [11] [51], and has been used for the simulations of various spectrographs, including the S800 spectrograph [60] at the National Superconducting Cyclotron Laboratory at Michigan State University and the spectrographs at Jefferson Laboratory. Figure 3.1 shows the measured field data of a bending magnet of the 8800 spectrograph[24] [23] [61] [62] [63], and these field data are used for the computation of the transfer map of the system. Figure 3.1: Measured field of a bending magnet of the S800 spectrograph. The picture shows the field at 65 X 74 points. 49 If there is reliable knowledge of the error bounds of the system, the method can be combined with the Remainder-enhanced Differential Algebraic (RDA) approach to obtain rigorous bounds for the induced errors in the transfer map. The topics related to the RDA are covered in the next chapters. 3.1 Wavelet Representations For the purpose of including measured midplane field data as a part of a particle optical system, a good interpolation method is required to obtain a field value from the limited information contained in the measured field data. The interpolation method has to fit to the DA technique, and in particular the result of the interpolation has to be differentiable as often as needed. Furthermore, the interpolation should be localized in the sense that for the interpolated fields at a certain point, only the fields at neighboring points contribute. Both of these requirements are met by the approach of wavelets [27]. Since our target data is more or less smooth as typically shown in Figure 3.1, we chose the method of Gaussian wavelets, which assures the required differentiability and locality. Assume a set of data Y, is given at N equidistant points :13,- for z' = 1, ..., N. Then, the interpolated value at a point :1: is expressed as N 1 (a: — 56,-)2 V613) = gYil/TTS- 9X1) ]__Aa:2—SZ’] , where A2: is the distance of two neighbouring points 2:,- and 23,-1.1, and S is the factor to control the width of Gaussian wavelets. Figure 3.2 shows how the Gaussian interpolation as a sum of Gaussian wavelets works for several one dimensional functions, including linear functions, a trigonomet- ric function and a Gaussian function. The data values Y, are supplied by taking their function values at each point :r, to simulate the function. 50 ] \ // \ l / I r,\\. ,{vi‘X/A: ‘{ N-rr \ ?\yxl{{)(y {[111] { ll'l’1l\ ‘ " V a? .w . ,1». [ (xii! [((‘xfix/Vv ><\\ .\ («“(KUY ()1,1,’,%\2 ill” l'ny H"? \ “16:51:81” (RA/.1}, _— -' C:~\;~._ .. ;.~ ‘:- -'.I .' .r x_ . ,(y Figure 3.2: Gaussian wavelets representation for f(:1:) = 1 (upper left), f (:1:) = :1: + 1 (upper right), f (:13) = cosa: + 1 (lower left) and f(:1:) = exp(-2:2) (lower right). 51 Table 3.1: Accuracy of the Gaussian wavelets representation for one dimensional functions. Test Number of Average Maximum Function Gaussian Error Error Wavelets N 1 10 2.6 x 10“14 2.6 x 10‘14 :1: +1 10 1.5 x 10‘14 2.6 x 10‘14 cos(:1:) +1 400 6.3 x 10-5 1.0 x 1041 exp(—:1:2) 600 4.6 x 10-5 1.6 x 10-4 There are several strategies for choosing the height of each Gaussian. One ap— proach is to set up a local least squares problem to determine the heights that fit the data optimally. While this approach usually yields the best approximation, it does not preserve the integral of the function, which for our cases is very desirable. Also, the method is not particularly well suited for the treatment of noisy data, which benefits from the smoothing effect achieved by wider Gaussians. Both of these criteria are important for our situation, and so we decided to choose the approach where the height of the data Y, directly determines the height of a Gaussian wavelet exp[—(:r — 36,-)2 / A1132S2] that is placed at 23,-. The resulting interpolated function is shown in each picture, which represents the original function very well. S is chosen to be 1.8 for all the cases. Table 3.1 summarizes the accuracy of the method for those functions, although local accuracy for our purposes is of a lesser concern compared to smoothing and preservation of area. Figure 3.3 shows the behavior of derivatives of interpolated functions up to third order. As an example, the Gaussian function f (11:) = exp(—:1:2) is chosen to be the original function shape. The advantage of the Gaussian function and many other wavelets is that it falls off quickly. Thus the potentially time consuming summation over all wavelets can be 52 replaced by the summation of only the neighboring Gaussian wavelets in the range of :l:8S, which is in the vein of other wavelet transforms and greatly improves efficiency. Figure 3.3: Derivatives of the function f(:1:) = exp(—:1:2) when represented by an ensemble of Gaussian wavelets. The interpolated function (upper left), and the first (upper right), the second (lower left) and the third (lower right) order derivatives of the interpolated function. 53 3.2 Elements Characterized by Measured Fields The Gaussian wavelets representation discussed above is utilized to represent fields which are specified by measured data. Because of the favorable features of the Gaus- sian function, the method allows the computation of the transfer map of such a magnetic field element as long as the out-of-plane expansion of the midplane data is in fact accurate enough to describe the field in the whole space. Similar to the procedure discussed in the previous section, the measured field data is given at equidistant grid points in two dimensional Cartesian coordinates. Figure 3.4 shows how the data grid is specified and the corresponding Cartesian coordinates to the data grid. Assume a set of data B(z'x, 12) is given at equidistant N1, x Nz points (:riz, 2,) for ix 2 1, ..., NI and 2'2 2 1, ..., NZ. Then, the interpolated value at a point (2:, z) is expressed as N; N; . . 1 (:1: __ xix? (Z __ Ziz)2 831(1), 2') = Z Z B(Z$,Zz)7r—S2' exp ———A$—28—2— — W , (3.1) 1'le iz=1 where A2: and Az are the grid spacing in :1: and 2: directions respectively, and S is the control factor of the width of Gaussian wavelets. A suitable choice for the control factor S depends on the behaviour of the origi- nal supplied data. If S is too small, the mountain structure of individual Gaussian wavelets is observed. On the other hand, if S is too large, the original value supplied by the data is washed out, which can also be used effectively for purposes of smoothing noisy data. For constant fields, the suitable S is about 1.8. For quickly varying fields, it should be about 1.0. Larger values of S usually provide more accurate evaluation of the derivatives. We implemented the Gaussian wavelet method in COSY INFINITY [11] [51] in the form of a general particle optical element to compute the transfer map from measured 54 ix x 1 data grid ] particle trajectory Sd> (55,52) 1234...N.—>z'. (0,0) Figure 3.4: Specification of measured field data for a particle optical element in COSY INFINITY. field data. For the purpose to specify the position of the particle in the element, the starting point (Sx, $2) and the direction S¢ of the trajectory of reference particle have to be given as shown in Figure 3.4. The rest of this chapter discusses the transfer maps computed using the new method with Gaussian wavelet representation to treat measured field data in comparison with those transfer maps obtained by default in COSY INFINITY. 3.3 Transfer Maps using Linear Field Data The first test is to check the performance of the new method using artificially created data. We begin with the two perhaps most common cases, namely the main fields of a homogeneous dipole, and an inhomogeneous dipole field with a linear inhomogeneity. We computed the high order transfer maps in four dimensions using the new method, and compared the result with the geometrically computed transfer maps which are obtained by the default procedures in COSY INFINITY. 55 3.3.1 The Homogeneous Dipole Field The first example is a set of data with a constant value, which corresponds to a homogeneous magnetic dipole field. A beam of 200MeV protons was chosen to be the reference, and a total arclength of 1m was demanded in the constant magnetic field of 1Tesla. The reference trajectory bends 26.65° with a radius of 2.15m. The Taylor transfer map computed geometrically by default in COSY INFINITY is listed in Table 3.2 in the standard output notation in COSY INFINITY. Because of the limitation of the space, we show the list of the transfer map only up to third order. The artificial data were prepared with the constant value of 1Tesla positioned at 200 x 200 equidistant points spaced by szAzzlcm, covering the area 2mx2m. The Taylor transfer map was computed using the Gaussian wavelet representation scheme, where the control factor S in (3.1) is set to be 1.8. The difference between the two transfer maps is listed in Table 3.3, where any errors smaller than 10‘10 are listed as 0. Very good agreement was found. 56 Table 3.2: Taylor transfer map of a homogeneous dipole of the radius 2.15m with a bending angle of 26.65° computed geometrically by the default procedure in COSY INFINITY. Transfer Map of a Homogeneous Dipole by default in COSY ( Radius: 2.15m, Angle: 26.65° ) Expansion coefficients of x,a,y,b depending on the exponents of xayb (x, (a, (y, (b, xayb 0.8937341 -0.2086851 0 0 1000 0.9643205 0.8937341 0 0 0100 0 0 1.000000 0 0010 0 0 1.000000 1.000000 0001 -0.4680778E-01 0 0 0 2000 0.4009265 0 0 0 1100 0.1020792 -0.2242986 0 0 0200 0 0 0.4485971 0 1001 0 0 0.2284330 0 0101 -0.1142165 -0.2242986 0 0 0002 -0.1006197 0 0 0 1200 0.4309231 0 0 0 0300 0 0 0.4821603 0 0201 -0.1006197 0 0 0 1002 0.4309231 0 0 0 0102 0 0 0.4821603 0 0003 57 Table 3.3: Difference between the Taylor transfer maps of a homogeneous dipole computed by the default procedure in COSY INFINITY and computed using the new element based on measured field data. Difference of Homogeneous Dipole Maps (Error > 10‘10 ) Expansion coefficients of x,a,y,b depending on the exponents of xayb (x, (a, (y, (b, xayb 0 0.1636561E-09 0 0 2000 0.1045300E-09 0.3457759E-09 0 0 1100 0 0.1274906E-09 0 0 0200 0 0 0 -0.3409426E-09 1010 0 0 -0.1101961E-09 -0.3730441E-09 0110 0 0 -0.1065634E-09 -0.3581052E-09 1001 0 0 0 -0.2630541E-09 0101 -0.3376048E-09 -0.5700778E-09 0 0 0020 -0.4243799E-09 -0.4655945E-09 0 0 0011 -0.4319352E-07 -0.4920250E-07 0 0 3000 -0.1845052E-07 -0.2518897E-07 0 0 2100 0.2875618E-08 -0.9350099E-08 0 0 1200 0.2609582E-08 -0.4923973E-08 0 0 0300 0 0 0.1326292E-06 0.1607858E-06 2010 0 0 0.3670618E-07 0.5794910E-07 1110 0 0 -0.3425433E-08 0.1218639E-07 0210 0 0 0.1880360E-07 0.28003608-07 2001 0 0 -0.5815043E-08 0.2099834E-07 1101 0 0 -0.7957903E-08 0.1552648E-07 0201 0.1161989E-06 0.1614849E-06 0 0 1020 0.1088566E-07 0.5549496E-07 0 0 0120 0.2288686E-07 0.1029208E-06 0 0 1011 -0.1375535E-07 0.9141796E-07 0 0 0111 -0.6363773E-08 0.4116490E-07 0 0 1002 -0.9467137E-08 0.5294305E-07 0 0 0102 0 0 -0.3965713E-07 -0.5944019E-07 0030 0 0 -0.1135347E-07 -0.5934921E-07 0021 0 0 0.6816748E-08 -0.4799466E-07 0012 0 O 0.3175650E-08 -0.1827340E-07 0003 58 3.3.2 Inhomogeneous Dipole Field The second example is a set of data which depend linearly on the radius r of the motion of the particles. The field strength is described as By(7‘)=Bo-(1—n-T;ZORO), where R0 is the radius of the reference trajectory, Bo is the field strength on the reference trajectory and n is the linear inhomogeneity constant. This is the field of an inhomogeneous dipole magnet with a linear inhomogeneity. The same setting was chosen, namely 200 MeV proton beam with the total arclength of 1m. By the choice of Bo = 1Tesla, the same bend is defined for the reference trajectory, namely an angle of 26.65° and a radius of 2.15m. The linear inhomogeneity constant 71 was set to be 0.5. The computations of the Taylor transfer maps were performed in a similar way as in the case of the homogeneous dipole. Table 3.4 shows the third order Taylor transfer map in six dimension computed geometrically by the default procedure in COSY INFINITY, and Table 3.5 shows the difference between the two transfer maps. 59 Table 3.4: Taylor transfer map of an inhomogeneous dipole of radius 2.15m and bending angle 26.65° with inhomogeneity 0.5 computed geometrically by the default procedure in COSY INFINITY. Transfer Map of a Inhomogeneous Dipole by default in COSY ( Radius: 2.15m, Angle: 26.65°, Inhomogeneity: 0.5 ) Expansion coefficients of x,a,y,b depending on the exponents of xayb (x, (a, (y, (b, xayb 0.9463845 -0.1062624 0 0 1000 0.9820634 0.9463845 0 0 0100 0 0 0.9463845 -0.1062624 0010 0 0 0.9820634 0.9463845 0001 0.2228769E-03 0.4766584E-01 0 0 2000 0.4486863 0.4810034E-01 0 0 1100 0.1131941 -0.2120960 0 0 0200 0 0 -0.4854609E-01 -0.4678238E-01 1010 0 0 -0.2449431E-01 -0.2360442E-01 0110 0 0 0.4323567 -0.4854609E-01 1001 0 0 0.2181490 -0.2449431E-01 0101 -0.1247084E-01 -0.2471637E-01 0 0 0020 -0. 1 152539 -0 . 2284255 0 0 0002 0.1103990E-01 0.4053979E-03 0 0 3000 0.1668532E-01 0.1144351E-01 0 0 2100 -0.1673401E-01 0.8570764E-02 0 0 1200 0.4663577 0.2705151E-01 0 0 0300 0 0 -0.1084069E-01 0.1201080E-02 2010 0 0 -0.1840065E-01 -0.9773169E-02 1110 0 0 -0.2976516E-01 -0.5771860E-03 0210 0 0 -0.9231508E-02 -0.1073923E-01 2001 0 0 0.3906643E-01 -0.2020779E-01 1101 0 0 0.4790956 -0.3070280E-01 0201 -0.5797469E-02 -0.6097624E-03 0 0 1020 -0.1991462E-02 0.5388325E-02 0 0 0120 0.3675865E-02 -0.1014555E-03 0 0 1011 -0.5031829E-01 0.1868718E-02 0 0 0111 -0.7733214E-01 -0.4731772E-02 0 0 1002 0.4297353 0.2417023E-01 0 0 0102 0 0 -0.1828454E-02 -0.3428371E-02 0030 0 0 0.1808343E-02 -0.5390550E-02 0021 0 0 -0.7450909E-01 0.9639483E-03 0012 0 0 0.4479615 -0.2512844E-01 0003 60 Table 3.5: Difference between the Taylor transfer maps of an inhomogeneous dipole computed by the default procedure in COSY INFINITY and computed using the measured field element. Difference of Inhomogeneous Dipole Maps (Error > 10‘10 ) Expansion coefficients of x,a,y,b depending on the exponents of xayb (b, (x, (a, (y, xayb 00000000 I 000 .9479292E-06 .3228052E-06 .4311228E-06 .3413067E-05 .9055087E-06 .20505803-06 .4247148E-06 .8749727E-06 .1363762E-06 .1537549E-06 .3419442E-06 .6032425E-06 .87479558-07 .2332158E-07 .2060690E-06 .7502058E-06 .1221911E-05 .1183294E-05 -0 -0 0 0 0 I O 0 000000 I I O O OOOOOOOOO .1894282E-05 .9815221E-06 .3774559E-06 .3804320E-06 .1825359E-05 .38082663-06 .1243242E-05 .14663163-05 .6688967E-07 .12272748-06 .5660030E-07 .4360596E-06 .3234010E-06 .2365737E-06 .3131300E-06 .3136700E-06 .3597991E-07 .3798253E-06 .9313577E-06 .3127151E-06 .3898681E-06 .4219984E-06 .3580534E-05 .1822555E-05 .4121860E-07 .2255669E-06 .49981288-06 .16133733-06 .6451915E-06 .3780246E-06 .5015356E-07 .2808341E-06 .1067846E-05 .8565441E-06 000000000 -0 .1828878E-05 .9313569E-06 .5115693E-06 .4020970E-07 .3969166E-06 .2028945E-06 .1766775E-06 .2357229E-06 .1999980E-07 .1114431E-06 .2972702E-06 .5160206E-06 .1325949E-06 .2199493E-06 .2708509E-06 .3690355E-06 1000 0100 0010 0001 2000 1100 0200 1010 0110 1001 0101 0020 0011 0002 3000 2100 1200 0300 2010 1110 0210 2001 1101 0201 1020 0120 1011 0111 1002 0102 0030 0021 0012 0003 61 3.4 Transfer Maps using Analytical Field Data In the previous section, we have discussed the performance of the method to compute Taylor transfer maps for typical main fields using the Gaussian wavelet representation. In the next example we apply the method to analytically created data of a typical fringe field profile. We prepared the data such that the computation can be checked with some existing procedures in COSY INFINITY. For the purpose to obtain the analytically described field values for the fringe field effects, we used an internal mechanism in COSY INFINITY for the treatment of the dipole field. It is based on the standardized description of the s-dependence of multipole strengths by an Enge function as in the program RAYTRACE [47]. The Enge function has the form 1 F(Z) : 1 + exp(a1+ 02 - (z/D) + + 05 - (z/D)5)’ where z is the distance perpendicular to the effective field boundary, and D is the full aperture of the dipole. a1 through as are the Enge coefficients, which depend on the details of the geometry of the element, including shimming and saturation effects. We used the same dipole magnet setup with the example in the previous section, namely a bend angle of 26.65° and a radius of 2.15m. For the purpose to have enough fringing region, we have drifts with the length of 0.4m before and after the dipole, which results in a total arclength of 1.8m. The aperture of the dipole magnet was set to be 10cm, and the default Enge coefficients for a dipole in COSY INFINITY were used. The analytical field distribution, which has a maximum field strength of 1 Tesla, is shown in Figure 3.5. In the upper picture, a beam enters from the center at the left edge rightward, then goes clockwise to the right down. We computed the first order transfer maps in four dimension using the default 62 Figure 3.5: Analytical dipole field distribution with fringe field effects. Bend: 26.65°. Radius: 2.15m. Aperture: 10cm. In the upper picture, a beam enters from the center at the left rightward then goes clockwise. The lower picture is viewed from a lower angle. 63 procedure in COSY INFINITY and using the analytically created field data with the Gaussian wavelet representation method with various values of the control factor S ranging from 1.8 to 5.0. The result in Table 3.6 shows the transfer map obtained by the default procedure in COSY INFINITY, and the difference between the default computation and the field data computation. The analytical data was prepared at N, x N, = 100 x 200 equidistant points spaced by szAzz-lcm, which covers the area 1m x2m. Because of the somewhat coarse spacing, the resulting transfer map agrees with the true map to only about 10’3. In the next computation, a finer spacing of the data grid points is being used, resulting in a higher overall accuracy of about 5 x 10“5 up to order five. Table 3.6: Transfer maps of a dipole with fringe field effects computed by default in 64 COSY INFINITY, and using analytically generated field data. Transfer Map of a Dipole with Fringe Field Effects by default in COSY ( Radius: 2.15m, Angle: 26.65°, Aperture: 10cm, two drifts of 0.4m ) Expansion coefficients of x,a,y,b depending on the exponents of xayb (x, (a, (y, (b, xayb 0.8102451 -0.2086852 0 0 1000 1.646121 0.8102227 0 0 0100 0 0 1.029143 0.3268649E-01 0010 0 0 1.820497 1.029503 0001 Errors of Transfer Map using Analytical Field Data (> 1040) 8:18 -0.2102107E-04 -0.1463520E-05 0 0 1000 -0.4156249E-05 0.1911708E-04 0 0 0100 0 0 -0.14957505-03 -0.1895449E-03 0010 0 0 -0.1400689E-03 -0.1901141E-03 0001 S=2.6 -0.2377301E-04 -0.1707500E-05 0 0 1000 -0.8942274E-05 0.2260579E-04 0 0 0100 0 0 -0.3316729E-03 -0.3975747E-03 0010 0 0 '0.2919843E-03 -0.3807626E-03 0001 S=3.4 -0.2754365E-04 -0.1973192E-05 0 0 1000 -0.1546368E-04 0.2751592E-04 0 0 0100 0 0 -0.5748396E-03 -0.6755873E-03 0010 0 0 -0.4952997E-03 -0.6357372E-03 0001 S242 -0.3237764E-04 -0.2361021E-05 0 0 1000 -0.2370848£-04 0.3368500E-04 0 0 0100 0 0 -0.8752840E-03 -0.1019264E-02 0010 0 0 -0.7472205E-03 -0.9510957E-03 0001 S=5.0 -0.3862635E-04 -0.4325084E-05 0 0 1000 -0.3407721E-04 0.3861323E-04 0 0 0100 0 0 -0.1228760E-02 -0.1422527E-02 0010 0 0 -0.1044417E-02 -0.1320222E-02 0001 65 The last computation using artificially created field data assesses the prospects of computing a Taylor transfer map from the measured field data of the S800 spec- trograph. The system setup consists of a dipole magnet, which has the same design parameters with those dipoles of the S800 spectrograph, and two drifts before and after the dipole. The parameters of the S800 system is listed in Table 3.8. The bend angle of the reference trajectory is 75° with a bending radius of 2.68m, and the dipole has an edge angle of 30°. The aperture is 7.62cm, and the length of the two drifts is 0.7m each, so the total arclength adds up to 4.91m. The default Enge coefficients for a dipole in COSY INFINITY were used to take the fringe field effects into account, and a beam of tritium 3H with the magnetic rigidity 2.704Tesla-m was used. The re- sulting analytical field distribution is shown in Figure 3.6, where the maximum field strength is 1.01Tesla. Despite of the bending angle of 75°, the angle of 45° seen in Figure 3.6 is due to the exit edge angle of 30°. In the upper picture, the beam enters from the upper left rightward, then goes clockwise and ends at the lower right. We computed the fifth order transfer maps in six dimensions using the default procedure in COSY INFINITY, and using the analytically created field data with the Gaussian wavelet representation method with the control factor S 1.8, 2.2 or 2.6. As a further check, it was calculated to what extent the transfer map breaks the symplectic symmetry, which is a general overall test for the accuracy of calculations. For reasons of space, we only list the relevant parts of the resulting maps in Table 3.7. The analytical data was prepared at N, X N, = 661 x 801 equidistant points spaced by Ax=Az=5mm, which covers the area 3.3mx4.0m. The result of the transfer maps using the data is very accurate, and shows the practical usefulness of the new method. 66 (u. e x ..1 Ln) .1. -1 u u :1 a. .d- .n 9 ..x l 1 m a. .u -: ) ,,, -:r m u. in- 1. '1 ~- 1 .1. (a) ...v «.1 J‘ 1 »u .2. in. .s w :2 .4 51 :a. 1m in. ..r u» :1 art a .m cm nm u; (a. 1:. 15. 5a .5. ”V ‘I Figure 3.6: Analytical S800-like dipole field distribution with fringe field effects. In the upper picture, a beam enters from the upper left rightward, then goes clockwise. The lower picture is viewed from a lower angle. 67 Table 3.7: A part of the fifth order Taylor transfer map of a S800-like dipole with the fringe field effects computed by the default procedure in COSY INFINITY, and the error of transfer map using analytical field data. Transfer Map of a S800-like Dipole with Fringe Field Effects and Two 0.7m Drifts by default in COSY ( Radius: 2.682m, Angle: 75°, Aperture: 7.62cm, Exit edge angle: 30° ) Symplectic Error: 2.07406 x 10‘7 Expansion coefficients of x,a,y,b,t depending on the exponents of xayb6 (x, (a, (y, (b, (t, xaybd 0.46777E-1 -0.30432 0 0 -0.49223 10000 3.19321 0.60316 0 0 -1.35783 01000 0 0 0.88841 -0.19817 0 00100 0 0 4.34429 0.15655 0 00010 1.50829 0.71012 0 0 0.93807 00001 0.46761 -0.51206 0 0 -0.90161 02000 -0.5951OE-1 -0.13825 0 0 0.66269 01001 -0.96599 -0.31226 0 0 -0.82018 00002 0.17337 0.71513E-1 0 0 -0.38834 03000 -0.11038 -0.21087 0 0 -0.24251 04000 -0.32441 0.46850E-1 0 0 -0.22583 05000 Errors of T1ansfer Map using Analytical Field Data (> 10‘”) S=1.8, Symplectic Error: 2.46947 x 10’6 -0.37943E-5 -0.40253E—6 0 0 -0.17484E-5 10000 0.31532E-5 0.95363E—6 0 0 -0.40188E-5 01000 0 0 -0.10045E-3 -0.42308E-4 0 00100 0 0 -0.18879E-3 -0.14706E-3 0 00010 0.12098E-4 0.32933E-5 0 0 0.387903-5 00001 -0.84006E-5 -0.49003£-5 0 0 -0.94809E-5 02000 -0.17281E-4 -0.29639£-5 0 0 -0.16631E-5 01001 -0.25206E-4 -0.34907E-5 0 0 -0.10033E-4 00002 0.16000E-4 0.31675E-5 0 0 -0.45384E-6 03000 -0.31202E-4 -0.80058E-5 0 0 -0.13688E-4 04000 0.70371E-3 0.99258E-4 0 0 -0.27143E-3 05000 [ S=2.2, Symplectic Error: 1.38298 x 10‘6 ] | S=2.6, Symplectic Error: 8.55948 x 10“7 ] 68 3.5 Transfer Maps of the 8800 Spectrograph The S800 spectrograph [60] at the National Superconducting Cyclotron Laboratory at Michigan State University is a modern high resolution spectrograph for nuclear physics, which offers large phase space acceptances with A6 = i60mr, Ad) : 21:90mr, and AE = :l:5%. Because of the large phase space acceptances, higher order aberra- tions severely affect the resolution. A method to reconstruct the particle trajectories in spectrographs was presented earlier [16], where an energy resolution of 1 / 9,100 was recovered with fifth order reconstructive correction from an uncorrected nonlinear res- olution of 1 / 100, almost reaching the linear resolution of 1 / 9,400. The limitation of the reconstructive correction method comes from various inaccuracies including the transfer map of the system. The two huge dipole magnets, which are the primary parts of the S800 spectrograph as a beam optical device, play a very important role for the transfer map, and a rough estimate shows that a knowledge of the overall fields with a relative accuracy of 10‘4 is necessary to achieve the design energy resolution of 1/10,000. In the following we utilize the method to include measured field data to transfer maps developed in this chapter, and assess to what extent the midplane field data are indeed sufficient to describe the field in the entire space. The laboratory layout and the system parameters of the S800 spectrograph are shown in Figure 3.7 and Table 3.8. 69 Figure 3.7: Layout of the 5800 spectrograph of the National Superconducting Cy- clotron Laboratory at Michigan State University. Courtesy Daniel Bazin. Table 3.8: Design parameters of the S800 spectrograph. Drift l = 60 cm Quad [2 40 cm Gm“ = 21 T/m r = 0.1 m Drift l : 20 cm Quad 1: 40 cm GM, 2 6.8 T/m r = 0.2 m Drift l = 50 cm Dipole p = 2.675 m Bma, : 1.5 T d) = 75° 61 = 0° 62 = 30° Drift l = 140 cm Dipole p = 2.675 m Bma, : 1.5 T (15 = 75° 61 = 30° 62 2 0° Drift l = 257.5 cm 3.5.1 Measured Field Data of the 8800 Dipole The fields of the magnets of the S800 have been mapped [24] [23] [61] [62] [63]; for con— venience, the data were taken in cylindrical coordinates whose z axis is perpendicular to the plane of Figure 3.7. The raw data was converted to the equidistant Cartesian grid data after a smoothing process to attempt the control of noise. The measured field data are prepared in a rectangular area. Figure 3.8 shows the two rectangular areas to cover the two dipole magnets. Figure 3.9 shows the field distribution of the 70 first S800 dipole magnet, which has a maximum field strength 0.965Tesla. In the up- per picture, a beam enters from the lower left rightward and travels counterclockwise. The central area of the magnet, which is seen to have a constant field in Figure 3.9, has a structure as shown in Figure 3.10 by cutting the field values to only those above 0.962Tesla. The measured field data have been supplied for use with the Gaussian wavelet representation method in the following settings. 1. N, x N, = 328 x 370 equidistant point data spaced by szAzz-lcm, covering the area 3.27mx3.69m. 2. N, x N, = 158 x 184 equidistant point data spaced by AzzAzz2cm, covering the area 3.14mx3.66m. 3. N, x N, = 631 x 736 equidistant point data spaced by szAz=5mm, covering the area 3.15mx3.675m. The second and the third are based on the same data with the different spacing. The directions of :1:- and z axes follow the notation in Figure 3.4. ‘- For the third data set, which is the finest, further details of the data are studied. Figure 3.11 covers a small rectangular area 34.5cmx34.5cm with grid points ranging z',=82,..,151, which is around 1 / 6 from the bottom vertically, and i,=331,..,400, which is around the center horizontally in the upper picture of Figure 3.9. In this region, the field values are ranging from 0.964Tesla to 0.965Tesla. The step structure seen in the top picture in Figure 3.11 is due to the fact that the data are available only to an accuracy of 6 digits. While any obvious noise from the raw data was removed, this step structure in the data will be a source of a numerical problem in the detailed high order Taylor transfer map computation. To study this aspect, the data is smoothed 71 to double precision accuracy on computers using the Gaussian interpolation method, which is shown in the middle picture of Figure 3.11. Then the original data (upper) and the smoothed data (middle) were compared, and the result is shown in the bottom picture. The difference ranges between —6.79 x 10‘6 and 7.44 x 10‘6, as expected. This subtle difference is sufficient to induce enough noise in the out-of-plane ex- pansion to become noticeable in the transfer map. The difference in the transfer maps range from about 1% in the linear part, and about 10% in the second order part, to as much as 100% for certain third order terms. This seems to be an indication that the common practice of midplane-only measurements has to be used with great care. Fortunately for the S800, out of plane data are available, and it is planned that they be used in the future for a more accurate representation of the overall field. On the other hand, the S800 dipoles are rather extreme regarding the total amount of deflection of the orbits, which is one of the main reasons why the higher order part of the transfer map is so sensitive. For other systems with smaller deflection angles, the situation is much less dramatic; as an example, for certain dipoles in the Fermilab recycler ring, which have a main field strength 0.13Tesla and a bend of 2°, several digits of accuracy in the matrix elements were found [65]. In any case, for the system like the S800 spectrograph, it turned out that it is doubtful whether in-plane data are sufficient at the current level of noise to describe the out-of-plane fields sufficiently accurately without further processing. Figure 3.8: Two rectangular areas of the measured field data of the two S800 dipoles. Courtesy Daniel Bazin. 73 Figure 3.9: Field distribution of the S800 dipole measured data. Maximum: 0.965 Tesla. In the upper picture, a beam enters from the lower left rightward and goes counterclockwise. The lower picture is viewed from a lower angle. 74 J Figure 3.10: Top part of the field distribution of the S800 dipole measured data. The picture shows the field values ranging from 0.962Tesla to the maximum 0.965Tesla. 75 Figure 3.11: Further detail of the S800 dipole measured data in a small area 34.5cm x 34.5cm on the top region shows the step structure due to the limited data digits (top). The data is smoothed (middle), and compared (bottom). Chapter 4 Remainder-enhanced Differential Algebra (RDA) Methods In this chapter it is shown how, in parallel to the accumulation of derivatives, Taylor remainder bounds for all functional dependencies can be carried along the computa- tion. Chapter 6 addresses the necessary Differential Algebraic (DA) formalism that allows the solutions of ODEs and the determination of flows. The additional com- putational effort is limited, and the resulting bounds are usually rather sharp, in particular at higher orders. Compared to other rigorous bounding methods, signifi- cant advantages arise for modest to complicated functional dependencies, especially in higher dimensions. The next chapter discusses an implementation of the method in the code COSY INFINITY [6] [7] [11] [15]. 4. 1 Introduction The significant advances in computer hardware that we have experienced particularly in the past decade allow the study of ever more complex problems. In many practi- cal problems, one must locally model nonlinear functional dependencies, for example to study parameter sensitivity or to perform optimization. This is typically done through the computation of derivatives, and the methods of computational differen- 76 V 'L... 77 tiation [13] [31], and for the solutions of ODEs and PDEs the Differential Algebraic (DA) enhancements, have excelled in providing such derivatives accurately and inex- pensively. While the derivatives themselves are accurate except for computational errors that are typically very small, rigor is lost when the derivatives are used to model a functional dependence, because of the lack of information about the size of remainder terms. The method of interval arithmetic (see, for example, [48]) often provides a means for keeping the mathematical rigor in the computation of model functions. However, the naive use of interval methods in large problems is prone to blow-up, which at times limits the practical usefulness of such methods. The new technique, the method of Remainder-enhanced Differential Algebras (RDA), combines the DA approach to obtain Taylor polynomials with interval tech- niques to determine bounds for the Taylor remainder. The new method uses the advantages and diminishes the disadvantages of either of the two methods. In beam physics and weakly nonlinear dynamics in general, the DA technique [3] [4] [5] [8] has offered a robust way to study the nonlinear behavior of beams and has evolved into one of the essential tools. However, many difficult yet important questions remained, including the long-term stability of beams in circular accelerators or other dynamical systems, such as the planets in the solar system. In light of modern stability theories, this problem can be cast in the form of an optimization problem [14] [15], which from the computational point of view is very complex. 4.2 Differential Algebra and Interval Arithmetic While DA methods [6] [7] [11] [15] can provide the derivatives of functional depen- dencies and solutions of ODEs to high orders, in a rigorous sense they fail to provide 78 information about the range of the function. A simple example that dramatically illustrates this phenomenon is the function shown in Figure 4.1 0 if:1: = 0 f(:r) : { exp(—1/:r2) else. (4'1) The value of the function and all the derivatives at :1: = 0 are 0. Thus the Taylor polynomial at the reference point :1: = 0 is just the constant 0. In particular, this also implies that the Taylor expansion of f converges everywhere, but it fails to agree with f (:1:) everywhere but at :1: = 0. 1.2 T MIT T I I I I I T _0.2 1 1 1 1 l 1 1 1 1 -10 -8 -6 -4 -2 0 2 4 6 8 10 Figure 4.1: Function f (:1:) = exp(—1/:1:2) if :1: # 0 ; 0 else, and its Taylor polynomial, which vanishes identically. For the purpose of bounding functional dependencies, the methods of interval arithmetic provides a conceptual contrast. Both extended domains of numbers as well as individual real numbers are represented via rigorous inclusions of floating point intervals. Arithmetic operations are introduced on intervals such that for any numbers in the intervals, a real arithmetic operation on the two numbers always 79 Table 4.1: Elementary properties of interval arithmetic; 11 = [(11,191], I2 = [(12, b2]. 11+12 1: [a1+a2,b1+bg] —11=[—b1,“611] 11-12 = [min(a1a2,a1b2, biaz. blbz),maX(a1a2,a1b2.b1a2,b1b2)] 110 gr 1,, 1/11 = [1/b1,1/a1] leads to a result that is contained in the interval obtained from the corresponding arithmetic operation on the intervals. Table 4.1 lists some elementary properties of interval arithmetic. By evaluating a function in interval arithmetic, it is thus possible to carry rigorous bounds information through the operations, and in the end obtain rigorous bounds of the function. However, while reasonably fast in practice, interval methods have some severe disadvantages, which limits their applicability for complicated functions. First, the width of resulting intervals scales with the width of the original intervals; and second, artificial blow-up usually occurs in extended calculations. Another practical limitation arises if scanning with small intervals is needed in the case of multiple dimensions because of the fast increase of the computational expense. To illustrate the blow-up phenomenon with a trivial example, we consider the interval I = [a, b], which has the width b — a. We compute the addition of I to itself and its subtraction from itself: 1+1 = [(1,1)] + [a,b] = [a+a,b+b] = [20,20] I—I = [(1,1)] —[a,b] = [a,b]+[—b,-a] = [a—b,b—a]. In both cases the resulting width is 2(b—a), which is twice the original Width, although 80 we know that regardless of what unknown quantity :6 is characterized by I , certainly :1: — a: should equal zero. Polynomials under their conventional operations form a commutative algebra with unity. However, interval arithmetic does not have even a group structure for either addition or multiplication, since intervals with nonzero width have no inverses. Fur- thermore, instead of distributivity, we have the sub-distributivity, 11'(12+13)§_11'12+11'I3. (4.2) The concepts of interval methods are discussed at length and in depth in several sources, including [48]. 4.3 Remainder-enhanced Differential Algebraic Operations In this section, we develop a new method that combines the advantage of rigor of the interval approach, while largely avoiding the blow-up problem through the use of DA techniques. The key idea is to describe the bulk of the functional dependence through a Taylor polynomial, and bound the deviation of the original function from the Taylor polynomial by an interval. In this endeavor, the Taylor theorem plays an important role. Theorem (Taylor): Suppose that a function f : [(1', 5] C R” —> R is (71 +1) times continuously partially differentiable 0n [6, 6] Assume :30 6 [613,5]. Then for each :6 E [(3, 5], there is 0 E R with 0 < 0 < 1 such that Mia 61.0 1 O(.._.,,..,n+l,,.,..._...,, V20 1" 81 -o —o k where the partial differential operator (h - V) operates as 8k a,?...a,3‘ k! - . 21 z . 'h’l ...hv" ill'ulv. (W)? 2 Depending on the situation at hand, the remainder term also can be cast into a variety of well-known other forms. Taylor’s theorem allows a quantitative estimate of the error that is to be expected when approximating a function by its Taylor polynomial. Furthermore, it even offers a way to obtain bounds for the error in practice, based on bounding the (n + 1)st derivative, a method that has sometimes been employed in interval calculations. Roughly speaking, Taylor’s theorem suggests that in many cases the error de- creases with the order as the width of the interval raised to the order being consid- ered, and its practical use is often connected to this observation. However, certain examples illustrate that this behavior does not have to occur; one such example is (4.1) in the previous section. For notational convenience, we introduce a parameter oz to describe the details of a given Taylor expansion, namely, the order of the Taylor polynomial n, and the reference point of expansion :30. For the purpose to derive bounds for the remainder, it is also necessary to include the domain interval [513, 5] on which the function is to be considered; altogether, we have -9 a = (n, (30, [6, 0]). (4.3) We now write an (n+1) times continuously partially differentiable function f : [6, l1] C R” —> R as a sum of its Taylor polynomial P0,; of nth order and a remainder 50,; as f(f) = POJCE — £0) + Ea,f(f— f0), 82 where sa,f(:i:' — 2'50) is continuous (even continuously differentiable) on the domain interval and thus bounded. Let the interval [0,, be such that —o vs 6 [6, b], 5,,(6 — 65,) e 1,, Then V55 6 [5, 13'], f(:1‘:') e P0,,(5' — :80) + 1,, (4.4) Because of the special form of the Taylor remainder term End) in practice the re- mainder usually decreases as If — :E'OI’H’I. Hence, if [:15 —— 550] is chosen to be small, the interval Imf, which from now on we refer to as the interval remainder bound, can become so small that even the effect of considerable blow-up is not detrimental. The set Pa,f(:1:' —— :50) + 16,} containing f consists of the Taylor polynomial Pa,f(:7:' — :30) and the interval remainder bounds IaJ. We say a pair (Pay, 10,!) of a Taylor polynomial Pa,f(§:' — .50) and an interval remainder bounds [0,, is a Taylor model of f if and only if (4.4) is satisfied. In this case, we denote the Taylor model by T0,, = (POJJOJ). We call n the order of the Taylor model, 550 the reference point of the Taylor model, [(1', 5] the domain interval of the Taylor model, and a the parameter of the Taylor model. In the following, we deve10p tools that allow us to efficiently calculate Taylor models for all functions representable on a computer. The key is to begin with the Taylor model for the identity function, which is trivial, and then successively build up Taylor models for the total function from its pieces. This requires methods to determine Taylor models for sums and products from those of the summands or factors, as well as from intrinsics applied to functions with known Taylor model. 83 4.3.1 Addition and Multiplication In this subsection, we discuss how a Taylor model of a sum or product of two func- tions can be obtained from the Taylor models of the two individual functions. This represents the first step toward the computation of Taylor models for any function that can be represented on a computer. Let the functions f, g : [5, h] C R” ——> R have Taylor models T0,] : (Pa,f)Ia,f) and Ta,g : (Pa,gaIa,g)a which entails that Vf€[a',li], f(:I:’) e Faye—sonny and -o 9(5) 6 Pa,g(f— 0) + [0,9' Then it is straightforward to obtain a Taylor model for f +9; in fact, for any :3 E [6, 5], f(;i:') + g(:z‘:') e (Pa,f(.’1—3‘ — £0) + 10,) + (Pa,g(:i:’ — 2‘0) + 1,9) = (Pm/(53' - 550) + Pa.g(i" - 550)) + (16,; + 16.9). so that a Taylor model T0,“, for f + 9 can be obtained via P0,,” = 0,; + P0,, and Ia,f+g = 0,; + 10,9. (4.5) Thus we define Ta,f + Ta,g : (Pa,f + Pa,ga 10,] + 10,9), and we obtain that T0,; + T0,, = (Pa,f+g,Ia,f+g) is a Taylor model for f + 9. Note that the above addition of Taylor models is both commutative and associative. The goal in defining a multiplication of Taylor models is to determine a Taylor model for f - g from the knowledge of the Taylor models T a, f and To, for f and g. 84 Observe that for any :1? E [6, ], ft?) 9b?) 6 (Ruff/E ‘“ £0) + Imf) ' (Push-3. _ 550) + 10.9) +Pa,f(:1':' — *0) ~10, + Pa,g(:1":' — f0) - I,” + 10,, - 10,9. Note that Pa’f - P0,, is a polynomial of (2n)th order. We split it into the part of up to nth order, which agrees with the Taylor polynomial P0,,” of order n of f - g, and the extra polynomial P,, so that we have Pa,f(:i" — f0) - Pa,g(:1‘:' — :30) = P0,].g(:1“:' — :50) + P,(:1:' — 550). (4.6) A Taylor model for f - 9 can now be obtained by finding a bound interval for all the terms except P0,}, For this purpose, let B(P) be bounds of the polynomial P : [(1', 6] C R” -—) R, namely, Apparently the efficient practical determination of B(P) is not completely trivial; depending on the order and number of variables, different strategies may be employed, ranging from analytical estimates to interval evaluations. However, thanks to the specific circumstances, the occurring contributions are very small, and even moderate overestimation is not critical. Various methods for determining B(P) will be discussed below. Altogether, interval remainder bounds for f - 9 can be found via I,” = B(Pe) + B(POJ) - Ia, + B(ng) ~10“; + 10,, '10,. (4.7) Thus we define TQJ - T0,, = (190,”, 10,”), and obtain that T0,; -Ta,g is a Taylor model for f - 9. Note that commutativity of multiplication holds, T0,; To, = To, -T,,,f, while 85 multiplication is not generally associative, and also distributivity does not generally hold. While the idea of Taylor models of constant functions is almost trivial, we mention it for the sake of completeness. For a constant function f(.1':') E t, the Taylor model of f is Ta,f E Ta,t : (Pa,t1[a,t) Z (t) [070]) Having introduced addition and multiplication as well as scalar multiplication, we can compute any polynomial of a Taylor model. Let Q( f) be a polynomial of a function f, that is, Q(f) : to +t1f + t2f2 + ---+ tkfk. In practice it is useful to evaluate Q( f ) via Horner’s scheme, Q(f):t0+f'(tl'l'f'(t2+f'('°'(tk—1+f’tk)"')))1 (4.8) in order to minimize operations. Furthermore, Horner’s scheme is often of advan- tage for interval related arithmetic because of the sub-distributivity (4.2) of interval arithmetic. Assume that we have already found the Taylor model of the function f to be T a, f = (POJ, 10,1). Then, using additions and multiplications of Taylor models described above, we can compute a Taylor model for the function Q( f ) via Tape) = (Porno), 10.00)) - 4.3.2 Intrinsic Functions In the preceding subsection, we showed how Taylor models for sums and products of functions can be obtained from those of the individual functions. The computation led to the definition of addition and multiplication of Taylor models. Here we study the computation of Taylor models for intrinsic functions, including the reciprocal applied to a given function f from the Taylor model of f. 86 The key idea is to employ Taylor’s theorem of the function under consideration: However, in order to ensure that the resulting remainder term yields a small remainder interval and does not contribute anything to the Taylor polynomial, several additional steps are necessary. Let us begin the study with the exponential function. Assume that we have already found the Taylor model of the function f to be To“ 2 (RM, IaJ). Write the constant part of the function f around :1:}, as emf, which agrees with the constant part of the Taylor polynomial Pa,f2 and write the remaining part as f; that is, 1" f(i") = 66.1 + f0?)- A Taylor model of f— is then T0,; = (me, [0]), where Pa,f(f - f0) = Pa‘f(.’i" — f0) — Ca,f and IQ")? = IQJ. Now we can write exp(f(§:')) = exp (ca,,+ f(1:*))=exp(c..,;)~exp (f'(:i')) : BXP(Ca,f) ' {1+ R5) + 2170(5))? + ' ° ' + filqfink 1 — .. Ic+l ’ -' +0. + 06“”) exp (9'f($))l’ where 0 < 6 < 1. Taking It 2 n, where n is the order of Taylor model, the part we...» - {1+ fa) + 2, (116))? + - - - + $661)") is a polynomial of f, of which we can obtain the Taylor model as outlined in the preceding subsection. The remainder part of exp( f (27)), 1 €XP(Ca,f) ' {NUQDHI + ' ' ' + 1 (1+1)! mar“ exp (0 - 1(6)}, (49) will be bounded by an interval. Since Pad-(if — f0) does not have a constant part, (RM-(:1? — 50))m starts from mth order. Thus, in the Taylor model computation, the 87 remainder part (4.9) has vanishing polynomial part. The remainder bound interval for the Lagrange remainder term 1 "’ ~ +1 - _. CXp(Ca,I) (k +1),(f(113))’c exp (9 - f(:r)) can be estimated because, for any :1? E [5, f1], PaJ-(f - 2'50) 6 B(Pafl, and 0 < 6 < 1, and so mar“ exp (6 - 1(6)) 6 (B(P.,;) + 11;)” exp ([0, 11 - (801,) + I..,,—)). (4.10) Since the exponential function is monotonically increasing, the estimation of the in- terval bounds of the part exp ([0,1] - (B(Paf) + 1a,,» is achieved by inserting the upper and lower bounds of the argument in the exponential. A Taylor model for the logarithm of a function f can be computed in a similar manner from the Taylor model of the function. In this case, there is the limitation that it has to be ensured that the range of the function f lies entirely within the range of definition of the logarithm, which will be the case if V6 H1135], 190,06 — 60) + 1,, c (0, 00). For the actual computation, we again split the constant part of the function f around 550 from the rest f(:1':') : coy + f(:i:'). Then we obtain 106063)) = 108(CaJ + 13(5)) =108{Ca,f' (1 + Hi)» can! R537) 2 logcaJ + log (I + C0,; _ , ff?) 1 (1755))2 110(5))" ._ logCn‘f'l‘EZ- g—C—if— +(-1)k+ )5 CE; +(_ )k+2 1 (f(f))k+1 1 88 where 0 < 6 < 1. Taking k 2 n, the part ft?) 16(5))? .. 110(6))" ..., 2737+W+Hl+aw log c, ,- +— is again treated as a polynomial of f in the Taylor model computation. The Lagrange remainder part of log( f (27)) becomes part of the remainder bound interval of the Taylor model of log( f ) The remainder term can be estimated as 1 (B(P, ,—)+ I, —)‘°+l 1 (“llk+2 k k +1 +1' 1 +1 c...) (1 + [0,1] - (B(p,,—) + Ian/Ger) In a rather similar fashion, it is possible to determine Taylor models of reciprocals, square roots, trigonometric functions and so on as soon as a Taylor model for the argument is known. In the following, we give short recipes to obtain Taylor models for other intrinsic functions. Multiplicative inverse: Under the condition V1.7. 6 [6,5], 0 ¢ Pa,f(f_ £0) + [0J7 1 = L. _f(i’) (f0???)2 _... _ 1(f(f))k} f(i") 6a,; {1 00.! + 031,; +( 1) 65,! +(__1)k+1 (f(f))k+l 1 (4.11) k — k ' 6,372 (1+ 6 - flea/cap) +2 Square root: Under the condition v5 6 [6, 13'], 19,,(5' — :50) + 1,, c (0, 66), +1f(:i:‘ 1(fri" 2 1-1 216—3” m: mlHZ L212? (2)) +'”+(_1) ( k12k) (68.23)} 2Ca’f Ca)! k (Qk—l)” (f(zic‘))"+1 1 +(-—1)\/537'(k+1)!2k+1 65.31 (1+9.f(§.~')/ca,; )k+1/2' 89 Multiplicative inverse of square root: Under the condition we E[(_1°,6], P,,(6 — 6,) + 1,, c (0, 66), 1 _ 1 1R6) 3!! ((6))? (2k -1)!!( (60* f(f) — M.{1_ 2 CQJ +2l22 Ci", + O . . + (—1)k 11312,“ of,” } 1 . (2k +1)!! (f(f))k+1 1 +(_1)k+1 fr...) (k+1)!2‘=+1 c’gi,‘ (1+6-f(:i:’)/ca,f)k+3/2. Sine: sin(f(:1":')) = sin(ca,f) + COS(CQJ) - flf) — 51;Sin(0a,f)‘ (Rf)? _, pose.» - (f(6’))3+ - - . + (k j 1),(f(6))k+1 J, where J _ —Jo if mod(k,4) = 1,2 _ Jo else J _ cos(c,,,f + 6 - f(:2')) if k is even 0 _ sin(ca,f + 6 - f-(a'?)) else. Cosine: ces(f(6)) = ees(c.,,) — enema-1(6) — %cos(cp.f) -(f(6))2 +§ sine...» - (f(6‘))3+ - - - + (,j1),(f(6))k+l 6 where J _ {—Jo if mod(k,4)=0,1 J0 else J _ sin(c,,,f + 6 - f_(:i:')) if k is even 0 — cos(c,,f +6 - f(2’:')) else. 90 Hyperbolic sine: sinh(f(i")) = sinh(ca,f) + cosh(ca,f) - it?) + l sinh(ca,f) - (f—(:E))2 2! 1 (f(f))"+‘ - J, (k+1)! +§ COSh(Ca,I) - (f(:z'))3+ - - - + where J _ COSh(CQ’f + 6 - f(:i:‘)) if k is even _ sinh(ca,f + 6 - f(5:')) else. Hyperbolic cosine: cosh(f(:i’)) = COSh(Ca,f) + sinh(ca,f) - f-(f) + i COSh(Ca,f) ' 6(5))2 2! +1 sinh(ca,f) - (f(f))3+ - - - + 1 (f(i‘))"“ - J, 3! (k +1)! where _ sinh(ca,f + 6 - f:(:i:')) if k is even _ COSh(Ca’f + 6 - f(:i:')) else. Arcsine: Under the condition Vi? e [5, 13'], P0,,(r — so) + 10,, c (—1, 1), using an addition formula for the arcsine, we have arcsin(f(:i:')) = arcsin(ca,f) + arcsin (ft?) - ‘/1 — of” — ca,f-\/1— (f(i"))2) . Utilizing that 9(5) 2 1(1) (/1— c3, - ca) - (/1 — (1(8))? does not have a constant part, we have arcsin(g(a-s)) = g(f)+§(g(i~'))3+-:—,(g(i))5+37,5 (9(5))7+--- (k j 1)! (9(i’))"+1 - arcsinw 1(5)), + 91 where arcsin'(a) : 1/\/1 — a2, arcsin"(a) ——— a/(l—a2)3/2, arcsinl3)(a) = (1+2a2)/(1— a2)5/2, A recursive formula for the higher order derivatives of arcsin arcsinlk+2)(a) = {(2k + 1)a arcsin(k+l)(a) + k2 arcsinlk)(a)} 1 — a2 is useful [58]. Arccosine: Use arccos(f(:ic’)) = 71/2 — arcsin(f(:i:')). Arctangent: Using an addition formula for the arctangent, we have f(i") -Ca.f ). 1+CaJ f(f) arctan(f(:i’)) = arctan(ca,f) + arctan ( Utilizing that f(f)_ca,f : f—(f) 1+ca,f-f(:i:’) 1+ca,f-f(:1':') does not have a constant part, we obtain 9(5) arctan(g(a:*)) -—— 9(a)) — gee)? + goo-2')? — ;(g(:z~'))7 + - . - + fi—lwnk“ -cosk+1(arctan(6-g(:i:')) - sin ((19 +1) - (arctan(6 - g(§:’)) + g)) . Altogether, it is now possible to compute Taylor models along for any function that can be represented in a computer environment along with the mere evaluation of the function by simple operator overloading, in much the same way as the mere computation of derivatives, Taylor polynomials, or interval bounds, along with the mere evaluation of the function. 4.3.3 Derivations and Antiderivations In the spirit of the idea of embedding the elementary operations of addition, multi- plication, and differentiation and their inverses that are defined on the class of C°° functions onto the structure of Taylor Models, we now come to the mapping of the derivation Operation 6 as well as its inverse 8‘1. Similar to the case of the DA, and 92 following one of the main thrusts of the history of differential algebra, we will use these for the solution of the initial value problem (1 - ~ .. are) = F(r(t).t), where F is continuous and bounded. We are interested in both the case of a specific initial condition F0, as well as the case in which the initial condition F0 is a variable, in which case our interest is in the flow of the differential equation 7.3“.) = M(F0, t). As in the case of the conventional DA method, in order to prevent loss of order in the differentiation process, the derivation 6 can be evaluated only in the context of a Lie derivative L9 = g - 3, where g(§:’0) = 0. However, in the case of Taylor models, an additional complication is connected to the fact that from the Taylor model alone, it is impossible to determine bounds for the derivative, since nothing is known about the rate of change of the function (f —Pa,f) within the remainder bounds I f. The situation can be remedied by a further extension of the Taylor model concept to contain not only bounds for the remainder, but also a low-parameter bounding sequence for all the higher derivatives that can occur. In contrast to the derivation 8, Taylor models of its inverse 8‘1 are readily available. Given an n—th order Taylor model (Pn,f,1n,f) of a function f : [6,5] C R” —> R around the reference point 50, we can determine a Taylor model for the indefinite integral 8,?" 1f = f f dz, with respect to the variable 16,-. The Taylor polynomial part is obviously just given by f0“ Pn_1,f(:f:')d:r,-. Since the part of the Taylor polynomial Pm; that is of precise order n is Pm, — Pn._1, f, remainder bounds can be obtained as (B(P,” — Pn_1, f) + Imf) -B(;r,-), where B (22,) is obtained from the range of definition of :1:,- as b,- — a,. We thus define the operator 6,-— 1on the space of Taylor models as 0:1(Pn,faln,f) = (Pn,6-1faIn,8'1f) 93 = ([0 Pn_1,,(5)dx,, (3(me — P,_1,,)+I,,,,)-B(a:,~)) . (4.12) With this definition, bounds for a definite integral over variable 2:,- from 112,) to :1:“, both in [a,-, b,], the domain of validity of the Taylor model of a function, can be obtained as A f (ijda‘i E (Pn,3‘1f(f|$i=$su-$io) _ Pnaa_1f(fl$i:xil_$i0)l In,6‘1f)' (4'13) 1'! 4.4 Examples RDA have many applications, including global optimization, quadrature, and solution of differential equations. We begin our discussion with the determination of sharp bounds for a simple example function using RDA. The sharpness of the resulting bounds will be compared with the results that can be obtained in other ways. Sec- ondly, we show schematically how the method compares with the interval method to obtain bound enclosures of functions in one and two dimensional cases. 4.4.1 A Simple Function The function under consideration is f(:L) = :-+ x. (4.14) For an actual computation, we set the parameter a of (4.3) to a = (71,270, [a,b]) = (3, 2, [1.9, 21]). As in the case of DA, the evaluation begins with the representation of the identity function, expressed in terms of a Taylor polynomial expanded at the reference point. This identity function i has the form i($)=x=xo+(:z:—:ro)=2+(3:—2). 94 Since this representation is exact, the remainder bound interval is [0,0]. Hence, a Taylor model of the identity function i is T0,,- = (2:0 + (a: — 2:0), [0,0]) = (2 + (a: — 2), [0, 0]). The constant part of 2' around 1:0 2 2 is cm,- = 2:0 2 2, and the nonconstant part of z' is 2(23) 2 :r — x0 = :1: — 2. The Taylor model off is T0,? = ((3: — 330)) [0,0]) = ((3: — 2)) [020]) ° The computation of the inverse requires the knowledge of bounds of P03, which here is readily obtained: B(Paj) = B(a: — 2:0) 2 [a — 2:0,b — 2:0] = [——0.1,0.1]. We have furthermore B(Paj) + 10,; = [—0.1,0.1] + [0,0] = [—0.1,0.1]. Using (4.6) and (4.7), we have for the Taylor model of (f)2 TOW = ((2: — 2)2, [0, 0]) . The Taylor model of (f)3 is computed similarly: Toma = ((1: — 2)3, [0, 0]) As can be seen, so far all remainder intervals are of zero size. The first nonzero remainder interval comes from the evaluation of the Taylor remainder term, which is (f(.’1‘))4 1 E (B(Pafi) + [0,?)4 0‘3.) (1 + 9 - Eco/C41)" x8 . (1 + [0,1] . (B(P,,;) + row/1:0)5 (4.15) [0,0.0001] _6 c 4. . - 25-([0.95,1.05])5-[0’ 038x“) 1 As expected, this remainder term is “small of order four”. According to (4.11), the Taylor model of 1/2' is then 1 1 1 .2 1 3 _6 T0,. = (§—¥(x—2)+§(x—2) —?($—2) ,[0,4.038x 10 1), and the remainder interval is indeed still very sharp. Using (4.5), we obtain as the final Taylor model of 1 / z' + 2' T0,?”- = T04, + TC“- = (Pa’%+,~, 10%“) (4.16) 95 Table 4.2: The remainder bound interval 11 ”+1, for various orders; 2:0 = 2, [a, b] = [1.9, 2.1]. Order The remainder bound interval : 0 , 1.4579384 x 10-3 j Z—7.6733603 x 10-5, 7.6733603 x 10-5 j I 0 ,4.0386107 x 10-6 j :—2.1255845 x 10-7 ,2.1255845 x 10—7 j [ 0 ,1.1187287 x 10-8 j j —5.8880459 x 10-10 , 5.8880459 x 10-10j : 0 , 3.0989715 x 10-11 j :—1.6310376 x 10-12 , 1.6310376 x 10-1 j 1 0 , 8.5844087 x 10-1‘11 j :—4.5181098 x 10-15 , 4.5181098 x 10-15j : 0 , 2.3779525 x 10-16j :—1.2515539 x 10-17 , 1.2515539 x 10-17 j i 0 , 6.5871262 x 10-19j :-3.4669085 x 10-2’0 , 3.4669085 x 10-20j ' 0 , 1.8246887 x 10-“1 l—‘F—‘l—‘Hi—‘H L = ((2 + %) + (1 -217) (a: — 2) + .2130; — 2)2 — in: — 2)3, {0,4038 x 10-61) = (2.5 + 075(2- — 2) + 0-125(-’r - 2)? - 0062503 — 2):], [0) 4-038 X 10‘6” Since the polynomial Pa, 1. +17 is monotonically increasing in the domain [a,b] = [1.9, 2.1], the bound interval of the polynomial is B (Pa,%+,.) = [Pa,%+,-(—O.1),P 1 (0.1)] = [2.42631,2.57618]. 0,744: The width of the bound interval of the Taylor polynomial is 0.14987, and the width of the remainder bound interval is 4.038 x 10‘6 in the third-order Taylor model eval- uation; thus the remainder part is just a minor addition. The size of the remainder bounds depends strongly on the order and decreases quickly with order. Table 4.2 shows the remainder bound interval for various orders in the Taylor model computa- tion. 96 Table 4.3: The width of the bound interval of f (:1:) = 1/2: + :1: by various methods; 2:0 2 2, [a,b] = [19,21]. Method Width of Bound Interval Intervals 71,, = 1 025012531 nd 2 10 0._15993589 nd = 102 0._15088206 nd = 103 0.1_49_97543 nd = 104 014988476 nd = 105 0.14987569 71,; = 106 014987478 nd = 107 014987469 nd = 108 014987468 Taylor models lst order 045145793 2nd order 0._15015346 3rd order 0.14987903 4th order 0.14987542 5th order 014987469 6th order 0.14987468 Exact 0.14987468 The Taylor model computation is assessed by noting the bound interval B of the original function (4.14), which is 1 1 1 B (E + :1:) = [E + a, E + b] = [2.42631, 2.57619]. It is illuminating to compare the sharpness of the bounding of the function with the sharpness that can be obtained from conventional interval methods. Evaluating the function with just one interval yields ;+[a,b] [a, bl The width of the bound interval obtained by interval arithmetic is 0.25012, and so this 1 = [79—217 + [19,21] <_: [2.37619,2.62631]. simple example already shows a noticeable blow-up. By dividing the domain interval into many subintervals, the blow-up can be suppressed substantially. However, to achieve the sharpness of the third-order Taylor model, the domain has to be split 97 into about 24, 000 subintervals. Table 4.3 shows a comparison Of the widths of bound interval for the exact value, the method of Taylor models, and the divided interval method, where 72,, indicates the number Of division of the domain interval. Of course, sophisticated interval optimization methods [33] [34] [41] [43] can find sharp bounds for the function using substantially fewer interval evaluations. Practically more important are optimization problems in several variables, and in this case, the situation becomes more dramatic. We wish first to illustrate the computational effort necessary for an accurate calculation Of the result by estimating the required number of floating-point operations. We use a simple example function of six variables such as f (:1?) = 2,-(1 / :1:,- +5121) to get a rough idea Of the computational expense in the case of functions Of many variables. In the one-dimensional case, one interval calculation 1/[a, b] + [a, b] requires two additions and two divisions. TO com- pare with the third-order Taylor model computation, we divide the domain into 104 subintervals, on which additions and divisions total ~105 floating-point operations. Thus, in the multidimensional case with six independent variables, the number of floating-point Operations explodes to (104)6 x (~10) 2 ~10”. Again, sophisticated interval Optimization methods will be more favorable than these numbers suggest, but typically there is still a very noticeable growth of complexity. TO estimate the performance of the Taylor model approach, we note that the one- dimensional Taylor model in the third-order computation involves a total of about 35 additions, multiplications, and divisions, as counted in (4.15) and (4.16). As we use more variables, however, the total number Of terms in the polynomial grows only modestly. For example, order three in six variables requires only a total of 84 terms. Thus in total, the number of floating-point Operations Of the third-order Taylor model is ~104. A summary of the number of floating-point Operations is given in Table 4.4. 98 Table 4.4: The total number of PP Operations required to bound a simple function like f(f) = 21(1/1133' + 233'). One Dimensional Six Dimensional Interval ~ 10 ~ 10 104 divided intervals ~ 105 ~ 1025 3rd order Taylor model ~ 10 ~ 104 4.4.2 Bound Enclosures of Functions In this subsection, we use some simple functions in one dimension and two dimensions to show schematically how the RDA method bounds functions in comparison with the interval method. The first function is a one dimensional function f(:z:) = :r(:1: —1.1)(:1: + 2)(:1: + 2.2)(2: + 2.5)(1‘ + 3) -sin(1.7:z: + 0.5). Figure 4.2 shows the function and its bound enclosures in the domain [—0.5,1.0]. The interval method is applied to bound the function using smaller domain intervals divided into 25 subintervals and 50 subintervals. The method of Taylor models com- putes the polynomial part Of the function and the remainder bound interval. The pictures show the bands of the enclosures Of the function around its polynomial parts by the remainder bounds using 7th order and 8th order Taylor models. A similar schematical comparison between the interval method and the method Of Taylor models is made for a two dimensional function. We worked on a function f(.r, y) = sin(1.7:1: + 0.5) - (y + 2) -sin(1.5y) in the domain {—1, 1] x [—1, 1]. Figure 4.3 shows the function enclosures by the interval method with the smaller domain intervals divided into 10 x 10, 20 x 20, 40 x 40 and 80 x 80. Figure 4.4 shows the function enclosures around its polynomial parts by the remainder bounds using 7th, 8th, 9th and 10th order Taylor models. Even in this 99 modest case of only two dimensions, the Taylor model approach requires much less effort to provide a similar level Of sharpness; the 1600 subintervals used to include the function are in contrast to only 66 expansion coefficients, plus one remainder bound interval. As dimension is increased, the number of subintervals necessary to provide an accurate modeling increases dramatically, while the number Of Taylor coefficients grows much more slowly. 100 Figure 4.2: One dimensional function and its bound enclosures. From top to bottom right: the function, the bounds by the interval method with the 25 and 50 divided domain intervals, the bounds by the 7th and 8th order Taylor models. 101 V ..“"’.”.‘3. v v c; . ~ 0 .‘o‘.¢‘.¢‘. ‘d‘vwfl . a“. “W- -.‘. ‘. .3. ‘ 'éamergeego - ‘ ' ¢ 0‘ - .o o ‘o-.H‘.- ‘ ’ .. «e»... b - .o.-.¢-:o‘:‘“ 6"” «.039»... .* - " ‘.‘¢"

n By denoting the remainder bound intervals [if and 1,139 by I?“ and 1;“, we can simplify the above expression (5.1) as R n+1 j k In,” g 2 1,1,. (5.2) j,k=0 j+k>n Practically, considering the semi-distributibity of interval arithmetic (4.2), we use the following estimate R n+1 i n+1 j n+1 i j In.” 9 230 If' .20 I, r) :30 19' Z If - (5-3) 2: j = 1,: i+j >11- i+j>n 114 The square Of a Taylor model TnJ = (Pn,f,1,fff) with order bound intervals 1?, I},..., I? can be computed as follows. Again, we denote the remainder bound interval [if by I?“ for the notational convenience. The order bound intervals of the resulting squared Taylor model Tng is D1» If2g2. i. i+J ,- ‘ 1M2 2 ifkis even else H k- H ”I ‘ 0 i V <1 and the resulting remainder bound interval is n+1 . 2—1 . n+1 . 2 H 2 J 2 1,,fzg2- £30 1,- go I, + £0 (1,) . i+j>n 2i>n Here, we can take advantage of explicit expressions for the square in interval arith- metic; for an interval 1 = [a, b] 12— 0 , max(a2,b2) ] ifabgO 54 _ [ min(a2,b2) , max(aZ,b2) ] else - (-) 5.3.3 Intrinsic Functions The actual computation of intrinsic functions Of Taylor models is performed in two steps. Let us take the exponential function as an example. Firstly, we prepare a Taylor model of flu?) :- f (if) — cm, from the Taylor model of f (f) This can be done trivially according to subsection 5.3.1, and thus we have a Taylor model of 1(6), Tm!- = (me, 153]) with order bound intervals I}, I},..., 1}}. As described in subsection 4.3.2, the exponential of a function f (63') can be expressed as exp(f(i’)) = €XP(Cn,f)°{1+f(6)+—21—!(f(f))2+m+ (n j 1)! (f(f))"+1 exp (6 - 1(8)). + exp(c,,,f) - ‘l" ' "'7 115 Then, we compute the polynomial Of f—(2'7') exp(c.,)) - {1+ f0?) + 2, (11(5))? + - - - + n, 0(a)") (5.5) according to the recipe for addition and multiplication Of Taylor models. Practically, we use Horner’s scheme (4.8) to make the computation efficient. The polynomial part of the resulting Taylor model, PM,p f, is computed in this process, and we denote the remainder bound interval of the resulting Taylor model of the polynomial (5.5) by Ianoly‘ As the second process, we compute the remainder bound interval of the resulting final Taylor model, Ifexpf, using (4.10), as 1 n+1 171:0”)!- Q 1fipoly+exp(cn,f)-m(B(Pn,f-)+ 1:11) exp ([0,1] ' (B(PnJ) + 1:») , where B (Pu, f) is readily Obtained as 2&0 I}. The other intrinsic functions Of a Taylor model including the multiplicative reciprocal are estimated in a similar way. 5.3.4 Integral Assume we have a Taylor model Of a function f (if), Tm, = (Pn,f,1,ff) with order bound intervals 1?, 1},..., 1?. To obtain a Taylor model for the indefinite integral 6,71 f = f f 622,- with respect to variable x,- on the Taylor model T n, f, we just follow the description in subsection 4.3.3. The polynomial part Of the resulting Taylor model Pn,3"1f is Obtained in the framework of DA. The remainder bound interval Of the resulting Taylor model 156.” is estimated as 1384, g (B(P” — Pn_1,,)+ 1,15,) - 3(a) = (1}: + 1,15,.) - (b,- — 6,). In the process of performing the indefinite integral of polynomial part, the information on the order bound intervals are lost. To reconstruct the order bound intervals Of the resulting Taylor model, 13-1f, 1,1,_1f,..., 13-1], we use a method to estimate order bound intervals discussed in the next section. 2" I rn r-.‘. J . v 116 5.4 Methods to Tighten Order Bound Intervals The preceding section discussed the algorithms to compute Taylor models by com- bining DA techniques and interval arithmetic. While the polynomial part of the Taylor model computation is performed precisely using DA, the width Of the remain- der bound intervals depends on the method to estimate the bounds. For example, we saw that the estimate Of the remainder bound interval of a Taylor model as a product of two Taylor models discussed in subsection 5.3.2 could be done in three ways, (5.1), (5.2) and (5.3), and the last way gives the narrowest width for the resulting remainder bound interval 15,”. However, still a naive application may not result in optimally sharp bounds, and it is useful to study the question of polynomial bounding. In this section, we discuss several already implemented and working algorithms to bound higher order multidimensional polynomials as well as some future ideas to narrow the bounds further. 5.4.1 Implemented Methods Currently, we have five ways to tighten order bound intervals of a Taylor model, beside an option to do nothing. The number of algorithm is switched freely using the intrinsic service procedure RDITIG in section 5.1. These algorithms are 0: NO further tightening is done. 1: Default method. Each monomial Of the polynomial part Pm; is estimated us— ing interval arithmetic including an efficient estimate of interval powers, and summed to each order bound interval. 2: For each polynomial Of the exact order as a part Of the polynomial me, utilize the existing multidimensional polynomial evaluator in interval arithmetic. 117 3: Same as 2, but using Horner’s scheme in the multidimensional polynomial eval- uator. 4: Each monomial of the polynomial part Pm, is estimated using a real number scanning technique. Obviously this method is not rigorous and thus used for test purposes only. It furthermore yields widths that are narrower than the sharpest true bounds. 5: Same as 4, but for the purpose to compensate the underestimation, a correction factor is imposed. The default algorithm 1 is used in all the Taylor model computations except for this section, and it is showing the practical strength of the RDA technique. The method to evaluate the powers of intervals is the same with the squaring method (5.4). It is most efficient when the expansion point 60 is positioned at the center of the domain [6, 0], since [6 —- :60, b — 60] is used as an argument for interval arithmetic. The existing evaluator Of multidimensional polynomials in COSY INFINITY, the intrinsic procedure POLVAL, accepts all arithmetic data types as arguments. The current interval evaluation in POLVAL does not use optimally efficient way of com- puting powers of intervals, and this presently limits the practical strength of the method of the algorithms 2 and 3. The method would be more useful in the future when another kind of order bounds information is needed, for example 15", bounds of polynomial Of the order k. The interval based algorithms 1, 2 and 3 can give a mathematically rigorous estimate, but the blow-up problem, a crucial problem of interval arithmetic, again limits the effort to tighten the bounds of a polynomial. In practice this is Often not so significant because the remainder is by nature much smaller than the original function. In addition, there is another mechanism to get an estimate on bounds of polynomials 118 by a scanning in real numbers. Scanning can give a good feeling on the behaviour of a function, however the estimate is not rigorous. And it is even dangerous when the function has many local minima and maxima in the domain and when we cannot have enough sample points to scan especially in high dimensional cases. Knowing this mathematical limit, we can use this approach to check the efficiency Of Taylor model arithmetic, and it can supply a good measure of the size of order bound intervals. For the real number scanning algorithms 4 and 5, the whole domain [6, 6] is covered by the total of NM,n equidistant points. Let m be the number Of points in each dimension. Then the total number of points Nscan is calculated as m”, where v is the number of dimensions. The number m is determined to be the maximum integer which satisfies the condition m” g Ntota), where NM“) can be specified using an intrinsic service procedure RDNPNT in section 5.1. The default value of Ntom; is 1000. To increase the performance on computers regarding CPU time, we utilize an array of VE data type in COSY INFINITY. The 2th array element Of the VE array is filled with Nscan coordinate values Of the 2th independent variable sci, which are the m“—1 times iteration Of m equidistant values in [a,, 0,] including a, and 0,, and as a total, the VE array with 2) elements can represent the whole equidistant points NW,n in the whole domain. We use a correction factor 1 + l/m in the algorithm 5. After getting the estimate Of order bound intervals using the VE array described above, shift the centers of the bound intervals to 0, then multiply with the correction factor, and shift back to the original centers. 5.4.2 Examples Of Performance In this subsection, some examples are presented to compare the performance Of the al- gorithms in the previous subsection. The first example is a three dimensional example .i1_‘ i 119 function treated in 5th order Taylor model arithmetic. The function 4 tan(3:cg) f(.’I?1,.’L'2,IL'3) = —120—2£E1—7$3(1+2III2) 6131 3$1+$1 m _ 1_ 6 3 +132 332 )+(332 ) — ' h (0.5 s1n + 82:2 + 7 31:3 51:1 tanh(0.9:1:3) V 5232 is computed around the reference point :60 = (2, 1, 1) in the domain [61,01] x [a2, b2] x —20:z:3(2:1:3 — 5) + —- 202:2 sin(3:r3) (5.6) [(23, b3] 2 [1.9,2.1]x[0.9,1.1]x[0.9,1.1]. The following is the resulting order bound and remainder bound intervals, and the resulting total bound intervals for each tightening algorithm. 0: No tightening. Order Bound interval 0 [-.3928616701165386 ,-.3928616701165386 I 1 [-32.85040954846235 , 32.85040954846235 ] 2 [*2.611562619634850 , 2.641562619634851 ] 3 [-.2545082014874667 ,0.2545082014874666 ] 4 [-.2634629627970872E-01,0.2634629627970871E-01] 5 [-.2739252867613695E-02,0.2739252867613695E-02] :Renuunder [-.3160819763199970E-03,0.5077650001011438E-03] jknal [-36.13874367082485 , 35.38321201361556 ] 1: Default. Order Bound interval 0 [-.3928616701165386 ,-.3928616701165386 ] 1 [-4.036832719935209 , 4.036832719935209 ] 2 [-.1517612694112098 ,0.7550980336507371 ] 3 [-.2288204979419630 ,0.2288204979419625 1 4 [-.1303457714927084E-01,0.2049495829307956E-01] 5 I‘.9783281805866122E-03,0.9783281805866111E-03] :Renunnder [-.1772503727079576E-03,0.3571133959329409E-03] Tbtal [‘4.824466313107486 , 4.649719981280969 ] 120 2, 3: Using the polynomial evaluator POLVAL. Order Bound interval 0 ['.3928616701165386 ,-.3928616701165386 ] 1 [‘4.036832719935209 , 4.036832719935209 I 2 [-.7553220623800865 ,0.7553220623800874 ] 3 [-.2288204979419630 ,0.2288204979419629 ] 4 [-.20883689423264833-01,0.2088368942326483E-01] 5 [-.9783281805866133E-03,0.97832818058661305-03] Renufinder [-.1932338291952928E-03,0.3848962108177413E-03] Tbtal [’5.435892201806844 , 4.650360523955390 ] 4: Real number scanning. ! 1 Order Bound interval 0 ['.3928616701165386 ,-.3928616701165386 ] 1 [-4.036832719935209 , 4.036832719935209 ] 2 I“.7674533702027161E-03,0.7352759591802416 ] 3 [-.2275570744837629 ,0.2275570744837623 ] 4 [-.5020266876434484E-02,0.20043055989984293-01] 5 [-.9628726613536728E-03,0.9628726613536709E-03] .Renufinder [-.1770676431478265E-03,0.3568468875725592E-03] 'Tbtal [-4.664179125086649 , 4.628166859021585 ] 5: Real number scanning with a correction factor. Order Bound interval 0 [-.3928616701165386 ,-.3928616701165386 ] 1 [-4.440515991928730 , 4.440515991928730 ] 2 [-.3756962399772495E-01,0.7720781298077638 ] 3 [-.2503127819321392 ,0.2503127819321386 ] 4 [-.6273433019755423E-02,0.2129622213330522E-01] 5 [-.10591599274890405-02,0.1059159927489038E-02] :Renufinder [-.2095690927896201E~03,0.3882500914646420E-03] Tkmal ['5.128802230015166 , 5.092788865704353 I The size of the remainder bound intervals are in the same range in all the algo- rithms even without any tightening (0). However, the total bounds estimate without tightening shows a big blow up. The function evaluation in interval arithmetic gives the resulting bounds 121 [-31.84061451917199 , 33 . 73372787735954 I , so the necessity of any tightening mechanism is clear. Other than that, all the tight- ening algorithms are working more or less well. The default algorithm 1 is working well compared to the uncorrected scanning algorithm 4. TO see the performance Of the total bounds, the function evaluation in real number scanning using the same scanning condition with the algorithm 4, gives a good idea, which is [-4.129818215892230 , 4.397287227019497 ]. The estimate Of the total bounds in Taylor model arithmetic is comparable to the bounds estimate Of the function in real number scanning. The second example shows the order bounds estimates of a multidimensional polynomial with randomly created coefficients. The polynomial is 5th order in four dimensions, and about half Of its coefficients are nonzero and range from —1 to 1; P5(:1:1,:1:2,:1:3,:r4) = 0.4007583772763610271 + 0.96514073479920634:2 —0.05361263593658805x4 — 0.984300413983874021222 +0.213129849638789933324 + 0.939666689839214123 —0.29482657765038312::1‘ + - -- -0.03892374457791448a:§x§ — 0.9986851064022630233232. The estimate on order bound intervals is made around the reference point :60 = 0 in the domain [6, 5] = [411,071]. The following are the resulting order bounds and total bounds for each tightening algorithm. 122 1: Default. Order Bound interval 1 [-.1419511748012155 ,0.1419511748012155 ] 2 [-.11974302636226648-01,0.2137096953461878E-01] 3 [-.7059430076740683E-02,0.7059430076740683E-02] 4 [-.8311970126233067E-03,0.78475300251739138-03] 5 [-.1136391405879840E-03,0.1136391405879840E-03] Tkmal [-.1619297436673942 ,0.1712799665556804 ] 2: Using the polynomial evaluator POLVAL. Order Bound interval ,, 1 [-.1419511748012155 ,0.1419511748012155 1 2 [-.2137096953461878E-01,0.2137096953461878E-01] 3 [-.7059430076740683E-02,0.7059430076740683E-02] 4 [-.8483902763924564E-03,0.8483902763924564E-03] 5 [-.1136391405879840E-03,0.11363914058798403-03] Tbtal [-.1713436038295554 ,0.1713436038295554 ] .2. in '7' h“)- 3: Using the polynomial evaluator POLVAL with Horner’s scheme. Order Bound interval 1 [-.1419511748012155 ,0.1419511748012155 ] 2 [-.2137096953461878E-01,0.2137096953461878E-01] 3 [-.7059430076740683E-02,0.7059430076740683E-02] 4 [-.8483902763924563E-03,0.8483902763924563E-03] 5 [-.1136391405879840E-03,0.1136391405879840E-03] Tbtal [-.1713436038295554 ,0.1713436038295554 J 4: Real number scanning. Order Bound interval 1 [-.1419511748012155 ,0.1419511748012155 1 2 [-.9843004139838742E-02,0.2137096953461878E-01] 3 [-.3913808809127659E-02,0.3913808809127659E-02] 4 [-.2922135995118879E-03,0.4383278226363474E-03] 5 [-.3743960746774976E-O4,0.3743960746774976E-04] Tbtal [-.1560376409571615 ,0.1677117205750661 ] 123 5: Real number scanning with a correction factor. Order Bound interval 1 I-. 1703414097614586 , O . 1703414097614586 J 2 I- . 1296440150728449E-01 , O . 244923669020645413-01] 3 [- . 4696570570953191E-02 , O . 4696570570953191E-02] 4 I:- . 36526774172671 15E-03 , 0 . 51 1381964851 171015-03] 5 I- . 4492752896129971E-04 , O . 4492752896129971E-04] Total I- . 1884125771 103843 , 0.2000866567282888 ] Again, the default algorithm 1 is working well and yields a result that is comparable to the uncorrected scanning algorithm 4. 5.4.3 Further Tightening Methods As discussed above, to maximize the sharpness of the remainder bounds it is useful to employ efficient polynomial bounders, and in this subsection, some algorithms for this purpose are presented. First we observe that if d is the width of the domain interval, then very roughly we expect the widths Of the 1’“ to scale like width(1") z d". Hence the dominating part of the total bounds of a polynomial comes from the linear part and the successive lower order parts. The linear part is easily bounded, and so it is useful to derive methods to bound the second and third order parts. Let us begin the discussion with a second order bounder, which can give an exact bounds for 152. At first, let us bound a one dimensional function F(CL‘) = (311332 + 20113 “I" C0 in the domain [a, 0]; while it is trivial, it is the basics Of the whole effort to bound a multidimensional quadratic function. If cm 2 0, the function is linear, so B(F) = [min(F(a), F(b)), max(F(a), F(b))]. (5.7) 124 If 611 76 0, by writing the function as 2 2 F(;r) =611 (2+2) — 3“+00, 011 011 if—cl/c11€[a,b], c2 (:2 B(F) = [min (Fm), F(b), —C—1- + co) ,max _ 2J( 1132 > + 2 ( C2 — 0, (58) where * 011 012 “ 2 J = < , and [J] = C11C22 — 612. 012 6'22 If |.1| 35 0, there is only one stationary point 331.9 __ _J‘_1 Ci _ __1_ C2201“ C1202 $23 C2 |J| 61102 — 61261 125 If [J] = 0, then 011022 = 0%, thus cu 7f 0 or 022 aé 0 must hold for the function to be Of second order. The equations of stationary points (5.8) have solutions if 02201 — 01202 = 0 and 01102 — 01201 = 0. Now, if (311 = 0, then it follows that 022 75 0, 012 = 0 and 01 = 0, so 0 2 (:2 0613111172) = 6223761" 202172 + Co = C22 (332 + ‘3) — —2 + Co; C22 C22 hence the stationary points lay on 2:23 2 —02/022. Similarly, if 022 = 0, the stationary points lay on 2:1, 2 —cl/cll. If cu # 0 and C22 76 0, the two equations in (5.8) are identical, and the stationary points lay on them, namely 6 :1: , + c 02 02 $25 = ——13L—2—, and G(:1:1s,a:23) = ——1 + co = ——2— + co. C22 611 C22 These stationary points are inside the domain if . ( Ci2ai+C2 C12bi+C2) ( 612014-02 Cl2b1+02)] min ———,——— ,max —— —— 7150. [0’21 b2] fl 7 C22 C22 622 C22 Using induction, the bounding of quadratic functions in more than two variables can be reduced to the above two variable case. Consider a 2) dimensional function _ 2 2 2 H(.2:1,a:2, - - -,:1:,,) — 611231 + 2612;171:122 + ' - 3 + 201051315131, + c2222 + - ' - + 6,,va +2c1231+~-+ 2cvxv + co in the domain [(11,111] x x [a,,bv]. If cu = 2 cm, 2 0, the bounds of H is determined by the maximum and the minimum of the set {H($lba ' ' 'avaHIBlb E {01,b1},‘ ° 'ava E {avabv}}a which has 2” elements. If at least one Of cij is nonzero, the maximum and the minimum are on the 22) boundaries, 271 = 61, 2:1 2 bl, ~ - -, 23,, 2 av, 18,, = by, or at the stationary r’J 126 points inside the domain. To find the extrema on each boundary reduces tO the bounding Of quadratic functions in v — 1 dimensions. On the other hand, the inner stationary points satisfy “31 C1 011 ‘ ‘ ' Clv VH=2J +2 3 =0, where]: , xv C‘U Cl?) . . . va so the stationary points 63', are found as follows. If If] 75 0, then 12“,, = —|/1,~] / [J | where A,- is Obtained by replacing the 2th column of .1 by (01,: - -,c,,). If it happens that [1] = 0, the task to find if, is quite complicated for v 2 3. and usually not possible in closed form. The implementation Of the multidimensional quadratic function bounder requires some careful consideration because of the increase of the effort with dimension. As- sume we have an explicit k-dimensional bounder. Then the generation Of a k + 1 dimensional bounder requires 216 applications Of the k dimensional bounder for the boundaries in addition to the study of internal stationary points. To Obtain a 2) di- mensional bounder from a k dimensional one thus requires 2) — k internal searches, and 212 - 2(k + 1) - - 22) applications of the k-dimensional bounder. Table 5.5 shows the number of applications Of a two, three or four dimensional bounder at the boundaries to bound a higher dimensional quadratic function. The implementation of the three dimensional bounder looks most useful. “hung-“2.11: * _ ' I | . 127 Table 5.5: Number Of applications Of the lower dimensional quadratic function bound- ers to bound a higher dimensional quadratic function at the boundaries. .fi um—Wfl l , Target by by by dimension two dimensional three dimensional four dimensional v bounder bounder bounder 3 6 1 4 48 8 1 5 480 80 10 6 5 , 760 960 120 7 8O , 640 13 , 440 1 , 680 8 1,290,240 215,040 26,880 Formula 2""32)! 2”“4211/3 2”‘7vl/3 128 +b_3 +b 2 Co+a_1 +a _2 +21 _3 Domain [a,b] Figure 5.1: The order bounds 12 = [a _2, b _2] and 13 = [a _3, b_3] around the linear line Co + c_12: E [Co + a_1,c0 + b_1]. We can pursue the effort to tighten a polynomial bounds further. Another method is to refine bounds of a third order polynomial 153 iteratively using the order bounds 1°, 11, 12, and 13. The method uses the fact that the linear part dominates the behavior of the function. Let us bound a one dimensional third order polynomial P(:r) = Co + c_12: + c322 + C_3.’L‘3 (5.9) in the domain [a, b]. Assume that we have the order bounds 10 = [c0,c0], 11=[a_1,b_1],12 = [a_2,b_2] and 13 = [a_3,b_3]. From these, we can estimate 153 simply I53 C_: [co+a_1+a_2+a_3, co+b_1+b_2+b_3], but a careful look at Figure 5.1, which shows how the order bounds 12 and 13 add the width of the band around the linear line co + c_1:r E [Co + a_1,co + b_1], gives 129 / / I /,’ / / /,’ /’I 2’ / ,' ,I I Co-l-b_1+a_2+a_3 /” ’I// / I - , ,I , / Co+b_l+a_2+a_3 -(b _2 +b_3) / I ” / /” I / / ’ ” 4’ //’ I I / I ’ / //I ’I/ / I ,’ ’ r’ ’ / I’ / , / I / II, I / z x“) Figure 5.2: An upper bound of 1S3 is above on + b_1 + a._2 + a_3. Find a point 23(1) to specify a new domain. an insight to tighten 153. Namely, an upper bound of 153 certainly must lie above co + b_1 + 0._2 + a_3. Now, find a point mm in the domain such that 004—642(1)=co+b_1+a_2+a_3-(b_2+b_3). (5.10) For the sake of simplicity of the argument, assume c_1 > 0 in the following. The new domain [23(1), b] contains the point which gives the maximum of the polynomial (5.9). Now, compute the order bounds 131) = [692,632] and [(31) = [69,59] in the new domain. Since the new domain is smaller than the original domain, the new order bounds [(21) and [(31) are narrower than the original order bounds I2 and 13 respectively. Thus, new bounds [(51]; computed from [(21) and [(31), has a smaller upper bound co + b_1 + 5‘,” + 58,). This procedure can be applied iteratively until the desired sharpness Of the upper bound of 15:)3 is achieved. The lower bound Of 153 is refined in the same way. The key in multidimensional case is how to define a new domain. Pick up the 130 boundary corner point A which gives the upper bound of 11, search the points which satisfy a condition similar to (5.10) at each neighbouring boundary of A, and these points and the point A define a new domain. As an example, let us bound P(a:) = :1: — 2:2 — 51:3 in the domain [—0.1,0.1]. The order bounds are 1° = [0,0],11 = [—0.1,0.1], 12 = [—0.01,0] and 13 = [—0.001,0.001]. From these, 153 Q [—0.111,0.101]. Find mm from the condition (5.10) to be 0.088, so the first new domain is [0.088, 0.1]. In this new domain, compute the order bounds 1,2,, = [~0.01,-7.744 x 10-3] and 13,, = [~0.001,—6.81472 x 10-4] to get 15,? g [0.077, 0.091574528]. Finding 112(2) gives the second new domain [0.0974, 0.1]. The new estimate Of the order bounds [(22) and 132) gives 15?. Altogether, after two iterations, the upper bound of 1S3 is found to be 0.089589229, which is much smaller than 0.101. The methods discussed in this subsection will help to get exact bounds of a second order polynomial [S2 and sharper bounds of third and also fourth order polynomials 153 and 154, and it will tighten the total bounds of a polynomial B (P) significantly. However, as discussed in subsection 5.3.2, the carriage of the order intervals 1°, 11, 12, ~ - ~, I” as a part Of Taylor model data on computers is indispensable, unless we have a direct mechanism to compute the remainder bound interval of a product of two Taylor models. One also should try a path in this direction, where the problem is cast into the following form; find bounds of a function F(:1:, 3,2) = P,,,.,(:3 — :80) + Pn,f(6 — :30) - t + P,,,(* — :30) . s + s .2 v.3 e [51, 13'], V8 6 1,5, and Vt 6 1,539. 5.5 Examples of Computation In this section, we present some examples and applications to show its practical power. First, we make an inspection of the method using the multidimensional function (5.6) 131 on page 119. Second, we apply the method to compute multidimensional definite integrals. The last example in this section bounds a normal form deviation function discussed in chapter 1, which is a six dimensional polynomial of roughly 200th order, and hence has many local extrema. 5.5.1 A Small Multidimensional Function A three dimensional function (5.6) on page 119 is computed in Taylor model arith- metic around the same reference point as in subsection 5.4.2, 550 = (2,1,1). The do- main is chosen centered at :30 such that the width is 0.1 in all the dimensions, namely [(11,121] = [195,205], [(12,122] = [095,105], [a3,b3] = [095,105]. The bounds estimate of the function by a scanning in real numbers with the same scanning condition with subsection 5.4.2 is [ -2.31165715, 1.78168226]. When the function is evaluated with the domain intervals in interval arithmetic, it gives the bounds [- 16 . 36393303 , 16 . 09747985] , which is almost ten times wider than the bounds evaluated by the scanning. Table 5.6 is the summary of the remainder bound intervals of the function via Taylor model computation in various orders. The table also lists the maximum number of terms of polynomial in each order. While the number of terms of polynomial increases moderately in order, the width of remainder bound interval drops down as expected from the study in the previous chapter 4. 132 Table 5.6: Remainder bound intervals of a small multidimensional function for various orders. [ Order] terms [ Remainder bound interval ] 1 4 [-.3914054034075695 ,0.7252479186770013 1 2 10 ['.3395018823172723E-01,0.3394051630619547E-01] 3 20 [-.1020280049382976E-02,0.16096620942791573-02] 4 35 [‘.8413202994873543E-04,0.8402878906072845E-04] 5 56 [-.2410738165297321E-05,0.4383393685666119E-05] 6 84 [-.33555368946081235-06,0.33431706247521083—06] 7 120 [’.16319419092290073-07,0.2051811636562382E-07] 8 165 [-.2424624748948692E-08,0.24107813420281415-08] 9 220 ['.1721943529722954E-09,0.1736796027250248E-09] 10 286 [‘.2313819358227192E-10,0.2298654603065415E-10] 11 364 [‘.1928098381883883E-11,0.1821047609084028E-11] 12 455 [-.2424312972351917E-12,0.2407755156049230E“12J 13 560 [-.2163485755166328E-13,0.2012630346355465E-13] 14 680 [-.2614793046930413E'14,0.2596691498167008E-14] 15 816 [-.2417223411345871E-15,0.2242821131576220E-15] 133 5.5.2 Multidimensional Integrals We can use antiderivations of Taylor models to compute bounds for multidimensional definite integrals that cannot be treated analytically. The interval method allows us to get a verified estimate of a definite integral. Figure 5.3 shows schematically how to obtain a definite integral of a one dimensional function f in the domain [a, b]. Interval arithmetic can give a bounds estimate of the function in the domain as [fL, fly], so the bounds of the definite integral is estimated as [fL, fU] ~ B (3:) = [fL, fU] - (b— a). On the other hand, Taylor model arithmetic would give a Taylor model of the the function (Pf,1 f) in the domain, from which we can get a Taylor model for the indefinite integral (Pa—1f, 15.1!) as discussed in subsection 4.3.3. And the bounds of the definite integral can be computed using (4.13) on page 93 as P5—1f(b) — Pa_1f(a) + 1511,. fu f1. Figure 5.3: Bounds estimates of a definite integral of a function f by the interval method (left) and the RDA method (right). In the following example calculation based on the RDA method, we study the following analytical double definite integral, which is found in section 4.621 of Grad- shteyn and Ryzhik [30]. .1 1 sin 1 — k2 sin2 :z:sin2 [’2] y‘/ ydzr dy= 7‘ (5.11) O 0 1— kzsin2y 2\/1— 18' The definite integral for k2 = 0.1 is 1.655764710966017. 134 The interval method gives an estimate of the bounds covering the domain with one interval box [0, 7r/2] X [0,7r/2] as [ 0.1433315723E-15, 2.741556778 ], which is too wide to be useful. To increase the sharpness of the bounds, we divided the whole domain into many smaller interval boxes and we obtained the following result shown in 5.7, where the case of 100 X 100 smaller domain interval boxes gives a reasonable estimate. Table 5.7: Bounds estimates of a two dimensional integral with the interval method. Interval Method domain interval boxes bounds estimate 1 [ 0.1433315723E-15 , 2 . 741556778 ] 10 X 10 [ 1.511886201 , 1 .793726593 ] 100 X 100 [ 1 . 641637834 , 1 . 669832424 ] The bounds estimates with Taylor model arithmetic are much sharper, even with- out any division of the domain. Table 5.8 shows a clear superiority of the Taylor model approach. Table 5.8: Bounds estimates of a two dimensional integral with the RDA method. RDA Method order domain interval boxes bounds estimate 1 [ 1 . 445463096 , 1 . 865215766 ] 5 2 X 2 [ 1 . 651327785 , 1 . 660154065 ] 4 X 4 [ 1 . 655670149 , 1 . 655858472 ] 1 I: 1 . 649379060 , 1 . 662771976 ] 10 2 X 2 [ 1 . 655760180 , 1 . 655769308 1 4 X 4 I: 1 . 655764708 , 1 . 655764714 1 135 The situation becomes more dramatic in higher dimensions. Since it is very dif- ficult to find any suitable formula for such an example in books, we extended the definite integral (5.11) to the four dimensional case. /§/§/-§/% (sin y\/1 — k2 sin2xsin2y + sin w\/1 — k2 sin2zsin2w o o o o 1— k2 sin2y 1— k2 sin2w ) drr dy dz dw 71.3 4\/1 -k2' The definite integral for k2 = 0.1 is 8.170871339259325. Similar to before, computa- tions are made to obtain the bounds with the interval method and the RDA method as shown in Table 5.9 and 5.10. Table 5.9: Bounds estimates of a four dimensional integral with the interval method. Interval Method domain interval boxes bounds estimate 1 [ o 7073129582E-15, 13.52904042 J 44 [ 6.344461569 , 9.814754327 J 164 [ 7.730435550 , 8.599902240 J Table 5.10: Bounds estimates of a four dimensional integral with the RDA method. RDA Method order domain interval boxes bounds estimate 1 [ 7.133074468 , 9.204470869 J 5 24 [ 8.148975985 , 8.192531932 J 44 [ 8.170404693 , 8.171334032 J 1 E 8.139359412 , 8.205450807 J 10 2‘r [ 8.170848978 , 8.170894025 J 44 [ 8.170871325 , 8.170871354 J '4 mil-.1 \ 136 5.5.3 Bounds of Normal Form Deviation Function The last example in this section is a bounds estimate of a deviation function from a normal form invariance in six dimensions discussed in chapter 1. The function is a six dimensional polynomial up to roughly 200th order which involves a large number of local minima and maxima and has many cancellations, thus representing a substantial challenge for interval methods. We made the estimate in the domain around the reference point 0'50 2 0.05 with the width in each dimension 0.02 [004,006] X [004,006] X [004,006] X [004,006] X [004,006] X [004,006]. The value of the function at the reference point is f(§:‘0) = 0.6976700784514303 x 10-5, and a bounds estimate is obtained by scanning in real numbers at the total of 1729 points in the whole domain, which consist of the 36 equidistant points including boundary points, and 1000 random points, as [-0.31211856E-05, 0.42124293E-04]. While this bounds estimate can give some idea, it is not very reliable considering the small number of sampling points and the specific feature of the function. The interval arithmetic covering the whole domain by one interval box gives a mathematically rigorous bounds [ -4.47134 , 4.80774 1, which dramatically shows the typical blow-up problem in interval arithmetic. By dividing the domain in question into smaller interval boxes, we expect to have a narrower bounds estimate. Table 5.11 shows bounds estimates in successively smaller don 9m 1110 Ta ll'lf domain interval boxes at various locations. Only the smallest boxes yield bounds of a size comparable to those obtained by the scanning estimate. However, to cover the entire domain in this fashion requires 1024 small interval boxes, showing the practical 137 limitations of the interval approach for this problem. Table 5.11: Bounds estimates of a normal form deviation function with the interval method. Interval Method domain interval box [ bounds estimate ] width [[0.040000, 0.060000]6[ E —4.47134 , 4.80774 J [9.27908 ] [0.040000, 0.042000]6 [-O.281964E-02, 0.424588E-02] 0.70655E-02 [0.040000, 0.040200]6 {-0.3113035-03, 0.327498E-03] 0.638808—03 [0.040000, 0.040020]6 [-0.304618E-04, 0.329529E-04] 0.63414E-04 [[0.040000, 0.040002]6 [-O.199253E-05, 0.434435E-05] 0.63368E-05 [0.049000, 0.051000]6 [-0.544365E-02, 0.1003323—01] 0.15476E-01 [0.049900, 0.050100]6 [-0.697752E-03, 0.762484E-03] 0.14602E-O2 [0.049990, 0.050010]6 [-0.657704E-O4, 0.802314E-04] 0.14600E-03 [0.049999, 0.050001]6 [-O.320844E-06, 0.1427933-041 0.14600E-04 [0.058000, 0.060000]6 [-0.13307OE-01, 0.265293E-01] 0.39836E-01 [0.059800, 0.060000]6 [-0.160977E-02, 0.188511E-02] 0.34948E-O2 [0.059980, 0.060000]6 [-O.144312E-03, 0.207897E-03] 0.352202—03 [0.059998, 0.060000]6 [ 0.1312909-04, 0.483782E-04] 0.352492—04 On the other hand, the Taylor model computation gives very small remainder bounds as shown in Table 5.12 computed in various orders, and the total bounds estimates are comparable to the bounds by the scanning. Table 5.12 also shows the number of terms of polynomials, which is very moderate compared to the division number necessary for a comparable interval evaluation. 138 Table 5.12: Bounds estimates of a normal form deviation function with the RDA method. RDA Method order terms remainder bounds interval total bounds 6 924 [-0 . 53585E-05 , O . 53588E-05] [-O . 346615-04 , O . 535813-04] 7 1716 [-0 . 83873E-06 , 0 . 83884E-06] ['0 . 30163-04 , O . 4902E-04] 8 3003 [-O . 1232113-06 , 0. 12321E-06] [~O . 29458-04 , 0 . 483 1121-04] Chapter 6 Verified Integration of ODES and Flows In this chapter, Remainder-enhanced Differential Algebraic (RDA) methods are ap- plied for the development of verified integration algorithms for ODEs and flows of ODES. We will use anti-derivations of Taylor models for the solution of the initial value problem $70) = 1777(1), t), (61) where F(t0) 2 F0 and F" is continuous and bounded. We are interested in both the case of a Specific initial condition F0, as well as the case in which the initial condition F0 is a variable and our interest is in the flow of the differential equation 77(1) = M(7‘~'0, t). 6.1 Verified Integration with Taylor Models The goal is to establish a Taylor model for M(f’o, t), and in particular rigorous bounds for the remainder term of the flow of the differential equation over a domain [7:01, F02] X [t0, t1]. In particular, F0 itself may be a Taylor model, as long as its range is known to lie in [F01, F02]. We have to deviate from the direct use of conventional numerical integrators, because they don’t provide rigorous estimates for the integration error, 139 140 but only approximate estimates. Rather, we have to start from scratch from the foundations of the theory of differential equations [18]. 6.1.1 Schauder’s Fixed Point Theorem As is common for the application of functional analysis tools to the study of differential equations, we re-write the differential equation (6.1) as an integral equation t_. MQ=$+ Fwijfl noting that the initial value problem has a (unique) solution if and only if the corre- sponding integral equation has a (unique) solution. Now we introduce the operator A : (30m, 11] —> 6000,11] on the space of continuous functions from [t0, t1] to R" via t _. _. A (f) (t) = 7'0 + F(f (t’) .t’) dt'; (62) to so a general function f in 50 [130, t1] is transformed into a new function in 6°[to, t1] via the insertion into F" and subsequent integration. Having introduced the operator A, the problem of finding a solution to the differential equation is reduced to a fixed-point problem 7‘: A(f‘). It is common fare in the theory of differential equations to establish that Schauder’s fixed point theorem asserts the existence of a solution of an ODE over [t0, t1] in case F is continuous on [t0, t1] X R" and bounded there. If F is even Lipschitz with respect to the first argument f, then Banach’s fixed point theorem even asserts a locally unique solution. We will now apply Schauder’s fixed point theorem in a different way to rigorously obtain a Taylor Model for the flow. 141 Theorem (Schauder): Let A be a continuous operator on the Banach Space X. Let M C X be compact and convex, and let A(M) C M. Then A has a fixed point in A1, i.e. there is an F6 A1 such that A(f) = 7". One should be reminded that the fixed point is not necessarily unique, for ex- ample, the identity map on M has every element of M as fixed points; furthermore compactness and convexity of M are essential, as simple counter-examples show. 6.1.2 Strategy to Satisfy the Requirements In our specific case, X = €O[t0, t1], the space of continuous vector functions on the interval, equipped with the usual maximum norm, and A is the integral operator in (6.2). From continuity of F, it follows easily that A is continuous on X. The process of our application of Schauder’s theorem now has three major steps: 1. Determine a sufficiently large family Y of subsets of X from which to draw candidates for the set M. To satisfy the requirements of Schauder’s theorem, the sets in Y have to be compact and convex; and to fit within our computational framework, it should be possible to contain them in suitable Taylor models. 2. Using the Differential Algebraic structure on Taylor models, construct an initial set Mo 6 Y that satisfies the inclusion property A(Mo) C M0. Once this set has been determined, all requirements of the fixed point theorem are satisfied, and the existence of a solution in M0 has been established. Since the sets in Y were chosen in such a way that they can be contained in Taylor models, a Taylor model inclusion of a solution of the ODE has been found. 3. Finally, the set M0 is iteratively reduced in size in order to obtain bounds that are as sharp as possible. For this purpose, for 2' = 1, 2,3, we construct the 142 sequence M,- = A(M,_1). We have the chain A11 3 A12 3..., and we continue to iterate until no significant further reduction in size is possible. 6.1.3 Schauder Candidate Sets For the first step, it is necessary to establish a family of sets Y from which to draw candidates for Mo. We define Y in the following way. Let (13 + 1) be a Taylor model depending on time t as well as the initial condition F0. Then we define the associated set Mp+f as follows: 111,; ~ c 6000,11]; and forr'e MP4,; F(t0) 2 F0 7(1) 6 [5+1 Vt€[t0,t1] V770 [f’(t')—F(t”)| < klt’—t”[ Vt’,t”€[to,t1] Vf'o. In the last condition, It is bounds for 13", which exists because 13“ is continuous and the solutions can cover only finite range over interval [t0, 151]. The last condition means that all 1" E M ,3 + f are uniformly Lipschitz with constant k. Define the family of candidate sets Y as 6.1.4 Convexity, Compactness, and Invariance First let us check the convexity of the Schauder candidate sets defined above. Definition (Convexity): Let X be a real vector space. M C X is called convex z'f V131,.T2 E M, Va 6 [0,1], (1231+(1— a)x2 E M. Let M C Y be a Schauder candidate set. Then M is convex, because _' _. $1,332 6 M=> 143 X1 I /\__/ (XX|+(1'(1)X2 Figure 6.1: Convexity; a set of numbers M (left) and a set of functions inside a Taylor model (right). af1+(1—a)x'2 E M‘v’aE[0,1], as any such linear combination of two k-Lipschitz functions is k-Lipschitz, is in the same Taylor models as 51 and 52, and assumes the value F0 at to. It is a little more involved to show that M is compact, and it needs a help of the Ascoli-Arzela Theorem. Definition (Compactness): Let (X, (1) be a metric space. M C X is called compact if every sequence in M has at least one cluster point in M. Theorem (Ascoli-Arzela): Let (xn) be a uniformly equicontinuous sequence of functions on [t0, t1], and let (xn) be uniformly bounded, i.e. 3K e R: |:z:,,(t)| < K Vt e [t0,t1], Vn e N. Then, there exists a subsequence of (x,,) that converges uniformly on [t0, t1]. To see that M is compact, let (in) be a sequence of functions in M. Then all in are k-Lipschitz and hence uniformly equicontinuous; since they are in the same Taylor model, they are uniformly bounded. Thus according to the Ascoli-Arzela theorem, (in) has a uniformly convergent subsequence. Let x" be the limit of this subsequence. Since the a, are continuous, so is 5:“, and we obviously have 55*(t0) = F0. Since the elements of the subsequence converging to x“ are k-uniformly Lipschitz, so is 53'“ itself, 144 to \ >§ . 9d /\ t1 \_/\ -K Figure 6.2: Functions on [t0, t1] are uniformly equicontinuous and uniformly bounded. (The Ascoli-Arzela Theorem) as a simple indirect proof reveals. Similarly, since the subsequence converging to a?" is in 13 + 1: so is x“. The interesting to notice is that the schematic explanation of the Ascoli-Arzela theorem in Figure 6.2 reminds us of Taylor models. Finally, the operator A maps any set in Y into another set in Y. Indeed, the image functions of A go through F0 and are continuous because they are integrals, and they are k-Lipschitz because F" is bounded by 11:. Finally, since A is continuous, all images of functions inside a Taylor model are bounded and hence themselves in a Taylor model. Hence the entire problem is reduced to finding a Taylor model I3 + L such that -O —O MP+BCP+, a condition which can be checked computationally using the differential algebraic operations on the set of Taylor models. 145 6.1.5 Satisfying the Inclusion Requirement with Differential Algebraic Methods For practical purposes it is of course in addition desirable to have 1 small. For this purpose it turns out to be important to determine a starting candidate that is on the one hand sufficiently small in width, but on the other hand shaped in such a way as to contain the true solution. This thought leads to attempt sets M * of the form * _ -O M — MMn(f’,t)+I*’ where Mn(f', t) is n-th order Taylor expansion of the solution. If n is high enough, we may expect that the true solution of the ODE and hence the fixed point problem is sufficiently close to the n—th order expansion, and hence that it may be possible to choose 1"" rather small. This approach requires the knowledge of the solution Mn(r’, t), and contrary to the usual situation in which we are only interested in Man: t) at the final value of t, here the explicit dependence on t is required. This quantity can be obtained by iterating (6.2) within the framework of DA. To begin with, one chooses an initial function MS’MF, t) = I, where I is the identity function, and then iteratively sets Mgk+1) =11 A(ng) This process converges to the exact result Mn as a Truncated Taylor Series in n + 1 steps. As the next step, we try to find 1" such that in fact the inclusion property necessary for Schauder’s theorem is satisfied, namely —0 -¢ * A(Mn(r‘, t) + 1*) C Mum t) + 146 Pu =Pn + [0,0] Figure 6.3: Finding an inclusion set as a Taylor model satisfying Schauder require— ment. The suitable choice of 1 requires a little experimenting, it is however greatly simplified by the observation that it is necessary that computationally, PDFW where 110) is defined as Mn(F, t) + in” = A(Mn(F, t) + [0, 0]). (6.3) 110) is a good benchmark for the size of intervals that is to be expected; and so we iteratively try the sequence 11k) = 2k-11"), (6.4) until a computational inclusion can be found, which means that we have established A(Mn(f’, t) + 11“) c Mn(F, t) + it“. (6.5) Once this computational inclusion has been determined, a solution of the ODE is certainly contained in the Taylor model Mn(f', t) + 11"), which satisfies our demand. 6.1.6 Iterative Refinement of the Inclusion For practical purposes it is useful that the sharpness of this solution can be improved. Denoting 11 = 11"), the first obtained interval satisfying our requirement, we itera- 147 I1=Im 12 E [7: Figure 6.4: Iterative refinement of the inclusions as Taylor models. tively define a sequence of Taylor models Mn(f', t) + 1‘; = A(Mn(i’, t) + 111.1). We then must have 1;c C 11-1 to get refinement for all k = 2,3,... To see this, we observe that this is the case for It 2 2 by definition of 1;, and then we infer inductively Mum t) + 1,. C Mn(1‘~‘,t) +124 => A(Mnb", t) + 1}.) C A(M,,(r‘,t) + 17.4) => 14,4731)“;+1 c Mn(r',t)+1;,. Furthermore, the fixed point function F must actually be contained in each of the elements of the sequence of Taylor models Mn(F, t) + 1;. Again by definition of 11, the fixed point function is contained in Mn(7"', t) + 1;, as mentioned in the previous subsection. And by induction, we see 1* e M,(r,t)+i;=> A(r‘) e A(Mnult) + 17,) i» F e Mum t) + 17.1.1. So this provides a mechanism to iteratively refine the inclusion until no further worthwhile decrease in size can be obtained. ‘g‘- ' 'l 148 6.2 Example: Remainder Bounds for a Dipole of the S800 Spectrograph In this section, we will provide an example for the practical use and performance of the method discussed above. In this example, we analyze the motion of a charged particle in a dipole of the S800 spectrograph [60] with a deflection radius of 2.67m and a deflection angle of 75 degrees (Refer to Table 3.8). The Taylor transfer map of the system of differential equations with remainder bounds for a region of initial conditions is determined. The motion is described by the coupled differential equa- tions (2.68), (2.69), (2.70) and (2.71) derived in chapter 2, where h = 1/R with R denoting the deflection radius and By 2 pO/eR, and p 2 p0. The integration was carried out via many steps through the dipole, and the Taylor polynomials describing the dependence of the four final coordinate values 0:, a, y, b on the initial coordinate values and their error bounds were determined. At each step, Taylor models as solutions of the fixed point problem were sought to satisfy the computational inclusion requirement. The computation was performed with various conditions for the comparison. The initial conditions of four phase space variables were chosen to be within the domain interval box {—01, .01], or [—.02, 02]. The order in time and initial conditions was chosen to be 5, 10 or 12, and the fixed step size of 1 degree or 0.5 degree was used. The error bounds estimates under those computational conditions are listed in Table 6.1, 6.2 and 6.3. The overall accuracy goes below 10‘8 using the twelfth order Taylor models with a fixed step size of 0.5 degree with the initial condition within the domain interval box [—.01,.01]. Since no automatic step size control was utilized, the estimates proved conservative and the actual resulting errors were somewhat lower. Table 6.1: Error bounds of the Taylor transfer map of a S800-like dipole with initial 149 condition within [—0.02, 0.02]4 and step size 1°. Error Bounds of the Map with Initial Condition within [—0.02, 0.02]4 step order error bounds estimate size x [- . 1 149558375574990E-03 , 0 . 9296078153993057E-04] 5 a [-.4666650061895221E-04,0.4946294642968968E-04] y [- . 4526632143133069E-05 , 0 . 4526632143133069E-05] b [O . 0000000000000000E+00 , 0 . 0000000000000000E+OO] x [-.4736194115913304E-06,0.7045006167034195E-06] 1° 10 a [- . 4457265344693460E-07 , 0.401 1287904726302E-07] y [- . 1399063567696992E-07 , 0. 1399063567696992E-07] b [0.0000000000000000E+00,0.0000000000000000E+00] x [- . 64018281747543475-07 , 0 . 4534782588950649E-07] 12 a [- . 9293828606969869E-08 , 0. 1002164821 121 11013-07] y [- . 2715296293606005E-08 , 0 . 2715296293606005E-08] b [0 . 0000000000000000E+00 , 0 . 0000000000000000E+00] Table 6.2: Error bounds of the Taylor transfer map of a S800-like dipole with initial condition within [—0.01,0.01]4 and step size 1°. Error Bounds of the Map with Initial Condition within [—0.01,0.01]4 step order error bounds estimate size as [- . 4774670173984800E-04 , 0 . 4224639157994040E-04] 5 a [-.2325461997662418E-04,0.2395429647130834E-04] y [-.1118914478174159E-05,0.1118914478174159E-05] b [0.0000000000000000E+00,0.0000000000000000E+00] x [-.1148378507775330E-06,0.1721148653001733E-06] 1o 10 a [-.7449115285435551E-08,0.6534666433615125E-08] y [-.1680086507721105E-08,0.1680086507721105E-08] b [0.0000000000000000E+00,0.0000000000000000E+00] x [- . 1295756988293815E-07 , 0 . 8486627448846389E-08] 12 a [-.6603760216278913E-09,0.7494712313226484E-09] y [-.1777086692946516E-09,0.1777086692946516E-09] b [0.0000000000000000E+00,0.0000000000000000E+OO] Table 6.3: Error bounds of the Taylor transfer map of a S800-like dipole with initial 150 condition within [—0.01,0.01]4 and step size 0.5°. Error Bounds of the Map with Initial Condition within [—0.01,0.01]4 step order error bounds estimate size :1: [- . 2309808273859498E-04 , 0 . 2031288462862392E-04] 5 a [- . 1 148663798183059E-04 , 0. 1 181920378129873E-04] y [- . 5527454999463349E-06 , 0 . 5527454999463349E-O6] b [0 . 0000000000000000E+00 , 0 . 0000000000000000E+00] x [- . 5756202685083773E-07 , 0 . 853862265475821 1152-07] 0.50 10 a [- . 3131363953715242E-08 , 0 . 28681609104831 14E-08] y [- . 84074155986697513-09 , 0 . 8407415598669751E-09] b [O . 0000000000000000E+00 , O . 0000000000000000E+00] x [- . 6429678314475644E-08 , 0 . 42562248959972198-08] 12 a [- . 30037410547976433-09 , 0 . 32987188012486723-09] y [- . 891375 1 143938564E-10 , 0.8913751 143938564E-10] b [0 . 0000000000000000EZ+00 , 0 . 0000000000000000E+OOJ Table 6.4 shows an example of the typical computational behavior in the process to find Taylor models as solutions of the fixed point problem in a step of the ver- ified integration of the 8800-like dipole. We are considering the first 3 step of the integration using fifth order Taylor model computation with a step size of 1 degree with the initial condition within the domain interval box {—01, .01]. The benchmark intervals 11°) and ilk) in (6.3), (6.4) and their mapped intervals via the operation A are estimated as follows until the inclusion requirement of (6.5) is met, which here was achieved in two steps. The resulting mapped Taylor models are the first solution of the fixed point problem to satisfy the inclusion requirement in this 3 step. The list 151 shows only the error bound intervals of Taylor models. Table 6.4: Computational behavior of the iteration to find the first solution of the inclusion requirement. Iteration to Find the First Solution of Inclusion Requirement iteration step error bounds candkhne x [-.851609168266E-10,0.117999365650E-09] set a [-.107333446719E-11,0.110180867568E-11] y [-.549821641117E-10,0.549821641117E-IO] k=:() 1“” b [0.000000000000E+00,0.000000000000E+00] Inapped x [-.426306772250E-10,0.590512338784E-10] set a [-.130800710095E-11,0.110758538377E-11] y [-.274910820558E-10,0.274910820558E-10] b [0.000000000000E+00,0.000000000000E+00] candkkne x [-.170321833653E-09,0.235998731300E-09] set a (“.214666893438E-11,0.220361735136E-11] y [-.1099643282233-09,0.109964328223E-09] k==l. in) b [0.000000000000E+00,0.000000000000E+00] Inapped x [-.426808960367E-10,0.591027849318E-10] set a [-.207934696832E-11,0.166426642970E-11] y [-.274910820558E-10,0.274910820558E-10] b [0.000000000000E+00,0.000000000000E+00] The trial to refine the solution discussed in subsection 6.1.6 was made, but because the intervals are already so small, in this case no significant refinement resulted. 152 After the integration through the whole dipole via many 3 steps was done, the resulting Taylor polynomials describing the dependence of final on initial coordinates were compared with those obtained by DA computation, and agreement was found. A check about the validity of the remainder bound intervals was done by launching a large collection of rays through the dipole. Because of the homogeneity of the field, orbits can be calculated from purely geometric arguments, and these results were compared by the prediction of the twelfth order map obtained from the verified integrator. For all the 34 rays studied, originating from the interval box of the initial conditions, the difference between the final coordinates determined geometrically and those predicted by the Taylor polynomials via the verified integrator were within the calculated remainder bounds. BIBLIOGRAPHY Bibliography [1] G. Alefeld and J. Herzberger. Introduction to Interval Computations. Academic Press, 1983. [2] A. O. Barut. Electrodynamics and Classical Theory of Fields and Particles. Dover, 1964, 1980. [3] M. Berz. Differential algebraic description of beam dynamics to very high orders. Particle Accelerators, 242109, 1989. [4] M. Berz. Arbitrary order description of arbitrary particle optical systems. Nu- clear Instruments and Methods, A298z426, 1990. [5] M. Berz. High-Order Computation and Normal Form Analysis of Repetitive Systems, in. M. Month (Ed), Physics of Particle Accelerators, volume AIP 249, page 456. American Institute of Physics, 1991. [6] M. Berz. COSY INFINITY Version 6. In M. Berz, S. Martin and K. Ziegler (Eds. ), Proc. Nonlinear Effects in Accelerators, page 125. IOP Publishing, 1992. [7] M. Berz. New features in COSY INFINITY. In Third Computational Accelerator Physics Conference, page 267. AIP Conference Proceedings 297, 1993. [8] M. Berz. Modern map methods for charged particle Optics. Nuclear Instruments and Methods, 3632100, 1995. [9] M. Berz. Differential algebras with remainder and rigorous proofs of long-term stability. In Fourth Computational Accelerator Physics Conference, volume 391, page 221. AIP Conference Proceedings, 1996. [10] M. Berz. COSY INFINITY version 7. In 1997 Particle Accelerator Conference. APS, 1997. [11] M. Berz. COSY INFINITY Version 8 reference manual. Technical Report MSUCL-1088, National Superconducting Cyclotron Laboratory, Michigan State University, East Lansing, MI 48824, 1997. [12] M. Berz. From Taylor series to Taylor models. AIP, 405:1—25, 1997. [13] M. Berz, C. Bischof, A. Griewank, G. Corliss, and Eds. Computational Differ— entiation: Techniques, Applications, and Tools. SIAM, Philadelphia, 1996. 153 [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] 154 M. Berz and G. Hoffstatter. Exact bounds of the long term stability of weakly nonlinear systems applied to the design of large storage rings. Interval Compu- tations, 2:68—89, 1994. M. Berz, G. Hoffstatter, W. Wan, K. Shamseddine, and K. Makino. COSY IN- FINITY and its applications to nonlinear dynamics. in. Computational Differ- entiation: Techniques, Applications, and Tools, M. Berz, C. Bischof, G. Corliss, A. Criewank (Eds), SIAM, 1996. M. Berz, K. Joh, J. A. Nolen, B. M. Sherrill, and A. F. Zeller. Reconstructive correction of aberrations in nuclear particle spectrographs. Physical Review C, 47,2:537, 1993. M. Berz and K. Makino. Arbitrary order maps, remainder terms, and long term stability in particle accelerators. In Optical Science, Engineering and Instrumen- tation ’97. SPIE, 1997. M. Berz and K. Makino. Verified integration of ODES and flows with differential algebraic methods on Taylor models. Reliable Computing, in print. M. Berz and K. Makino. Perturbative determination of first integrals and Lya- punov functions for dynamical systems near fixed points. Journal of Symbolic Computation, submitted. M. Berz, K. Makino, K. Shamseddine, and W. Wan. Applications of Modern Map Methods in Particle Beam Physics. Academic Press, Orlando, Florida, in print, 1998. M. Berz, S. Martin, and K. Ziegler (Eds). Proceedings International Workshop on Nonlinear Problems in Accelerators Berlin. Institute of Physics Publishing, Bristol, 1993. I. M. Bomze, T. Csendes, R. Horst, and P. M. Pardalos, editors. Developments in Global Optimization. KLUWER, 1997. J. A. Caggiano, D. Bazin, B. S. Davids, R. Foutus, D. Karnes, P. Johnson, B. Sherrill, and A. Zeller. S800 spectrograph dipole mapping. Annual Report of the Michigan State University National Superconducting Cyclotron Laboratory, 1996. J. A. Caggiano and B. M. Sherrill. Data analysis techniques for S800 dipole magnetic field maps. Annual Report of the Michigan State University National Superconducting Cyclotron Laboratory, 1995. D. C. Carey. The Optics of Charged Particle Beams. Harwood, 1987. A. W. Chao. Physics of Collective Beam Instabilities in High Energy Accelerators. Wiley, 1993. Ingrid Daubechies. Ten Lectures on Wavelets. SIAM, Philadelphia, 1992. D. A. Edwards and M. J. Syphers. An Introduction to the Physics of High Energy Accelerators. Wiley, 1993. 155 [29] H. Goldstein. Classical Mechanics. Addison-Wesley, Reading, MA, 1980. [30] I. S. Gradshteyn and I. M. Ryzhik. Table of Integrals, Series, and Products. Academic Press, New York, 1980. [31] A. Griewank, G. F. Corliss, and Eds. Automatic Differentiation of Algorithms. SIAM, Philadelphia, 1991. [32] E. Hansen. Global Optimization using Interval Analysis. Marcel Dekker, 1992. [33] E. R. Hansen. Global optimization using interval analysis — the one—dimensional case. J. Optim. Theor. and Appl., 29:331—334, 1979. [34] E. R. Hansen. An Overview of Global Optimization Using Interval Analysis, pages 289—307. Academic Press, New York, 1988. [35] P. W. Hawkes and E. Kasper. Principles of Electron Optics. Academic Press, London,1989. [36] G. Hoffstatter and M. Berz. Efficient computation of fringe fields using symplectic scaling. In Third Computational Accelerator Physics Conference, page 467. AIP Conference Proceedings 297, 1993. [37] G. Hoffstatter and M. Berz. Accurate and fast computaton of high—order fringe field maps via symplectic scaling. Nuclear Instruments and Methods, 363:124, 1995. [38] G. Hoffstatter and M. Berz. Symplectic scaling of transfer maps including fringe fields. Physical Review E, 54,4, 1996. [39] G. H. Hoffstatter. Rigorous bounds on survival times in circular accelerators and efi‘lcient computation of fringe-field transfer maps. PhD thesis, Michigan State University, East Lansing, Michigan, USA, 1994. also DESY 94-242. [40] S. Humphries. Principles of Charged Particle Acceleration. Wiley, New York, 1986. [41] K. Ichida and Y. Fujii. An interval arithmetic method for global optimization. Computing, 23:85—97, 1979. [42] J. D. Jackson. Classical Electrodynamics. Wiley, New York, 1975. [43] C. Jansson. A global optimization method using interval arithmetic. [MACS Annals of Computing and Applied Mathematics, 1992. [44] X. Jiye. Aberration Theory in Electron and Ion Optics. Advances in Electronics and Electron Physics, Supplement 17. Academic Press, Orlando, Florida, 1986. [45] R. Baker Kearfott. Rigorous Global Search: Continuous Problems. KLUWER, 1996. [46] R. Baker Kearfott and Vladik Kreinovich, editors. Applications of Interval Com— putations, volume 3. Kluwer, 1996. 156 [47] S. Kowalski and H. Enge. RAYTRACE. Technical report, MIT, Cambridge, Massachussetts, 1985. [48] U. W. Kulisch and W. F. Miranker. Computer Arithmetic in Theory and Practice. Academic Press, New York, 1981. [49] R. S. MacKay and J. D. Meiss. Hamiltonian Dynamical Systems. Adam Hilger, 1987. [50] K. Makino. Rigorous integration of maps and long-term stability. In 1997 Particle Accelerator Conference. APS, 1997. [51] K. Makino and M. Berz. COSY INFINITY Version 7. In Fourth Computa- tional Accelerator Physics Conference, volume 391, page 253. AIP Conference Proceedings, 1996. [52] K. Makino and M. Berz. Remainder differential algebras and their applica- tions. in. Computational Differentiation: Techniques, Applications, and Tools, M. Berz, C. Bischof, G. Corliss, A. Griewank (Eds. ), SIAM, 1996. [53] K. Makino and M. Berz. Arbitrary order aberrations for elements characterized by measured fields. In Optical Science, Engineering and Instrumentation ’97. SPIE, 1997. [54] K. Makino and M. Berz. Implementation and applications of Taylor model methods. Reliable Computing, submitted, 1997. [55] K. Makino and M. Berz. Verified quadrature and ODE integration involving elementary functions of high complexity. Journal of Symbolic Computation, sub- mitted. [56] L. Michelotti. Intermediate Classical Dynamics with Applications to Beam Physics. Wiley, 1995. [57] Ramon E. Moore. Methods and Applications of Interval Analysis. SIAM, 1979. [58] S. Moriguchi, K. Udagawa, and N. Ichimatsu. Mathematics Formulas I. Iwanami Zensho. Iwanami Shoten, Tokyo, 1956. [59] A. Neumaier. Interval Methods for Systems of Equations. Cambridge, 1990. [60] J. Nolen, A.F. Zeller, B. Sherrill, J. C. DeKamp, and J. Yurkon. A proposal for construction of the S800 spectrograph. Technical Report MSUCL-694, National Superconducting Cyclotron Laboratory, 1989. [61] Annual Report of the Michigan State University National Superconducting Cy- clotron Laboratory, 1994. [62] Annual Report of the Michigan State University National Superconducting Cy- clotron Laboratory, 1995. [63] Annual Report of the Michigan State University National Superconducting Cy- clotron Laboratory, 1996. I 71:.in if? “AV; 5.. I h—h‘um .. 157 [64] M. Reiser. Theory and Design of Charged Particle Beams. Wiley, 1994. [65] W. Wan. Private communication, 1997. [66] R. L. Warnock and R. D. Ruth. Long—term bounds on nonlinear Hamiltonian motion. Physica D, 56(14):188—215, 1992. also SLAC—PUB—5267. [67] H. Wiedemann. Particle Accelerators Physics. Springer-Verlag, 1994. [68] H. Wiedemann. Particle Accelerators Physics 11. Springer-Verlag, 1995. [69] H. Wollnik. Charged Particle Optics. Academic Press, Orlando, Florida, 1987. I- :— “E“‘T‘l "‘111111111111111ES