On the parallel computation of invariant tori a, ` Enric Castell` Angel Jorba and Estrella Olmedo Universitat de Barcelona December 4, 2008 1 / 45
Setting We focus on dynamical systems of the form ¯ x = f ( x , θ ) , ¯ θ = θ + ω, where x ∈ R n , θ ∈ T r . The frequency vector ω ∈ R r is supposed to be irrational. The autonomous case ¯ x = f ( x ) is included in this setting. 2 / 45
Setting An invariant torus can be represented by a map x : ( θ, ϕ ) ∈ T r × T s �→ R n , and it must satisfy the invariance condition f ( x ( θ, ϕ )) = x ( θ + ω, ϕ + ν ) . For simplicity, we will explain the methods for tori that do not depend on the “inner” angles ϕ , although in some of the examples used later on we will compute tori depending on inner angles. The parametrization satisfies the equation f ( x ( θ )) = x ( θ + ω ) . 3 / 45
Setting Suppose that the map has an invariant curve with rotation number ω . The curve is given (in parametric form) by a map x : T 1 → R n . Let us write x ( θ ) as a real Fourier series, � x ( θ ) = a 0 + a k cos( k θ ) + b k sin( k θ ) , k > 0 where a k , b k ∈ R n , k ∈ N . As it is usual in numerical methods, we look for a truncation of this series. So, let us fix in advance a truncation value N (the selection of N will be discussed later on), and let us try to determine an approximation to the 2 N + 1 unknown coefficients a 0 , a k and b k , 0 < k ≤ N . 4 / 45
Newton method We consider the map x ( θ ) �→ F ( x ( θ )) = f ( x ( θ )) − x ( θ + ω ) , where x ( θ ) denotes the parametrization of a torus. The main idea is to apply a Newton method to find x ( θ ) such that F ( x ( θ )) ≡ 0. We note that F acts on a space of periodic functions. First, let us define a mesh of 2 N + 1 points on T 1 , 2 π j θ j = 2 N + 1 , 0 ≤ j ≤ 2 N . It is not difficult to compute x ( θ j ), f ( θ j ), f ( θ j + ω ) and, hence, F ( x ( θ j )). From these values, we can derive the Fourier coefficents of F ( x ( θ )). 5 / 45
Newton method Therefore, we have a procedure to compute the map F . As this procedure can be easily differentiated, we can also obtain DF . Then, a Newton method can be applied: x m +1 = x m − ( DF ( x m )) − 1 F ( x m ) . 6 / 45
Error estimates A natural question is about the size of the error of the obtained curve. To measure such error we use E ( x , ω ) = max θ ∈ T 1 � f ( x ( θ )) − x ( θ + ω ) � . We estimate E ( x , ω ) using a much finer mesh than the one used in the previous computations. If this error is too big, we increase N and we apply the Newton process again. 7 / 45
Linearized normal behaviour Let h represent a small displacement with respect to an arbitrary point x ( θ ) on the invariant curve. Then, f ( x ( θ ) + h ) = f ( x ( θ )) + D x f ( x ( θ )) h + O ( � h � 2 ) . As f ( x ( θ )) = x ( θ + ω ), it follows that the linear normal behaviour is described by the following dynamical system, � ¯ x = A ( θ ) x , ¯ = θ + ω θ where A ( θ ) = D x f ( x ( θ )). 8 / 45
Linearized normal behaviour Definition The previous linear system is called reducible iff there exists a (may be complex) change of variables x = C ( θ ) y such that it becomes ¯ = � y By , ¯ θ = θ + ω where the matrix B ≡ C − 1 ( θ + ω ) A ( θ ) C ( θ ) does not depend on θ . 9 / 45
Linearized normal behaviour We define the operator T ω : ψ ( θ ) ∈ C ( T 1 , C n ) �→ ψ ( θ + ω ) ∈ C ( T 1 , C n ) , and let us consider now the following generalized eigenvalue problem: to look for couples ( λ, ψ ) ∈ C × C ( T 1 , C n ) such that A ( θ ) ψ ( θ ) = λ T ω ψ ( θ ) . 10 / 45
Linearized normal behaviour Definition Two eigenvalues λ 1 and λ 2 are said to be unrelated iff λ 1 � = exp( i k ω ) λ 2 , ∀ k ∈ Z . Otherwise, we refer to them as related. Lemma Assume that there exist n unrelated eigenvalues λ 1 , . . . , λ n for the previous eigenproblem. Then, the linear skew product can be reduced to constant coefficients, where B = diag ( λ 1 , . . . , λ n ) . 11 / 45
y L 4 + The Bicircular problem E L L x It is a model for the study of the dynamics of a small particle in the 2 1 L 3 t y + + + Earth-Moon-Sun system. M � + y L 5 S 12 / 45
The Bicircular problem The BCP can be described by the Hamiltonian system, 1 + yp x − xp y − 1 − µ µ − m S p 2 x + p 2 y + p 2 � � = − H BCP z 2 r PE r PM r PS − m S ( y sin θ − x cos θ ) , a 2 S where ( x − µ ) 2 + y 2 + z 2 , r 2 = PE ( x − µ + 1) 2 + y 2 + z 2 , r 2 = PM ( x − x S ) 2 + ( y − y S ) 2 + z 2 , r 2 = PS being x S = a S cos θ , y S = − a S sin θ and θ = ω S t . 13 / 45
The Bicircular problem 0.55 0.545 VF2 VF2 0.54 0.535 VF1 VF1 0.53 VF3 VF3 0.525 0.52 0.515 0.51 0.505 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 14 / 45
The Bicircular problem 1.5 1 0.5 0 -0.5 -1 -1.5 -1.5 -1 -0.5 0 0.5 1 1.5 N = 16 (total dimension: 198). 15 / 45
The Bicircular problem Modulus Argument λ 1 1.091942641437887 0.000000000000000 λ 2 0.915799019152856 0.000000000000000 λ 3 0.999999999999985 2.035517841801725 λ 4 0.999999999999985 -2.035517841801725 Normal eigenvalues around an unstable invariant curve of the family VF1. The rotation number is ω = 0 . 535033339385478, and the value of the ˙ z coordinate when z = 0 is ˙ z = 0 . 080508698608030. We can check that | λ 1 λ 2 − 1 | ≈ 4 × 10 − 15 . 16 / 45
The Bicircular problem 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04 Motion of one of the couples of eigenvalues in the complex plane, near the change of stability in the families VF1 and VF2. 17 / 45
Computing the unstable manifold We call ψ j the eigenfunction corresponding to λ j , j = 1 , . . . , 4, and we focus on the couple ( λ 1 , ψ 1 ) ∈ R × C ( T 1 , R n ). The linearized unstable manifold is given by x ( θ ) + h ψ 1 ( θ ). To estimate a suitable value for h , we note that f ( x ( θ )) + hD x f ( x ( θ )) ψ 1 ( θ ) + O ( h 2 ) f ( x ( θ ) + h ψ 1 ( θ )) = x ( θ + ω ) + h λ 1 ψ 1 ( θ + ω ) + O ( h 2 ) . = Hence, the size of the term O ( h 2 ) can be estimated by a numerical evaluation of E ( h ) = max θ ∈ T 1 � f ( x ( θ ) + h ψ 1 ( θ )) − x ( θ + ω ) − h λ 1 ψ 1 ( θ + ω ) � 2 . It follows that h = 10 − 7 is enough to have E ( h ) < 10 − 13 . We define the curve C 1 ⊂ R n as the image of the map θ �→ x ( θ ) + h ψ 1 ( θ ) and, for j > 1, C j = f ( C j − 1 ). 18 / 45
0.8705 1 0.87 0.98 0.8695 0.96 0.869 0.94 0.8685 0.868 0.92 0.8675 0.9 0.867 0.88 0.8665 0.866 0.86 -0.493 -0.4925 -0.492 -0.4915 -0.491 -0.4905 -0.49 -0.4895 -0.489 -0.4885 -0.488 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 19 / 45
0.878 1 0.876 0.874 0.95 0.872 0.87 0.9 0.868 0.866 0.85 0.864 0.862 0.8 0.86 0.858 0.75 -0.52 -0.51 -0.5 -0.49 -0.48 -0.47 -0.46 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 20 / 45
Tori of higher dimensions Now we assume that x : ( θ, ϕ ) ∈ T r × T s �→ R n , with r > 1. We use multidimensional truncated Fourier series to approximate the torus. DF is a full matrix of large dimension and, hence, it requires a large amount of memory, solving the linear system ( DF ) h m = f m requieres a lot of CPU time. 21 / 45
Parallelization (I) We assume we have a system with distributed memory (for instance, a cluster). The matrix can be distributed in different machines The linear system is solved in parallel For this scheme, we have coded a (parallel) QR factorization using Householder reflections. Advantages: We can deal with Fourier expansions with many terms, the computations are distributed Inconveniences: When the number of machines increases, the communications become the new bottleneck. 22 / 45
y L Example: The EBCP 4 + The Elliptic Bicircular Problem E L L x 2 1 L 3 t y + + + M � + y L 5 S Here we assume that the motion of Earth and Moon is elliptical. 23 / 45
Example: The EBCP The equations of motion can be written in Hamiltonian form, H = H RTBP ( x , y ) + ˆ H ( x , y , θ ) + � I θ , ω � , where x ∈ R 3 , y ∈ R 3 and θ ∈ T 2 . We are interested in computing the 1-parametric family of 3-D tori corresponding to the vertical family of periodic orbits of the RTBP near L 5 . These tori are parametrized by 3 angles, the two of the perturbation and an inner angle coming from the vertical periodic orbit considered. 24 / 45
Example: The EBCP We use the Poincar´ e section θ 1 = 0 (mod 2 π ). When applying the Newton method, we obtained linear systems of dimension up to 20,000. We used up to 16 nodes to solve the problem. 25 / 45
Example: The EBCP 13 14 15 12 0.541 16 16 0.54 11 17 11 10 9 18 8 7 19 17 1 21 2 345 6 20 10 22 0.53 0.538 23 9 8 0.52 18 4 24 0.535 7 19 6 2 3 5 0.51 0 0.2 0.4 0.6 0.8 0.15 0.2 0.25 26 / 45
Example: The EBCP Torus number 2 0.16 0.905 0.08 0.88 0 0.855 -0.08 0.83 -0.16 -0.515 -0.5 -0.485 -0.47 -0.515 -0.5 -0.485 -0.47 0.16 0.08 0 -0.08 -0.16 -0.16 -0.08 0 0.08 0.16 27 / 45
Example: The EBCP Torus number 6 0.2 0.88 0.1 0.85 0 0.82 -0.1 0.79 -0.2 -0.58 -0.56 -0.54 -0.52 -0.58 -0.56 -0.54 -0.52 0.2 0.1 0 -0.1 -0.2 -0.2 -0.1 0 0.1 0.2 28 / 45
Recommend
More recommend
