Continuation and bifurcations of periodic orbits and symbolic dynamics in the Swift-Hohenberg equation
Jakub Czwórnóg Daniel Wilczak

*
    Axel Klawonn0000-0003-4765-7387 
 Martin Lanser0000-0002-4232-9395
    
======================================================================

§ ABSTRACT 
  Steady states of the Swift–Hohenberg <cit.> equation are studied. For the associated four–dimensional ODE we prove that on the energy level E=0 two smooth branches of even periodic solutions are created through the saddle-node bifurcation. We also show that these orbits satisfy certain geometric properties, which implies that the system has positive topological entropy for an explicit and wide range of parameter values of the system. 
  
  The proof is computer-assisted and it uses rigorous computation of bounds on certain Poincaré map and its higher order derivatives.




§ INTRODUCTION

The Swift-Hohenberg equation <cit.> is a fundamental partial differential equation that plays a crucial role in the study of pattern formation and models various phenomena in physics and biology <cit.>. Originally the one-dimensional equation has the form

    ∂ U/∂ t = 
    	-  ( ∂^2/∂ X^2 + 1  )^2 U + α U - U^3.

The existence of some types of stationary solutions such as periodic and homoclinic  <cit.> in one-dimensional case is studied analytically. There are also many numerical simulations that demonstrate the existence of complicated structures and behaviour <cit.>. However, many observed phenomena are not proved by means of mathematical rigour. 

The aim of this article is to reproduce and extend results from <cit.> about one-dimensional stationary Swift-Hohenberg <cit.> equation


    -U”” - 2 U” + (α - 1)U - U^3 = 0.

This equation has conserved energy

    E = U”' U' - 1/2 (U”)^2 + U'^2 - α - 1/2 U^2 + 1/4 U^4 + (α - 1)^2/4.

The authors of <cit.> proved the following theorem.
 <cit.> 
	The dynamics of the ODE (<ref>) on the energy level E=0 is chaotic for all α≥ 2 in the sense that certain Poincaré map is semi-conjugated to a subshift of finite type with positive topological entropy.

The proof in <cit.> splits into two parts. First, the authors prove <cit.> that the existence of a periodic orbits of (<ref>) with the parameter value α>3/2 on the energy level E=0 and satisfying certain geometric properties (see Section <ref> for details) implies the existence of symbolic dynamics for (<ref>). Then in <cit.> the authors show, that the assumptions of <cit.> are checked for all α≥ 2. The proof of <cit.> is computer-assisted and it is based on so-called radii polynomial approach. The authors reformulate the problem of the existence of a periodic orbit as a zero-finding problem for the Fourier coefficients of the orbit to be found. Then, by means of interval arithmetics <cit.> and Newton-Kantorowitch type argument, it is shown that this infinite-dimensional system of algebraic equations has a branch of isolated solution parameterized by  α≥ 2. 

The aim of this paper is to prove the following result.

	There exists α^*∈ 1.9690842080_101989^293001 and two smooth curves U_±:[α^*,∞)→ℝ^4 such that 
	
		
  * E(U_±(α))=0 for α≥α^*, 
		
  * the solution to (<ref>) with the parameter value α and with the initial condition (U,U',U”,U”')=U_±(α) is an even and periodic function satisfying geometric properties from <cit.>,
		
  * U_-(α)≠ U_+(α) for α>α^* and U_-(α^*)=U_+(α^*),
		
  * the fold bifurcation occurs at α^*.
	

For the proof of Theorem <ref> we reformulate the problem of the existence of an even periodic orbit as a zero-finding problem for an univariate map (certain Poincaré map). In that sense, the proposed approach is much simpler and geometric as we work directly in the phase space of the ODE rather than an infinite dimensional space of Fourier coefficients. Moreover, computer-assisted verification of the result is significantly faster than the one presented in <cit.> – less than 13 minutes for two branches of periodic orbits and the fold bifurcation versus 12 hours for one branch.

The paper is organised as follows. In Section <ref> we recall the required geometric properties of these periodic orbits from <cit.>. In Section <ref> and derive a simple scalar equation for even periodic orbits of (<ref>). In Section <ref> we present details of the computer-assisted proof of Theorem <ref>.




§ SYMBOLIC DYNAMICS FROM A SINGLE PERIODIC ORBIT


For self-consistency of the article we recall here the forcing theorem <cit.> and required geometric conditions for periodic orbits.

The system (<ref>) has two equilibria ±√(α-1) on the energy level E=0. For 1<α≤3/2 they are stable foci (pure imaginary eigenvalues) and for all α>3/2 they are of saddle-focus type. Thus, for α>3/2 near these equilibria families of hyperbolic (on isolated energy level) periodic orbits appear giving rise to more complicated dynamics. 

In order to make the interval of parameters bounded, and thus much easier for computer analysis, following the ideas from <cit.> we perform a change of coordinates 

    y = X/√(α - 1),    u(y) = U(X)/√(α - 1),   ξ = 2/√(α - 1).

Now the parameter range α≥3/2 corresponds to 0<ξ≤√(8) and (<ref>) becomes 

    -u”” - ξ u” + u - u^3 = 0

with the energy

    E = u”' u' - 1/2(u”)^2+ ξ/2(u')^2 + 1/4(u^2 - 1)^2.


 <cit.> 
	Let ξ∈ [0,√(8)) and suppose that there exists a periodic solution ũ of (<ref>) at the energy level E=0, satisfying the following geometric conditions
	
		
  * ũ has exactly four monotone laps in one period and extrema ũ_1,ũ_2,ũ_3,ũ_4;
		
  * ũ_1, ũ_3 are minima, and ũ_2,ũ_4 are maxima;
		
  * ũ_1 < -1 < ũ_3 < 1 < ũ_2, ũ_4;
		
  * ũ is symmetric at its minima.
	
	Then the system is chaotic in the sense that there exists a two-dimensional Poincaré return
	map which has a compact invariant set on which the topological entropy is positive.


Geometric conditions about periodic orbits are illustrated in Fig.<ref>. The solutions that belong to the chaotic invariant set resulting from Theorem <ref> are encoded by their extrema. While the maximum is always bigger than 1 we are free to choose minimum to be either less than -1 or between -1 and 1. This leads to symbolic dynamics on three symbols, where each symbol is a building block of solution – as illustrated in Fig. <ref>.







§ EVEN PERIODIC ORBITS AS A ZERO-FINDING PROBLEM.


Equation (<ref>) can rewritten as a system of first order ODEs

    x'(t)    = y(t) 
    
    		y'(t)    = z(t) 
    
    		z'(t)    = w(t) 
    
    		w'(t)    = -ξ z(t) + x(t) - x^3(t).

Define a Poincaré section

    Π = {(x,y,z,w)  :  y = 0 }.

and let 𝒫_ξ:Π→Π be the associated Poincaré map for (<ref>) with the fixed parameter value ξ. The choice y=u'=0 of Poincaré  section is very natural when we are looking for periodic orbits with fixed number of extrema. Then, periodic points for 𝒫_ξ of principal period n correspond to functions with exactly n local extrema in one period.

Another property of this section is that Fix(R)⊂Π, where 

    R(x,y,z,w) = (x,-y,z,-w).

is a reversing symmetry of (<ref>). It is well known <cit.> that in such case the Poincaré map is also R-reversible. In order to find an even (R-symmetric) periodic orbit of (<ref>) it suffices to find a point u=(x,0,z,0)∈Fix(R) such that 𝒫_ξ^k(u)∈Fix(R). Then 𝒫_ξ^2k(u)=u and the trajectory of u is periodic and R-symmetric.

Hence, even periodic solutions of (<ref>) correspond to solutions of the scalar equation

    π_w 𝒫^2_ξ(x,0,z,0) = 0,

where π_w is the projection onto w variable. Observe that the intersection of Poincaré section with the energy level E=0

    E(x,0,z,w) = -1/2z^2+1/4(x^2-1)=0
 
gives the following relation

    u”=z = ±1/√(2)(x^2 - 1).

This means, that a solution u at the energy level E=0 has always proper extrema provided they are taken at x≠±1. Shifting minumum ũ_1<-1 to t=0 (see geometric conditions in Theorem <ref> and Fig. <ref>) we must take z(0)=1/√(2)(x(0)^2-1)>0. Substituting (<ref>) to (<ref>), we eventually transform the question of finding even periodic solutions of (<ref>) to a zero finding problem of the following equation

    G(ξ,x) =  (π_w ∘𝒫_ξ^2  )  (x,0, 1/√(2)(x^2 - 1), 0 ) = 0.


Finally, we have to check that an orbit corresponding to a solution of G(ξ,x_0)=0 satisfies the geometric conditions from Theorem <ref>. It is sufficient to check if

    (x_0 < -1) ∧ (x_1 > 1)∧ (-1<x_2<1),

where 
x_i= π_x𝒫_ξ^i(x_0,0,1/√(2)(x_0^2 - 1),0), i=1,2. Indeed, by (<ref>) the function has at most four extrema in one period. From (<ref>) the minima x_0≠ x_2 are different. Hence, the solution has exactly four extrema in period. All of them are proper because x_i≠±1, i=0,1,2. 







§ COMPUTER-ASSISTED PROOF OF THEOREM <REF>

As already mentioned (see (<ref>)) the parameter range α≥3/2 corresponds to 0<ξ≤√(8). Numerical simulation shows (see Fig. <ref>, Fig. <ref> and <cit.>) that a saddle-node bifurcation occurs at ξ≈ 2.0316516135713902 creating two branches of even periodic orbits. These branches continue to exist until ξ=0 (in fact they continue until ξ≈ -2.3 but we restrict here to the range ξ≥ 0 due to change of variables (<ref>)).


Computer-assisted proof of Theorem <ref> is split into the following steps.

	
  Step 1. We fix a threshold value =266 291 · 2^-17≈2.0316390991210938 ( is IEEE-754 <cit.> representable number) and prove that there are two smooth branches x_±:[0,]→ℝ such that for all ξ∈[0,] there holds x_-(ξ)<x_+(ξ), G(ξ,x_±(ξ))≡ 0 and x_±(ξ) satisfies (<ref>).
	
  Step 2. We prove that there is a smooth and concave function ξ̃:[x_-(),x_+()]→ℝ such that for x∈ [x_-(),x_+()] there holds G(ξ̃(x),x)≡ 0 and x satisfies (<ref>) in the system with parameter ξ̃(x).
	
  Step 3. We prove that ξ̃ has unique maximum ξ^*=ξ̃(x^*) for some x^*∈[x_-(),x_+()]. 
	
  Step 4. We prove that graphs of x_± and ξ̃ glue into a smooth curve in the (ξ,x) plane, as shown in Fig. <ref> and Fig. <ref>.


Finally, from these steps we will conclude that the functions x_± can be extended beyond  to ξ^*, at which value they are equal and the saddle-node bifurcation occurs at (ξ^*,x^*). Via change of variables (<ref>) we obtain the assertion of Theorem <ref>.



 §.§ Interval Newton Operator

In a computer-assisted verification of Steps 1–4 we will use the standard Interval Newton Operator for proving the existence and uniqueness of zeros of maps <cit.>.

	Let f:X⊂ℝ^n→ℝ^n be 𝒞^1 smooth, where X is compact and convex. Fix x_0∈int X and let [A] be an interval matrix such that 
	
    {Df(x) : x∈ X}   ⊂    [A].

	If the interval Newton operator
	
    N := x_0 - [A]^-1f(x_0)

	is well defined and N⊂intX then f has unique zero x in X and x∈ N.

Using implicit function theorem and Theorem <ref> it is easy to prove the following extension of Theorem <ref> to parameterised functions.


	Let f:Z× X⊂ℝ^k×ℝ^n→ℝ^n be 𝒞^1 smooth, where Z,X are compact and convex. Fix x_0∈int X and let [A] be an interval matrix and [e] be an interval vector such that 
	
    {D_xf(z,x) : (z,x)∈ Z× X}   ⊂    [A],
    {f(z,x_0) : z∈ Z}   ⊂    [e].

	If the interval Newton operator
	
    N := x_0 - [A]^-1[e]

	is well defined and N⊂intX then the set of zeroes of f in Z× X is a graph of a 𝒞^1 smooth function x:Z→ N⊂ X. That is, for (z,x)∈ Z× X
	
    f(z,x)=0⟺ x=x(z).





 §.§ Details of verification of Steps 1–4.


Step 1. To prove the existence of two smooth curves  x_±:[0,]→ℝ solving (<ref>) we perform an adaptive subdivision of the parameter range [0,]=Ξ^1_±∪⋯∪Ξ_±^N_±, where Ξ^j_± are closed intervals satisfying minΞ_±^1=0, maxΞ_±^N_±= and maxΞ_±^j=minΞ_±^j+1 for j=1,…,N_±-1. The sizes of subdivisions N_-=24157, N_+=142821 are returned by our adaptive algorithm. Then for each subinterval Ξ^j_± of parameters we construct (using continuation algorithms) a closed interval X_±^j and using Theorem <ref> we prove that the solution set of (<ref>) in Ξ^j_±× X_±^j is a smooth curve x_±^j:Ξ_±^j→ X_±^j. Moreover (<ref>) holds true for every point on this curve. Finally we show, that 

    X_±^j∩ X_±^j+1   ≠   ∅,   j=1…,N_±-1.

Because maxΞ^j_±=minΞ_±^j+1 and by the uniqueness property of the Interval Newton Operator we conclude that x_±^j(maxΞ_±^j)=x_±^j+1(minΞ_±^j+1) for j=1,…,N_±-1. Hence, the functions x_±^j, defined on subintervals Ξ_±^j of the parameter range join into continuous functions x_± defined on [0,]. Sine at every ξ∈[0,] we have D_xG(ξ,x_±(ξ))≠ 0 (the interval Newton operator is defined), by the implicit function theorem these functions are smooth. 

Finally, using Theorem <ref> we validate that

    x_+()∈ -1.5824440318_2912^327613,       
    	x_-()∈ -1.58253506275_30319^63035.

From these bounds it is clear, that x_-()< x_+() and thus, by the uniqueness property of the interval Newton operator this relation extends to x_-(ξ) < x_+(ξ) for all ξ∈[0,].

Step 2.
To prove the existence of a concave curve ξ̃ : [x_-(), x_+()] → solving (<ref>) we need to compute second order derivatives of G with respect to the parameter ξ. Because of the interface of the CAPD library <cit.>, in order to compute derivatives with respect to a parameter we have to treat it as a variable. Therefore we extend the system (<ref>) to

    x'(t)    = y(t) 
    
    		y'(t)    = z(t) 
    
    		z'(t)    = w(t) 
    
    		w'(t)    = -ξ(t) z(t) + x(t) - x^3(t) 
    ξ'(t)    = 0.

We are going to seek the zeros of the modified function

    G̅(ξ,x) =  (π_w ∘𝒫^2  )  (x,0, 1/√(2)(x^2 - 1), 0, ξ),

where 𝒫 : Π̅→Π̅ is the Poincaré map for the section

    Π̅ = {(x,y,z,w,ξ) : y = 0 }.

Clearly zeroes of G are in on-to-one correspondence with zeroes of G̅. 




The rest of this step is similar to the Step 1 with the addition of verification of concavity. Denote by x^_- the lower bound for x_-() and by x^_+ the upper bound for x_+() – see (<ref>). To find a curve solving (<ref>) we perform an adaptive subdivision of the range X_*:=[x^_-, x^_+] = ⋃_j=1^N X^j into N = 6903 closed intervals with max X^j=min X^j+1, j=1,…, N-1. Then for each subinterval X^j we construct an interval Ξ^j and using Theorem <ref> we prove, that the solution set of (<ref>) in Ξ^j× X^j is the graph of a smooth function ξ̃^j:Ξ^j→ X^j. We also check Ξ^j∩Ξ^j + 1≠∅ for j = 1,...,N - 1. By the uniqueness property of the interval Newton operator the functions ξ̃^j join into a smooth function ξ̃:X_*→ such that for x∈ X_* there holds G̅(ξ̃(x),x)=0 and x satisfies (<ref>). 

Additionally, on each subinterval X^j we compute a bound on the second order derivative of the implicit function ξ̃. For this purpose we use the 𝒞^r-Lohner algorithm <cit.> for integration of second order variational equations for (<ref>) needed for the second derivatives of Poincaré map 𝒫. From these computation we obtain a bound

    ξ̃”(x) ∈  [-74010.849232287583, -12744.872650106316],

for x∈ X_*. Hence the function ξ̃ is concave on the interval [x^_-,x^_+].

Step 3. We apply Theorem <ref> to the following equation

    H(ξ,x) = (G(ξ,x),∂ G/∂ x(ξ,x))=0.

Using standard Newton method we have found an approximate solution

    (ξ_0,x_0) = (2.0316516135713902,-1.5824941113082425)

to H(ξ,x)=0. Then, using Theorem <ref> we proved that there is a unique zero (ξ^*,x^*) of H in the set (ξ_0,x_0)+[-1,1]^2· 10^-10 that belongs to 

    ξ^*   ∈    2.0316516135_613893^814116, 
    
    	x^*   ∈    -1.582494111^3301776_2863635.

Using bound (<ref>) and change of variables (<ref>) we obtain a bound for α^* in Theorem <ref>. Finally using bounds (<ref>)-(<ref>) we checked, that (ξ^*,x^*) belongs to one of the sets Ξ^j× X^j used to verify the existence of ξ̃ in Step 2. Hence, by the uniqueness property of the interval Newton operator, we obtain ξ^*=ξ̃(x^*). From the implicit function theorem we have ξ̃'(x^*)=-∂ G/∂ x(ξ^*,x^*)/∂ G/∂ξ(ξ^*,x^*)=0. Given that ξ̃ is concave we conclude, that ξ̃ has unique maximum at ξ^*.

Step 4.
By the construction of ξ̃ we know that x_±()∈ X_*.
To show, that the graphs x_± and ξ̃ glue into a smooth curve in the (ξ,x) plane it suffices to check that 

    ξ̃(x_±())=.

Indeed, because the interval Newton operator is well defined in Steps 1–2, both partial derivatives ∂ G/∂ x(,x_±()) and ∂ G/∂ξ(,x_±()) are nonzero. Hence the curves x_±() can be extended in a smooth and unique way in a neighborhood of  and therefore graphs of ξ̃ and x_± must locally coincide near (,x_±()).

In order to check (<ref>), using bounds (<ref>) we verify that 

    (,x_-())∈Ξ^1× X^1  and   (,x_+())∈Ξ^N× X^N.

Then, from construction of ξ̃ in Step 2 we obtain (<ref>).



 §.§ Implementation notes.

In the computer-assisted verification of Steps 1–4 we need to compute bounds on Poincaré map and its derivatives up to the second order. For this purpose we used interval arithmetics <cit.>, algorithms for rigorous integration of ODEs and associated (higher order) variational equations and algorithms computation of Poincaré maps <cit.> all implemented in the CAPD library <cit.>.

The C++ program that performs verification of Steps 1–4 is a supplement to this article and available at <cit.>. 

The computation of the bound for the curve x_-(ξ) takes approximately 109 seconds. This curve corresponds to the solution from <cit.>. The curve x_+(ξ) takes less than 10 minutes (537 seconds) to compute, and the curve ξ̃(x) 1.5 minutes (100 seconds). Overall, the total computation takes less than 13 minutes.



§ CONCLUSIONS AND FUTURE WORKS

In this paper we reproduced and extended some results from <cit.> about dynamics of the stationary Swift-Hohenberg equation on the energy level E=0. After change of variables (<ref>) we obtain a system with a parameter ξ related to original parameter α. Numerical simulation shows that the two curves resulting from Theorem <ref> continue to exist for negative values of ξ, which apparently undergo further bifurcation near ξ≈ -2.3. 

Our method is computationally less expensive than the one proposed in <cit.>. This gives a hope to obtain results as in Theorem <ref> for a range of energy levels and, perhaps, study codimension two bifurcations of this family of periodic orbits. 

99
	
SwiftHohenberg J Swift and PC Hohenberg. Hydrodynamic fluctuations at the convective instability. Physical Review A, 15(1):319, 1977.

	
TLIDI19941475 M. Tlidi and Paul Mandel. Spatial patterns in nascent optical bistability. Chaos, Solitons & Fractals, 4(8):1475–1486, 1994.

PhysRevLett.73.2978 J. Lega, J. V. Moloney, and A. C. Newell. Swift-Hohenberg equation for lasers. Phys. Rev. Lett., 73:2978–2981, Nov 1994.

MERON201270 Ehud Meron. Pattern-formation approach to modelling spatially extended
ecosystems. Ecological Modelling, 234:70–82, 2012.

Glebsky Lerman LM Glebsky LY. On small stationary localized solutions for the generalized 1-D Swift-Hohenberg equation, 1995.

BurkeKnoblach John Burke and Edgar Knobloch. Localized states in the generalized Swift-Hohenberg equation. Phys. Rev. E, 73:056211, May 2006.

Deng Shengfu Deng. Periodic solutions and homoclinic solutions for a Swift-Hohenberg	equation with dispersion. Discrete and Continuous Dynamical Systems - S, 9(6):1647–1662,
2016.

Yang Junxiang Yang and Junseok Kim. Numerical simulation and analysis of the Swift–Hohenberg equation by the stabilized Lagrange multiplier approach, 2022.

Su Jian Su, Weiwei Fang, Qian Yu, and Yibao Li. Numerical simulation of Swift–Hohenberg equation by the fourth-order compact scheme, 2019.

sanchez2013numerical S Sánchez Pérez-Moreno, S Ruiz Chavarría, and
G Ruiz Chavarría. Numerical solution of the Swift–Hohenberg equation. In Experimental and computational fluid mechanics, pages 409–416. Springer, 2013.

chaos Jan Bouwe Van Den Berg and Jean-Philippe Lessard. Chaotic braided solutions via rigorous numerics: Chaos in the Swift–Hohenberg equation. SIAM Journal on Applied Dynamical Systems, 7(3):988–1031, 2008.

Moore1966 Ramon E. Moore. Interval analysis. Prentice-Hall, Inc., Englewood Cliffs, N.J., 1966.

IEEE1788-2015 IEEE Computer Society. 1788-2015 - IEEE Standard for Interval Arithmetic, 2015.

lamb1992reversing Jeroen SW Lamb. Reversing symmetries in dynamical systems. Journal of Physics A: Mathematical and General, 25(4):925, 1992.

wilczak2003chaos Daniel Wilczak. Chaos in the Kuramoto–Sivashinsky equations—a
computer-assisted proof. Journal of Differential Equations, 194(2):433–459, 2003.

IEEE754-2019 IEEE Computer Society. IEEE Standard for Floating-Point Arithmetic. IEEE Std 754-2019, July 2019.

Neumaier1990 Arnold Neumaier. Interval methods for systems of equations, volume 37 of Encyclopedia of Mathematics and its Applications. Cambridge University Press, Cambridge, 1990.

capd Tomasz Kapela, Marian Mrozek, Daniel Wilczak, and Piotr Zgliczyński. CAPD::DynSys: a flexible C++ toolbox for rigorous numerical analysis of dynamical systems. Communications in nonlinear science and numerical simulation, 101:105578, 2021.

WilczakZgliczynski2011 Daniel Wilczak and Piotr Zgliczyński. 𝒞^r-Lohner algorithm. Schedae Informaticae, 20:9–46, 2011.

poincare Tomasz Kapela, Daniel Wilczak, and Piotr Zgliczyński. Recent advances in a rigorous computation of Poincaré maps. Communications in Nonlinear Science and Numerical Simulation, 110:106366, 2022.

repo <https://github.com/Jacob4leg/SH_bifurcation>.
	


