Geometry and numerical continuation of multiscale orbits in a nonconvex variational problem
Annalisa Iuorio Christian Kuehn Peter Szmolyan




Annalisa Iuorio]Annalisa Iuorio

Institute for Analysis and Scientific Computing, Vienna University of Technology

Wiedner Hauptstraße 8-10,

1040, Vienna

Austria
annalisa.iuorio@tuwien.ac.at


Faculty of Mathematics, Technical University of MunichBoltzmannstraße 385748 Garching bei München,Germany
ckuehn@ma.tum.de




Institute for Analysis and Scientific Computing, Vienna University of Technology

Wiedner Hauptstraße 8-10,

1040, Vienna

Austria
szmolyan@tuwien.ac.at


Primary 70K70; Secondary 37G15 











We investigate a singularly perturbed, non-convex variational problem arising in materials science with a
combination of geometrical and numerical methods.
Our starting point is a work by Stefan Müller, where it is proven that the solutions of the variational 
problem are periodic and exhibit a complicated multi-scale structure. In order to get more insight into the
rich solution structure, we transform the corresponding Euler-Lagrange equation 
into a Hamiltonian system of first order ODEs and then use geometric singular perturbation theory to study 
its periodic solutions. Based on the geometric analysis we construct an initial periodic orbit to start numerical 
continuation of periodic orbits with respect to the key parameters. This allows us to observe the influence of the 
parameters on the behavior of the orbits and to study their interplay in the minimization process. Our results 
confirm previous analytical results such as the asymptotics of the period of minimizers predicted by Müller. 
Furthermore, we find several new structures in the entire space of admissible periodic orbits.



Geometry and numerical continuation of multiscale orbits in a nonconvex variational problem
    Peter Szmolyan
    January 27, 2017
===========================================================================================





§ INTRODUCTION



The minimization problem we consider is to find

    min_u ∈ U{ℐ^ε(u):=∫_0^1( ε^2 u_XX^2 + 
    W(u_X) + u^2 ) dX },

where U is a space containing all sufficiently regular functions u:[0,1] 
→ℝ of the spatial variable X ∈[0,1],
0 < ε≪ 1 is a small parameter, 
u_X=∂ u/∂ X, u_XX=∂^2 u/∂ X^2, and the 
function W is a symmetric, double well potential; in particular, here W is chosen as

    W(u_X)=1/4(u_X^2-1)^2.


This model arises in the context of coherent solid-solid phase transformation to describe 
the occurrence of simple laminate microstructures in one-space dimension. Simple laminates 
are defined as particular structures where two phases of the same material (e.g., 
austenite/martensite) simultaneously appear in an alternating pattern <cit.>. This 
situation is shown schematically in Figure <ref>(a). These and related structures 
have been intensively studied both in the context of geometrically linear 
elasticity <cit.> and in the one of fully nonlinear 
elasticity <cit.>. A comparison 
between these two approaches is given by Bhattacharya <cit.>. We focus here on the 
one-dimensional case starting from the work of Müller <cit.>, but a 2D approach has also 
been proposed <cit.>. 
An alternative choice of the functional W which sensibly simplifies energy calculations
for equilibria has been recently adopted by Yip <cit.>. The same functional with more general boundary
conditions has been treated by Vainchtein et al <cit.>.
In all these cases, very significant theoretical and experimental advances have been reached.
Nevertheless, many interesting features concerning the asymptotics and dynamics of these problems can still be explored.

We start from the one-dimensional model (<ref>)-(<ref>) analyzed 
by Müller and introduce a different approach based on geometric singular perturbation 
theory <cit.> which allows us to better understand the critical points of the 
functional ℐ^ε and to obtain an alternative method eventually able 
to handle more general functionals.

       

In <cit.> minimizers are proven to exhibit a periodic multi-scale structure (Figure <ref>): 
a fast scale of order (ε) describes the “jumps” between the two values of the 
derivative u_X, and a slow scale of order (ε^1/3) represents the distance between 
two points with equal value of u_X. From a physical viewpoint, the two values of the derivative 
u_X=± 1 model the two different phases of the material. The jumps describe the transition 
between the phases and the regions with almost constant values of u_X correspond to parts of 
the material occupied by the same phase.
One of the key results in <cit.> consists in an
asymptotic formula for the period of minimizing solutions, when the solution space U is chosen
as the set of all u ∈ H^2(0,1) subject to Dirichlet boundary conditions. 
For ε→ 0, the period P^ε behaves as

    P^ε=2(6 A_0 ε)^1/3 + 𝒪(ε^2/3),

where A_0=2 ∫_-1^1  W^1/2(w)  dw. 

The approach based on fast-slow analysis of the Euler-Lagrange equation applied here
allows us to identify geometrically certain classes
of periodic orbits. These orbits are used as starting solutions for numerical continuation using the software 
package  <cit.>. This powerful tool has been adopted for example by
Grinfeld and Lord <cit.> in their numerical analysis of small amplitude periodic solutions of (<ref>).
We provide here a detailed study of periodic solutions based upon one-parameter continuation in the parameters 
ε and μ. It turns out that several fold bifurcations of periodic orbits  
structure the parameter space. A numerical comparison with the law (<ref>) will be presented, 
by means of a minimization process of the functional ℐ^ε along certain families 
of periodic orbits. Our work also leads to new insights into the dependence of the period 
on the parameters ε and μ for non-minimizing sequences of periodic orbits.

The paper is structured as follows: Section <ref> introduces the approach based on geometric singular 
perturbation theory using the intrinsic multi-scale structure of the problem. We describe the 
transformation of the Euler-Lagrange equation associated to the functional ℐ^ε 
into a multiscale ODE system, along with the decomposition of periodic orbits into slow and fast pieces 
using the Hamiltonian function. 
We identify a family of large amplitude singular periodic orbits and prove their persistence for ϵ small.
A crucial point is the construction of an initial periodic 
orbit for ε≠ 0 in order to start numerical continuation: the strategy we use is 
illustrated in Section <ref>, where the continuation of the orbits with respect to the main 
parameters is also performed. This section includes also the comparison between the analytical expression 
of the period given by Müller and our numerical results as well as the general parameter study of 
periodic solutions. Section <ref> is devoted to conclusions and an outline for future 
work.



§ THE EULER-LAGRANGE EQUATION AS A FAST-SLOW SYSTEM



In this section, the critical points (not only the minimizers) of the functional 
ℐ^ε are analyzed. A necessary condition they 
have to satisfy is the Euler-Lagrange equation <cit.>. The Euler-Lagrange
equation associated to ℐ^ε is the singularly perturbed, 
fourth order equation

    ε^2u_XXXX-1/2σ(u_X)_X+u=0,

where σ(u_X)=W'(u_X)=u_X^3-u_X. Equation (<ref>) can be rewritten 
via 


    w    := u_X, 
    
      v    := -ε^2 w_XX+1/2σ(w), 
    
      z    := ε w_X,


as an equivalent system of first order ODEs

    u̇   =w,
    v̇   =u,
    εẇ   = z,
    εż   = 1/2 (w^3-w)-v,

where d/dX= .
Equations (<ref>) exhibit the structure of a (2,2)-fast-slow system, 
with u, v as slow variables and w, z as fast variables. We recall that a 
system is called (m,n)-fast-slow <cit.> when it has the form

    εẋ   = f(x,y,ε),
    ẏ   = g(x,y,ε),

where x ∈ℝ^m are the fast variables and y ∈ℝ^n 
are the slow variables. The re-formulation of system (<ref>)
on the fast scale is obtained by using the change of variable ξ=X/ε, i.e.

    d x/dξ = x'    = f(x,y,ε),
    d y/dξ = y'    = ε g(x,y,ε).

On the fast scale, system (<ref>) has the form

    u'    =ε w,
    
    v'    =ε u,
    
    w'    = z,
    
    z'    = 1/2 (w^3-w)-v,

which for ε>0 is equivalent to (<ref>).

The system possesses the unique equilibrium 

    p_0= ( 0,0,0,0 ),

which is a center, since the eigenvalues of the Jacobian are all purely imaginary. 
An important property of system (<ref>) is stated in the following result:



Equations (<ref>) and (<ref>) are singularly perturbed Hamiltonian 
systems
=0.7pt1.4

    [        u̇         = - ∂ H/∂ v;        v̇         =   ∂ H/∂ u;        εẇ         =  -∂ H/∂ z;        εż         =  ∂ H/∂ w, ]    ξ=X/ε⇔    [        u'         = -ε∂ H/∂ v;        v'         =  ε∂ H/∂ u;        w'         =  -∂ H/∂ z;        z'         =  ∂ H/∂ w, ]


i.e., they are Hamiltonian systems with respect to the symplectic form 
dz d w + 1/εdv du and with Hamiltonian function 

    H(u,v,w,z)=1/8(4 u^2 - 8 v w - 2 w^2 + w^4 - 4 z^2).




The result follows by differentiating (<ref>) with respect to 
the four variables. For more background on fast-slow Hamiltonian systems 
of this form see <cit.>.


Since the Hamiltonian (<ref>) is a first integral of the system, the dynamics 
take place on level sets, defined by fixing H(u,v,w,z) to a constant value 
μ∈ℝ. This allows us to reduce the dimension of the system by one, which 
we use both in analytical and numerical considerations.

One main advantage in the use of geometric singular perturbation theory is that the 
original problem can be split into two subsystems by analyzing the singular limit 
ε→ 0 on the slow scale (<ref>) and on the fast 
scale (<ref>). The subsystems are usually easier to handle. Under suitable
conditions the combination of both subsystems allows us to obtain information for 
the full system when 0 < ε≪ 1. In particular, if one can construct 
a singular periodic orbit by combining pieces of slow and fast orbits, then the existence of 
a periodic orbit 𝒪(ε)-close to the singular one for small 
ε≠ 0 can frequently be proven under suitable technical 
conditions by tools from geometric singular perturbation 
theory <cit.>.

The slow singular parts of an orbit are derived from the reduced problem (or 
slow subsystem), obtained by letting ε→ 0 in (<ref>)

    0    = f(x,y,0),
    ẏ   = g(x,y,0),

which describes the slow dynamics on the critical manifold 

    _0:={(x,y)∈ℝ^m×ℝ^n:f(x,y,0)=0.}

Considering ε→ 0 on 
the fast scale (<ref>) yields the layer problem (or 
fast subsystem)

    x'    = f(x,y,0),
    
    y'    = 0,

where the fast dynamics is studied on “layers” with constant values of the 
slow variables. Note that _0 can also be viewed as consisting of equilibrium
points for the layer problem. _0 is called normally hyperbolic if 
the eigenvalues of the matrix Df_x(p,0)∈ℝ^m× m 
do not have zero real parts for p∈_0. For normally hyperbolic invariant
manifolds, Fenichel's Theorem applies and yields the existence of a slow
manifold _ε. The slow manifold lies at a distance (ε)
from _0 and the dynamics on _ε is well-approximated by the reduced
problem; for the detailed technical statements of Fenichel's Theorem we refer 
to <cit.>.

In our Hamiltonian fast-slow context, we focus on the analysis of families of 
periodic orbits for system (<ref>) which are parametrized by the 
level set parameter μ. The first goal is to geometrically construct  
periodic orbits in the singular limit ε=0. The reduced 
problem is given by

    u̇   = w,
    v̇   = u,

on the critical manifold (see Fig. <ref>)

    𝒞_0={ (u, v, w, z) ∈ℝ^4 : z=0, v=1/2(w^3-w) }.

The equations of the layer problem are 

    w'    = z
    
    z'    = 1/2(w^3-w)-v̅,

on “layers” where the slow variables are constant (u=u̅, v=v̅).
Note that for Hamiltonian fast-slow systems such as (<ref>), both reduced and layer
problems are Hamiltonian systems with one degree of freedom.



 §.§ The Reduced Problem



Equations (<ref>) describe the reduced problem on _0, if w is considered
as a function of (u, v) on _0.



𝒞_0 is normally hyperbolic except for two fold lines 

    _-    ={(u,(w_-^3-w_-)/2 ,w_-,0)∈ℝ^4},
    _+    ={(u,(w_+^3-w_+)/2 ,w_+,0)∈ℝ^4},

where w_± are defined by σ'(w_±)=0, i.e., w_±=±1/√(3). 
For p∈_±, the matrix D_x f(p,0) has a double zero 
eigenvalue.

      
The lines _± naturally divide 𝒞_0 into three parts

    𝒞_0,l=𝒞_0 ∩{ w < w_- }, 𝒞_0,m=𝒞_0 ∩{w_- ≤ w ≤ w_+}, 𝒞_0,r=𝒞_0 ∩{ w > w_+ },

as shown in Figure <ref>. The submanifolds involved in our analysis are 
only 𝒞_0,l and 𝒞_0,r, which are normally hyperbolic. More
precisely, 𝒞_0,l and 𝒞_0,r are of saddle-type, since 
the matrix D_x f(p,0) along them always has two real eigenvalues of opposite
sign. We remark that saddle-type critical manifolds have played an important role in the
history of fast-slow systems in the context of the travelling wave problem for the 
FitzHugh-Nagumo equation, see for 
example <cit.>.

       
      

On _0-_±, the flow of the reduced system is, up to a time rescaling, given by

    u̇   = (3w^2-1)w,
    ẇ   = 2u.




We differentiate v=1/2(w^3-w) with respect to X, re-write the equation in 
(u,w)-variables and apply the time rescaling corresponding to the multiplication
of the vector field by the factor (3w^2-1) (cf. <cit.>). On 
𝒞_0,m this procedure changes the direction of the flow, but it does not affect
the parts of the critical manifold involved in our analysis.


The Hamiltonian function allows us to restrict our attention to two subsets of 
𝒞_0,l^μ and 𝒞_0,r^μ by fixing the value of μ. 
Analyzing the slow flow on these two normally hyperbolic branches, we see that
u decreases along _0,l and increases along _0,r as shown in Figure <ref>. 

      
      



 §.§ The Layer Problem



The layer problem is obtained by setting ε=0 in (<ref>). We 
obtain a two-dimensional Hamiltonian vector field on “layers” where the slow variables 
are constant (u=u̅,v=v̅)

    w'    = z,
    
    z'    = 1/2(w^3-w)-v̅.

The two branches 𝒞_0,l^μ and 𝒞_0,r^μ are 
hyperbolic saddle equilibria for the system (<ref>) for every value 
of u̅, v̅. To construct a singular limit periodic orbit we are 
particularly interested in connecting orbits between equilibria of the 
layer problem.

      
      


The layer problem (<ref>) has a double heteroclinic connection if and only 
if v̅=0. These are the only possible heteroclinic connections of the layer 
problem (<ref>).



System (<ref>) is Hamiltonian, with v̅ as a parameter and Hamiltonian
function
 
    H_f(w,z)=-z^2/2 + w^4/8 - w^2/4 - v̅ w.

The lemma follows easily by discussing the level curves of the Hamiltonian;
for the convenience of the reader we outline the argument.

Indexing the level set value of H_f as θ, the solutions 
of (<ref>) are level curves {H_f(w,z)=θ}. The equilibria 
of (<ref>) are {z=0,w=w_l,w_m,w_r}; here w_l,w_m,w_r are the 
three solutions of 

    2v̅-w^3+w=0

which depend upon v̅. We only have to
consider the case where there are at least two real equilibria w_l and w_r 
which occurs for v̅∈[-1/(3√(3)),1/(3√(3))]. Let 

    H_f(w_l,0)=:θ_l,     H_f(w_r,0)=:θ_r

and note that since (<ref>) is cubic we can calculate θ_l,r
explicitly. To get a heteroclinic connection we must have θ_l=θ_r
and by an explicit calculation this yields the condition v̅=0. Hence, 
heteroclinic connections of (<ref>) can occur only if v̅=0. 
For v̅=0 one easily finds that the relevant equilibria are located at
w_l=-1 and w_r=1 so that θ_l=-1/8=θ_r. The double heteroclinic
connection is then explicitly given by the curves {z=±1/2(1-w^2)} 
(see also Figure <ref>).

  
The next step is to check where the relevant equilibria of the layer problem 
are located on the critical manifold _0^μ for a fixed value of the 
parameter μ since we have a level set constraint for the full system. Using 
Lemma <ref> one must require w=±1,v=0 while z=0 is the critical
manifold constraint, hence

    H(u,0,±1,0)=1/2 u^2-1/8!=μ.

Therefore, the transition points where fast jumps from _0,l to _0,r
and from _0,r to _0,l are possible are given by

    _0^μ∩{v=0,w=± 1}={u=±√(2μ+1/4),v=0,w=±1,z=0}.

Observe that fast orbits corresponding to positive values of u connect 
𝒞_0,r^μ to 𝒞_0,l^μ, while the symmetric orbits 
with respect to the u=0 plane connect 
𝒞_0,l^μ to 𝒞_0,r^μ.

Recall that w = ± 1 represent the two phases of the material.
Hence, the heteroclinic orbits of the layer problem can be interpreted as
instantaneous transitions between these phases.









 §.§ Singular Fast-Slow Periodic Orbits



The next step is to define singular periodic orbits by combining pieces of orbits of
the reduced and layer problem. Figure <ref> illustrates the situation. The entire singular orbit 
γ_0^μ is obtained connecting two pieces of orbits of the reduced problem with
heteroclinic orbits of the fast subsystem for a 
fixed value of μ, see Figure <ref>(a). The configuration of the two-dimensional
critical manifold and the singular periodic orbit is indicated in Figure <ref>(b).

      

Here we are only interested in singular periodic orbits which have nontrivial slow
and fast segments. Therefore, we do need transitions points from the fast 
subsystem to the slow subsystem. This requirement implies, by using the 
result (<ref>), the lower bound μ>-1/8. A second requirement
we impose is that the slow subsystem orbits lie inside the normally hyperbolic 
parts C_0,l^μ and C_0,r^μ. The u-coordinate of the slow segment
closest to the lines _± is located at u=0. Hence, we calculate 
the value of the Hamiltonian under the condition that the slow trajectory is tangent to
_± which yields

    H(0,1/2(w_±^3-w_±),w_±,0)=1/24.

Combining these considerations with the results from 
Sections <ref>-<ref> gives the following result on the existence 
of singular periodic orbits (ε=0):



For ε=0, the fast-slow system (<ref>),(<ref>) 
has a family of periodic orbits {γ_0^μ}_μ consisting 
of precisely two fast and two slow subsystem trajectories with slow parts
lying entirely in _0,l and _0,r if and only if 

    μ∈ I_μ,  I_μ := 
    (-1/8,1/24).





The persistence of these periodic orbits for 0 < ε≪ 1 on each individual surface level
of the Hamiltonian can be proven by using an argument based on the theorem introduced
by Soto-Treviño in <cit.>.


For every μ∈ I_μ and for ε > 0 sufficiently small,
there exists a locally unique periodic orbit of the fast-slow system (<ref>),(<ref>) that is
𝒪(ε) close to the corresponding singular .



The Hamiltonian structure of the system suggests to study the individual levels as parametrized
families by directly applying the Hamiltonian function as a first integral
to reduce the dimension of the system. At first sight, a convenient choice is to express v 
as a function of the variables (u, w, z) and μ

    v=4u^2-8μ-2w^2+w^4-4z^2/8w.

Consequently, equations (<ref>) transform into a (2,1)-fast-slow system

    u̇   = w,
    εẇ   = z,
    εż   = 1/2 (w^3-w)-4u^2-8μ-2w^2+w^4-4z^2/8w.

Theorem 1 in <cit.> for a C^r (r ≥ 1) (2,1)-fast-slow system guarantees the persistence of periodic orbits
consisting of two slow pieces connected by heteroclinic orbits for 0 < ε≪ 1 when the following conditions hold:

 
  * The critical manifolds are one-dimensional and normally hyperbolic (given by Lemma <ref>).
 
  * The intersection between W^u(_0,l) (resp. W^u(_0,r)) and W^s(_0,r) (resp. W^s(_0,l)) is transversal (confirmed by Lemma <ref>, see Fig. <ref>).
 
  * The full system possesses a singular periodic orbit and the slow flow on the critical manifolds is transverse to touch-down and take-off sets,
	which reduce to 0-dimensional objects in this case as we explicitly obtained in (<ref>).


      
      
However, system (<ref>) appears to be nonsmooth at w=0 and the fast orbits  necessarily cross w=0. To overcome this (apparent) difficulty we use other charts 
for the manifold H(u,v,w,z) = μ for parts of the singular orbit close to w=0. Instead of (<ref>) we now express u as a function of the other variables, i.e.

    u = ±1/2√(8 v w + 2 w^2 - w^4 + 4 z^2 + 8 μ).

This leads to the following description of the dynamics:

 
  * System (<ref>) describes the dynamics on the slow pieces away from w=0;
 
  * The heteroclinic connection corresponding to u>0 is expressed by
      
    v'    = + ε/2√(8 v w + 2 w^2 - w^4 + 4 z^2 + 8 μ),
    
    	w'    = z,
    
    	z'    = 1/2 (w^3-w)-v.

 
  * The heteroclinic connection corresponding to u<0 is expressed by
      
    v'    = - ε/2√(8 v w + 2 w^2 - w^4 + 4 z^2 + 8 μ),
    
    	w'    = z,
    
    	z'    = 1/2 (w^3-w)-v.


If we consider system (<ref>) as a smooth dynamical system on the manifold defined by H=μ, the proof given in <cit.> (based on proving the transversal intersection of two manifolds obtained by flowing suitably chosen initial conditions forward and backward in time) goes through without being affected by the fact that we have to work with several coordinate systems, as described above.






























System (<ref>) has two parameters μ,ε, which naturally
leads to the question how periodic orbits deform and bifurcate when the two 
parameters are varied. Furthermore, the fast-slow structure with orbits consisting
of two fast jumps and two slow segments as shown in Figure <ref> and the 
three-dimensional form (<ref>) provide analogies to the travelling wave
frame system obtained from the partial differential equation version of the 
FitzHugh-Nagumo <cit.> (FHN) equation. The three-dimensional fast-slow 
FHN system has been studied in great detail using various fast-slow systems 
techniques (see e.g. <cit.>).
One particular approach to investigate the FHN parameter space efficiently is to 
employ numerical continuation methods <cit.>. In fact,
numerical approaches to FHN have frequently provided interesting conjectures and 
thereby paved the way for further analytical studies. Adopting this approach, we are 
going to investigate the problem (<ref>) considered here using numerical 
continuation to gain better insight into the structure of periodic orbits.



§ NUMERICAL CONTINUATION



This section is devoted to the numerical investigation of the the critical points 
of the functional ℐ^ε via the Euler-Lagrange 
formulation (<ref>). A powerful tool for such computations 
is  <cit.>.  is able to numerically
track periodic orbits depending upon parameters using a combination of a boundary
value problem (BVP) solver with a numerical continuation algorithm. Using
such a framework for fast-slow systems often yields a wide variety
of interesting numerical and visualization results; for a few recent examples we refer 
to <cit.>.

The first task one has to deal with is the construction of a starting orbit for fixed 
ε≠ 0. For (<ref>) this is actually a less trivial task than
for the FHN equation as we are going to explain in Section <ref>. 
In Section <ref>, we are also going to construct a starting periodic 
orbit based upon the geometric insights of Section <ref>.

Once the starting periodic orbit is constructed, we use  to perform 
numerical continuation in both parameters μ and ε. This yields 
bifurcation diagrams and the solutions corresponding to some interesting points 
on the bifurcation branches. Then, the connection between the parameters in the 
minimization process is investigated, in order to numerically determine the 
correspondence that leads to the functional minimum. Finally, a comparison with 
the period law (<ref>) predicted by Müller <cit.> is performed. 



 §.§ Construction of the starting orbit



As indicated already, the construction of a starting periodic orbit is not trivial:



  * The singular orbit itself, obtained by matching slow and fast subsystem orbits 
for a fixed value of μ and for ε=0, cannot be used owing to re-scaling 
problems (the fast pieces would all correspond to x=0).

  * The computation of a full periodic orbit using a direct initial value solver 
approach for  is hard to perform since the slow manifolds are 
of saddle type and an orbit computed numerically would diverge from them exponentially 
fast <cit.>.

  * Matching slow segments obtained with a saddle-type algorithm <cit.> 
and fast parts computed with an initial value solver may cause problems at the points 
where the four pieces should match.

  * In contrast to the FitzHugh-Nagumo case <cit.>, the periodic 
orbits we are looking for cannot be detected as Hopf bifurcations from the zero equilibrium. 
In that case, we could use  to locate such bifurcations and then find a periodic 
orbit for 0 < ε≪ 1 fixed by branch-switching at the Hopf bifurcation point. In 
our case, however, the origin p_0 is a center equilibrium, and an infinite number of periodic 
orbits exist around it in the formulation (<ref>). 

  * Starting continuation close to the equilibrium p_0 is difficult due to its degenerate nature (w=0).


 


Our strategy is to use the geometric insight from Section <ref> in combination with
a slow manifolds of saddle-type (SMST) algorithm and a homotopy approach. We construct an 
approximate starting periodic orbit using a value of μ which leads to a short 
“time” spent on the slow parts of the orbits, so that the saddle-type branches 
do not lead to numerical complications. Then we use an SMST algorithm to find a suitable 
pair of starting points lying extremely close to the left and right parts of the slow manifolds 
^μ_ε,l and ^μ_ε,r. In the last step we employ 
numerical continuation to study the values of μ we are actually interested in; this 
is the homotopy step.

The value μ=-1/8 is peculiar, since in this case the singular slow segments for (<ref>) reduce 
to two points

    𝒞_0,l^-1/8 = { (0,-1,0) },    and    𝒞_0,r^-1/8 = { (0,1,0) },

i.e., touch-down and take-off sets for the fast dynamics coincide in this case.
The range (<ref>) we are considering does not include μ=-1/8; 
however, this property makes it an excellent candidate for the first step of our strategy. 
Indeed, we know already from the geometric analysis in Section <ref> that the time
spent near slow manifolds is expected to be very short in this case.

Although it is still not possible to compute the full orbit using forward/backward integration, 
we can compute two halves, provided we choose the correct initial condition. We aim to find 
a point on the slow manifolds ^μ_ε,l and ^μ_ε,r as an
initial value. The SMST algorithm <cit.> helps to solve this problem. The 
procedure is based on a BVP method to compute slow manifolds of saddle-type in fast-slow 
systems. Fixing ε and μ, we select manifolds B_l and B_r, which are 
transverse to the stable and unstable eigenspaces of 𝒞_0,l^μ and 
𝒞_0,r^μ, respectively (Figure <ref>). The plane B_l and the line
B_r provide the boundary conditions for the SMST algorithm.

Implementing the algorithm for μ=-1/8 and ε=0.001 for (<ref>)
shows that there are actually two points (0,w_L,0) and (0,w_R,0) which 
are contained in the slow manifold even for ε≠ 0 as well as in the critical 
manifold _0. From the geometric analysis in Section <ref> we know that at 
μ=-1/8 the take-off and touch-down points coincide and a singular 
double-heteroclinic loop exists for v=0. This motivates the choice of (u,v,w,z)=(0,0,w_L,0) 
and (u,v,w,z)=(0,0,w_R,0) in the following algorithm: a numerical integration 
of the full four-dimensional problem (<ref>) forward and backward in x is performed, 
imposing the Hamiltonian constraint using a projective algorithm. The computation is stopped 
once the hyperplane {w=0} is reached. The full periodic orbit is then constructed by matching 
two symmetric pieces together.

In principle, there are different ways how one may arrive at a useful construction of a
highly accurate starting periodic orbit. In our context, the geometric analysis guided the
way to identify the simplest numerical procedure, which is an approach that is likely to be
successful for many other non-trivial fast-slow numerical continuation problems.



 §.§ Continuation in μ



A detailed analysis of the critical points' dependence on the Hamiltonian is performed. 
The value of μ can be arbitrarily chosen only in the interval I_μ, while 
ε is fixed to 0.001. Continuation is performed on system (<ref>) 
using the initial orbit obtained numerically in Section <ref>. Starting at
μ=-1/8,  is able to compute the variation of the orbits up to 
μ=1/24. The bifurcation diagram of the period P with respect to the 
parameter μ is shown in Figure <ref>(a).

      
      
The first/upper branch of the continuation displays fast-slow orbits corresponding perturbations 
of the singular ones {γ_0^μ}_μ∈ I_μ for fixed 
ε≠ 0. As predicted by the geometric analysis we observe that decreasing 
μ reduces the length of the slow parts, so that the orbits almost correspond to the double
heteroclinic one analytically constructed at μ = -1/8; see Figure <ref>. Near 
μ=-1/8 the bifurcation branch has a fold in (μ,P)-space leading to 
the second/lower bifurcation branch. The difference between the orbits on the two 
branches for a fixed value of μ is shown in Figure <ref>(b). Along the 
second branch, periodic solutions around the center equilibrium appear, which collapse 
into it with increasing μ (Figure <ref>(b)).

Furthermore, numerical continuation robustly indicates that the upper branch has another 
fold when continued from μ=0 to higher values of μ as shown in 
Figure <ref>(a). The orbits obtained by fixing a value of μ on the upper branch
and its continuation after the fold differ only because of the appearance of two new fast 
parts near the plane {u=0} as shown in Figure <ref>(a). We conjecture that 
these parts arise due to the loss of normal hyperbolicity at _±; see also 
Section <ref>.

      
      



 §.§ Continuation in ε



We perform numerical continuation in ε by fixing three values of 
μ in order to capture the behavior of the solutions for the range I_μ 
from Proposition <ref>. We consider μ_l ≈ -1/8 with
μ_l>-1/8, μ_c ≈ 0, and μ_r ≈1/24 with 
μ_r<1/24; or more precisely μ_l = -0.12489619925, 
μ_c = 1.5378905702· 10^-5, and μ_r = 0.04100005066. 
For each of these values, we find two bifurcation branches connected via a fold
in (μ,P)-space; see Figure <ref>.

The bifurcation diagrams and the associated solutions shown in Figure <ref>
nicely illustrate the dependence of the period on the singular perturbation parameter
ε. When ε→ 0 there are two very distinct limits for 
the period P=P(ε) (Figure <ref>(a)-(b)) depending whether we are 
on the upper and lower parts of the main branch of solutions. In the case with 
μ_r≈1/24 when orbits come close to non-hyperbolic singularities on
_0, we actually seem to observe that P(0) 
seems to be independent on whether we consider the upper or lower part of the branch
(see Figure <ref>(c)). Furthermore, functional forms of P(ε) are 
clearly different for small ε so the natural conjecture is that there 
is no universal periodic scaling law if we drop the functional minimization constraint.

The deformation under variation of ε of the periodic orbits in (w,z,u)-space 
is also interesting. For μ=μ_l (Figure <ref>(a)), we observe that the upper 
branch corresponds to the solutions that we expect analytically from 
Proposition <ref> consisting of two fast and two slow segments when 
approaching ε=0. A similar scenario occurs also 
for the other values of μ (Figure <ref>(b) and Figure <ref>(c)). 
When the ε value is too large, or when we are on 
a different part of the branch of solutions, the orbits closed to the equilibrium of 
the full system or additional pieces resembling new fast contributions appear.

      
      


 §.§ Period scaling



So far, no boundary conditions have been imposed; moreover, all the computed solutions 
are not necessarily minimizers of the functional, but only critical points. Our 
conjecture is that the interaction between the two main parameters of the system 
μ and ε should allow us to obtain the true minimizers via a 
double-limit. In other words, for every value of ε there is a 
corresponding orbit which minimizes the functional ℐ^ε, 
and since along this orbit the Hamiltonian has to constantly assume a certain 
value μ̅, the minimization process should imply a direct connection 
between the parameters. Consequently, it is interesting to investigate this 
ansatz from the numerical viewpoint.

      
      
A first possibility is to establish a connection between the two parameters 
ε and μ via a direct continuation in both parameters, starting 
from certain special points, such as the fold points detected in 
Sections <ref>-<ref>. However, it turns out that this
process does not lead to the correct scaling law for minimizers of 
_ε as shown in Figure <ref>.

Another option is instead to check if among the critical points of the Euler-Lagrange 
equation (<ref>) we have numerically obtained there are also the minimizers of 
the functional ℐ^ε respecting the power law (<ref>). 
In <cit.>, boundary conditions on the interval [ 0,1 ] are also included 
in the variational formulation, and from the results obtained from the continuation in 
ε, one may expect that high values of μ would not be able to fit them, 
since the period is always too high. Lower values of μ, instead, seem to have sufficiently 
small period. Hence, one could fix one of those (for example, μ_l) and look at what happens 
as ε→ 0. The hope is that the (ε^1/3) leading-order scaling 
for the period naturally emerges. Unfortunately, this does not happen, as we can see in 
Figure <ref>; the lower branch seems to give a linear dependence on ε, 
while the upper branch gives a quadratic one.


      

The reason why from this naive approach the (ε^1/3) leading-order scaling does 
not emerge lies in the lack of connection with the minimization process. However, 
Figure <ref> demonstrates that there are several nontrivial scalings of natural 
families of periodic orbits as ε→ 0.


So far, we have just assumed that the Hamiltonian value of the minimizers should be “low”, 
but indeed there is a strict connection between the values of ε one is considering 
and the value of μ of the minimizers. In other words, there is not a unique value of 
μ given by the minimizers for every ε small but minimizers move over different
Hamiltonian energy levels as ε→ 0. Starting from this consideration, 
another option, which turns out to be the correct one to recover the scaling (<ref>), 
is to use the periodic orbits from numerical continuation to compute the numerical value of 
the functional ℐ^ε as a function of the period P fixing different 
values of ε in a suitable range, such as:

    I_ε= [ 10^-7, 10^-1].

			
Then, we obtain different parabola-shaped diagrams where we can extract the value of the period 
minimizing the functional (Figure <ref>). When plotting these values related to the 
value of ε for which they have been computed, one obtains the results shown in 
Figure <ref>. The values numerically extracted from our solutions match the analytical 
results on the period proven by Müller (<ref>) when the value of ε is 
sufficiently small. As ε increases, the period law is less accurate, as one 
would expect. 

      

      



§ CONCLUSION & OUTLOOK



In summary, we have shown that geometric singular perturbation theory and numerical continuation 
methods can be very effective tools to understand nonconvex multiscale variational
problems via the Euler-Lagrange formulation. We have proven the existence of a class of singular 
periodic orbits based upon a fast-slow decomposition approach and we have shown that these orbits persist
for ε small.

The geometric insight was used
to determine a starting solution for numerical continuation in the context of a reduced 
three-dimensional fast-slow system. Then we studied the dependence of periodic orbits on the 
singular perturbation parameter as well as the Hamiltonian energy level set parameter arising
in the reduction from a four- to a three-dimensional system. The parameter space is structured
by several fold points. Furthermore, we were able to study the shape of non-minimizing periodic
orbits for very broad classes of parameters. Finally, we showed that several natural scaling 
laws for non-minimizing sequences of periodic solutions exist and also confirmed numerically
the leading-order scaling predicted by Müller for minimizing sequences.

Based upon this work, there are several open problems as well as generalizations one might
consider. In particular, it would be desirable to extend the persistence result 
to the general class of singularly-perturbed Hamiltonian fast-slow systems (<ref>); this
is the subject of ongoing work.

Another important observation of our numerical study are the intricate orbits that seem to
arise when parts of the slow segments start to interact with the singularities _± where
the critical manifold is not normally hyperbolic. The natural conjecture is that the 
additional small fast loops that we observe numerically could correspond to homoclinic “excursions”
in the fast subsystem anchored at points close to _±. The blow-up method <cit.> is 
likely to provide an excellent tool to resolve the non-normally hyperbolic singularities; see 
e.g. <cit.> where the existence of complicated fast-slow 
periodic orbits involving loss of normal hyperbolicity is proven.

The construction of an initial orbit has been one of the hardest problems to tackle. It was solved 
using analytical and numerical tools, after discarding several other plausible approaches. 
The SMST algorithm <cit.> has been a helpful tool to determine good starting points on the slow 
manifolds and then use an initial value solver to obtain segments of a complete whole orbit. Although
our approach works well in practical computations, there are interesting deep numerical analysis
questions still to be answered regarding the interplay between certain classes of fast-slow “initial
guess” starting orbits and the success or failure of Newton-type methods for the associated BVPs. In
particular, can one prove certain geometric conditions or restrictions on ε to guarantee
the convergence for the first solution?

Another highly relevant direction would be to extend our approach to more general classes of 
functionals. There are many different singularly-perturbed variational problems, arising 
e.g. in materials science, to which one may apply the techniques presented here. In this context,
it is important to emphasize that we expect that particularly other non-convex functionals could 
be excellent candidates for future work.

From the viewpoint of applications, it would be interesting to study the practical relevance of 
non-minimizing sequences of periodic solutions. Although we expect the long-term behavior to be 
governed by minimizers, it is evident that non-minimizing periodic orbits can have a high impact 
on time-dependent dynamics, e.g., either via transient behavior, via noise-induced phase 
transitions, or as dynamical boundaries between different regimes.

Acknowledgements: AI and PS would like to thank the Fonds zur Förderung der 
wissenschaftlichen Forschung (FWF) for support via a doctoral school (project W1245). 
CK would like to thank the VolkswagenStiftung for
support via a Lichtenberg professorship. CK and PS also acknowledge partial support 
of the European Commission (EC/REA) via a Marie-Curie International Reintegration Grant (MC-IRG). 

plain


