Geometry and numerical continuation of multiscale orbits in a nonconvex variational problem

TL;DR

This paper combines geometric and numerical methods to analyze multiscale periodic solutions in a complex non-convex variational problem relevant to materials science.

Contribution

It transforms the Euler-Lagrange equation into a Hamiltonian system and uses geometric singular perturbation theory for numerical continuation of solutions.

Findings

Confirmed asymptotic behavior of minimizer periods

Constructed initial periodic orbits for numerical continuation

Discovered new structures in the space of periodic orbits

Abstract

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\"uller, 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\"uller. Furthermore, we find several new structures in the entire space of admissible periodic orbits.

Equations

$$\underset{u\in U}{\min }\left\{{\mathcal{I}}^{\epsilon }(u):={\int }_0^1\left({\epsilon }^2{u}_{XX}^2+W(u_X)+u^2\right)\mathrm{d}X\right\},$$

$$W(u_X)=\frac{1}{4}(u_X^2-1{)}^2.$$

$$P^{\epsilon }=2(6A_0\epsilon {)}^{1/3}+\mathcal{O}({\epsilon }^{2/3}),$$

$${\epsilon }^2{u}_{XXXX}-\frac{1}{2}\sigma (u_X{)}_X+u=0,$$

$$w$$

$$v$$

$$z$$

$$\dot{u}$$

$$\dot{v}$$

$$\epsilon \dot{w}$$

$$\epsilon \dot{z}$$

$$\epsilon \dot{x}$$

$$\dot{y}$$

$$\frac{\mathrm{d}x}{\mathrm{d}\xi }=x^{\prime }$$

$$\frac{\mathrm{d}y}{\mathrm{d}\xi }=y^{\prime }$$

$$u^{\prime }$$

$$v^{\prime }$$

$$w^{\prime }$$

$$z^{\prime }$$

$$p_0=\left(0,0,0,0\right),$$

$$\begin{array}[]{rcr}\dot{u}&=&-\frac{\partial H}{\partial v}\\ \dot{v}&=&\frac{\partial H}{\partial u}\\ \varepsilon\dot{w}&=&-\frac{\partial H}{\partial z}\\ \varepsilon\dot{z}&=&\frac{\partial H}{\partial w},\end{array}\qquad\overset{\textstyle{\xi=\frac{X}{\varepsilon}}}{\Leftrightarrow}\qquad\begin{array}[]{rcr}u^{\prime}&=&-\varepsilon\frac{\partial H}{\partial v}\\ v^{\prime}&=&\varepsilon\frac{\partial H}{\partial u}\\ w^{\prime}&=&-\frac{\partial H}{\partial z}\\ z^{\prime}&=&\frac{\partial H}{\partial w},\end{array}$$

$$H(u,v,w,z)=\frac{1}{8}(4u^2-8vw-2w^2+w^4-4z^2).$$

$$0$$

$$\dot{y}$$

$${\mathcal{C}}_0:=\{(x,y)\in {\mathbb{R}}^m\times {\mathbb{R}}^n:f(x,y,0)=0.\}$$

$$x^{\prime }$$

$$y^{\prime }$$

$$\dot{u}$$

$$\dot{v}$$

$${\mathcal{C}}_0=\left\{(u,v,w,z)\in {\mathbb{R}}^4:z=0,v=\frac{1}{2}(w^3-w)\right\}.$$

$$w^{\prime }$$

$$z^{\prime }$$

$${\mathcal{L}}_{-}$$

$${\mathcal{L}}_{+}$$

$${\mathcal{C}}_{0,l}={\mathcal{C}}_0\cap \left\{w<w_{-}\right\},{\mathcal{C}}_{0,m}={\mathcal{C}}_0\cap \left\{w_{-}\le w\le w_{+}\right\},{\mathcal{C}}_{0,r}={\mathcal{C}}_0\cap \left\{w>w_{+}\right\},$$

$$\dot{u}$$

$$\dot{w}$$

$$w^{\prime }$$

Full text

1. Introduction

The minimization problem we consider is to find

[TABLE]

where $U$ is a space containing all sufficiently regular functions $u:\left[0,1\right]\rightarrow\mathbb{R}$ of the spatial variable $X\in\left[0,1\right]$, $0<\varepsilon\ll 1$ is a small parameter, $u_{X}=\frac{\partial u}{\partial X}$, $u_{XX}=\frac{\partial^{2}u}{\partial X^{2}}$, and the function $W$ is a symmetric, double well potential; in particular, here $W$ is chosen as

[TABLE]

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 [37]. This situation is shown schematically in Figure 1(a). These and related structures have been intensively studied both in the context of geometrically linear elasticity [28, 29, 41] and in the one of fully nonlinear elasticity [1, 2, 3, 4, 6, 13, 39, 40, 43, 44]. A comparison between these two approaches is given by Bhattacharya [5]. We focus here on the one-dimensional case starting from the work of Müller [36], but a 2D approach has also been proposed [30, 18, 25]. An alternative choice of the functional $W$ which sensibly simplifies energy calculations for equilibria has been recently adopted by Yip [47]. The same functional with more general boundary conditions has been treated by Vainchtein et al [46]. 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 (1)-(2) analyzed by Müller and introduce a different approach based on geometric singular perturbation theory [26, 34] which allows us to better understand the critical points of the functional $\mathcal{I}^{\varepsilon}$ and to obtain an alternative method eventually able to handle more general functionals.

In [36] minimizers are proven to exhibit a periodic multi-scale structure (Figure 1): a fast scale of order ${\mathcal{O}}(\varepsilon)$ describes the “jumps” between the two values of the derivative $u_{X}$, and a slow scale of order ${\mathcal{O}}(\varepsilon^{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}=\pm 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 [36] consists in an asymptotic formula for the period of minimizing solutions, when the solution space $U$ is chosen as the set of all ${u\in H^{2}(0,1)}$ subject to Dirichlet boundary conditions. For $\varepsilon\rightarrow 0$, the period $P^{\varepsilon}$ behaves as

[TABLE]

where $A_{0}=2\int_{-1}^{1}\!W^{1/2}(w)\,\mathrm{d}w$.

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 AUTO [12]. This powerful tool has been adopted for example by Grinfeld and Lord [19] in their numerical analysis of small amplitude periodic solutions of (1). We provide here a detailed study of periodic solutions based upon one-parameter continuation in the parameters $\varepsilon$ and $\mu$. It turns out that several fold bifurcations of periodic orbits structure the parameter space. A numerical comparison with the law (3) will be presented, by means of a minimization process of the functional $\mathcal{I}^{\varepsilon}$ along certain families of periodic orbits. Our work also leads to new insights into the dependence of the period on the parameters $\varepsilon$ and $\mu$ for non-minimizing sequences of periodic orbits.

The paper is structured as follows: Section 2 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 $\mathcal{I}^{\varepsilon}$ 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 $\epsilon$ small. A crucial point is the construction of an initial periodic orbit for $\varepsilon\neq 0$ in order to start numerical continuation: the strategy we use is illustrated in Section 3, 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 4 is devoted to conclusions and an outline for future work.

2. The Euler-Lagrange equation as a fast-slow system

In this section, the critical points (not only the minimizers) of the functional $\mathcal{I}^{\varepsilon}$ are analyzed. A necessary condition they have to satisfy is the Euler-Lagrange equation [9]. The Euler-Lagrange equation associated to $\mathcal{I}^{\varepsilon}$ is the singularly perturbed, fourth order equation

[TABLE]

where $\sigma(u_{X})=W^{\prime}(u_{X})=u_{X}^{3}-u_{X}$. Equation (4) can be rewritten via

[TABLE]

as an equivalent system of first order ODEs

[TABLE]

where $\frac{\mathrm{d}}{\mathrm{d}X}=\dot{}$ . Equations (5) 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 [15, 26, 34] when it has the form

[TABLE]

where $x\in\mathbb{R}^{m}$ are the fast variables and $y\in\mathbb{R}^{n}$ are the slow variables. The re-formulation of system (6) on the fast scale is obtained by using the change of variable $\xi=\frac{X}{\varepsilon}$, i.e.

[TABLE]

On the fast scale, system (5) has the form

[TABLE]

which for $\varepsilon>0$ is equivalent to (5).

The system possesses the unique equilibrium

[TABLE]

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

Lemma 2.1**.**

Equations (5) and (8) are singularly perturbed Hamiltonian systems

[TABLE]

i.e., they are Hamiltonian systems with respect to the symplectic form $\mathrm{d}z\land\mathrm{d}w+\frac{1}{\varepsilon}\mathrm{d}v\land\mathrm{d}u$ and with Hamiltonian function

[TABLE]

Proof.

The result follows by differentiating (11) with respect to the four variables. For more background on fast-slow Hamiltonian systems of this form see [17]. ∎

Since the Hamiltonian (11) 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 $\mu\in\mathbb{R}$. 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 $\varepsilon\rightarrow 0$ on the slow scale (6) and on the fast scale (7). 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<\varepsilon\ll 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 $\mathcal{O}(\varepsilon)$-close to the singular one for small $\varepsilon\neq 0$ can frequently be proven under suitable technical conditions by tools from geometric singular perturbation theory [15, 26, 34, 42].

The slow singular parts of an orbit are derived from the reduced problem (or slow subsystem), obtained by letting $\varepsilon\rightarrow 0$ in (6)

[TABLE]

which describes the slow dynamics on the critical manifold

[TABLE]

Considering $\varepsilon\rightarrow 0$ on the fast scale (7) yields the layer problem (or fast subsystem)

[TABLE]

where the fast dynamics is studied on “layers” with constant values of the slow variables. Note that ${\mathcal{C}}_{0}$ can also be viewed as consisting of equilibrium points for the layer problem. ${\mathcal{C}}_{0}$ is called normally hyperbolic if the eigenvalues of the matrix $\textnormal{D}f_{x}(p,0)\in\mathbb{R}^{m\times m}$ do not have zero real parts for $p\in{\mathcal{C}}_{0}$. For normally hyperbolic invariant manifolds, Fenichel’s Theorem applies and yields the existence of a slow manifold ${\mathcal{C}}_{\varepsilon}$. The slow manifold lies at a distance ${\mathcal{O}}(\varepsilon)$ from ${\mathcal{C}}_{0}$ and the dynamics on ${\mathcal{C}}_{\varepsilon}$ is well-approximated by the reduced problem; for the detailed technical statements of Fenichel’s Theorem we refer to [15, 26, 34].

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

[TABLE]

on the critical manifold (see Fig. 2)

[TABLE]

The equations of the layer problem are

[TABLE]

on “layers” where the slow variables are constant ($u=\bar{u},\ v=\bar{v}$). Note that for Hamiltonian fast-slow systems such as (10), both reduced and layer problems are Hamiltonian systems with one degree of freedom.

2.1. The Reduced Problem

Equations (15) describe the reduced problem on ${\mathcal{C}}_{0}$, if $w$ is considered as a function of $(u,v)$ on ${\mathcal{C}}_{0}$.

Lemma 2.2**.**

$\mathcal{C}_{0}$* is normally hyperbolic except for two fold lines*

[TABLE]

where $w_{\pm}$ are defined by $\sigma^{\prime}(w_{\pm})=0$, i.e., $w_{\pm}=\pm 1/\sqrt{3}$. For $p\in{\mathcal{L}}_{\pm}$, the matrix $\mathrm{D}_{x}f(p,0)$ has a double zero eigenvalue.

The lines ${\mathcal{L}}_{\pm}$ naturally divide $\mathcal{C}_{0}$ into three parts

[TABLE]

as shown in Figure 2. The submanifolds involved in our analysis are only $\mathcal{C}_{0,l}$ and $\mathcal{C}_{0,r}$, which are normally hyperbolic. More precisely, $\mathcal{C}_{0,l}$ and $\mathcal{C}_{0,r}$ are of saddle-type, since the matrix $\mathrm{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 [27, 33, 21].

Lemma 2.3**.**

On ${\mathcal{C}}_{0}-{\mathcal{L}}_{\pm}$, the flow of the reduced system is, up to a time rescaling, given by

[TABLE]

Proof.

We differentiate $v=\frac{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. [34, Sec.7.7]). On $\mathcal{C}_{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 $\mathcal{C}_{0,l}^{\mu}$ and $\mathcal{C}_{0,r}^{\mu}$ by fixing the value of $\mu$. Analyzing the slow flow on these two normally hyperbolic branches, we see that $u$ decreases along ${\mathcal{C}}_{0,l}$ and increases along ${\mathcal{C}}_{0,r}$ as shown in Figure 3.

2.2. The Layer Problem

The layer problem is obtained by setting $\varepsilon=0$ in (8). We obtain a two-dimensional Hamiltonian vector field on “layers” where the slow variables are constant ($u=\bar{u},v=\bar{v}$)

[TABLE]

The two branches $\mathcal{C}_{0,l}^{\mu}$ and $\mathcal{C}_{0,r}^{\mu}$ are hyperbolic saddle equilibria for the system (20) for every value of $\bar{u},\bar{v}$. To construct a singular limit periodic orbit we are particularly interested in connecting orbits between equilibria of the layer problem.

Lemma 2.4**.**

The layer problem (20) has a double heteroclinic connection if and only if $\bar{v}=0$. These are the only possible heteroclinic connections of the layer problem (20).

Proof.

System (20) is Hamiltonian, with $\bar{v}$ as a parameter and Hamiltonian function

[TABLE]

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 $\theta$, the solutions of (20) are level curves $\{H_{f}(w,z)=\theta\}$. The equilibria of (20) are $\{z=0,w=w_{l},w_{m},w_{r}\}$; here $w_{l},w_{m},w_{r}$ are the three solutions of

[TABLE]

which depend upon $\bar{v}$. We only have to consider the case where there are at least two real equilibria $w_{l}$ and $w_{r}$ which occurs for $\bar{v}\in[-1/(3\sqrt{3}),1/(3\sqrt{3})]$. Let

[TABLE]

and note that since (21) is cubic we can calculate $\theta_{l,r}$ explicitly. To get a heteroclinic connection we must have $\theta_{l}=\theta_{r}$ and by an explicit calculation this yields the condition $\bar{v}=0$. Hence, heteroclinic connections of (20) can occur only if $\bar{v}=0$. For $\bar{v}=0$ one easily finds that the relevant equilibria are located at $w_{l}=-1$ and $w_{r}=1$ so that $\theta_{l}=-1/8=\theta_{r}$. The double heteroclinic connection is then explicitly given by the curves $\{z=\pm\frac{1}{2}(1-w^{2})\}$ (see also Figure 4). ∎

The next step is to check where the relevant equilibria of the layer problem are located on the critical manifold ${\mathcal{C}}_{0}^{\mu}$ for a fixed value of the parameter $\mu$ since we have a level set constraint for the full system. Using Lemma 2.4 one must require $w=\pm 1,v=0$ while $z=0$ is the critical manifold constraint, hence

[TABLE]

Therefore, the transition points where fast jumps from ${\mathcal{C}}_{0,l}$ to ${\mathcal{C}}_{0,r}$ and from ${\mathcal{C}}_{0,r}$ to ${\mathcal{C}}_{0,l}$ are possible are given by

[TABLE]

Observe that fast orbits corresponding to positive values of $u$ connect $\mathcal{C}_{0,r}^{\mu}$ to $\mathcal{C}_{0,l}^{\mu}$, while the symmetric orbits with respect to the $u=0$ plane connect $\mathcal{C}_{0,l}^{\mu}$ to $\mathcal{C}_{0,r}^{\mu}$.

Recall that $w=\pm 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.

2.3. 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 5 illustrates the situation. The entire singular orbit $\gamma_{0}^{\mu}$ is obtained connecting two pieces of orbits of the reduced problem with heteroclinic orbits of the fast subsystem for a fixed value of $\mu$, see Figure 5(a). The configuration of the two-dimensional critical manifold and the singular periodic orbit is indicated in Figure 5(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 (22), the lower bound $\mu>-\frac{1}{8}$. A second requirement we impose is that the slow subsystem orbits lie inside the normally hyperbolic parts $C_{0,l}^{\mu}$ and $C_{0,r}^{\mu}$. The $u$-coordinate of the slow segment closest to the lines ${\mathcal{L}}_{\pm}$ is located at $u=0$. Hence, we calculate the value of the Hamiltonian under the condition that the slow trajectory is tangent to ${\mathcal{L}}_{\pm}$ which yields

[TABLE]

Combining these considerations with the results from Sections 2.1-2.2 gives the following result on the existence of singular periodic orbits ($\varepsilon=0$):

Proposition 2.5**.**

For $\varepsilon=0$, the fast-slow system (5),(8) has a family of periodic orbits $\{\gamma_{0}^{\mu}\}_{\mu}$ consisting of precisely two fast and two slow subsystem trajectories with slow parts lying entirely in ${\mathcal{C}}_{0,l}$ and ${\mathcal{C}}_{0,r}$ if and only if

[TABLE]

The persistence of these periodic orbits for $0<\varepsilon\ll 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 [42].

Theorem 2.6**.**

For every $\mu\in I_{\mu}$ and for $\varepsilon>0$ sufficiently small, there exists a locally unique periodic orbit of the fast-slow system (5),(8) that is $\mathcal{O}(\varepsilon)$ close to the corresponding singular orbit $\gamma_{0}^{\mu}$.

Proof.

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 $\mu$

[TABLE]

Consequently, equations (5) transform into a $(2,1)$-fast-slow system

[TABLE]

Theorem $1$ in [42] for a $C^{r}$ $(r\geq 1)$ $(2,1)$-fast-slow system guarantees the persistence of periodic orbits consisting of two slow pieces connected by heteroclinic orbits for $0<\varepsilon\ll 1$ when the following conditions hold:

• •

The critical manifolds are one-dimensional and normally hyperbolic (given by Lemma 2.2).

• •

The intersection between $W^{u}({\mathcal{C}}_{0,l})$ (resp. $W^{u}({\mathcal{C}}_{0,r})$) and $W^{s}({\mathcal{C}}_{0,r})$ (resp. $W^{s}({\mathcal{C}}_{0,l})$) is transversal (confirmed by Lemma 2.4, see Fig. 6).

• •

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 (22).

However, system (26) 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)=\mu$ for parts of the singular orbit close to $w=0$. Instead of (25) we now express $u$ as a function of the other variables, i.e.

[TABLE]

This leads to the following description of the dynamics:

• •

System (26) describes the dynamics on the slow pieces away from $w=0$;

• •

The heteroclinic connection corresponding to $u>0$ is expressed by

[TABLE]

• •

The heteroclinic connection corresponding to $u<0$ is expressed by

[TABLE]

If we consider system (5) as a smooth dynamical system on the manifold defined by $H=\mu$, the proof given in [42] (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 (26) has two parameters $\mu,\varepsilon$, 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 5 and the three-dimensional form (26) provide analogies to the travelling wave frame system obtained from the partial differential equation version of the FitzHugh-Nagumo [16, 38] (FHN) equation. The three-dimensional fast-slow FHN system has been studied in great detail using various fast-slow systems techniques (see e.g. [7, 21, 27, 33]). One particular approach to investigate the FHN parameter space efficiently is to employ numerical continuation methods [22, 8]. 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 (26) considered here using numerical continuation to gain better insight into the structure of periodic orbits.

3. Numerical Continuation

This section is devoted to the numerical investigation of the the critical points of the functional $\mathcal{I}^{\varepsilon}$ via the Euler-Lagrange formulation (26). A powerful tool for such computations is AUTO [12]. AUTO 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 [10, 11, 20, 23, 24, 45].

The first task one has to deal with is the construction of a starting orbit for fixed $\varepsilon\neq 0$. For (26) this is actually a less trivial task than for the FHN equation as we are going to explain in Section 3.1. In Section 3.1, we are also going to construct a starting periodic orbit based upon the geometric insights of Section 2.

Once the starting periodic orbit is constructed, we use AUTO to perform numerical continuation in both parameters $\mu$ and $\varepsilon$. 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 (3) predicted by Müller [36] is performed.

3.1. 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 $\mu$ and for $\varepsilon=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 $0<\varepsilon\ll 1$ is hard to perform since the slow manifolds are of saddle type and an orbit computed numerically would diverge from them exponentially fast [20].

• •

Matching slow segments obtained with a saddle-type algorithm [20, 32] 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 [22, 8], the periodic orbits we are looking for cannot be detected as Hopf bifurcations from the zero equilibrium. In that case, we could use AUTO to locate such bifurcations and then find a periodic orbit for $0<\varepsilon\ll 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 (5).

• •

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 2 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 $\mu$ 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 ${\mathcal{C}}^{\mu}_{\varepsilon,l}$ and ${\mathcal{C}}^{\mu}_{\varepsilon,r}$. In the last step we employ numerical continuation to study the values of $\mu$ we are actually interested in; this is the homotopy step.

The value $\mu=-\frac{1}{8}$ is peculiar, since in this case the singular slow segments for (26) reduce to two points

[TABLE]

i.e., touch-down and take-off sets for the fast dynamics coincide in this case. The range (24) we are considering does not include $\mu=-\frac{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 2 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 ${\mathcal{C}}^{\mu}_{\varepsilon,l}$ and ${\mathcal{C}}^{\mu}_{\varepsilon,r}$ as an initial value. The SMST algorithm [20] 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 $\varepsilon$ and $\mu$, we select manifolds $B_{l}$ and $B_{r}$, which are transverse to the stable and unstable eigenspaces of $\mathcal{C}_{0,l}^{\mu}$ and $\mathcal{C}_{0,r}^{\mu}$, respectively (Figure 7). The plane $B_{l}$ and the line $B_{r}$ provide the boundary conditions for the SMST algorithm.

Implementing the algorithm for $\mu=-\frac{1}{8}$ and $\varepsilon=0.001$ for (26) 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 $\varepsilon\neq 0$ as well as in the critical manifold ${\mathcal{C}}_{0}$. From the geometric analysis in Section 2 we know that at $\mu=-\frac{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 (5) 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.

3.2. Continuation in $\mu$

A detailed analysis of the critical points’ dependence on the Hamiltonian is performed. The value of $\mu$ can be arbitrarily chosen only in the interval $I_{\mu}$, while $\varepsilon$ is fixed to $0.001$. Continuation is performed on system (26) using the initial orbit obtained numerically in Section 3.1. Starting at $\mu=-\frac{1}{8}$, AUTO is able to compute the variation of the orbits up to $\mu=\frac{1}{24}$. The bifurcation diagram of the period $P$ with respect to the parameter $\mu$ is shown in Figure 8(a).

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

Furthermore, numerical continuation robustly indicates that the upper branch has another fold when continued from $\mu=0$ to higher values of $\mu$ as shown in Figure 9(a). The orbits obtained by fixing a value of $\mu$ 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 9(a). We conjecture that these parts arise due to the loss of normal hyperbolicity at ${\mathcal{L}}_{\pm}$; see also Section 4.

3.3. Continuation in $\varepsilon$

We perform numerical continuation in $\varepsilon$ by fixing three values of $\mu$ in order to capture the behavior of the solutions for the range $I_{\mu}$ from Proposition 2.5. We consider $\mu_{l}\approx-\frac{1}{8}$ with $\mu_{l}>-\frac{1}{8}$, $\mu_{c}\approx 0$, and $\mu_{r}\approx\frac{1}{24}$ with $\mu_{r}<\frac{1}{24}$; or more precisely $\mu_{l}=-0.12489619925$, $\mu_{c}=1.5378905702\cdot 10^{-5}$, and $\mu_{r}=0.04100005066$. For each of these values, we find two bifurcation branches connected via a fold in $(\mu,P)$-space; see Figure 10.

The bifurcation diagrams and the associated solutions shown in Figure 10 nicely illustrate the dependence of the period on the singular perturbation parameter $\varepsilon$. When $\varepsilon\rightarrow 0$ there are two very distinct limits for the period $P=P(\varepsilon)$ (Figure 10(a)-(b)) depending whether we are on the upper and lower parts of the main branch of solutions. In the case with $\mu_{r}\approx\frac{1}{24}$ when orbits come close to non-hyperbolic singularities on ${\mathcal{C}}_{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 10(c)). Furthermore, functional forms of $P(\varepsilon)$ are clearly different for small $\varepsilon$ 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 $\varepsilon$ of the periodic orbits in $(w,z,u)$-space is also interesting. For $\mu=\mu_{l}$ (Figure 10(a)), we observe that the upper branch corresponds to the solutions that we expect analytically from Proposition 2.5 consisting of two fast and two slow segments when approaching $\varepsilon=0$. A similar scenario occurs also for the other values of $\mu$ (Figure 10(b) and Figure 10(c)). When the $\varepsilon$ 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.

3.4. 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 $\mu$ and $\varepsilon$ should allow us to obtain the true minimizers via a double-limit. In other words, for every value of $\varepsilon$ there is a corresponding orbit which minimizes the functional $\mathcal{I}^{\varepsilon}$, and since along this orbit the Hamiltonian has to constantly assume a certain value $\bar{\mu}$, 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 $\varepsilon$ and $\mu$ via a direct continuation in both parameters, starting from certain special points, such as the fold points detected in Sections 3.2-3.3. However, it turns out that this process does not lead to the correct scaling law for minimizers of ${\mathcal{I}}_{\varepsilon}$ as shown in Figure 11.

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

The reason why from this naive approach the ${\mathcal{O}}(\varepsilon^{1/3})$ leading-order scaling does not emerge lies in the lack of connection with the minimization process. However, Figure 12 demonstrates that there are several nontrivial scalings of natural families of periodic orbits as $\varepsilon\rightarrow 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 $\varepsilon$ one is considering and the value of $\mu$ of the minimizers. In other words, there is not a unique value of $\mu$ given by the minimizers for every $\varepsilon$ small but minimizers move over different Hamiltonian energy levels as $\varepsilon\rightarrow 0$. Starting from this consideration, another option, which turns out to be the correct one to recover the scaling (3), is to use the periodic orbits from numerical continuation to compute the numerical value of the functional $\mathcal{I}^{\varepsilon}$ as a function of the period $P$ fixing different values of $\varepsilon$ in a suitable range, such as:

[TABLE]

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

4. 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 $\varepsilon$ 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 (10); 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 ${\mathcal{L}}_{\pm}$ 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 ${\mathcal{L}}_{\pm}$. The blow-up method [14] is likely to provide an excellent tool to resolve the non-normally hyperbolic singularities; see e.g. [31, 35] 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 [20] 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 $\varepsilon$ 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).