An Introduction to Numerical Continuation Methods with Applications - - PowerPoint PPT Presentation

An Introduction to Numerical Continuation Methods with Applications Eusebius Doedel IIMAS-UNAM

July 28 - August 1, 2014

Persistence of Solutions

• Newton’s method for solving a nonlinear equation [B83]1

G(u) = 0 , G(·) , u ∈ Rn , may not converge if the “initial guess” is not close to a solution.

• However, one can put a homotopy parameter in the equation.
• Actually, most equations already have parameters.
• We will discuss persistence of solutions to such equations.

1 See Page 83+ of the Background Notes on Elementary Numerical Methods.

1

The Implicit Function Theorem

Let G : Rn × R → Rn satisfy (i) G(u0, λ0) = 0 , u0 ∈ Rn , λ0 ∈ R . (ii) Gu(u0, λ0) is nonsingular (i.e., u0 is an isolated solution) , (iii) G and Gu are smooth near u0 . Then there exists a unique, smooth solution family u(λ) such that

• G(u(λ), λ) = 0 ,

for all λ near λ0 ,

• u(λ0) = u0 .

NOTE : The IFT also holds in more general spaces · · · 2

EXAMPLE : A Simple Homotopy . ( Course demo : Simple-Homotopy2) Let g(u, λ) = (u2 − 1) (u2 − 4) + λ u2 e

1 10 u .

When λ = 0 the equation g(u, 0) = 0 , has four solutions , namely, u = ± 1 , and u = ± 2 . We have gu(u, λ)

• λ=0

≡ d du(u, λ)

• λ=0

= 4u3 − 10u .

2 http://users.encs.concordia.ca/ doedel/

3

Since gu(u, 0) = 4u3 − 10u , we have gu(−1, 0) = 6 , gu( 1, 0) = −6 , gu(−2, 0) = −12 , gu( 2, 0) = 12 , which are all nonzero . Thus each of the four solutions when λ = 0 is isolated . Hence each of these solutions persists as λ becomes nonzero, ( at least for “small” values of | λ | · · · ). 4

Solution families of g(u, λ) = 0 . Note the fold . 5

NOTE :

• Each of the four solutions at λ = 0 is isolated .
• Thus each of these solutions persists as λ becomes nonzero.
• Only two of the four homotopies reach λ = 1 .
• The other two homotopies meet at a fold .
• IFT condition (ii) is not satisfied at the fold.

(Why not? ) 6

In the equation G(u, λ) = 0 , u , G(·, ·) ∈ Rn , λ ∈ R , let x ≡

• u

λ

• .

Then the equation can be written G(x) = 0 , G : Rn+1 → Rn . DEFINITION : A solution x0 of G(x) = 0 is regular if the matrix G0

x ≡ Gx(x0) ,

(with n rows and n + 1 columns) has maximal rank , i.e., if Rank(G0

x) = n .

7

In the parameter formulation , G(u, λ) = 0 , we have Rank(G0

x) = Rank(G0 u | G0 λ) = n

⇐ ⇒                (i) G0

u is nonsingular,

• r

(ii)    dim N(G0

u) = 1 ,

and G0

λ ∈ R(G0 u) .

Here N(G0

u) denotes the null space of G0 u ,

and R(G0

u) denotes the range of G0 u ,

i.e., R(G0

u) is the linear space spanned by the n columns of G0 u .

8

COROLLARY (to the IFT) : Let x0 ≡ ( u0 , λ0 ) be a regular solution of G(x) = 0 . Then, near x0 , there exists a unique one-dimensional solution family x(s) with x(0) = x0 . PROOF : Since Rank( G0

x ) = Rank( G0 u | G0 λ ) = n , we have that

(i) either G0

u is nonsingular and by the IFT we have

u = u(λ) near x0 , (ii)

• r else we can interchange colums

in the Jacobian G0

x to see that the

solution can locally be parametrized by one of the components of u . Thus a (locally) unique solution family passes through x0 . QED ! 9

NOTE :

• Such a solution family is sometimes called a solution branch .
• Case (i) is where the IFT applies directly .
• Case (ii) is that of a simple fold .
• Thus even near a simple fold there is a unique solution family .
• However, near such a fold, the family cannot be parametrized by λ.

10

More Examples of IFT Application

• We give examples where the IFT shows that a given solution persists

(at least locally ) when a problem parameter is changed.

• We also consider cases where the conditions of the IFT are not satisfied .

11

EXAMPLE : The A → B → C Reaction . ( Course demo : Chemical-Reactions/ABC-Reaction/Stationary ) u′

1

= −u1 + D(1 − u1)eu3 , u′

2

= −u2 + D(1 − u1)eu3 − Dσu2eu3 , u′

3

= −u3 − βu3 + DB(1 − u1)eu3 + DBασu2eu3 , where 1 − u1 is the concentration of A , u2 is the concentration of B , u3 is the temperature, α = 1 , σ = 0.04 , B = 8 , D is the Damkohler number , β > 0 is the heat transfer coefficient . NOTE : The zero stationary solution at D = 0 persists (locally), because the Jacobian is nonsingular there, having eigenvalues −1, −1, and −(1 + β) . 12

Families of stationary solutions of the A → B → C reaction. (From left to right : β = 1.1 , 1.3 , 1.5 , 1.6 , 1.7 , 1.8 . ) 13

NOTE : In the preceding bifurcation diagram:

• u

=

• u2

1 + u2 2 + u2 3 .

• Solid/dashed curves denote stable/unstable solutions.
• The red squares are Hopf bifurcations .

From the basic theory of ODEs:

• u0 is a stationary solution of u′(t) = f(u(t)) if f(u0) = 0 .
• u0 is stable if all eigenvalues of fu(u0) are in the negative half-plane.
• u0 is unstable if one or more eigenvalues are in the positive half-plane.
• At a fold there is zero eigenvalue.
• At a Hopf bifurcation there is a pair of purely imaginary eigenvalues.

14

EXAMPLE (of IFT application) : The Gelfand-Bratu Problem . ( Course demo : Gelfand-Bratu/Original ) The boundary value problem    u′′(x) + λ eu(x) = 0 , ∀x ∈ [0, 1] , u(0) = u(1) = 0 , defines the stationary states of a solid fuel ignition model. If λ = 0 then u(x) ≡ 0 is a solution. This problem can be thought of as an operator equation G(u; λ) = 0 . We can use (a generalized) IFT to prove that there is a solution family u = u(λ) , for |λ| small . 15

The linearization of G(u; λ) acting on v , i.e., Gu(u; λ)v , leads to the homogeneous equation v′′(x) + λ eu(x)v = 0 , v(0) = v(1) = 0 , which for the solution u(x) ≡ 0 at λ = 0 becomes v′′(x) = 0 , v(0) = v(1) = 0 . Since this equation only has the zero solution v(x) ≡ 0 , the IFT applies. Thus (locally) a unique solution family passes through u(x) ≡ 0 , λ = 0 . 16

In Course demo : Gelfand-Bratu/Original the BVP is implemented as a first order system : u′

1(t) = u2(t) ,

u′

2(t) =

− λ eu1(t) , with boundary conditions u1(0) = 0 , u1(1) = 0 . A convenient solution measure in the bifurcation diagram is the value of 1 u1(x) dx . 17

Bifurcation diagram of the Gelfand-Bratu equation. 18

Some solutions of the Gelfand-Bratu equation. (The solution at the fold is colored red ). 19

EXAMPLE : A Boundary Value Problem with Bifurcations . ( Course demo : Basic-BVP/Nonlinear-Eigenvalue ) u′′ + ˆ λ u(1 − u) = 0 , u(0) = u(1) = 0 , has u(x) ≡ 0 as a solution for all ˆ λ . QUESTION : Are there more solutions ? Again, this problem corresponds to an operator equation G(u; ˆ λ) = 0 . Its linearization acting on v leads to the equation Gu(u; ˆ λ)v = 0 , i.e., v′′ + ˆ λ (1 − 2u)v = 0 , v(0) = v(1) = 0 . 20

In particular, the linearization about the zero solution family u ≡ 0 is v′′ + ˆ λ v = 0 , v(0) = v(1) = 0 , which for most values of ˆ λ only has the zero solution v(x) ≡ 0 . However, when ˆ λ = ˆ λk ≡ k2π2 , then there are nonzero solutions , namely, v(x) = sin(kπx) , Thus the IFT does not apply at ˆ λk = k2π2 . (We will see that these solutions are bifurcation points .) 21

In the implementation we write the BVP as a first order system . We also use a scaled version of λ. The equations are then u′

1 = u2 ,

u′

2 = λ2π2 u1 (1 − u1) ,

with ˆ λ = λ2π2 . A convenient solution measure in the bifurcation diagram is γ ≡ u2(0) = u′

1(0) .

22

Solution families to the nonlinear eigenvalue problem. 23

Some solutions to the nonlinear eigenvalue problem. 24

Hopf Bifurcation

THEOREM : Suppose that along a stationary solution family (u(λ), λ) , of u′ = f(u, λ) , a complex conjugate pair of eigenvalues α(λ) ± i β(λ) ,

• f fu(u(λ), λ) crosses the imaginary axis transversally , i.e., for some λ0 ,

α(λ0) = 0 , β(λ0) = 0 , and ˙ α(λ0) = 0 . Also assume that there are no other eigenvalues on the imaginary axis . Then there is a Hopf bifurcation, that is, a family of periodic solutions bifurcates from the stationary solution at (u0, λ0) . NOTE : The assumptions imply that f 0

u is nonsingular, so that the stationary

solution family is indeed (locally) a function of λ . 25

EXAMPLE : The A → B → C reaction . ( Course demo : Chemical-Reactions/ABC-Reaction/Homoclinic ) A stationary (blue) and a periodic (red) family of the A → B → C reaction for β = 1.2 . The periodic orbits are stable and terminate in a homoclinic orbit . 26

The periodic family orbit family approaching a homoclinic orbit (black). The red dot is the Hopf point; the blue dot is the saddle point on the homoclinic. 27

Course demo : Chemical-Reactions/ABC-Reaction/Periodic Bifurcation diagram for β = 1.1, 1.3, 1.5, 1.6, 1.7, 1.8 . (For periodic solutions u = 1

T

T

• u2

1 + u2 2 + u2 3 dt , where T is the period.)

28

EXAMPLE : A Predator-Prey Model . ( Course demo : Predator-Prey/ODE/2D )    u′

1 =

3u1(1 − u1) − u1u2 − λ(1 − e

−5u1) ,

u′

2 =

−u2 + 3u1u2 . Here u1 may be thought of as “fish ” and u2 as “sharks ”, while the term λ (1 − e

−5u1) ,

represents “fishing”, with “fishing-quota ” λ . When λ = 0 the stationary solutions are 3u1(1 − u1) − u1u2 = 0 −u2 + 3u1u2 = 0    ⇒ (u1, u2) = (0, 0) , (1, 0) , (1 3, 2) . 29

The Jacobian matrix is Gu(u1, u2 ; λ) =

• 3 − 6u1 − u2 − 5λe

−5u1

−u1 3u2 −1 + 3u1

• so that

Gu(0, 0 ; 0) = 3 −1

• ;

real eigenvalues 3 , -1 (unstable) Gu(1, 0 ; 0) = −3 −1 2

• ;

real eigenvalues -3 , 2 (unstable) Gu(1 3, 2 ; 0) = −1 −1

3

6

• ;

complex eigenvalues − 1 2 ± 1 2 √ 7 i (stable) All three Jacobians at λ = 0 are nonsingular . Thus, by the IFT, all three stationary points persist for (small) λ = 0 . 30

In this problem we can explicitly find all solutions: Family 1 : (u1, u2) = (0 , 0) Family 2 : u2 = 0 , λ = 3u1(1 − u1) 1 − e−5u1 ( Note that lim

u1 → 0 λ =

lim

u1 → 0

3(1 − 2u1) 5e−5u1 = 3 5 ) Family 3 : u1 = 1 3 , 2 3 − 1 3 u2 − λ(1−e−5/3) = 0 ⇒ u2 = 2−3λ(1−e−5/3) These solution families intersect at two bifurcation points , one of which is (u1, u2, λ) = (0 , 0 , 3/5) . 31

quota

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

sharks

0.0 0.5 1.0 1.5

fish

0.00 0.25 0.50 0.75

Stationary solution families of the predator-prey model. Solid/dashed curves denote stable/unstable solutions. Note the bifurcations and Hopf bifurcation (red square). 32

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

quota

0.0 0.2 0.4 0.6 0.8

fish

Stationary solution families, showing fish versus quota. Solid/dashed curves denote stable/unstable solutions. 33

Stability of Family 1 : Gu(0, 0 ; λ) = 3 − 5λ −1

• ;

eigenvalues 3 − 5λ, − 1 . Hence the zero solution is : unstable if λ < 3/5 , and stable if λ > 3/5 . Stability of Family 2 : This family has no stable positive solutions. 34

• Stability of Family 3 :

At λH ≈ 0.67 the complex eigenvalues cross the imaginary axis: − This crossing is a Hopf bifurcation , − Beyond λH there are stable periodic solutions . − Their period T increases as λ increases. − The period becomes infinite at λ = λ∞ ≈ 0.70 . − This final orbit is a heteroclinic cycle . 35

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

quota

−0.25 0.00 0.25 0.50 0.75 1.00

min fish, max fish

Stationary (blue) and periodic (red) solution families of the predator-prey model. ( For the periodic solution family both the maximum and minimum are shown. ) 36

Periodic solutions of the predator-prey model. The largest orbits are close to a heteroclinic cycle. 37

The bifurcation diagram shows the solution behavior for (slowly) increasing λ :

• Family 3 is followed until λH ≈ 0.67 .
• Periodic solutions of increasing period until λ = λ∞ ≈ 0.70 .
• Collapse to trivial solution (Family 1).

38

Continuation of Solutions

Parameter Continuation

Suppose we have a solution (u0, λ0) of G(u, λ) = 0 , as well as the derivative ˙ u0 . Here ˙ u ≡ du dλ . We want to compute the solution u1 at λ1 ≡ λ0 + ∆λ . 39

• "u"

u u λ λ λ u 1

(0) du d λ at λ0)

u

1 1

∆λ

(=

Graphical interpretation of parameter-continuation. 40

To solve the equation G(u1 , λ1) = 0 , for u1 (with λ = λ1 fixed) we use Newton’s method Gu(u(ν)

1 , λ1) ∆u(ν) 1

= − G(u(ν)

1 , λ1) ,

u(ν+1)

1

= u(ν)

1

+ ∆u(ν)

1

. ν = 0, 1, 2, · · · . As initial approximation use u(0)

1

= u0 + ∆λ ˙ u0 . If Gu(u1, λ1) is nonsingular , and ∆λ sufficiently small then this iteration will converge [B55]. 41

After convergence, the new derivative ˙ u1 is computed by solving Gu(u1, λ1) ˙ u1 = − Gλ(u1, λ1) . This equation is obtained by differentiating G(u(λ), λ) = 0 , with respect to λ at λ = λ1 . Repeat the procedure to find u2 , u3 , · · · . NOTE :

• ˙

u1 can be computed without another LU-factorization of Gu(u1, λ1) .

• Thus the extra work to compute

˙ u1 is negligible . 42

EXAMPLE : The Gelfand-Bratu Problem . u′′(x) + λ eu(x) = 0 for x ∈ [0, 1] , u(0) = 0 , u(1) = 0 . We know that if λ = 0 then u(x) ≡ 0 is an isolated solution . Discretize by introducing a mesh , 0 = x0 < x1 < · · · < xN = 1 , xj − xj−1 = h , (1 ≤ j ≤ N) , h = 1/N . The discrete equations are : uj+1 − 2uj + uj−1 h2 + λ euj = 0 , j = 1, · · · , N − 1 , with u0 = uN = 0 . 43

Let u ≡     u1 u2 · uN−1     . Then we can write the discrete equations as G( u , λ ) = 0 , where G : RN−1 × R → RN−1 . 44

Parameter-continuation : Suppose we have λ0 , u0 , and ˙ u0 . Set λ1 = λ0 + ∆λ . Newton’s method : Gu(u(ν)

1 , λ1) ∆u(ν) 1

= − G(u(ν)

1 , λ1) ,

u(ν+1)

1

= u(ν)

1

+ ∆u(ν)

1

, for ν = 0, 1, 2, · · · , with u(0)

1

= u0 + ∆λ ˙ u0 . After convergence compute ˙ u1 from Gu(u1, λ1) ˙ u1 = − Gλ(u1, λ1) . Repeat the procedure to find u2 , u3 , · · · . 45

Here Gu( u , λ ) =       − 2

h2 + λeu1 1 h2 1 h2

− 2

h2 + λeu2 1 h2

. . . . . .

1 h2

− 2

h2 + λeuN−1

      . Thus we must solve a tridiagonal system for each Newton iteration. NOTE :

• The solution family has a fold where parameter-continuation fails !
• A better continuation method is “pseudo-arclength continuation ”.
• There are also better discretizations, namely collocation , as used in AUTO .

46

Pseudo-Arclength Continuation

This method allows continuation of a solution family past a fold . It was introduced by H. B. Keller (1925-2008) in 1977. Suppose we have a solution (u0, λ0) of G( u , λ ) = 0 , as well as the normalized direction vector ( ˙ u0, ˙ λ0) of the solution family. Pseudo-arclength continuation consists of solving these equations for (u1, λ1) : G(u1, λ1) = 0 , u1 − u0 , ˙ u0 + (λ1 − λ0) ˙ λ0 − ∆s = 0 . 47

• u 0

u 0

∆ s

• λ 0
• "u"

λ λ λ u

1 1

u ( ) , λ 0

Graphical interpretation of pseudo-arclength continuation. 48

Solve the equations G(u1, λ1) = 0 , u1 − u0 , ˙ u0 + (λ1 − λ0) ˙ λ0 − ∆s = 0 . for (u1, λ1) by Newton’s method :   (G1

u)(ν)

(G1

λ)(ν)

˙ u∗ ˙ λ0  

• ∆u(ν)

1

∆λ(ν)

1

• = −

  G(u(ν)

1 , λ(ν) 1 )

u(ν)

1

− u0 , ˙ u0 + (λ(ν)

1

− λ0)˙ λ0 − ∆s   . Compute the next direction vector by solving   G1

u

G1

λ

˙ u∗ ˙ λ0   ˙ u1 ˙ λ1

• =

  1   , and normalize it. 49

NOTE :

• We can compute ( ˙

u1, ˙ λ1) with only one extra backsubstitution .

• The orientation of the family is preserved if ∆s is sufficiently small.
• Rescale the direction vector so that indeed ˙

u1 2 + ˙ λ2

1 = 1 .

50

FACT : The Jacobian is nonsingular at a regular solution. PROOF : Let x ≡

• u

λ

• ∈ Rn+1 .

Then pseudo-arclength continuation can be simply written as G(x1) = 0 , x1 − x0 , ˙ x0 − ∆s = 0 , ( ˙ x0 = 1 ) .

∆ s

• 1

x x x

Pseudo-arclength continuation. 51

The pseudo-arclength equations are G(x1) = 0 , x1 − x0 , ˙ x0 − ∆s = 0 , ( ˙ x0 = 1 ) . The Jacobian matrix in Newton’s method at ∆s = 0 is

• G0

x

˙ x∗

• .

At a regular solution N(G0

x) = Span{ ˙

x0} . We must show that

• G0

x

˙ x∗

• is nonsingular at a regular solution.

52

If on the contrary

• G0

x

˙ x∗

• is singular then for some vector z = 0 we have :

G0

x z

= 0 , ˙ x0 , z = 0 , Since by assumption N(G0

x) = Span{ ˙

x0} , we have z = c ˙ x0 , for some constant c . But then = ˙ x0 , z = c ˙ x0 , ˙ x0 = c ˙ x0 2 = c , so that z = 0 , which is a contradiction . QED ! 53

EXAMPLE : The Gelfand-Bratu Problem . Use pseudo-arclength continuation for the discretized Gelfand-Bratu problem. Then the matrix

• Gx

˙ x∗

• =
• Gu

Gλ ˙ u∗ ˙ λ

• ,

in Newton’s method is a bordered tridiagonal matrix :               ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆               . which can be decomposed very efficiently . 54

Following Folds and Hopf Bifurcations

• At a fold the the behavior of a system can change drastically.
• How does the fold location change when a second parameter varies ?
• Thus we want the compute a locus of folds in 2 parameters.
• We also want to compute loci of Hopf bifurcations in 2 parameters.

55

Following Folds

Treat both parameters λ and µ as unknowns , and compute a solution family X(s) ≡ ( u(s) , φ(s) , λ(s) , µ(s) ) , to F(X) ≡            G(u, λ, µ) = 0 , Gu(u, λ, µ) φ = 0 , φ , φ0 − 1 = 0 , and the added continuation equation u − u0 , ˙ u0 + φ − φ0 , ˙ φ0 + (λ − λ0)˙ λ0 + (µ − µ0) ˙ µ0 − ∆s = 0 . As before, ( ˙ u0 , ˙ φ0 , ˙ λ0 , ˙ µ0 ) , is the direction of the family at the current solution point ( u0 , φ0 , λ0 , µ0) . 56

EXAMPLE : The A → B → C Reaction . ( Course demo : Chemical-Reactions/ABC-Reaction/Folds-SS ) The equations are u′

1

= −u1 + D(1 − u1)eu3 , u′

2

= −u2 + D(1 − u1)eu3 − Dσu2eu3 , u′

3

= −u3 − βu3 + DB(1 − u1)eu3 + DBασu2eu3 , where 1 − u1 is the concentration of A , u2 is the concentration of B , u3 is the temperature , α = 1 , σ = 0.04 , B = 8 , D is the Damkohler number , β is the heat transfer coefficient . 57

A stationary solution family for β = 1.20. Note the two folds and the Hopf bifurcation . 58

0.00 0.05 0.10 0.15 0.20

D

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

β

0.14 0.15 0.16 0.17 0.18 0.19 0.20

D

1.20 1.25 1.30 1.35 1.40

β

A locus of folds (with blow-up) for the A → B → C reaction. Notice the two cusp singularities along the 2-parameter locus. ( There is a swallowtail singularity in nearby 3-parameter space. ) 59

0.12 0.14 0.16 0.18 0.20 0.22

D

1 2 3 4 5 6 7

||u||

Stationary solution families for β = 1.20, 1.21, · · · , 1.42. ( Open diamonds mark folds, solid red squares mark Hopf points. ) 60

Following Hopf Bifurcations

The extended system is F(u, φ, β, λ; µ) ≡            f(u, λ, µ) = 0 , fu(u, λ, µ) φ − i β φ = 0 , φ , φ0 − 1 = 0 , where F : Rn × Cn × R2 × R → Rn × Cn × C , and to which we want to compute a solution family ( u , φ , β , λ , µ ) , with u ∈ Rn, φ ∈ Cn, β, λ, µ ∈ R . Above φ0 belongs to a “reference solution ” ( u0 , φ0 , β0 , λ0 , µ0 ) , which normally is the latest computed solution along a family. 61

EXAMPLE : The A → B → C Reaction . ( Course demo : Chemical-Reactions/ABC-Reaction/Hopf ) The stationary family with Hopf bifurcation for β = 1.20 . 62

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

D

0.5 1.0 1.5 2.0

β

The locus of Hopf bifurcations for the A → B → C reaction. 63

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

D

1 2 3 4 5 6 7

||u|| Stationary solution families for β = 1.20, 1.20, 1.25, 1.30, · · · , 2.30, with Hopf bifurcations (the red squares) . 64

Boundary Value Problems

Consider the first order system of ordinary differential equations u′(t) − f( u(t) , µ , λ ) = 0 , t ∈ [0, 1] , where u(·) , f(·) ∈ Rn , λ ∈ R, µ ∈ Rnµ , subject to boundary conditions b( u(0) , u(1) , µ , λ ) = 0 , b(·) ∈ Rnb , and integral constraints 1 q( u(s) , µ , λ ) ds = 0 , q(·) ∈ Rnq . 65

This boundary value problem (BVP) is of the form F( X ) = 0 , where X = ( u , µ , λ ) , to which we add the continuation equation X − X0 , ˙ X0 − ∆s = 0 , where X0 represents the latest solution computed along the family. In detail , the continuation equation is 1 u(t) − u0(t) , ˙ u0(t) dt + µ − µ0 , ˙ µ0 + (λ − λ0)˙ λ0 − ∆s = 0 . 66

NOTE :

• In the context of continuation we solve this BVP for (u(·) , λ , µ) .
• In order for problem to be formally well-posed we must have

nµ = nb + nq − n ≥ 0 .

• A simple case is

nq = 0 , nb = n , for which nµ = 0 . 67

Discretization: Orthogonal Collocation

Introduce a mesh { 0 = t0 < t1 < · · · < tN = 1 } , where hj ≡ tj − tj−1 , (1 ≤ j ≤ N) , Define the space of (vector) piecewise polynomials Pm

h

as Pm

h

≡ { ph ∈ C[0, 1] : ph

• [tj−1,tj] ∈ Pm } ,

where Pm is the space of (vector) polynomials of degree ≤ m . 68

The collocation method consists of finding ph ∈ Pm

h ,

µ ∈ Rnµ , such that the following collocation equations are satisfied : p′

h(zj,i) = f( ph(zj,i) , µ, λ ) ,

j = 1, · · · , N , i = 1, · · · , m , and such that ph satisfies the boundary and integral conditions . The collocation points zj,i in each subinterval [ tj−1 , tj ] , are the (scaled) roots of the mth-degree orthogonal polynomial (Gauss points3).

3 See Pages 261, 287 of the Background Notes on Elementary Numerical Methods.

69

• 1

t t t t t

1 2 N

zj,1 zj,2 zj,3 t t t j-1

j-1 j j

l l

j,3 j,1

(t) (t)

j-1/3

t

j-2/3

t

The mesh {0 = t0 < t1 < · · · < tN = 1} , with collocation points and extended-mesh points shown for m = 3 . Also shown are two of the four local Lagrange basis polynomials . 70

Since each local polynomial is determined by (m + 1) n , coefficients, the total number of unknowns (considering λ as fixed) is (m + 1) n N + nµ . This is matched by the total number of equations : collocation : m n N , continuity : (N − 1) n , constraints : nb + nq ( = n + nµ ) . 71

Assume that the solution u(t) of the BVP is sufficiently smooth. Then the order of accuracy of the orthogonal collocation method is m , i.e., ph − u ∞ = O(hm) . At the main meshpoints tj we have superconvergence : maxj | ph(tj) − u(tj) | = O(h2m) . The scalar variables λ and µ are also superconvergent . 72

Implementation

For each subinterval [ tj−1 , tj ] , introduce the Lagrange basis polynomials { ℓj,i(t) } , j = 1, · · · , N , i = 0, 1, · · · , m , defined by ℓj,i(t) =

m

• k=0,k=i

t − tj− k

m

tj− i

m − tj− k m

, where tj− i

m ≡ tj −

i m hj . The local polynomials can then be written pj(t) =

m

• i=0

ℓj,i(t) uj− i

m .

With the above choice of basis uj ∼ u(tj) and uj− i

m ∼ u(tj− i m) ,

where u(t) is the solution of the continuous problem. 73

The collocation equations are p

′

j(zj,i) = f( pj(zj,i) , µ , λ ) ,

i = 1, · · · , m, j = 1, · · · , N . The boundary conditions are bi( u0 , uN , µ , λ ) = 0 , i = 1, · · · , nb . The integral constraints can be discretized as

N

• j=1

m

• i=0

ωj,i qk( uj− i

m , µ , λ) = 0 ,

k = 1, · · · , nq , where the ωj,i are the Lagrange quadrature weights . 74

The continuation equation is 1 u(t) − u0(t) , ˙ u0(t) dt + µ − µ0 , ˙ µ0 + (λ − λ0) ˙ λ0 − ∆s = 0 , where ( u0 , µ0 , λ0 ) , is the previous solution along the solution family, and ( ˙ u0 , ˙ µ0 , ˙ λ0 ) , is the normalized direction of the family at the previous solution . The discretized continuation equation is of the form

N

• j=1

m

• i=0

ωj,i uj− i

m − (u0)j− i m ,

( ˙ u0)j− i

m

+ µ − µ0 , ˙ µ0 + (λ − λ0) ˙ λ0 − ∆s = 0 . 75

Numerical Linear Algebra

The complete discretization consists of m n N + nb + nq + 1 , nonlinear equations , with unknowns {uj− i

m} ∈ RmnN+n ,

µ ∈ Rnµ , λ ∈ R . These equations are solved by a Newton-Chord iteration . 76

We illustrate the numerical linear algebra for the case n = 2 ODEs , N = 4 mesh intervals , m = 3 collocation points , nb = 2 boundary conditions , nq = 1 integral constraint , and the continuation equation.

• The operations are also done on the right hand side , which is not shown.
• Entries marked “◦” have been eliminated by Gauss elimination.
• Entries marked “·” denote fill-in due to pivoting .
• Most of the operations can be done in parallel .

77

u0 u 1

3

u 2

3

u1 u2 u3 uN µ λ

• The structure of the Jacobian .

78

u0 u 1

3

u 2

3

u1 u2 u3 uN µ λ

• The system after condensation of parameters, which can be done in parallel .

79

u0 u 1

3

u 2

3

u1 u2 u3 uN µ λ

• The preceding matrix, showing the decoupled • subsystem .

80

u0 u 1

3

u 2

3

u1 u2 u3 uN µ λ

• ·

·

• ·

·

• ·

·

• ·

·

• Stage 1 of the nested dissection to solve the decoupled • subsystem.

81

u0 u 1

3

u 2

3

u1 u2 u3 uN µ λ

• ·

·

• ·

·

• ·

·

• ·

·

• ·

·

• ·

·

• Stage 2 of the nested dissection to solve the decoupled • subsystem.

82

u0 u 1

3

u 2

3

u1 u2 u3 uN µ λ

• ·

·

• ·

·

• ·

·

• ·

·

• ·

·

• ·

·

• The preceding matrix showing the final decoupled • subsystem .

83

u0 u 1

3

u 2

3

u1 u2 u3 uN µ λ

• ·

·

• ·

·

• ·

·

• ·

·

• ·

·

• ·

·

• A

A

• B

B

• A

A

• B

B

• The approximate Floquet multipliers are the eigenvalues of M ≡ −B−1A .

84

Accuracy Test

The Table shows the location of the fold in the Gelfand-Bratu problem, for 4 Gauss collocation points per mesh interval, and N mesh intervals . N Fold location 2 3.5137897550 4 3.5138308601 8 3.5138307211 16 3.5138307191 32 3.5138307191 85

Periodic Solutions

• Periodic solutions can be computed efficiently using a BVP approach.
• This method also determines the period very accurately.
• Moreover, the technique can compute unstable periodic orbits.

86

Consider u′(t) = f( u(t) , λ ) , u(·) , f(·) ∈ Rn , λ ∈ R . Fix the interval of periodicity by the transformation t → t T . Then the equation becomes u′(t) = T f( u(t) , λ ) , u(·) , f(·) ∈ Rn , T , λ ∈ R . and we seek solutions of period 1 , i.e., u(0) = u(1) . Note that the period T is one of the unknowns . 87

The above equations do not uniquely specify u and T : Assume that we have computed ( uk−1(·) , Tk−1 , λk−1 ) , and we want to compute the next solution ( uk(·) , Tk , λk ) . Then uk(t) can be translated freely in time : If uk(t) is a periodic solution, then so is uk(t + σ) , for any σ . Thus, a phase condition is needed. 88

An example is the Poincar´ e phase condition uk(0) − uk−1(0) , u

′

k−1(0) = 0 .

( But we will derive a numerically more suitable integral phase condition . ) u k-1 (0)

• u k-1 (0)

u (0)

k

Graphical interpretation of the Poincar´ e phase condition. 89

An Integral Phase Condition If ˜ uk(t) is a solution then so is ˜ uk(t + σ) , for any σ . We want the solution that minimizes D(σ) ≡ 1 ˜ uk(t + σ) − uk−1(t) 2

2 dt .

The optimal solution ˜ uk(t + ˆ σ) , must satisfy the necessary condition D′(ˆ σ) = 0 . 90

Differentiation gives the necessary condition 1 ˜ uk(t + ˆ σ) − uk−1(t) , ˜ u′

k(t + ˆ

σ dt = 0 . Writing uk(t) ≡ ˜ uk(t + ˆ σ) , gives 1 uk(t) − uk−1(t) , u′

k(t) dt = 0 .

Integration by parts, using periodicity, gives 1

0 uk(t) , u

′

k−1(t) dt = 0 .

This is the integral phase condition. 91

Continuation of Periodic Solutions

• Pseudo-arclength continuation is used to follow periodic solutions .
• It allows computation past folds along a family of periodic solutions.
• It also allows calculation of a “vertical family ” of periodic solutions.

For periodic solutions the continuation equation is 1 uk(t)−uk−1(t) , ˙ uk−1(t) dt + (Tk −Tk−1) ˙ Tk−1 + (λk −λk−1)˙ λk−1 = ∆s . 92

SUMMARY : We have the following equations for periodic solutions : u′

k(t) = T f( uk(t) , λk ) ,

uk(0) = uk(1) , 1 uk(t) , u

′

k−1(t) dt = 0 ,

with continuation equation 1 uk(t)−uk−1(t) , ˙ uk−1(t) dt + (Tk −Tk−1) ˙ Tk−1 + (λk −λk−1)˙ λk−1 = ∆s , where u(·) , f(·) ∈ Rn , λ , T ∈ R . 93

Stability of Periodic Solutions

In our continuation context, a periodic solution of period T satisfies u′(t) = T f( u(t)) , for t ∈ [0, 1] , u(0) = u(1) , (for given value of the continuation parameter λ). A small perturbation in the initial condition u(0) + ǫ v(0) , ǫ small , leads to the linearized equation v′(t) = T fu( u(t) ) v(t) , for t ∈ [0, 1] , which induces a linear map v(0) → v(1) , represented by v(1) = M v(0) . 94

v(1) = M v(0) The eigenvalues of M are the Floquet multipliers that determine stability. M always has a multiplier µ = 1 , since differentiating u′(t) = T f( u(t)) , gives v′(t) = T fu( u(t) ) v(t) , where v(t) = u′(t) , with v(0) = v(1) . 95

v(1) = M v(0)

• If M has a Floquet multiplier µ with | µ |

> 1 then u(t) is unstable .

• If all multipliers (other than µ = 1) satisfy | µ |

< 1 then u(t) is stable .

• At folds and branch points there are two multipliers µ = 1 .
• At a period-doubling bifurcation there is a real multiplier µ = −1 .
• At a torus bifurcation there is a complex pair on the unit circle.

96

EXAMPLE : The Lorenz Equations . ( Course demo : Lorenz ) These equations were introduced in 1963 by Edward Lorenz (1917-2008) as a simple model of atmospheric convection : x′ = σ (y − x) , y′ = ρ x − y − x z , z′ = x y − β z , where (often) σ = 10 , β = 8/3 , ρ = 28 . 97

Course demo : Lorenz/Basic Bifurcation diagram of the Lorenz equations for σ = 10 and β = 8/3 . 98

Course demo : Lorenz/Basic Unstable periodic orbits of the Lorenz equations. 99

In the Lorenz Equations :

• The zero stationary solution is unstable for ρ > 1 .
• Two nonzero stationary families bifurcate at ρ = 1 .
• The nonzero stationary solutions are unstable for ρ > ρH .
• At ρH ≈ 24.7 there are Hopf bifurcations .
• Unstable periodic solution families emanate from the Hopf bifurcations.
• These families end in homoclinic orbits (infinite period) at ρ ≈ 13.9 .
• At ρ = 28 (and a range of other values) there is the Lorenz attractor .

100

EXAMPLE : The A → B → C Reaction . ( Course demo : Chemical-Reactions/ABC-Reaction/Periodic ) Stationary and periodic solution families of the A → B → C reaction: β = 1.55 . Note the coexistence of stable solutions, for example, solutions 1 and 2 . 101

0.20 0.25 0.30 0.35

D

2 3 4 5 6 7 8 9

max u3

0.20 0.25 0.30 0.35

D

2 3 4 5 6 7 8 9

max u3

0.20 0.25 0.30 0.35

D

2 3 4 5 6 7 8 9

max u3

0.20 0.25 0.30 0.35

D

2 3 4 5 6 7 8 9

max u3

Top left: β = 1.55, right: β = 1.56, Bottom left: β = 1.57, right: β = 1.58. ( QUESTION : Is something missing somewhere ? ) 102

Following Folds for Periodic Solutions

Recall that periodic orbits families can be computed using the equations u′(t) − T f( u(t) , λ ) = 0 , u(0) − u(1) = 0 , 1 u(t) , u

′

0(t) dt = 0 ,

where u0 is a reference orbit , typically the latest computed orbit. The above boundary value problem is of the form F( X , λ ) = 0 , where X = ( u , T ) . 103

At a fold with respect to λ we have FX( X , λ ) Φ = 0 , Φ , Φ = 1 , where X = ( u , T ) , Φ = ( v , S ) , i.e., the linearized equations about X have a nonzero solution Φ . In detail : v′(t) − T fu( u(t) , λ ) v − S f( u(t) , λ ) = 0 , v(0) − v(1) = 0 , 1 v(t) , u

′

0(t) dt = 0 ,

1 v(t) , v(t) dt + S2 = 1 . 104

The complete extended system to follow a fold is F( X , λ , µ ) = 0 , FX( X , λ , µ ) Φ = 0 , Φ , Φ − 1 = 0 , with two free problem parameters λ and µ . To the above we add the continuation equation X−X0 , ˙ X0 + Φ−Φ0 , ˙ Φ0 + (λ−λ0) ˙ λ0 + (µ−µ0) ˙ µ0 − ∆s = 0 . 105

In detail : u′(t) − T f( u(t) , λ , µ ) = 0 , u(0) − u(1) = 0 , 1 u(t) , u

′

0(t) dt = 0 ,

v′(t) − T fu( u(t) , λ , µ ) v − S f( u(t) , λ , µ ) = 0 , v(0) − v(1) = 0 , 1 v(t) , u

′

0(t) dt = 0 ,

with normalization 1 v(t) , v(t) dt + S2 − 1 = 0 , and continuation equation 1 u(t)−u0(t) , ˙ u0(t) dt + 1 v(t)−v0(t) , ˙ v0(t) dt + + (T0 − T) ˙ T0 + (S0 − S) ˙ S0 + (λ − λ0)˙ λ0 + (µ − µ0) ˙ µ0 − ∆s = 0 . 106

EXAMPLE : The A → B → C Reaction . ( Course demo : Chemical-Reactions/ABC-Reaction/Folds-PS ) Stationary and periodic solution families of the A → B → C reaction . (with blow-up) for β = 1.55 . Note the three folds , labeled 1, 2, 3 . 107

0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28

D

1.54 1.55 1.56 1.57 1.58 1.59 1.60 1.61 1.62

β

Loci of folds along periodic solution families for the A → B → C reaction. 108

Stationary and periodic solution families of the A → B → C reaction: β = 1.56 . 109

0.15 0.20 0.25 0.30 0.35 0.40

D

1 2 3 4 5 6 7 8 9

max u3 Stationary and periodic solution families of the A → B → C reaction: β = 1.57 . 110

Stationary and periodic solution families of the A → B → C reaction: β = 1.58 . 111

Stationary and periodic solution families of the A → B → C reaction: β = 1.61 . 112

Stationary and periodic solution families of the A → B → C reaction: β = 1.62 . 113

Periodic solutions along the isola for β = 1.58 . (Stable solutions are blue, unstable solutions are red.) 114

Following Period-doubling Bifurcations

Let ( u(t) , T ) be a periodic solution , i.e., a solution of u′(t) − T f( u(t) , λ ) = 0 , u(0) − u(1) = 0 , 1 u(t) , u

′

0(t) dt = 0 ,

where u0 is a reference orbit . A necessary condition for a period-doubling bifurcation is that the following linearized system have a nonzero solution v(t) : v′(t) − T fu( u(t) , λ ) v(t) = 0 , v(0) + v(1) = 0 , 1 v(t) , v(t) dt = 1 . 115

The complete extended system to follow a period-doubling bifurcation is u′(t) − T f( u(t) , λ , µ ) = 0 , u(0) − u(1) = 0 , 1 u(t) , u

′

0(t) dt = 0 ,

v′(t) − T fu( u(t) , λ ) v(t) = 0 , v(0) + v(1) = 0 , 1 v(t) , v(t) dt − 1 = 0 , and continuation equation 1 u(t)−u0(t) , ˙ u0(t) dt + 1 v(t)−v0(t) , ˙ v0(t) dt + + (T0 − T) ˙ T0 + (λ − λ0)˙ λ0 + (µ − µ0) ˙ µ0 − ∆s = 0 . 116

EXAMPLE : Period-Doubling Bifurcations in the Lorenz Equations . ( Course demo : Lorenz/Period-Doubling )

• The Lorenz equations also have period-doubling bifurcations .
• In fact, there is a period-doubling cascade for large ρ .
• We start from numerical data .
• (Such data may be from simulation , i.e., initial value integration .)
• We also want to compute loci of period-doubling bifurcations .

117

( Course demo : Lorenz/Period-Doubling ) Left panel : Solution families of the Lorenz equations. The open diamonds denote period-doubling bifurcations . Right panel : Solution 1 was found by initial value integration . 118

( Course demo : Lorenz/Period-Doubling ) Left panel : A primary period-doubled solution. Right panel : A secondary period-doubled solution. 119

( Course demo : Lorenz/Period-Doubling ) Loci of period-doubling bifurcations for the Lorenz equations (with blow-up) . Black: primary, Red: secondary, Blue: tertiary period-doublings . 120

Periodic Solutions of Conservative Systems

EXAMPLE : A Model Conservative System . ( Course demo : Vertical-HB ) u′

1 =

− u2 , u′

2 = u1 (1 − u1) .

PROBLEM :

• This system has a family of periodic solutions, but no parameter !
• The system has a constant of motion , namely the Hamiltonian

H(u1, u2) = − 1 2 u2

1 − 1

2 u2

2 + 1

3 u3

1 .

121

REMEDY : Introduce an unfolding term with unfolding parameter λ : u′

1 = λ u1 − u2 ,

u′

2 = u1 (1 − u1) .

Then there is a vertical Hopf bifurcation from the trivial solution at λ = 0 . 122

Bifurcation diagram of the vertical Hopf bifurcation problem. ( Course demo : Vertical-HB ) 123

NOTE :

• The family of periodic solutions is vertical .
• The parameter λ is solved for in each continuation step.
• Upon solving, λ is found to be zero , up to numerical precision.
• One can use standard BVP continuation and bifurcation algorithms.

124

A phase plot of some periodic solutions. 125

EXAMPLE : The Circular Restricted 3-Body Problem (CR3BP) . ( Course demo : Restricted-3Body/Earth-Moon/Orbits ) x′′ = 2y′ + x − (1 − µ) (x + µ) r3

1

− µ (x − 1 + µ) r3

2

, y′′ = −2x′ + y − (1 − µ) y r3

1

− µ y r3

2

, z′′ = −(1 − µ) z r3

1

− µ z r3

2

, where r1 =

• (x + µ)2 + y2 + z2 ,

r2 =

• (x − 1 + µ)2 + y2 + z2 .

and ( x , y , z ) denotes the position of the zero-mass body . NOTE : For the Earth-Moon system µ ≈ 0.01215 . 126

The CR3BP has one integral of motion , namely, the “Jacobi-constant” : J = x′2 + y′2 + z′2 2 − U(x, y, z) − µ 1 − µ 2 , where U = 1 2(x2 + y2) + 1 − µ r1 + µ r2 , and r1 =

• (x + µ)2 + y2 + z2 ,

r2 =

• (x − 1 + µ)2 + y2 + z2 .

127

Boundary value formulation : x′ = T vx y′ = T vy z′ = T vz v′

x

= T [ 2vy + x − (1 − µ)(x + µ)r−3

1

− µ(x − 1 + µ)r−3

2

+ λ vx ] v′

y

= T [ − 2vx + y − (1 − µ)yr−3

1

− µyr−3

2

+ λ vy ] v′

z

= T [ − (1 − µ)zr−3

1

− µzr−3

2

+ λ vz ] with periodicity boundary conditions x(1) = x(0) , y(1) = y(0) , z(1) = z(0) , vx(1) = vx(0) , vy(1) = vy(0) , vz(1) = vz(0) , + phase constraint + continuation equation . Here T is the period of the orbit. 128

NOTE :

• One can use BVP continuation and bifurcation algorithms.
• The unfolding term λ ∇v regularizes the continuation.
• λ will be zero , once solved for.
• Other unfolding terms are possible.

129

Schematic bifurcation diagram of periodic solution families of the Earth-Moon system . 130

The planar Lyapunov family L1. 131

The Halo family H1. 132

The Halo family H1. 133

The Vertical family V1. 134

The Axial family A1. 135

Stable and Unstable Manifolds

EXAMPLE : Phase-plane orbits: Fixed length . These can be computed by orbit continuation . Model equations are x′ = ǫx − y3 , y′ = y + x3 . where ǫ > 0 is small.

• There is only one equilibrium , namely, (x, y) = (0, 0) .
• This equilibrium has eigenvalues ǫ and 1 ; it is a source .

136

For the computations :

• The time variable t is scaled to [0, 1].
• The actual integration time T is then an explicit parameter :

x′ = T ( ǫx − y3 ) , y′ = T ( y + x3 ) . 137

These constraints are used :

• To put the initial point on a small circle around the origin :

x(0) − r cos(2πθ) = 0 , y(0) − r sin(2πθ) = 0 .

• To keep track of the end points :

x(1) − x1 = 0 , y(1) − y1 = 0 .

• To keep track of the length of the orbits

1

• x′(t)2 + y′(t)2 dt − L = 0 .

138

The computations are done in 3 stages :

• In the first run an orbit is grown by continuation :
• The free parameters are T, L, x1, y1 .
• The starting point is on the small circle of radius r.
• The starting point is in the strongly unstable direction .
• The value of ǫ is 0.5 in the first run.
• In the second run the value of ǫ is decreased to 0.01 :
• The free parameters are ǫ, T, x1, y1 .
• In the third run the initial point is free to move around the circle :
• The free parameters are θ, T, x1, y1 .

139

( Course demo : Basic-Manifolds/2D-ODE/Fixed-Length ) Unstable Manifolds in the Plane: Orbits of Fixed Length . (The right-hand panel is a blow-up, and also shows fewer orbits.) 140

EXAMPLE : Phase-plane orbits: Variable length . These can also be computed by orbit continuation . Model equations are x′ = ǫx − y2 , y′ = y + x2 .

• The origin (x, y) = (0, 0) is an equilibrium .
• The origin has eigenvalues ǫ and 1 ; it is a source .
• Thus the origin has a 2-dimensional unstable manifold .
• We compute this stable manifold using continuation .
• (The equations are 2D; so we actually compute a phase portrait.)

141

For the computations :

• The time variable t is scaled to [0, 1].
• The actual integration time T is then an explicit parameter :

x′ = T(ǫx − y2) , y′ = T(y + x2) . NOTE :

• There is also a nonzero equilibrium

(x, y) = (ǫ

1 3, −ǫ 2 3) .

• It is a saddle (1 positive, 1 negative eigenvalue).

142

These constraints are used :

• To put the initial point on a small circle at the origin :

x(0) − r cos(2πθ) = 0 , y(0) − r sin(2πθ) = 0 .

• To keep track of the end points :

x(1) − x1 = 0 , y(1) − y1 = 0 .

• To keep track of the length of the orbits we add an integral constraint :

1

• x′(t)2 + y′(t)2 dt − L = 0 .
• To allow the length L to contract :

(Tmax − T)(Lmax − L) − c = 0 . 143

Again the computations are done in 3 stages :

• In the first run an orbit is grown by continuation :
• The free parameters are T, L, x1, y1, c .
• The starting point is on a small circle of radius r.
• The starting point is in the strongly unstable direction .
• In this first run ǫ = 0.5 .
• In the second run the value of ǫ is decreased to 0.05 :
• The free parameters are ǫ, T, L, x1, y1 .
• In the third run the initial point is free to move around the circle :
• The free parameters are θ, T, L, x1, y1 .

144

( Course demo : Basic-Manifolds/2D-ODE/Variable-Length ) Unstable Manifolds in the Plane: Orbits of Variable Length . 145

( Course demo : Basic-Manifolds/2D-ODE/Variable-Length ) Unstable Manifolds in the Plane: Orbits of Variable Length (Blow-up). 146

EXAMPLE : A 2D unstable manifold in R3 . This can also be computed by orbit continuation. The model equations are x′ = ǫx − z3 , y′ = y − x3 , z′ = −z + x2 + y2 .

• We take ǫ = 0.05.
• The origin is a saddle with eigenvalues ǫ, 1, and −1.
• Thus the origin has a 2-dimensional unstable manifold .
• The initial point moves around a circle in the 2D unstable eigenspace .
• The equations are 3D; so we will compute a 2D manifold in R3 .
• There is also a nonzero saddle , so we use retraction .
• The set-up is similar to the 2D phase-portrait example.

147

( Course demo : Basic-Manifolds/3D-ODE/Variable-Length ) Unstable Manifolds in R3: Orbits of Variable Length . 148

EXAMPLE : Another 2D unstable manifold in R3 . The model equations are x′ = ǫx − y3 + z3 , y′ = y + x3 + z3 , z′ = −z − x2 + y2 .

• We take ǫ = 0.05.
• The origin is a saddle with eigenvalues ǫ, 1, and −1.
• Thus the origin has a 2-dimensional unstable manifold .
• The initial point moves around a circle in the 2D unstable eigenspace .
• The equations are 3D; so we will compute a 2D manifold in R3 .
• No retraction is needed, so we choose to compute orbits of fixed length .
• The set-up is similar to the 2D phase-portrait example.

149

( Course demo : Basic-Manifolds/3D-ODE/Fixed-Length ) Unstable Manifolds in R3: Orbits of Fixed Length . 150

The Lorenz Manifold

• For ρ > 1 the origin is a saddle point .
• The Jacobian has two negative and one positive eigenvalue .
• The two negative eigenvalues give rise to a 2D stable manifold .
• This manifold is known as as the Lorenz Manifold .
• The set-up is as for the earlier 3D model, using fixed length .

151

Course demo : Lorenz/Manifolds/Origin/Fixed-Length Part of the Lorenz Manifold (with blow-up). Orbits have fixed length L = 60 . 152

Course demo : Lorenz/Manifolds/Origin/Fixed-Length Part of the Lorenz Manifold. Orbits have fixed length L = 200 . 153

Heteroclinic Connections.

• The Lorenz Manifold helps understand the Lorenz attractor .
• Many orbits in the manifold depend sensitively on initial conditions .
• During the manifold computation one can locate heteroclinic orbits .
• These are also in the 2D unstable manifold of the nonzero equilibria.
• The heteroclinic orbits have a combinatorial structure 4.
• One can also continue heteroclinic orbits as ρ varies.

4 Nonlinearity 19, 2006, 2947-2972.

154

Course demo : Lorenz/Heteroclinics Four heteroclinic orbits with very close initial conditions 155

One can also determine the intersection of the Lorenz manifold with a sphere . The set-up is as follows : x′ = T σ (y − x) , y′ = T (ρ x − y − x z) , z′ = T (x y − β z) , which is of the form u′(t) = T f( u(t) ) , for 0 ≤ t ≤ 1 , where

• T is the actual integration time , which is negative !

To this we add boundary and integral constraints . 156

The complete set-up consists of the ODE u′(t) = T f( u(t) ) , for 0 ≤ t ≤ 1 , subject to the following constraints : u(0) − ǫ (cos(θ) v1 − sin(θ) v2) = 0 u(0) is on a small circle u(1) − u1 = 0 to keep track of the end point u(1) u1 − R = 0 distance of u1 to the origin u1/ u1 , f(u1)/ f(u1) − τ = 0 to locate tangencies, where τ = 0 T 1 f(u) ds − L = 0 to keep track of the orbit length (T − Tn) (L − Ln) − c = 0 . allows retraction into the sphere 157

The continuation system has the form F(Xk) = 0 , where X = ( u(·) , Λ ) . with continuation equation Xk − Xk−1 , ˙ Xk−1 − ∆s = 0 , ( ˙ Xk−1 = 1 ) . The computations are done in 2 stages :

• In the first run an orbit is grown by continuation :
• The starting point is on the small circle of radius ǫ.
• The starting point is in the strongly stable direction .
• The free parameters are Λ = (T , L , c , τ , R , u1) .
• In the second run the orbit sweeps the stable manifold.
• The initial point is free to move around the circle :
• The free parameters are Λ = (T , L , θ , τ , R , u1) .

158

Course demo : Lorenz/Manifolds/Origin/Sphere Intersection of the Lorenz Manifold with a sphere 159

NOTE :

• We do not just change the initial point (i.e., θ) and integrate !
• Every continuation step requires solving a boundary value problem .
• The continuation stepsize ∆s controls the change in X .
• X can only change a little in any continuation step.
• This way the entire manifold (up to a given length L) is computed.
• The retraction constraint allows the orbits to retract into the sphere.
• This is necessary when heteroclinic connections are encountered.

160

EXAMPLE : Unstable Manifolds of a Periodic Orbit . ( Course demo : Lorenz/Manifolds/Orbits/Rho24.3 ) Left: Bifurcation diagram of the Lorenz equations. Right: Labeled solutions. 161

Both sides of the unstable manifold of periodic orbit 3 at ρ = 24.3 . 162

EXAMPLE : Unstable Manifolds in the CR3BP . ( Course demo : Restricted-3Body/Earth-Moon/Manifolds/H1 )

• ”Small” Halo orbits have one real Floquet multiplier outside the unit circle.
• Such Halo orbits are unstable .
• They have a 2D unstable manifold .
• The unstable manifold can be computed by continuation .
• First compute a starting orbit in the manifold.
• Then continue the orbit keeping, for example, x(1) fixed .

163

Part of the unstable manifold of three Earth-Moon L1-Halo orbits. 164

• The initial orbit can be taken to be much longer

· · ·

• Continuation with x(1) fixed can lead to a Halo-to-torus connection!

165

• The Halo-to-torus connection can be continued as a solution to

F( Xk ) = 0 , < Xk − Xk−1 , ˙ Xk−1 > − ∆s = 0 . where X = ( Halo orbit , Floquet function , connecting orbit) . 166

In detail , the continuation system is du dτ − Tuf(u(τ), µ, l) = 0 , u(1) − u(0) = 0 , 1 u(τ) , ˙ u0(τ) dτ = 0 , dv dτ − TuDuf(u(τ), µ, l)v(τ) + λuv(τ) = 0 , v(1) − sv(0) = 0 (s = ±1) , v(0) , v(0) − 1 = 0 , dw dτ − Twf(w(τ), µ, 0) = 0 , w(0) − (u(0) + εv(0)) = 0 , w(1)x − xΣ = 0 . 167

The system has 18 ODEs , 20 boundary conditions , 1 integral constraint . We need 20 + 1 + 1 - 18 = 4 free parameters . Parameters :

• An orbit in the unstable manifold:

Tw , l , Tu , xΣ

• Compute the unstable manifold:

Tw , l , Tu , ε

• Follow a connecting orbit:

λu , l , Tu , ε 168

169

170

171

The Solar Sail Equations

The equations in Course demo : Solar-Sail/Equations/equations.f90 : x′′ = 2y′ + x − (1 − µ)(x + µ) d3

S

− µ(x − 1 + µ) d3

P

+ β(1 − µ)D2Nx d2

S

y′′ = − 2x′ + y − (1 − µ)y d3

S

− µy d3

P

+ β(1 − µ)D2Ny d2

S

z′′ = − (1 − µ)z d3

S

− µz d3

P

+ β(1 − µ)D2Nz d2

S

where dS =

• (x + µ)2 + y2 + z2 , dP =
• (x − 1 + µ)2 + y2 + z2 , r =
• (x + µ)2 + y2

Nx = [cos(α)(x + µ) − sin(α)y] [cos(δ) − sin(δ)z r ]/dS Ny = [cos(α)y + sin(α)(x + µ)] [cos(δ) − sin(δ)z r ]/dS Nz = [cos(δ)z + sin(δ)r]/dS , D = x + µ dS Nx + y dS Ny + z dS Nz 172

Course demo : Solar-Sail/Sun-Jupiter/Libration/Points Sun-Jupiter libration points, for β = 0, α = 0, δ = 0. 173

Course demo : Solar-Sail/Sun-Jupiter/Libration/Points Sun-Jupiter libration points, for β = 0.02, α = 0.02, δ = 0. 174

Course demo : Solar-Sail/Sun-Jupiter/Libration/Loci Sun-Jupiter libration points, with δ ∈ [−π

2, π 2], for various β, with α = 0.

175

Course demo : Solar-Sail/Sun-Jupiter/Libration/Loci Sun-Jupiter libration points, with δ ∈ [−π

2, π 2], for various α, with β = 0.15.

176

Course demo : Solar-Sail/Sun-Jupiter/Libration/Homoclinic Sun-Jupiter: detection of a homoclinic orbit at β = 0.050698, with α = 0, δ = 0. 177

Course demo : Solar-Sail/Sun-Jupiter/Libration/Manifolds Sun-Jupiter: unstable manifold orbits for δ ∈ [−π

2, π 2], with β = 0.05, α = 0.1.

178

Course demo : Solar-Sail/Sun-Jupiter/Libration/Manifolds

0.926 0.927 0.928 0.929 0.930 0.931 0.932 0.933

x

−0.0075 −0.0050 −0.0025 0.0000 0.0025 0.0050 0.0075

z

−1.05 −1.00 −0.95 −0.90 −0.85 −0.80 −0.75 −0.70 −0.65

y(T)

−0.04 −0.03 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04

z(T)

The libration points The end points 179

Course demo : Solar-Sail/Sun-Jupiter/Libration/Manifolds Some connecting orbits for α = 0.07 and varying β and δ. 180

Course demo : Solar-Sail/Sun-Jupiter/Orbits V1-orbits with β = 0.15, T = 6.27141, δ ∈ [0 , 0.6415] . 181