Symbolic and Numerical Modeling of Nonlinear Dynamics of Particles in Accelerators

This paper is devoted to one of the methods of symbolic computation, which based on perturbation theory and constructed in the framework of matrix formalism. This method can be used for solving optimization problems and preliminary modeling in accelerator physics. There are presented the main theoretical positions and demonstration the result of this method on example of one of classical problems.


Introduction
The main instruments of high-energy physics are accelerators, the effective work of which depends both on the mathematical models and the corresponding software.Nowadays there are several widely used packages for this kind of tasks, for example, MAD [1], COSY Infinity [2] and MaryLie [3], which allow getting consistent results with the data.The first two (MAD, COSY Infinity) use phase variables to describe the beam, and both are based on numerical methods.The main approach used in them is the trajectory analysis.But numerical methods do not have the functionality that is required in the modern world, for example, a feature of numerical methods is the need to repeat all calculations under system parameter changing, it is time-consuming.The MaryLie package, unfortunately, does not support the Lego-ideology, and therefore requires new calculations for the entire system when replacing any of its elements.
Symbolic calculations in turn are cheap in time; they can support the Lego-ideology, and also taking into account symplecticity.In this paper, we propose to use a matrix approach [4] to modeling of nonlinear dynamical systems.It allows the use of symbolic calculations, such that cumbersome calculations are performed just once to yield final formulas which can be stored and used later.Also, the matrix approach allows us consider trajectories beam as a single object.Particular this allows us to move from coordinates and impulses to quantities that can be measured experimentally.

Theoretical basis
In this section we describe shortly a schematic description of the theory of constructing a solution of a system of ordinary differential equations in explicit form using perturbation theory.Consider the equation of motion in the form dX/ds = F(X, s), F(0, s) = 0, where X is a vector of phase moments and the arbitrary analytic function F(X, s) is defined within a neighborhood of X = 0, X ∈ R n and is measurable at s ∈ R n .The following definition for the vector of phase moment is necessary: X is a vector of phase moments of the k-th order, representing the k-th Kronecker degree of the phase vector X.
, the elements of which are the derivatives of the components of the vectorvalued function F(0, s) = 0 of the k-th order.Note that P 1k (s) are matrices which include the control parameters of the system.Thus, taking into account the assumption that F(0, s) = 0, we can write: The solution of such a nonlinear system in the form of 2-dimensional matrices, which can be calculated according to the algorithm presented in [4], is sought in the form: Let us note that the matrices R 1k (s) contain coefficients for the corresponding orders of the elements of the vector X 0 , where X 0 = X(s 0 ).A finite cut-off of this infinite series can be defined which follows from the properties of the object under consideration (accelerator).
However, multi-rotation in cyclic machines and the necessity to make a huge number of integration steps leads to the need to preserve the integrals of motion, for example, the law of conservation of energy, and symplecticity.The requirements for modern cyclic systems lead to the necessity of using nonlinear control fields up to some k-th order, which in turn leads to a violation of the symplectic property.This is a critical moment for this kind of tasks.In terms of the matrix formalism, the symplecticity condition can be ensured using the Jacobi matrix M(X, s|s 0 ) of the transformation M(X, s|s 0 ; H) of the dynamical system as follows: where H is the Hamiltonian of the system.Traditionally, symplectic properties can be restored using various path-based methods, in particular, the methods of canonical transformations [5].This approach is used for each individual trajectory of particles (10 9 -10 12 ) and needs a lot of time.In this paper it is proposed to make corrections not in the trajectory, but in the matrices R 1k .The proposed approach is based on the correction of the elements of the transformation matrices R 1k (s) to the desired order of nonlinearity [6].As a result, the amount of computations performed is much less than on-trajectory analysis.Note that the simplification procedure must be performed sequentially for all orders (from 1 to the k-th), because the elements of each current transformation matrix have a relationship with the elements of the previous one.

Demonstration example
The general approach described above is demonstrated using the example of the dynamics of a pendulum with a cubic nonlinearity in comparison with the solution of this equation by the fourth-order Runge-Kutta method.The analytic solution presented in [7] was used to check the numerical results.The system of differential equations corresponding to the dynamics of such a 2-dimensional pendulum is Hamiltonian and has the form: where α is a variable parameter of the system, s is length along the reference trajectory, H is the Hamiltonian of the system.The solution of the system (4), describing dynamics of the 2-dimensional pendulum, i.e.X = (q, p) T was obtained under initial conditions q(0) = 0.8, p(0) = 0. Taking into account symplectic corrections, final solutions were obtained, the phase portraits of which are consistent with the ones presented [7]: Fig. 1 provides plots of the phase portraits under various α.The parameter α is responsible for the compression along the p axis: the larger it is, the stronger the compression.At α = 9.8 the Runge-Kutta method of the 4th order fails to generate the right solution, whereas within framework of the matrix formalism the true solution is found.Let's pay attention to fact that phase portraits of the numerical solution are closed for different values of s.The solutions obtained within framework of matrix formalism always close at s = 2π, due to the correction procedure performed.Fig. 2 shows the time dependence for numerical calculations of the integrals used to obtain the analytical solution with different accuracies.It can be seen up to an accuracy of ±10 −3 or less, the analytical method allows almost instantaneous derivation of the solution (0.14 seconds).To reach an accuracy of ±10 −4 , the numerical method already requires significant time costs -13 seconds.

Conclusion
This work is for the most part demonstrative of the main ideas of our approach.The features of the proposed approach can be summarized as follows: the correction of the mappings does not depend on the initial state of the system.The computational cost is decreased due to the possibility to perform only once and reuse a similar correction which is needed at each order of the step by step algorithm.In addition, these operations with matrices are performed independently of the initial dynamic system and corresponding results can be stored in a special database.
In the case described in this paper, the laboriousness of symbolic computations is not too high, however, in the case of dynamics of particle beam it is.In addition to all of the above, we would like to emphasize once again the intrinsic property of this method of being independent of the concrete representation of dynamic (Hamiltonian) system.

Figure 2 .
Figure 2. Comparing the time spent on obtaining a solution without symplectic corrections analytically and numerically.