Rippling rectangular waves for a modiﬁed Benney equation

One parameter family of rectangular periodic traveling wave solutions are known to exists in a perturbed system of the modiﬁed KdV equation. The rectangular periodic traveling wave consists basically of front and back transitions. It turns out that the rectangular traveling wave becomes unstable as its period becomes large. More precisely, torus bifurcation occurs successively along the branch of the rectangular traveling wave solutions. And, as a result, a “rippling rectangular wave” appears. It is roughly the rectangular traveling wave on which small pulse wave trains are superim-posed.Thebifurcationbranchisconstructedbyanumericaltoruscontinuationmethod.Theinstabilityisexplainedbyusingtheaccumulationofeigenvaluesontheessentialspectrumaroundthestationarysolutions.Moreover,thecriticaleigenfunctionswhichcorrespondtothetorusbifurcationcanbecharacterizedtheoretically.


Introduction
We study a perturbed system: of the modified KdV (m-KdV) equation: on the whole line. Here, u = u(t, x) is a real-valued function of time t and space x ∈ R. Moreover, the constant ε is supposed to be positive. The m-KdV equations are well-known integrable systems appearing in the study of the KdV equation. There are two different types depending on the coefficient of the nonlinear term. In fact, the m-KdV equation (2) has traveling front solutions u(z) = ± √ 3c tanh √ c/2(x + ct) for an arbitrary c > 0. They are sometimes called "kink" and "anti-kink". On the other hand, the different version of the m-KdV equation, which is obtained by replacing the term −(1/3)∂ x u 3 in (2) by (1/3)∂ x u 3 , has solitary wave solutions.
The m-KdV equation is also obtained from a long wave approximation of a certain fluid dynamics problem. By taking higher order approximation terms we can also obtain (1) instead of (2). In fact, Komatsu and Sasa [12] obtained both (2) and (1) from the so-called optimal velocity model which characterizes traffic congestion. The unknown variable u means the density of the traffic vehicles in their formulation. Therefore, sequences of kink and anti-kink can be considered as the bands of traffic jams. Now, the behavior of the solutions of the equations (2) and (1) are similar in the sense that both have traveling front solutions and by taking ε smaller the front solutions to (1) become close to that of (2) as we see later. However, the perturbed equation has a remarkable property in the following sense. In fact, only one traveling front solution can persist for ε > 0 among the uncountably many front solutions to (2). We can say that the velocity and amplitude of the front solution to (2) are chosen by the balance of dissipative perturbation terms.
The similar kind of observation and analysis can be found in [1,11,15]. They deal with the so-called Benney equation: which is obtained by adding dissipative perturbations to the KdV equation: Let us roughly explain the reason why the amplitude of the KdV soliton solution: where A is an arbitrary positive number, is selected to survive for the perturbations.
Consider the time development of the two conservations M = udx and W = u 2 dx for the KdV equation. We can obtain the following from (3): It is clear that both M and W are constants along the traveling wave solution, therefore we can determine the amplitude by plugging (5) to dW/dt = 0. We are going to discuss this amplitude selection later mathematically for (1) and even for the periodic traveling wave solutions. In fact, we can conclude there is a one-parameter family of periodic traveling wave solutions for each small ε parameterized by the wavelength. However, most significant point in this paper is not this amplitude selection of periodic solutions for each period but the stability change of the periodic traveling wave solutions. We can observe a new solution bifurcates by the torus bifurcation along the branch of the periodic traveling wave solutions. It is roughly the rectangular periodic traveling wave solution on which small pulse wave trains are superimposed. And, moreover, these superimposed pulses move at the velocity which is different from the background rectangular wave. Therefore, it is not a traveling wave after the bifurcation. We call this solution the rippling rectangular wave. See Fig. 1.
One might consider that it is reasonable for the equation (1) to have such type of rippling rectangular wave by the following heuristic argument. Let us linearize (1) at the lowest and highest flat levels ± u 0 since the periodic traveling wave becomes close to the kink and anti-kink as its spatial period becomes large. By plugging u =ũ ± u 0 into the equation and neglecting the higher order terms we obtain Therefore, we may observe a similar instability aboutũ = 0 as the Benney equation so that pulse wave trains may bifurcate depending on the length of the flat part, α, and u 0 as well.
To approach the bifurcation we study solutions that are L-periodic with respect to x. In other words, we study (1) on the interval (0, L) under the periodic boundary condition. We show that the rectangular waves are unstable for sufficiently large L provided that α is less than some threshold. Moreover, we prove that the distance between two successive torus bifurcation points tends to 2π/ 1 − αu 2 0 as L goes to the infinity. The proofs of these theorems are based on a topological argument, which is simpler than an analytic approach and easy to understand the connection with numerical simulations. In addition to the theoretical results, we also show numerical results for (1) with α = 0. We compute the torus bifurcation points, from which the rippling rectangular waves bifurcate, as well as critical eigenfunctions. We compute numerically the branches of the rippling rectangular waves by applying the so-called parametrization method and the predictor-corrector method. We also compute the torus bifurcation curves drawn on the (L , α)-plane to observe the dependency on the parameter α. It may be noticed that the parameter α = 2/3 in the case of the reduced equation from the optimal velocity model by Komatsu and Sasa [12]. Therefore, the rippling waves may not be observed in their case. That means the nonlinear diffusion term − α∂ 2 x u 3 plays a role of stabilizing the rectangular wave. We have studied (1) with α < 2/3 from the mathematical interest. We, however, believe that whole information of the bifurcation structure may be useful to really understand the model. The rest of this paper is organized as follows. We present our main theorems in the next section. We postpone the proofs of the theorems to Sect. 4 for the sake of readability. Section 3 is devoted to numerical results. The results for Sects. 2-4 are discussed in Sect. 5. The existence of the rectangular waves is proved in Sect. "Appendix A". Finally, we present our numerical methods in Sect. "Appendix B".

Traveling waves
In this section, we consider the existence and stability of traveling wave solutions for (1). Introducing a moving frame by ξ = x + ct, we can rewrite (1) as Traveling waves then are steady states of (9), that is, they satisfy the ordinary differential equation(ODE): (9) can be rewritten as the first-order system in R 4 : where = d dξ and U = (u, v, w, z) ∈ R 4 . The following theorem shows that (11) has a family of periodic orbits and two heteroclinic orbits.
There are small constant ε * and smooth functions c(ε, L) and c ∞ (ε), 0 ≤ ε < ε * , such that the followings are satisfied for all 0 ≤ ε < ε * ; 1. There exists L 0 ≥ 0 such that for any L ∈ (L 0 , ∞), (11) with c = c(ε, L) has a periodic solution h per (ξ ; c(ε, L), L). 2. If c = c ∞ (ε), then (11) has two heteroclinic solutions h f (ξ ) and h b (ξ ) satisfying where ±U 0 (c ∞ (ε)) = (±u 0 , 0, 0, 0) are equilibria for (11) with u 0 = √ 3c ∞ (ε). 3. c(ε, L) converges to c ∞ (ε) as L → ∞, and hence, We present a proof of this theorem, which is based on [15] in Appendix. In the following arguments, we fix ε ∈ (0, ε * ), and we put c * = c ∞ (ε). Next, we consider the linearized stability problem of the stationary solutions {u per (ξ, c(L); L} for (9). The linearized equation of (9) associated with u per (ξ, c(L); L) is given by where u = u per (ξ, c(L); L). We denote the linear operator in the right hand side by L per,L , that is, L per,L is an operator on L 2 (I L ) whose domain is the Sobolev space where I L = [0, L]. We consider the eigenvalue problem of L per,L : We deonte the spectrum of L per,L by σ (L per,L ). In the above setting, σ (L per,L ) consists of only discrete eigenvalues because I L is compact. One of our results induces the spectral instability of periodic solutions u per with long period, that is, at least one eigenvalue of L per is contained in the right-half plane.
Theorem 2 There exists L * > 0 such that u per (ξ, c(L); L) are unstable for any L > L * and sufficiently small α ≥ 0.
Theorem 3 Suppose α is sufficiently small and nonnegative. There exists L ∞ > 0 and a sequence {h n } ∞ n=1 such that the intersection of σ (L per,L n ) and iR \ {0} is nonempty for any n = 1, 2, . . . , where L n is given by and h n converges monotonically to zero as n tends to infinity, where u 0 = u 0 (c * ).
This theorem suggests that u per may undergo the Hopf bifurcation at L = L n for each n, and L n+1 − L n ≈ 2π/ 1 − αu 2 0 holds for sufficiently large n. Here, the Hopf bifurcation in the moving coordinates corresponds to the torus bifurcation in the original coordinates.
Remark 1 It should be noticed that α needs not to be nonnegative. In fact, we can prove Theorem 1 even for negative α under the condition that some quantities are positive. See the denominators in (37) and (39) in Appendix A in detail. We have, however, assumed that α is nonnegative for simplicity in Theorems 1-3.
We shall show Theorems 2 and 3 in the later section. Instead, we introduce the numerical results in the following section to see the meanings of these theorems.

Numerical results
We present numerical results for (1) with the periodic boundary condition on (0, L). Especially, we focus on the cases where L < 100. In this section, we set parameters (α, ε) = (0, 5 × 10 −3 ), unless otherwise mentioned. We regard a rectangular wave with the velocity c for (1) as a time-periodic solution with the period L/|c| because this view is useful for our numerical bifurcation analysis. The rest of this section is organized in the following way. First, we present a bifurcation diagram of the rectangular wave for (1). In the bifurcation diagram below, we find torus bifurcation points, which are nothing but the Hopf bifurcation points in the moving coordinates (9). Next, we focus on individual solutions appearing from the torus bifurcation points. Finally, we present a result of two-parameter continuation for the torus bifurcation points with respect to (L , α), keeping ε = 5 × 10 −3 fixed. Details of numerical methods are presented in the appendix of this paper.  "TR" in the diagram, on the branch of rectangular wave in 65 < L < 100. Table 1 presents the values of L at the torus bifurcation points. The solutions on the branch coming from TRn are also referred to as TRn (n = 1, 2, . . . 8). Figure 3 presents profiles of the critical eigenfunctions at the torus bifurcations. The number of maxima (and minima) increases as L increases, except first two bifurcations.
We classify the bifurcating solutions, which we call rippling rectangular waves, into two types: in-phase and anti-phase solutions, respectively. Let us illustrate it by Figs. 4 and 5. First, Fig. 4 illustrates snapshots of an in-phase solution for L = 86: 1. At t = 1.8, there are two humps on the lower plateau and two hollows on the higher plateau. 2. At t = 2.7, an additional hump and hollow appear simultaneously from the interface between lower and higher plateaus. 3. At t = 3.6, the leftmost hump and hollow reach the cliff and disappear simultaneously. 4. At t = 4.5, the solution recovers a similar profile as t = 1.8, while it is shifted to the left.
The solution repeats these process. We call this type of solutions in-phase because the humps and hollows are synchronous. Next, Fig. 5 illustrates snapshots of an anti-phase solution for L = 87.675611: 1. At t = 1.6, there are two humps and two hollows on the lower and higher plateaus, respectively. 2. At t = 2.4, an additional hump appears on the lower plateau. 3. At t = 3.2, the leftmost hump disappears at the right cliff. Again, there are two humps and two hollows on the lower and higher plateaus, respectively. 4. At t = 4.0, an additional hollow appears on the higher plateau.
The solution repeats these process. We call this type of solutions anti-phase because the humps and hollows appear and disappear alternately. Figure 6 presents profiles of the rippling rectangular waves at L = 100. Each of them has different number of humps. Figure 6a-d present the rippling rectangular waves TR1, TR4, TR6, and TR8, which are in-phase solutions, while (e-h), which appears from TR2, TR3, TR5, and TR7, are anti-phase ones.
Next, we consider the case where α > 0. Figure 7 presents profiles of traveling waves at L = 65 for some values of α. As α increases, the amplitude of the traveling wave decreases.

Proof of theorem
. Let L 0,L be the linearized operator about u 0 given by where ξ ∈ (0, L). Remark that the linearized operator about −u 0 coincides with L 0,L . Theorem 2 can be proved by the following two key facts. First, the uniform stationary solutions ±u 0 are unstable on the whole line R if the parameter α > 0 is sufficiently small since the essential spectrum intersects with the right half plane. Sec- ond, eigenvalues of the linearized operator L per,L accumulate to the essential spectrum of L 0,L when the period L of u per tends to infinity. This is because the periodic orbit corresponding to the traveling wave solution stays close to the uniform stationary solution for a longer and longer period as we take the period L larger and larger. These two are similar to the discussion on the stability analysis for the periodic orbit close to a homoclinic orbit by Gardner [6]. Sandstede and Scheel [16] also discussed about the similar problem with an extension. In this paper, we consider a family of periodic orbits close to a heteroclinic cycle; This is different from the previous studies. Moreover, we take a topological approach, which has been developed in [19] and [18]. This simplifies the proof of Theorems 2 and 3, and this allows us to understand the connection of theoretical and numerical results as we shall discuss in Sect. 5. First, we consider the eigenvalue problem of L 0,L , In a similar way to the derivation of (11), we rewrite this equation in the normal form: where The characteristic polynomial of A 0 (λ) is given by

Definition 1
The essential spectrum of L 0,L is defined by Remark that Σ 0 ess coincides with the essential spectrum of L 0,L even if the operator L 0,L is defined on L 2 (R). Therefore, the following lemma holds.

Lemma 2
The essential spectrum Σ 0 ess of L 0,L is given by Hence, if αu 2 0 < 1, then Proof It is immediate from the dispersion relation of L 0,L .
Proof By a simple calculation, we obtain , Therefore, the assertion holds.
Let us also rewrite the eigenvalue problem of L per,L L per,L p = λp, Here, and u = u per (ξ, c(L); L). We also write the eigenvalue problem for the front traveling wave solution u * (ξ ), * = f, b as follows: It clearly holds that Let us define the subsets of R 4 as We take γ > 0 sufficiently small so that the sets Π ± f and Π ± b are local cross sections for the heteroclinic orbits h f and h b . Then it is clear that Π ± f and Π ± b are also local cross sections for h per (ξ ; c(L), L) as well for sufficiently large L.
Now let ξ ± f,b (L) are the times when the periodic orbit crosses the local sections: We can assume ξ − f = 0 without loss of generality. We can also assume depend on u alone and is independent of other variables v, w, z. Also F(U, c, ε, α) has a symmetry F(−U, c, ε, α) = −F (U, c, ε, α). Therefore, we can conclude that ξ and u per (ξ + ξ + b ; L) = −u per (ξ ; L). We now have Lemma 4 There is a number ξ ∞ such that tends to infinity as L goes to infinity.
Proof It is clear if we take ξ ∞ as the time interval for the front heteroclinic orbit h f (ξ ) passes from Π − f through Π + f . Let us prepare the following lemma, which means that the eigenvalues of the linearized operator L per,L accumulate to the essential spectrum Σ ess as L → ∞. Theorem 2 is a immediate consequence of the following lemma because S per ∩ {λ | Re λ > 0} = ∅ by Lemmas 2 and 3 if α is sufficiently small. Let λ 0 ∈ S per and take a disk B(λ 0 ; δ) centered at λ 0 with radius δ. This lemma follows from Lemmas 6-10 below.
Let us consider the extended ODE system of (19) to treat the boundary condition as a linear subspace.

The accumulation of eigenvalues occurs on the intervals
Therefore, we consider the eigenvalue problem locally on the intervals ] to see the accumulation of eigenvalues. If we write A per (ξ ; λ, L) = A 0 (λ) + B L (ξ ; λ, c), B L is negligible by taking L large enough and γ small enough.
Therefore, (23) can be approximated by on the both intervals.

Lemma 7 λ is an eigenvalue of L per,L if and only ifΦ(L
is a 1 dimensional complex vector space inĒ c (λ). Then the differential equation which is induced on CP 7 by (23) is controlled by the following restricted equation on P(Ē c (λ)/E e (λ)) ∼ = CP 1 : Here, E e (λ)is the generalized eigenspace corresponding to 0, and P(Ē c (λ)/E e (λ)) is projectivization of the quotient vector spaceĒ c (λ)/E e (λ). The above equation (25) can be described as after an appropriate transformation. Here, η = q/ p is an inhomogeneous coordinate on N + = {[p : q] | p = 0} ⊂ CP 1 and It can be similarly described as by using an inhomogeneous coordinate η = p/q on N − = {[p : q] | q = 0}. If L is sufficiently large and γ is small then |β L | is small. Let P(ξ ; λ) be a solution of (25) with P(ξ − b ; λ) = G 3 ∩Ē c (λ) and P per =Ē c (λ) ∩ U per . Now we have, Lemma 8 If P(L; λ) = P per then λ is an eigenvalue of L per .
We consider an ordinary differential equation which is induced on CP 1 from (24): Let P 0 (x; λ) be a solution of (26) with P 0 (ξ − b ; λ) ∈ G 3 . By taking the above inhomogeneous coordinate on CP 1 , we can express P 0 (L; . Choose δ > 0 sufficiently small so that the following two conditions hold: 1. B(λ 0 ; δ)\ S per consists of a disjoint union of half disks B 1 and B 2 , and Re μ * (λ) > 0 for any λ ∈ B 1 and Re μ * (λ) < 0 for any λ ∈ B 2 .

Lemma 10
There are δ > 0 and L * > 0 such that the following holds: for any L ≥ L * there is a number m (L) ≥ m such that ψ L is also an m (L)-covering map.
This completes the proof of Theorem 2.

it follows then, for any vector
by the symmetry. This claim is a special case of [16,Eq. (5.2)]. We then define a 4-dimensional linear subspace From the above argument, we have two solutions Y e and Y o to (24) satisfying the following properties. That is, Y e (0; λ 0 ), Y e (ξ + b ; λ 0 ), and Y e (L n ; λ 0 ) are contained in U per for an even n. On the other hand, Y o (0; λ 0 ) and Y o (L n ; λ 0 ) are contained in U per while Y o (ξ + b ; λ 0 ) ∈ U * per for an odd n. Translating it into the terms of the projective space, we have two solutions P e and P o to (26) with the initial data P per at ξ = 0 satisfying the following properties. P e corresponds to Y e and satisfies P e (ξ + b ; λ 0 ) = P e (L n ; λ 0 ) = P per for an even n. P o corresponds to Y o and satisfies P o (L n ; λ 0 ) = P per and P o (ξ + b ; λ 0 ) = P * per for an odd n. Let P(ξ ; λ) be a solution of (25) satisfying P(0; λ) = P per . Although P(L n ; λ 0 ) is no longer equal to P per , we can clearly adjust λ and L so that P(L; λ) = P per holds for (L , λ) ≈ (L n , λ 0 ) with λ ∈ iR \ {0}.
Summarizing the above arguments, there are complex conjugate eigenvalues of L per,L in iR near by Adjusting n, we can take L ∞ and h n so that (15) holds. Thus we obtain Theorem 3.
Remark 2 By using the solutions P e and P o to (26), we can distinguish to which type of torus bifurcation the critical eigenfunction leads. Indeed, if the index n is even, then the bifurcating rippling rectangular wave is anti-phase. On the other hand, an in-phase solution bifurcates if n is odd. Figure 9 illustrates these situations.

Discussion
We have demonstrated both theoretically and numerically that the torus bifurcation occurs successively on the branch of rectangular waves as the system size L increases.
Since the curve of essential spectrum crosses the imaginary axis we can construct critical eigenfunctions which correspond to the "TRn" torus bifurcation. The profile of the critical eigenfunction consists basically of the sinusoidal part and the transition part along the front. This is because the rectangular wave stays close to the two steady states ±u 0 except the two transition regions. The sinusoidal part can be characterized by the eigenfunction corresponding to a pure imaginary eigenvalue for the eigenvalue problem about ±u 0 . As mentioned in Remark 2, there are two types of critical eigenfunctions depending on their periods: L or L/2. The type of bifurcating solution is in-phase when the period of the critical eigenfunction is exactly L. While on the other hand, anti-phase solutions bifurcate when the period is L/2. Suppose that α = 0. If the first torus bifurcation is supercritical, then the solution TR1 is stable. This is natural from the viewpoint of bifurcation theory. The other rippling rectangular waves are considered to be unstable if L is near the corresponding torus bifurcation point. However, when L is far from the bifurcation point, the rippling rectangular waves can be observed stably as a numerical solution of an initial value problem for (1). Indeed, Fig. 5 is obtained by solving an initial value problem numerically. This suggests that they are stable at such a parameter value. An unstable rippling rectangular wave should recover its stability by some bifurcations. When α is varied from zero, as shown in Fig. 8, two torus bifurcation curves intersect. This means that two torus bifurcations occur simultaneously. This may be responsible for the secondary bifurcation at which the unstable rippling rectangular wave becomes stable. Therefore further study of whole bifurcation structure is necessary.
First of all, let us reformulate (10) as spatial dynamics. Integrating both sides of (10) with respect to ξ , we obtain where C is an integral constant. We concentrate on solutions with the symmetry u(ξ + L * ) = −u(x) for some L * . Therefore we assume C = 0. We also assume c > 0 because we regard (10) as a perturbed system of the modified KdV equation for which the velocity c of the nontrivial traveling wave solution is positive. Introducing new variable U and τ by u = √ cU and ξ = τ/ √ c, we rewrite (28) with C = 0 as a system of first order equations: where the dot means differentiation with respect to τ . A kink or anti-kink solution corresponds to a heteroclinic solution to (29), and a periodic traveling wave corresponds to a periodic solution to (29). Thus the problem is to find a heteroclinic solution and a periodic solution to (29). If we put ε = 0 formally, then we get U 3 /3 − U − W = 0 from the third equation of (29). Substituting W = U 3 /3 − U into (29) with ε = 0, we obtain a Hamiltonian system (34) below. It is easy to see that (34) has a heteroclinic orbit and a periodic orbit. We would like to continue these orbits for ε = 0. According to Fenichel's singular perturbation theory [5], as we will see below, there exists a two-dimensional invariant manifold of (29) in a neighborhood of {(U, V, W ) | U 3 /3 − U − W = 0} for sufficiently small ε. This allows us to consider such a perturbation problem on this manifold.
In order to apply Fenichel's theorem, let us rewrite the equations further. Introducing a new variable Z = U 3 /3 − U − W we obtain the following suspended system: where we consider ε as a variable. Moreover, introducing a new independent variable by s = τ/ε, we obtain where I is some open interval with 0 ∈ I . By the invariance of the center manifold, we obtain By substituting this into (30), we get the following planar system: We can regard (33) as a perturbed system of the following Hamiltonian system: This system has a Hamiltonian function H given by The level set {H = 3/4} corresponds to two heteroclinic orbits between the equilibria ( √ 3, 0) and (− √ 3, 0). For any 0 < κ < 3/4, the level set {H = κ} corresponds to a periodic orbit.
To complete the proof of Theorem 1, we need to prove that the wavelength L of the periodic solution to the perturbed system is monotone with respect to κ, at least for sufficiently large L. We can discuss the monotonicity in a similar manner as [2,3], and [15]. Let us briefly explain this.
Let Q, R, and S be where T = T (κ) is the period of the periodic solution to (34). Keep in mind that c 0 = Q/(R + αS). For simplicity of the argument, we parametrize the periodic solutions by K = 2κ instead of κ. Note that K is contained in (0, 3/2). Let us denote the two roots of 3) by β(K ) and −β(K ).
We can rewrite Q, R, and S in the following way: It is convenient to represent Q, R, and S by the following integrals: when α > −2/3. Therefore we have lim K →3/2−0 X (K ) = +∞, which implies that L is monotone with respect to K if K is sufficiently close to 3/2. This completes the proof of Theorem 1.

Appendix B: Numerical methods
In this section, we note our numerical methods. We discretize (1) by the Fourier spectral method. That is, we approximate a solution to (1) by for any t ∈ R and n ∈ {0, ±1, . . . , ±N }. In particular, u 0 (t) ∈ R. Remark that u 0 (t), which is the spatial average of u(t, x) over [0, L], is independent of time. Since we consider the case where u 0 ≡ 0, we have only to find u n for n = 1, 2, . . . , N to compute an approximate solution of (1). The Fourier coefficients are governed by the following ODEs:u n = −a n u n + g n (n = 1, 2, . . . , N ), where a n is given by and g n are Fourier coefficients of the nonlinear term of (1): We evaluate the nonlinear terms in the physical space and compute their derivative in the Fourier space. We used Ooura's fast Fourier transform (FFT) library to implement it. As u n (t) are complex-valued, (42) defines a 2N -dimensional continuous dynamical system. We apply the second order exponential Runge-Kutta method ( [14], [8, (2.40)] with c 2 = 1/2) for solving an initial value problem of (42).
Since a rectangular wave for (1) corresponds to a limit cycle, we apply AUTO-07p [4] for its parameter-continuation. AUTO-07p is a standard software for numerical bifurcation analysis, which supports path-following of a solution branch, stability check, and detection of bifurcation points for equilibria and limit cycles. Since it does not, however, support branching off from a torus bifurcation point, we have to implement it by ourselves.
We characterize a rippling rectangular wave (RRW) as a quasi-periodic solution to (1), which can be approximated by that of (42) on an invariant torus. We apply a parametrization method [10,17] for computing an invariant circle on a Poincaré section corresponding to the invariant torus. We define a Poincaré section S by S = {u(t, x) | Re (u 1 ) = 0, Re (u 1 ) > 0} , and consider a Poincaré mapP defined on on S. S is real 2N − 1 dimensional. Let be a chart on S, and define P : R 2N −1 → R 2N −1 by P = ϕ −1P ϕ. The invariant torus corresponds to an invariant closed curve (ICC) of P. we characterize it as a rigid rotation on an invariant closed curve for the Poincaré map.
We recall the parametrization method for computing the ICC [10,17]. The basic idea is to construct a topological conjugation v : T → R 2N −1 , θ → (v 1 (θ ), v 2 (θ ), . . . , v 2N −1 (θ )), between the rigid rotation of a circle and the Poincaré map: v(θ + ω) = P(v(θ )), where θ ∈ T 1 [0, 2π ] is the coordinate on T 1 and ω ∈ R is the angular rotation number. If v(θ) solves (43), then the image v(T) is an invariant closed curve of P. We for some M. As v j (θ ) is a real-valued function, v j,m = v j,−m for all j and m. Notice that ω is also unknown and that there is a degree of freedom corresponding to phase shift. Therefore we need to impose a phase condition. We adopt dv k dθ (0) = 0, for the phase condition, where k is a natural number. We take k = 1 for computing an in-phase solutions and k = 3 for computing anti-phase solutions. We employ the The Newton method is available for solving this system of equations. We need an initial guess for Newton iteration. In particular, it is a nontrivial task to select an initial guess of an unstable ICC. Since we are interested in bifurcations from the rectangular wave, we make use of the critical eigenfunctions at the torus bifurcation. Notice that the torus bifurcation is a Neimark-Sacker bifurcation of a fixed point of the Poincaré map P. Let {v * (L)|L > 0} be a family of fixed points of P, which is parametrized by the parameter L. Suppose that there exists L = L * such that the Jacobian matrix of P at v = v * (L * ) has a pair of complex conjugate eigenvalues located on the unit circle. Let φ be an eigenvector corresponding to the critical eigenvalue. We define an initial guess by u = u * + We apply a standard predictor-corrector technique with the secant prediction and the natural-parameter correction. We refer readers to [7], [13,Chapter 10] for detail.
We solve the initial value problem for (1) by the second order exponential Runge-Kutta method with a fixed step size h = 10 −3 . In this paper we present the numerical results for N = 100 and M = 8. Figure 10 shows the absolute value of Fourier coefficients of (1) with L = 100 computed for N = 100 and N = 200. As the tail of Remark 3 Figures in this paper were generated by gnuplot [20] and Python matplotlib [9].