On the parallel computation of invariant tori a, ` Enric Castell` - - PowerPoint PPT Presentation

On the parallel computation of invariant tori

Enric Castell` a, ` 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 ∈ Rn, θ ∈ Tr. The frequency vector ω ∈ Rr 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 : (θ, ϕ) ∈ Tr × Ts → Rn, 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 : T1 → Rn. Let us write x(θ) as a real Fourier series, x(θ) = a0 +

• k>0

ak cos(kθ) + bk sin(kθ), where ak, bk ∈ Rn, 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 2N + 1 unknown coefficients a0, ak and bk, 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 2N + 1 points on T1, θj = 2πj 2N + 1, 0 ≤ j ≤ 2N. 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: xm+1 = xm − (DF(xm))−1F(xm).

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

θ∈T1 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(θ)) + Dxf (x(θ))h + O(h2). As f (x(θ)) = x(θ + ω), it follows that the linear normal behaviour is described by the following dynamical system, ¯ x = A(θ)x ¯ θ = θ + ω

• ,

where A(θ) = Dxf (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(T1, Cn) → ψ(θ + ω) ∈ C(T1, Cn), and let us consider now the following generalized eigenvalue problem: to look for couples (λ, ψ) ∈ C × C(T1, Cn) 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

The Bicircular problem

It is a model for the study of the dynamics of a small particle in the Earth-Moon-Sun system.

• x

12 / 45

The Bicircular problem

The BCP can be described by the Hamiltonian system, HBCP = 1 2

• p2

x + p2 y + p2 z

• + ypx − xpy − 1 − µ

rPE − µ rPM − mS rPS −mS a2

S

(y sin θ − x cos θ) , where r2

PE

= (x − µ)2 + y2 + z2, r2

PM

= (x − µ + 1)2 + y2 + z2, r2

PS

= (x − xS)2 + (y − yS)2 + z2, being xS = aS cos θ, yS = −aS sin θ and θ = ωSt.

13 / 45

The Bicircular problem

0.505 0.51 0.515 0.52 0.525 0.53 0.535 0.54 0.545 0.55 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 VF1 VF2 VF3 VF1 VF2 VF3 14 / 45

The Bicircular problem

• 1.5
• 1
• 0.5

0.5 1 1.5

• 1.5
• 1
• 0.5

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.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(T1, Rn). The linearized unstable manifold is given by x(θ) + hψ1(θ). To estimate a suitable value for h, we note that f (x(θ) + hψ1(θ)) = f (x(θ)) + hDxf (x(θ))ψ1(θ) + O(h2) = x(θ + ω) + hλ1ψ1(θ + ω) + O(h2). Hence, the size of the term O(h2) can be estimated by a numerical evaluation of E(h) = max

θ∈T1 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 C1 ⊂ Rn as the image of the map θ → x(θ) + hψ1(θ) and, for j > 1, Cj = f (Cj−1).

18 / 45

0.866 0.8665 0.867 0.8675 0.868 0.8685 0.869 0.8695 0.87 0.8705

• 0.493 -0.4925 -0.492 -0.4915 -0.491 -0.4905
• 0.49
• 0.4895 -0.489 -0.4885 -0.488

0.86 0.88 0.9 0.92 0.94 0.96 0.98 1

• 0.5
• 0.4
• 0.3
• 0.2
• 0.1

0.1

19 / 45

0.858 0.86 0.862 0.864 0.866 0.868 0.87 0.872 0.874 0.876 0.878

• 0.52
• 0.51
• 0.5
• 0.49
• 0.48
• 0.47
• 0.46

0.75 0.8 0.85 0.9 0.95 1

• 0.8
• 0.7
• 0.6
• 0.5
• 0.4
• 0.3
• 0.2
• 0.1

0.1

20 / 45

Tori of higher dimensions

Now we assume that x : (θ, ϕ) ∈ Tr × Ts → Rn, 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)hm = fm 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

Example: The EBCP

The Elliptic Bicircular Problem

• x

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 = HRTBP(x, y) + ˆ H(x, y, θ) + Iθ, ω , where x ∈ R3, y ∈ R3 and θ ∈ T2. 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 L5. 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

0.51 0.52 0.53 0.54 0.2 0.4 0.6 0.8 1 2 345 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 0.535 0.538 0.541 0.15 0.2 0.25 2 3 4 5 6 7 8 9 10 11 16 17 18 19

26 / 45

Example: The EBCP

Torus number 2

0.83 0.855 0.88 0.905

• 0.515
• 0.5
• 0.485
• 0.47
• 0.16
• 0.08

0.08 0.16

• 0.515
• 0.5
• 0.485
• 0.47
• 0.16
• 0.08

0.08 0.16

• 0.16
• 0.08

0.08 0.16

27 / 45

Example: The EBCP

Torus number 6

0.79 0.82 0.85 0.88

• 0.58
• 0.56
• 0.54
• 0.52
• 0.2
• 0.1

0.1 0.2

• 0.58
• 0.56
• 0.54
• 0.52
• 0.2
• 0.1

0.1 0.2

• 0.2
• 0.1

0.1 0.2

28 / 45

Example: The EBCP

Torus number 14

0.92 0.935 0.95 0.965

• 0.15
• 0.13
• 0.11
• 0.08
• 0.04

0.04 0.08

• 0.15
• 0.13
• 0.11
• 0.08
• 0.04

0.04 0.08

• 0.08
• 0.04

0.04 0.08

29 / 45

Example: The EBCP

Torus number 18

0.85 0.88 0.91 0.94

• 0.4
• 0.36
• 0.32
• 0.28
• 0.2
• 0.1

0.1 0.2

• 0.4
• 0.36
• 0.32
• 0.28
• 0.2
• 0.1

0.1 0.2

• 0.2
• 0.1

0.1 0.2

30 / 45

A different Newton scheme

The proof of several KAM-related results for lower dimensional tori involve Newton methods. Most of them are not based on inverting an operator, but on reducing (to constant coefficients) the linearized dynamics around the torus. This is done iteratively, as a part of the Newton process. These schemes are explicit and constructive, so they can be used as an efficient numerical method (see also recent papers by A. Haro and R. de la Llave).

31 / 45

A different Newton scheme

Initial approximations: y0 = ¯ x0 − f (x0, θ) y0 ≈ ε ¯ x = A0(θ)x ¯ θ = θ + ω

• x=C0(θ)y

− → ¯ y = (B0 + Q0(θ))y ¯ θ = θ + ω A0(θ) = Dxf (x0(θ), θ), Q0(θ) ≈ ε

32 / 45

Another Newton scheme

After one step: y1 = ¯ x1 − f (x1, θ), y1 ≈ ε2 ¯ x = A1(θ)x ¯ θ = θ + ω

• x=C1(θ)y

− → ¯ y = (B1 + Q1(θ))y ¯ θ = θ + ω A1(θ) = Dxf (x1(θ), θ), Q1(θ) ≈ ε2

33 / 45

Improving the solution

We want x1(θ) = x0(θ) + h(θ) with h(θ) small h(θ) must satisfy: h(θ + ω) = A0(θ)h(θ) − y0(θ) We apply: h(θ) = C0(θ)u(θ) We obtain: u(θ + ω) = B0u(θ) + r(θ) We solve: 2N + 1 independent linear systems n × n

34 / 45

Improving the Floquet matrix

We want C1(θ) = C0(θ)(I + H(θ)) with H(θ) small We want to reduce up to order ε2 ¯ x = A1(θ)x, ¯ θ = θ + ω We apply x(θ) = C0(θ)(I + H(θ))z(θ) We obtain H(θ + ω)B1 = B1H(θ) + R(θ) We solve N independent linear systems 2n2 × 2n2

35 / 45

Example II

“Toy” example: ˙ x = y ˙ y = −α sin x + ε d + 2 +

d

• i=0

cos θi with α = 0.8, d = 4, θi = ωit + θ0

i , θ0 i are the initial phases and

ω0 = 1, ω1 = √ 2, ω2 = √ 3, ω3 = √ 5, ω4 = √ 7.

36 / 45

Example II

There is a 5-D dimensional torus that branches off from the origin when the perturbation is added. This torus is of the saddle type. We take the section θ0 = 0 (mod 2π) and we apply the previous scheme, to compute a 4-D torus for the Poincar´ e map. The accuracy to compute the torus is 10−10.

37 / 45

Example II

Several 2-D slices of the torus

• 0.08
• 0.06
• 0.04
• 0.02

0.02 0.04 0.06 0.08

• 0.02

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

• 0.15
• 0.1
• 0.05

0.05 0.1 0.15

• 0.02

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

• 0.15
• 0.1
• 0.05

0.05 0.1 0.15

• 0.02

0.02 0.04 0.06 0.08 0.1 0.12

• 0.06
• 0.04
• 0.02

0.02 0.04 0.06

• 0.04
• 0.02

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

• 0.02
• 0.015
• 0.01
• 0.005

0.005 0.01 0.015 0.02

• 0.08
• 0.06
• 0.04
• 0.02

0.02 0.04 0.06 0.08

• 0.15
• 0.1
• 0.05

0.05 0.1 0.15

• 0.04
• 0.02

0.02 0.04 0.06 0.08 0.1 0.12

38 / 45

Example II

Number of Fourier coefficients: 1,975,467. p Total 1 109m 45.770s 2 58m 5.449s 4 32m 7.966s 6 23m 28.204s 8 19m 8.460s 10 16m 55.771s 12 16m 43.602s

39 / 45

Example III

We use a model by Gomez, Masdemont & Mondelo for the motion of a particle in the Earth-Moon system. The model includes effects from the Sun and the noncircular motion of the Moon. These models are written as a perturbation of the RTBP H = 1 2{p2

x + p2 y + p2 z} + ypx − xpy −

1 − µ {(x − µ)2 + y2 + z2}(1/2) − µ {(x − µ + 1)2 + y2 + z2}(1/2) + εgi(x, y, z, px, py, pz, θ). When ε = 0, the equations are the equations of RTBP and when ε = 1, we have model SSSM1 if i = 1 with θ ∈ T; SSSM2 if i = 2, with θ ∈ T2; and SSSM3 if i = 3, with θ ∈ T3.

40 / 45

Example III

We focus on the computation of the invariant tori that substitute L1,2. Due to the instability of this region, we use a parallel shooting method. This means that, in the process of reducing the flow to a map, we will increase the dimension of the phase space.

41 / 45

Example III

For instance, to compute the substitute of L1 we have used 4 intermediate Poincar´ e sections. This means that we look for a torus of the map R24 × Td − → R24 × Td (x, θ) → (P(x), θ + ρ), To compute the substitute of L2 we have used 3 sections: R18 × Td − → R18 × Td (x, θ) → (P(x), θ + ρ), The tori have been computed with an accuracy of 10−10.

42 / 45

Example III

43 / 45

Example III

44 / 45

Conclusions

It is possible to approximate invariant tori of moderate dimension (2,3, 4

• r even 5) in a phase space of moderate dimension.

Computing the Floquet transformation at the same time as the torus increases the degree of parallelism, but sometimes the number of Fourier modes to approximate the Fourier transformation becomes very large. So, if the transformation needed to reduce the torus is much more complex than the torus, this method could be a bad option.

45 / 45