- Email: [email protected]
Stationary travelling-wave solutions of an unstable KdV–Burgers equation
Physica D 137 (2000) 228–236
Stationary travelling-wave solutions of an unstable KdV–Burgers equation Bao-Feng Feng, Takuji Kawahara ∗ Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto 606-8501, Japan Received 9 March 1999; received in revised form 22 July 1999; accepted 24 August 1999 Communicated by Y. Kuramoto
Abstract Both “solitary” and “periodic” stationary travelling wave solutions are investigated numerically for an unstable Korteweg–de Vries–Burgers equation ut + uux + uxxx − η(u + uxx ) = 0 (η > 0). A family of stationary solitary wave solutions whose members are distinguished by the number of “humps” is found for a given η. Corresponding to each solitary wave thus found, a family of stationary periodic waves with the same number of “humps” exists under periodic condition and ends up in the infinite periodicity to the corresponding solitary wave. The numerical results are consistent with the theoretical estimates based on the conservation properties. ©2000 Elsevier Science B.V. All rights reserved. PACS: 03.40.Kf; 02.70.H Keywords: Unstable Korteweg–de Vries–Burgers equation; Solitary and periodic waves; Multi-hump solutions; Rational Chebyshev and Fourier pseudo-spectral method
1. Introduction Much attention has been paid to the roles of localized solutions in the initial value problem of the nonlinear dispersive equations involving both instability and dissipation. One prototype of such equations is the one-dimensional Benney equation: ut + uux + uxx + δuxxx + uxxxx = 0.
(1)
It has been found that in certain cases time evolutions of Eq. (1) are well approximated by an array of similar amplitude localized pulses (travelling wave solu∗ Corresponding author. E-mail address: [email protected] (T. Kawahara).
tions) [1,2]. The coexistence of instability and dissipation is responsible for determining both local and global behaviours. Locally the scale of pulse is determined by a balance between instability and dissipation terms of Eq. (1), which leads to a pulse with a definite amplitude, unlike the case of KdV solitons where arbitrary amplitudes are admitted. Globally the equilibrium states (either deterministic or chaotic) of the initial value problem are determined also by the balance between instability and dissipation. The balance of the rate of energy inflow by instability with that of energy outflow by dissipation determines an equilibrium state through energy transfer due to nonlinearity influenced by dispersion. For strongly dispersive cases, overall evolutions are deterministic but for relatively weakly
0167-2789/00/$ – see front matter ©2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 2 7 8 9 ( 9 9 ) 0 0 1 8 3 - 9
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
dispersive cases they are chaotic and the balance is in a statistical sense. On the other hand, recent researches have demonstrated the existence of complicated families of localized solutions for the Kuramoto–Sivashinsky equation [3–5] or the Benney equation [6,7]. These results pose a problem to be investigated, i.e., what relations exist between a variety of localized solutions and the overall evolutions (either deterministic or chaotic) of the initial value problem. This problem has not yet been explored in detail because of the complication and existence of too many families of solutions. In order to carry such an exploration a step further, it will be worthwhile to consider similar other nonlinear dispersive equations involving instability and dissipation. One of such equations is an unstable Korteweg–de Vries–Burgers equation: ut + uux + uxxx − η(u + uxx ) = 0,
(2)
where the parameter η is a positive constant, representing the importance of instability and dissipation relative to dispersion and nonlinearity. Eq. (2) can be regarded as a simple model for a so-called extrinsic instability, i.e., a longwave instability whose growth rate takes a maximum at zero wave number [8], or a physical system such as a current-driven ion acoustic instability in a collision-dominated plasma [9]. In this paper, both stationary solitary and periodic wave solutions to Eq. (2) are computed via a pseudo-spectral method. Existence of multiply humped solitary wave solutions to Eq. (2) in the range 0 < η ≤ 1.5 is discovered numerically. It is also found that every existing solitary wave solution corresponds to a family of periodic wave solutions. These periodic waves tend to the corresponding solitary wave provided that the periodicity length L → ∞. The remainder of this paper is organized as follows. In Section 2, a preliminary analysis is given to show some basic properties of Eq. (2). In Section 3, we describe the numerical results for the solitary waves and in Section 4 those for the periodic waves. Then, in Section 5, we give some discussions and comments. Finally, in Appendix A, the employed pseudo-spectral methods are introduced briefly.
229
2. Preliminary analysis Substitution of u ∝ exp(iκx + σ t) into Eq. (2) yields the linear dispersion relation σ = iκ 3 + η(1 − κ 2 ). Thereby the small-amplitude sinusoidal waves are unstable for long wavelengths given by |κ| < κc ≡ 1 and damped for short wavelengths |κ| > κc . This is similar to Eq. (1) but differs in that the maximum growth rate is given at zero wave number. We shall restrict our attention to stationary travelling-wave solutions of Eq. (2) over a spatially periodic domain x ∈ [0, L] or an unbounded domain x ∈ (−∞, ∞). These stationary waves are solutions of Eq. (2) which travel steadily at a constant speed c without changing their shapes. Therefore, it is appropriate to carry out a moving coordinate transformation ∂/∂t → ∂/∂t − c∂/∂x, which leads to the governing equation for stationary waves: −cux + uux + uxxx − η(u + uxx ) = 0,
(3)
where c is the wave speed. As preliminaries of the numerical study, we examine some basic properties of Eq. (3). Firstly, for the dynamical system corresponding to Eq. (3), only one fixed point O1 (0, 0, 0) is allowed in the three-dimensional phase space defined by coordinates (u, ux , uxx ). This is due to the existence of the unstable term proportional to u in Eq. (3). Therefore, only solitary waves, which correspond to homoclinic trajectories of O1 , are possible over the whole x domain −∞ < x < ∞. The boundary conditions are given by u(±∞) = 0,
ux (±∞) = 0,
uxx (±∞) = 0 (4a)
for solitary waves or by u(0) = u(L),
ux (0) = ux (L),
uxx (0) = uxx (L) (4b)
for periodic waves with periodicity length L. In order to investigate the asymptotic behaviour of u, we linearize Eq. (3) about O1 , which yields the characteristic equation for the characteristic exponent µ: µ3 − ηµ2 − cµ − η = 0.
(5)
230
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
We can distinguish the property of the roots to Eq. (5) by the criterion in terms of the following quantity: R = 4η4 + (18c − c2 )η2 − 4c3 .
(6)
A positive R corresponds to a real positive root and a pair of complex conjugate roots with negative real part, whereas a negative R corresponds to three real roots. Thereby, an oscillatory structure is expected in the tail of the solitary waves in the case of R > 0, whereas a monotonous structure in the case of R < 0. This property is confirmed in the numerical computations afterwards. Secondly, integration of Eq. (3) over x on the whole domain Ω (where Ω = (−∞, ∞) for a solitary wave or Ω = [0, L] for a periodic wave) with (4a) and (4b) leads to a necessary condition for the existence of a stationary wave to Eq. (2), i.e., hui = 0,
(7)
R∞ where h· · · i denotes −∞ · · · dx for a solitary wave or RL 0 · · · dx for a periodic wave with periodicity length L. Condition (7) is set to be satisfied throughout our computations of stationary wave solutions. Moreover, multiplying Eq. (3) by u and integrating it over x on the domain Ω, the relation hu2 i = hu2x i
Fig. 1. Numerical solitary wave solutions for Eq. (2) with η = 0.1: (a) c = 4.602758, (b) c = 4.281452, (c) c = 3.991830, (d) c = 3.728704, (e) c = 3.487565, (f) c = 3.235719.
(8)
can be derived after some manipulations with (4a) and (4b). Although no other conservation properties are possible because of the non-integrability of Eq. (2), the conservation properties (7) and (8) can provide us with a powerful tool for verifying the numerical results. As a matter of fact, the result of such a verification is very satisfactory as will be discussed in Section 5 in an attempt to obtain some qualitative explanations for the general features of the numerical results.
3. Numerical results for solitary waves The stationary solitary wave solutions of Eq. (2) were investigated by Kawahara and Toh [10] via a shooting method. They found that for a given η, one
particular single hump solitary wave solution exists. Moreover, by a perturbation analysis, such a solution is shown to approach the one soliton solution of the Korteweg–de Vries equation with a particular velocity c = 5.0 when η → 0. Here, we employ a rational Chebyshev pseudospectral method (see Appendix A) to compute the solitary wave solutions of Eq. (2). In addition to a single hump solitary wave found previously by the shooting method, we find a series of solitary wave solutions involving a multi-hump structure for a fixed value of η. Six solitary wave solutions (single peak to six peaks) for η = 0.1 are shown in Fig. 1 with eigenvelocities found numerically as c = 4.602758, 4.281452, 3.991830, 3.728704, 3.487565, 3.235719. For a relatively large η = 1.0, three numericals
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
231
4. Numerical results for periodic waves
Fig. 2. Numerical solitary wave solutions for Eq. (2) with η = 1.0: (a) c = 1.399297, (b) c = 0.243849, (c) c = −0.412683.
solutions are shown in Fig. 2 with c = 1.399297, 0.243849, −0.412683. It is worthwhile to note that similar to 1D Benney equation (1), multi-hump solitary wave solutions exist for Eq. (2). It can be seen that for η = 0.1, all of the six solutions have a negative plateau in their tail part, whereas for η = 1.0, an oscillatory structure arises. This result is consistent with the analysis of characteristic equation in Section 2. As the numerical computations become increasingly inaccurate and difficult as the number of “humps” increases or the value of η becomes large, only six solutions are identified for η = 0.1 and three for η = 1.0. We conjecture here that for every η ∈ (0, 1.5], there is a family of solitary wave solutions with a countable number of members. Members in each family are distinguished by the number of “humps”. Because of the numerical difficulty, it is not sure whether such solutions exist or not if η is beyond this range. For a given η, the eigenvelocity for the first member (single hump solitary wave) of the family is the largest, and it decreases as the number of “humps” increases. For the members with the same number of “humps”, the eigenvelocity decreases as η increases. Moreover, in contrast to the KdV solitons, the negative eigenvelocity is possible.
The stationary periodic wave solutions to Eq. (2) are computed by means of a Fourier pseudo-spectral method, which will be introduced briefly in Appendix A together with the rational Chebyshev pseudo-spectral method. We begin to seek a periodic solution with a relatively large wave number α. Here, α ≡ 2π/L denotes the minimum wave number allowed for the periodic system with periodicity length L. The initial profile is taken from the equilibrium state in the temporal evolution of Eq. (2) when it is integrated over the same periodicity length via an implicit spectral method developed by Wineberg et al. [11]. The Newton–Kantorovich iteration method is used in every time step. Then we initiate the continuation scheme towards the small wave number. Thus a family of periodic waves is constructed for a given η. One part of the numerical results with different wave number α is shown in Figs. 3 and 4 with η = 0.1, 1.0, respectively. It is found numerically that similar to the case of solitary waves, multi-hump periodic wave solutions exist. We denote n-hump periodic wave solution by P (n) . Moreover, we found that P (1) starts from one near-sinusoidal periodic wave with α = α ∗ ≈ 0.96664 and converges to the corresponding solitary wave solution as the wave number α → 0. For a fixed value of η, we conjecture that families of periodic wave solutions P (1) , P (2) , . . . , P (n) , . . . exist, which originate from one near-sinusoidal periodic wave with wave number α = α ∗ , α ∗ /2, . . . , α ∗ /n, . . . , respectively. Each family of periodic waves gradually tend to the corresponding solitary wave as the wave number α → 0. It is also interesting to note that the way of convergence is associated with the tail structure of the solitary waves. For the solitary wave with monotonous tail, the set of the eigenvelocities for the periodic waves converges to the velocity of the corresponding solitary wave monotonously from the lower side, whereas for the solitary wave with an oscillatory tail, it also converges to the velocity of the corresponding solitary wave in an oscillatory way.
232
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
Fig. 3. One-hump and two-hump periodic wave solutions for Eq. (2) with η = 0.1. One-hump: (a) L = 20, c = 3.020606, (b) L = 40, c = 3.966167, (c) L = 60, c = 4.201104. Two-hump: (d) L = 20, c = 1.189170, (e) L = 40, c = 3.029772, (f) L = 60, c = 3.572864. (Note the difference in horizontal scale.)
5. Discussion and some comments The stationary solitary wave solutions and the periodic wave solutions to an unstable Korteweg–de Vries–Burgers equation are investigated, respectively, by means of a rational Chebyshev pseudo-spectral method and a Fourier pseudo-spectral method. Existence of multiple stationary wave solutions for both solitary and periodic waves is discovered numerically, and the relationship between them is also found. Within the numerical computations carried out so far, the following general features have been found: 1. All the velocities for the solitary waves or periodic waves are less than 5.0 and negative velocity is possible.
Fig. 4. One-hump and two-hump periodic wave solutions for Eq. (2) with η = 1.0. One hump: (a) L = 15, c = 1.454888, (b) L = 20, c = 1.395870, (c) L = 25, c = 1.399301. Two hump: (d) L = 25, c = 0.297846, (e) L = 27.5, c = 0.222642, (f) L = 40, c = 0.244304.
2. For a fixed value of η, there exists one family of solitary wave solutions to Eq. (2), in which it is conjectured that a countable number of members should exist. The members are distinguished with each other by the number of “humps”. The velocity for the first member is the greatest, then it decreases as the number of “humps” increases. 3. The velocity of the first member in the family of solitary waves approaches to the value of 5.0 as η → 0. For the member with the same number of “humps”, the corresponding velocity decreases as η increases. 4. For each solitary wave, one family of periodic waves with the same number of “humps” exists. These periodic waves gradually tend to the corre-
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
sponding solitary wave as the wave number α → 0. However, the way of the convergence of the eigenvelocities is different, and depends on the tail structure of the solitary waves. As Eq. (2) is nonintegrable and an analytical solution of stationary wave is not available, it has to be investigated numerically. However, some qualitative and quantitative explanations for the numerical results may be given by reviewing the basic properties of Eq. (2) such as (7) and (8). For example, it has been well explained in [10] by a perturbation analysis why the velocity c is close to the value of 5.0 when η becomes very small. Here, we derive an expression to estimate the velocity c of a stationary wave. Integration of Eq. (3) from the starting point of Ω to an arbitrary position ξ ∈ Ω leads to Z ξ 1 2 u dx + ux (ξ ) −cu(ξ ) + 2 u (ξ ) + uxx (ξ ) − η −∞
=0
(9a)
for a solitary wave or Z −cu(ξ ) +
1 2 2 u (ξ ) + uxx (ξ ) − η
= cu(0) −
ξ
u dx + ux (ξ )
0 1 2 2 u (0) − uxx (0) + ηux (0)
(9b)
Rξ
for a periodic wave. We set v(ξ ) = −∞ u dx for a Rξ solitary wave or v(ξ ) = 0 u dx for a periodic wave, and then obviously, we have dv/dξ = u. Then multiplying Eq. (9a) or (9b) by u(ξ ) and integrating it over ξ on the whole domain Ω, we have −chu2 i + 21 hu3 i + huuxx i − η(huvi + huux i) = 0, (10) here, the condition (7) is used. Thereafter, it can be shown readily that huux i = 0,
huvi = 0,
huuxx i = −hu2x i.
(11)
Therefore, using (8), we finally have c=
hu3 i − 1, 2hu2 i
(12)
which gives an estimate of the eigenvelocity in terms of the integrals of u.
233
We can divide the integral hui into two parts, i.e., the integral I1 over the domain where u is positive and the integral I2 over the domain where u is negative. The same procedure can be done for the integrals hu3 i and hu2 i too. Apparently, from (7), we have I1 > 0 and I1 = −I2 . It should be noted that the integrals of hu3 i in each domain keep the same sign as hui, but those for hu2 i are always positive. Bearing this fact in mind, we consider the stationary waves with the same η. As the number of “humps” increases, the negative level of the wave tail has to be lowered in order to maintain the zero mean value of u. Because of the cancellation between positive and negative contributions for hu3 i, the rate of increase in hu3 i is expected to be smaller than that in hu2 i. Thus it is reasonable that the quantity hu3 i/hu2 i will decrease, i.e., the velocity c will decrease. From the numerical results and the preliminary analysis, we have the fact that for a small value of η, there is a long and negative plateau in the tail part, whereas for a large value of η, an oscillatory structure arises in the tail part. Thereby, for the stationary waves with the same number of “humps”, the “basin” following the last “hump” is relatively deep and narrow for a large value of η, while being flat and broad for a small value of η. Henceforth, it is also reasonable to think that the quantity hu3 i/hu2 i will decrease as the value of η increases. One more point is that when the quantity hu3 i/hu2 i is less than 2, the velocity will be negative. The numerical values of velocity c versus wave number α for typical values of η = 0.1, 1.0 are shown in Figs. 5 and 6. The theoretical estimates of c from Eq. (12) are also computed and compared. Here, P (1) and P (2) denote the branches for the one-hump and two-hump stationary waves, respectively. Numerical integrals hu3 i and hu2 i are computed by Simpson’s rule. The numerical results for the rates of increase in hu2 i and hu3 i are consistent with the above-mentioned theoretical expectation. Figs. 5 and 6 show that the numerical values of c coincide with the theoretical estimates very well. In other words, this is the verification of our numerical results. Moreover, for each branch P (i) (i = 1, 2), the family of stationary periodic waves starts from
234
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
Fig. 5. Comparison of the numerical results and theoretical estimates for velocity c with η = 0.1. Solid lines: numerical results, (×): theoretical estimates for one-hump solutions, (•): theoretical estimates for two-hump solutions.
certain point in the α–c plane and ends up with the corresponding solitary wave. For the same value of η and wave number α, the velocity c of the two-hump stationary wave is smaller than that of one-hump one. When comparing Fig. 5 with Fig. 6, it is obvious that the velocity c for η = 0.1 is always larger than that for η = 1.0 provided that the wave number α and the number of “humps” are the same. It has been pointed out that for a very large value of η, there is no possibility of solitary wave solutions with nonzero c [10]. As the numerical computa-
tions become increasingly inaccurate and difficult as the number of “humps” becomes large or η is over 1.5, an upper limit for such η cannot be obtained at present. Furthermore, the distribution of the eigenvelocities within a family of stationary solitary wave with many “humps” is still unclear. The initial value problem of Eq. (2) was investigated in [10]. It was shown that for strongly dispersive cases (small η), the overall evolution appears to be nonstationary and somewhat chaotic, whereas for weakly dispersive case (large η), saturated saw-tooth
Fig. 6. Comparison of the numerical results and theoretical estimates for velocity c with η = 1.0. Solid lines: numerical results, (×): theoretical estimates for one-hump solutions, (•): theoretical estimates for two-hump solutions.
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
shock structures are observed as the final stage of temporal evolutions. In addition, transient appearance of multi-hump solitary waves similar to those given in Fig. 1 has been observed in the temporal evolution for small η (cf. Fig. 4.3 for η = 0.1 in [10]). Also some pulse structures are observed at transient stage for large η (Figs. 4.6 and 4.7 in [10]). This result again reminds us of the possibility that the stationary waves play fundamental roles in the initial value problem of the nonlinear evolution equations. Thus the stability of the stationary waves in dynamics becomes essential. In examining the stability, the system size L is an important parameter to be compared with the scale of solitary pulses. The scale of pulse is roughly estimated as l ≡ 2π/κc (κc = 1 as given in Section 2), because it is determined by a balance between instability and damping. The long wave instability of Eq. (2) becomes predominant for large system, because the growth rate takes a maximum at zero wave number. This is a crucial difference of Eq. (2) from the Benney (or the Kuramoto–Sivashinsky) equation for which the maximum growth rate is given at a finite wave number. Stability of one-peak and two-peak periodic stationary waves shown in Fig. 3 (for η = 0.1) is examined numerically by observing whether such solutions travel steadily or not in the initial value problem of Eq. (2). The one-peak case is found to be stable if the system size is roughly within 6.5 ≤ L ≤ 11.5 and two-peak case is stable if 14 ≤ L ≤ 18. Thus the stationary periodic waves can be an equilibrium solution of Eq. (2) so far as the system size is not very large. For large system size (L l), however, the stationary waves of typical scale (∼ l) will be unstable because the long wave instability becomes effective. This point has not been investigated yet, but numerical results [10] suggest the following conjectures. Although the stationary solitary wave solutions obtained in this paper cannot be equilibrium solutions of Eq. (2) for large system size case, they can arise transiently in the dynamical evolution [10]. For small η, some of the multi-hump solitary waves can be observed in chaotic-like evolutions. For large η, localized pulses arise in the transient stage, although the final state is
235
shock-like structures with the longest possible spatial scales as observed numerically in [10]. Recent investigations on the Benney equation have clarified the relation between the stationary waves and dynamical behaviours. Amongst them, the following point differing from the present case is worth noting here. Chaotic behaviour (or the multiplicity of solutions) in the weakly dispersive case (including Kuramoto–Sivashinsky limit) is reduced according to an increase of dispersive effect and leads to an equilibrium consisting of an array of one-peak solitary waves [6]. This is thought to be due to the mechanism that the maximum growth rate at a finite wave number governs to arrange solitary waves in order. In the present case, the longest scale is most unstable and no such finite scale exists. This difference in the growth rate will be a cause for somewhat contrasting dynamical behaviours, i.e., chaotic-like behaviour for strongly dispersive case and shock-like equilibrium structure for weakly dispersive case. However, the overall connection between the stationary wave solutions and the time-dependent solutions to Eq. (2) for large system size (large L) is still an unsettled problem. Appendix A. Numerical method A rational Chebyshev pseudo-spectral method, which was introduced originally by Boyd [12], is used to compute the stationary solitary waves. The basis functions we used are the mapped Chebyshev polynomials: TBn (x) = cos(n cot−1 (x/M)),
n = 0, 1, . . . , (A.1)
where √ M is the mapped parameter. With the map t = x/ M 2 + x 2 , our basis function, as defined in (A.1), is equal to Tn (t), where Tn (t) is the usual Chebyshev polynomial. We expand the solution of Eq. (3) as u(x) =
N X
ri TBi−1 (x),
i=1
where N is a given positive integer.
(A.2)
236
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
The corresponding collocation points xj (j = 1, 2, . . . , N) are chosen as xj = cot(π(2j − 1)/2N).
(A.3)
For the pseudo-spectral method, Eq. (3) has to be satisfied at the collocation points xj (j = 1, 2, . . . , N ). This result, together with the phase condition ux |x=0 = 0, yields N + 1 nonlinear equations for the N + 1 unknowns rj (j = 1, 2, . . . , N) and c, which can be solved by the Newton–Kantorovich iteration method. Practically, one sech-type function or the superposition of several sech-type functions is used as the initial profile to obtain one-hump, two-hump, etc., solitary wave solutions by adjusting the mapped parameter M simultaneously. Once the solitary wave solution is obtained, then we can improve the accuracy by virtue of a relatively large N . The typical value of N we used is 100–200. For the stationary periodic wave solution to Eq. (2), we expand it as a Fourier series u(x) =
N/2 X
(ak cos(kαx) + bk sin(kαx)),
(A.4)
k=1
where α = 2π/L, L is the periodicity length and N is a chosen positive number. Elimination of the constant from the Fourier expansion implies that condition (7) is satisfied automatically. The above expansion is substituted into the differential equation (3), which is set to be satisfied at each of the N collocation points: xj = π(j − 1)/N,
j = 1, 2, . . . , N.
(A.5)
The above N equations together with the phase condition being chosen as uN ≡ u0N here give N + 1 nonlinear equations for the N + 1 unknowns ak , bk (k = 1, 2, . . . , N/2) and c. We start from
a certain initial profile with the Fourier spectra ak0 , bk0 (k = 1, 2, . . . , N/2), which is taken from the equilibrium state in the temporal evolution of Eq. (2), and then conduct a Newton–Kantorovich iteration method until satisfactory convergence is reached. Since FFT is employed here, the typical value of N in our computations is 128 or 256. References [1] T. Kawahara, Formation of saturated solitons in a nonlinear dispersive system with instability and dissipation, Phys. Rev. Lett. 51 (1983) 381–383. [2] T. Kawahara, S. Toh, Pulse interaction in an unstable dissipative–dispersive nonlinear system, Phys. Fluids 31 (1988) 2103–2111. [3] D. Michelson, Steady solutions of the Kuramoto–Sivashinsky equation, Physica D 19 (1986) 89–111. [4] A.P. Hooper, R. Grimshaw, Travelling wave solutions of the Kuramoto–Shivashinsky equation, Wave Motion 10 (1988) 405–420. [5] Y.A. Demekhin, G.Yu. Tokarev, V.Ya. Shkadov, Hierarchy of bifurcations of space-periodic structures in a nonlinear model of active dissipative media, Physica D 52 (1991) 338– 361. [6] H.-C. Chang, E.A. Demekhin, D.I. Kopelevich, Laminarizing effects of dispersion in an active-dissipative nonlinear medium, Physica D 63 (1993) 299–320. [7] H.-C. Chang, Wave evolution on a falling film, Ann. Rev. Fluid Mech. 26 (1994) 103–136. [8] A. Nitzan, The critical behavior of nonequilibrium transitions in reacting diffusing systems, in: H. Haken (Ed.), Dynamics of Synergetic Systems, Springer, Berlin, 1980, pp. 119–131. [9] E. Ott, W.M. Manheimer, D.L. Book, J. Boris, Model equations for mode coupling saturation in unstable plasmas, Phys. Fluids 16 (1973) 855–862. [10] T. Kawahara, S. Toh, On some properties of solutions to a nonlinear equation including long-wavelength instability, in: A. Jeffrey (Ed.), Contributions to Nonlinear Wave Motion, Longman, New York, 1988, pp. 95–117. [11] S.B. Wineberg, J.F. McGrath, E.F. Gabl, L.R. Scott, C.E. Southwell, Implicit spectral methods for wave propagation problems, J. Comput. Phys. 97 (1991) 311–336. [12] J.P. Boyd, Spectral methods using rational basis functions on an infinite interval, J. Comput. Phys. 25 (1987) 113–143.
Stationary travelling-wave solutions of an unstable KdV–Burgers equation Bao-Feng Feng, Takuji Kawahara ∗ Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto 606-8501, Japan Received 9 March 1999; received in revised form 22 July 1999; accepted 24 August 1999 Communicated by Y. Kuramoto
Abstract Both “solitary” and “periodic” stationary travelling wave solutions are investigated numerically for an unstable Korteweg–de Vries–Burgers equation ut + uux + uxxx − η(u + uxx ) = 0 (η > 0). A family of stationary solitary wave solutions whose members are distinguished by the number of “humps” is found for a given η. Corresponding to each solitary wave thus found, a family of stationary periodic waves with the same number of “humps” exists under periodic condition and ends up in the infinite periodicity to the corresponding solitary wave. The numerical results are consistent with the theoretical estimates based on the conservation properties. ©2000 Elsevier Science B.V. All rights reserved. PACS: 03.40.Kf; 02.70.H Keywords: Unstable Korteweg–de Vries–Burgers equation; Solitary and periodic waves; Multi-hump solutions; Rational Chebyshev and Fourier pseudo-spectral method
1. Introduction Much attention has been paid to the roles of localized solutions in the initial value problem of the nonlinear dispersive equations involving both instability and dissipation. One prototype of such equations is the one-dimensional Benney equation: ut + uux + uxx + δuxxx + uxxxx = 0.
(1)
It has been found that in certain cases time evolutions of Eq. (1) are well approximated by an array of similar amplitude localized pulses (travelling wave solu∗ Corresponding author. E-mail address: [email protected] (T. Kawahara).
tions) [1,2]. The coexistence of instability and dissipation is responsible for determining both local and global behaviours. Locally the scale of pulse is determined by a balance between instability and dissipation terms of Eq. (1), which leads to a pulse with a definite amplitude, unlike the case of KdV solitons where arbitrary amplitudes are admitted. Globally the equilibrium states (either deterministic or chaotic) of the initial value problem are determined also by the balance between instability and dissipation. The balance of the rate of energy inflow by instability with that of energy outflow by dissipation determines an equilibrium state through energy transfer due to nonlinearity influenced by dispersion. For strongly dispersive cases, overall evolutions are deterministic but for relatively weakly
0167-2789/00/$ – see front matter ©2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 2 7 8 9 ( 9 9 ) 0 0 1 8 3 - 9
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
dispersive cases they are chaotic and the balance is in a statistical sense. On the other hand, recent researches have demonstrated the existence of complicated families of localized solutions for the Kuramoto–Sivashinsky equation [3–5] or the Benney equation [6,7]. These results pose a problem to be investigated, i.e., what relations exist between a variety of localized solutions and the overall evolutions (either deterministic or chaotic) of the initial value problem. This problem has not yet been explored in detail because of the complication and existence of too many families of solutions. In order to carry such an exploration a step further, it will be worthwhile to consider similar other nonlinear dispersive equations involving instability and dissipation. One of such equations is an unstable Korteweg–de Vries–Burgers equation: ut + uux + uxxx − η(u + uxx ) = 0,
(2)
where the parameter η is a positive constant, representing the importance of instability and dissipation relative to dispersion and nonlinearity. Eq. (2) can be regarded as a simple model for a so-called extrinsic instability, i.e., a longwave instability whose growth rate takes a maximum at zero wave number [8], or a physical system such as a current-driven ion acoustic instability in a collision-dominated plasma [9]. In this paper, both stationary solitary and periodic wave solutions to Eq. (2) are computed via a pseudo-spectral method. Existence of multiply humped solitary wave solutions to Eq. (2) in the range 0 < η ≤ 1.5 is discovered numerically. It is also found that every existing solitary wave solution corresponds to a family of periodic wave solutions. These periodic waves tend to the corresponding solitary wave provided that the periodicity length L → ∞. The remainder of this paper is organized as follows. In Section 2, a preliminary analysis is given to show some basic properties of Eq. (2). In Section 3, we describe the numerical results for the solitary waves and in Section 4 those for the periodic waves. Then, in Section 5, we give some discussions and comments. Finally, in Appendix A, the employed pseudo-spectral methods are introduced briefly.
229
2. Preliminary analysis Substitution of u ∝ exp(iκx + σ t) into Eq. (2) yields the linear dispersion relation σ = iκ 3 + η(1 − κ 2 ). Thereby the small-amplitude sinusoidal waves are unstable for long wavelengths given by |κ| < κc ≡ 1 and damped for short wavelengths |κ| > κc . This is similar to Eq. (1) but differs in that the maximum growth rate is given at zero wave number. We shall restrict our attention to stationary travelling-wave solutions of Eq. (2) over a spatially periodic domain x ∈ [0, L] or an unbounded domain x ∈ (−∞, ∞). These stationary waves are solutions of Eq. (2) which travel steadily at a constant speed c without changing their shapes. Therefore, it is appropriate to carry out a moving coordinate transformation ∂/∂t → ∂/∂t − c∂/∂x, which leads to the governing equation for stationary waves: −cux + uux + uxxx − η(u + uxx ) = 0,
(3)
where c is the wave speed. As preliminaries of the numerical study, we examine some basic properties of Eq. (3). Firstly, for the dynamical system corresponding to Eq. (3), only one fixed point O1 (0, 0, 0) is allowed in the three-dimensional phase space defined by coordinates (u, ux , uxx ). This is due to the existence of the unstable term proportional to u in Eq. (3). Therefore, only solitary waves, which correspond to homoclinic trajectories of O1 , are possible over the whole x domain −∞ < x < ∞. The boundary conditions are given by u(±∞) = 0,
ux (±∞) = 0,
uxx (±∞) = 0 (4a)
for solitary waves or by u(0) = u(L),
ux (0) = ux (L),
uxx (0) = uxx (L) (4b)
for periodic waves with periodicity length L. In order to investigate the asymptotic behaviour of u, we linearize Eq. (3) about O1 , which yields the characteristic equation for the characteristic exponent µ: µ3 − ηµ2 − cµ − η = 0.
(5)
230
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
We can distinguish the property of the roots to Eq. (5) by the criterion in terms of the following quantity: R = 4η4 + (18c − c2 )η2 − 4c3 .
(6)
A positive R corresponds to a real positive root and a pair of complex conjugate roots with negative real part, whereas a negative R corresponds to three real roots. Thereby, an oscillatory structure is expected in the tail of the solitary waves in the case of R > 0, whereas a monotonous structure in the case of R < 0. This property is confirmed in the numerical computations afterwards. Secondly, integration of Eq. (3) over x on the whole domain Ω (where Ω = (−∞, ∞) for a solitary wave or Ω = [0, L] for a periodic wave) with (4a) and (4b) leads to a necessary condition for the existence of a stationary wave to Eq. (2), i.e., hui = 0,
(7)
R∞ where h· · · i denotes −∞ · · · dx for a solitary wave or RL 0 · · · dx for a periodic wave with periodicity length L. Condition (7) is set to be satisfied throughout our computations of stationary wave solutions. Moreover, multiplying Eq. (3) by u and integrating it over x on the domain Ω, the relation hu2 i = hu2x i
Fig. 1. Numerical solitary wave solutions for Eq. (2) with η = 0.1: (a) c = 4.602758, (b) c = 4.281452, (c) c = 3.991830, (d) c = 3.728704, (e) c = 3.487565, (f) c = 3.235719.
(8)
can be derived after some manipulations with (4a) and (4b). Although no other conservation properties are possible because of the non-integrability of Eq. (2), the conservation properties (7) and (8) can provide us with a powerful tool for verifying the numerical results. As a matter of fact, the result of such a verification is very satisfactory as will be discussed in Section 5 in an attempt to obtain some qualitative explanations for the general features of the numerical results.
3. Numerical results for solitary waves The stationary solitary wave solutions of Eq. (2) were investigated by Kawahara and Toh [10] via a shooting method. They found that for a given η, one
particular single hump solitary wave solution exists. Moreover, by a perturbation analysis, such a solution is shown to approach the one soliton solution of the Korteweg–de Vries equation with a particular velocity c = 5.0 when η → 0. Here, we employ a rational Chebyshev pseudospectral method (see Appendix A) to compute the solitary wave solutions of Eq. (2). In addition to a single hump solitary wave found previously by the shooting method, we find a series of solitary wave solutions involving a multi-hump structure for a fixed value of η. Six solitary wave solutions (single peak to six peaks) for η = 0.1 are shown in Fig. 1 with eigenvelocities found numerically as c = 4.602758, 4.281452, 3.991830, 3.728704, 3.487565, 3.235719. For a relatively large η = 1.0, three numericals
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
231
4. Numerical results for periodic waves
Fig. 2. Numerical solitary wave solutions for Eq. (2) with η = 1.0: (a) c = 1.399297, (b) c = 0.243849, (c) c = −0.412683.
solutions are shown in Fig. 2 with c = 1.399297, 0.243849, −0.412683. It is worthwhile to note that similar to 1D Benney equation (1), multi-hump solitary wave solutions exist for Eq. (2). It can be seen that for η = 0.1, all of the six solutions have a negative plateau in their tail part, whereas for η = 1.0, an oscillatory structure arises. This result is consistent with the analysis of characteristic equation in Section 2. As the numerical computations become increasingly inaccurate and difficult as the number of “humps” increases or the value of η becomes large, only six solutions are identified for η = 0.1 and three for η = 1.0. We conjecture here that for every η ∈ (0, 1.5], there is a family of solitary wave solutions with a countable number of members. Members in each family are distinguished by the number of “humps”. Because of the numerical difficulty, it is not sure whether such solutions exist or not if η is beyond this range. For a given η, the eigenvelocity for the first member (single hump solitary wave) of the family is the largest, and it decreases as the number of “humps” increases. For the members with the same number of “humps”, the eigenvelocity decreases as η increases. Moreover, in contrast to the KdV solitons, the negative eigenvelocity is possible.
The stationary periodic wave solutions to Eq. (2) are computed by means of a Fourier pseudo-spectral method, which will be introduced briefly in Appendix A together with the rational Chebyshev pseudo-spectral method. We begin to seek a periodic solution with a relatively large wave number α. Here, α ≡ 2π/L denotes the minimum wave number allowed for the periodic system with periodicity length L. The initial profile is taken from the equilibrium state in the temporal evolution of Eq. (2) when it is integrated over the same periodicity length via an implicit spectral method developed by Wineberg et al. [11]. The Newton–Kantorovich iteration method is used in every time step. Then we initiate the continuation scheme towards the small wave number. Thus a family of periodic waves is constructed for a given η. One part of the numerical results with different wave number α is shown in Figs. 3 and 4 with η = 0.1, 1.0, respectively. It is found numerically that similar to the case of solitary waves, multi-hump periodic wave solutions exist. We denote n-hump periodic wave solution by P (n) . Moreover, we found that P (1) starts from one near-sinusoidal periodic wave with α = α ∗ ≈ 0.96664 and converges to the corresponding solitary wave solution as the wave number α → 0. For a fixed value of η, we conjecture that families of periodic wave solutions P (1) , P (2) , . . . , P (n) , . . . exist, which originate from one near-sinusoidal periodic wave with wave number α = α ∗ , α ∗ /2, . . . , α ∗ /n, . . . , respectively. Each family of periodic waves gradually tend to the corresponding solitary wave as the wave number α → 0. It is also interesting to note that the way of convergence is associated with the tail structure of the solitary waves. For the solitary wave with monotonous tail, the set of the eigenvelocities for the periodic waves converges to the velocity of the corresponding solitary wave monotonously from the lower side, whereas for the solitary wave with an oscillatory tail, it also converges to the velocity of the corresponding solitary wave in an oscillatory way.
232
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
Fig. 3. One-hump and two-hump periodic wave solutions for Eq. (2) with η = 0.1. One-hump: (a) L = 20, c = 3.020606, (b) L = 40, c = 3.966167, (c) L = 60, c = 4.201104. Two-hump: (d) L = 20, c = 1.189170, (e) L = 40, c = 3.029772, (f) L = 60, c = 3.572864. (Note the difference in horizontal scale.)
5. Discussion and some comments The stationary solitary wave solutions and the periodic wave solutions to an unstable Korteweg–de Vries–Burgers equation are investigated, respectively, by means of a rational Chebyshev pseudo-spectral method and a Fourier pseudo-spectral method. Existence of multiple stationary wave solutions for both solitary and periodic waves is discovered numerically, and the relationship between them is also found. Within the numerical computations carried out so far, the following general features have been found: 1. All the velocities for the solitary waves or periodic waves are less than 5.0 and negative velocity is possible.
Fig. 4. One-hump and two-hump periodic wave solutions for Eq. (2) with η = 1.0. One hump: (a) L = 15, c = 1.454888, (b) L = 20, c = 1.395870, (c) L = 25, c = 1.399301. Two hump: (d) L = 25, c = 0.297846, (e) L = 27.5, c = 0.222642, (f) L = 40, c = 0.244304.
2. For a fixed value of η, there exists one family of solitary wave solutions to Eq. (2), in which it is conjectured that a countable number of members should exist. The members are distinguished with each other by the number of “humps”. The velocity for the first member is the greatest, then it decreases as the number of “humps” increases. 3. The velocity of the first member in the family of solitary waves approaches to the value of 5.0 as η → 0. For the member with the same number of “humps”, the corresponding velocity decreases as η increases. 4. For each solitary wave, one family of periodic waves with the same number of “humps” exists. These periodic waves gradually tend to the corre-
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
sponding solitary wave as the wave number α → 0. However, the way of the convergence of the eigenvelocities is different, and depends on the tail structure of the solitary waves. As Eq. (2) is nonintegrable and an analytical solution of stationary wave is not available, it has to be investigated numerically. However, some qualitative and quantitative explanations for the numerical results may be given by reviewing the basic properties of Eq. (2) such as (7) and (8). For example, it has been well explained in [10] by a perturbation analysis why the velocity c is close to the value of 5.0 when η becomes very small. Here, we derive an expression to estimate the velocity c of a stationary wave. Integration of Eq. (3) from the starting point of Ω to an arbitrary position ξ ∈ Ω leads to Z ξ 1 2 u dx + ux (ξ ) −cu(ξ ) + 2 u (ξ ) + uxx (ξ ) − η −∞
=0
(9a)
for a solitary wave or Z −cu(ξ ) +
1 2 2 u (ξ ) + uxx (ξ ) − η
= cu(0) −
ξ
u dx + ux (ξ )
0 1 2 2 u (0) − uxx (0) + ηux (0)
(9b)
Rξ
for a periodic wave. We set v(ξ ) = −∞ u dx for a Rξ solitary wave or v(ξ ) = 0 u dx for a periodic wave, and then obviously, we have dv/dξ = u. Then multiplying Eq. (9a) or (9b) by u(ξ ) and integrating it over ξ on the whole domain Ω, we have −chu2 i + 21 hu3 i + huuxx i − η(huvi + huux i) = 0, (10) here, the condition (7) is used. Thereafter, it can be shown readily that huux i = 0,
huvi = 0,
huuxx i = −hu2x i.
(11)
Therefore, using (8), we finally have c=
hu3 i − 1, 2hu2 i
(12)
which gives an estimate of the eigenvelocity in terms of the integrals of u.
233
We can divide the integral hui into two parts, i.e., the integral I1 over the domain where u is positive and the integral I2 over the domain where u is negative. The same procedure can be done for the integrals hu3 i and hu2 i too. Apparently, from (7), we have I1 > 0 and I1 = −I2 . It should be noted that the integrals of hu3 i in each domain keep the same sign as hui, but those for hu2 i are always positive. Bearing this fact in mind, we consider the stationary waves with the same η. As the number of “humps” increases, the negative level of the wave tail has to be lowered in order to maintain the zero mean value of u. Because of the cancellation between positive and negative contributions for hu3 i, the rate of increase in hu3 i is expected to be smaller than that in hu2 i. Thus it is reasonable that the quantity hu3 i/hu2 i will decrease, i.e., the velocity c will decrease. From the numerical results and the preliminary analysis, we have the fact that for a small value of η, there is a long and negative plateau in the tail part, whereas for a large value of η, an oscillatory structure arises in the tail part. Thereby, for the stationary waves with the same number of “humps”, the “basin” following the last “hump” is relatively deep and narrow for a large value of η, while being flat and broad for a small value of η. Henceforth, it is also reasonable to think that the quantity hu3 i/hu2 i will decrease as the value of η increases. One more point is that when the quantity hu3 i/hu2 i is less than 2, the velocity will be negative. The numerical values of velocity c versus wave number α for typical values of η = 0.1, 1.0 are shown in Figs. 5 and 6. The theoretical estimates of c from Eq. (12) are also computed and compared. Here, P (1) and P (2) denote the branches for the one-hump and two-hump stationary waves, respectively. Numerical integrals hu3 i and hu2 i are computed by Simpson’s rule. The numerical results for the rates of increase in hu2 i and hu3 i are consistent with the above-mentioned theoretical expectation. Figs. 5 and 6 show that the numerical values of c coincide with the theoretical estimates very well. In other words, this is the verification of our numerical results. Moreover, for each branch P (i) (i = 1, 2), the family of stationary periodic waves starts from
234
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
Fig. 5. Comparison of the numerical results and theoretical estimates for velocity c with η = 0.1. Solid lines: numerical results, (×): theoretical estimates for one-hump solutions, (•): theoretical estimates for two-hump solutions.
certain point in the α–c plane and ends up with the corresponding solitary wave. For the same value of η and wave number α, the velocity c of the two-hump stationary wave is smaller than that of one-hump one. When comparing Fig. 5 with Fig. 6, it is obvious that the velocity c for η = 0.1 is always larger than that for η = 1.0 provided that the wave number α and the number of “humps” are the same. It has been pointed out that for a very large value of η, there is no possibility of solitary wave solutions with nonzero c [10]. As the numerical computa-
tions become increasingly inaccurate and difficult as the number of “humps” becomes large or η is over 1.5, an upper limit for such η cannot be obtained at present. Furthermore, the distribution of the eigenvelocities within a family of stationary solitary wave with many “humps” is still unclear. The initial value problem of Eq. (2) was investigated in [10]. It was shown that for strongly dispersive cases (small η), the overall evolution appears to be nonstationary and somewhat chaotic, whereas for weakly dispersive case (large η), saturated saw-tooth
Fig. 6. Comparison of the numerical results and theoretical estimates for velocity c with η = 1.0. Solid lines: numerical results, (×): theoretical estimates for one-hump solutions, (•): theoretical estimates for two-hump solutions.
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
shock structures are observed as the final stage of temporal evolutions. In addition, transient appearance of multi-hump solitary waves similar to those given in Fig. 1 has been observed in the temporal evolution for small η (cf. Fig. 4.3 for η = 0.1 in [10]). Also some pulse structures are observed at transient stage for large η (Figs. 4.6 and 4.7 in [10]). This result again reminds us of the possibility that the stationary waves play fundamental roles in the initial value problem of the nonlinear evolution equations. Thus the stability of the stationary waves in dynamics becomes essential. In examining the stability, the system size L is an important parameter to be compared with the scale of solitary pulses. The scale of pulse is roughly estimated as l ≡ 2π/κc (κc = 1 as given in Section 2), because it is determined by a balance between instability and damping. The long wave instability of Eq. (2) becomes predominant for large system, because the growth rate takes a maximum at zero wave number. This is a crucial difference of Eq. (2) from the Benney (or the Kuramoto–Sivashinsky) equation for which the maximum growth rate is given at a finite wave number. Stability of one-peak and two-peak periodic stationary waves shown in Fig. 3 (for η = 0.1) is examined numerically by observing whether such solutions travel steadily or not in the initial value problem of Eq. (2). The one-peak case is found to be stable if the system size is roughly within 6.5 ≤ L ≤ 11.5 and two-peak case is stable if 14 ≤ L ≤ 18. Thus the stationary periodic waves can be an equilibrium solution of Eq. (2) so far as the system size is not very large. For large system size (L l), however, the stationary waves of typical scale (∼ l) will be unstable because the long wave instability becomes effective. This point has not been investigated yet, but numerical results [10] suggest the following conjectures. Although the stationary solitary wave solutions obtained in this paper cannot be equilibrium solutions of Eq. (2) for large system size case, they can arise transiently in the dynamical evolution [10]. For small η, some of the multi-hump solitary waves can be observed in chaotic-like evolutions. For large η, localized pulses arise in the transient stage, although the final state is
235
shock-like structures with the longest possible spatial scales as observed numerically in [10]. Recent investigations on the Benney equation have clarified the relation between the stationary waves and dynamical behaviours. Amongst them, the following point differing from the present case is worth noting here. Chaotic behaviour (or the multiplicity of solutions) in the weakly dispersive case (including Kuramoto–Sivashinsky limit) is reduced according to an increase of dispersive effect and leads to an equilibrium consisting of an array of one-peak solitary waves [6]. This is thought to be due to the mechanism that the maximum growth rate at a finite wave number governs to arrange solitary waves in order. In the present case, the longest scale is most unstable and no such finite scale exists. This difference in the growth rate will be a cause for somewhat contrasting dynamical behaviours, i.e., chaotic-like behaviour for strongly dispersive case and shock-like equilibrium structure for weakly dispersive case. However, the overall connection between the stationary wave solutions and the time-dependent solutions to Eq. (2) for large system size (large L) is still an unsettled problem. Appendix A. Numerical method A rational Chebyshev pseudo-spectral method, which was introduced originally by Boyd [12], is used to compute the stationary solitary waves. The basis functions we used are the mapped Chebyshev polynomials: TBn (x) = cos(n cot−1 (x/M)),
n = 0, 1, . . . , (A.1)
where √ M is the mapped parameter. With the map t = x/ M 2 + x 2 , our basis function, as defined in (A.1), is equal to Tn (t), where Tn (t) is the usual Chebyshev polynomial. We expand the solution of Eq. (3) as u(x) =
N X
ri TBi−1 (x),
i=1
where N is a given positive integer.
(A.2)
236
B.-F. Feng, T. Kawahara / Physica D 137 (2000) 228–236
The corresponding collocation points xj (j = 1, 2, . . . , N) are chosen as xj = cot(π(2j − 1)/2N).
(A.3)
For the pseudo-spectral method, Eq. (3) has to be satisfied at the collocation points xj (j = 1, 2, . . . , N ). This result, together with the phase condition ux |x=0 = 0, yields N + 1 nonlinear equations for the N + 1 unknowns rj (j = 1, 2, . . . , N) and c, which can be solved by the Newton–Kantorovich iteration method. Practically, one sech-type function or the superposition of several sech-type functions is used as the initial profile to obtain one-hump, two-hump, etc., solitary wave solutions by adjusting the mapped parameter M simultaneously. Once the solitary wave solution is obtained, then we can improve the accuracy by virtue of a relatively large N . The typical value of N we used is 100–200. For the stationary periodic wave solution to Eq. (2), we expand it as a Fourier series u(x) =
N/2 X
(ak cos(kαx) + bk sin(kαx)),
(A.4)
k=1
where α = 2π/L, L is the periodicity length and N is a chosen positive number. Elimination of the constant from the Fourier expansion implies that condition (7) is satisfied automatically. The above expansion is substituted into the differential equation (3), which is set to be satisfied at each of the N collocation points: xj = π(j − 1)/N,
j = 1, 2, . . . , N.
(A.5)
The above N equations together with the phase condition being chosen as uN ≡ u0N here give N + 1 nonlinear equations for the N + 1 unknowns ak , bk (k = 1, 2, . . . , N/2) and c. We start from
a certain initial profile with the Fourier spectra ak0 , bk0 (k = 1, 2, . . . , N/2), which is taken from the equilibrium state in the temporal evolution of Eq. (2), and then conduct a Newton–Kantorovich iteration method until satisfactory convergence is reached. Since FFT is employed here, the typical value of N in our computations is 128 or 256. References [1] T. Kawahara, Formation of saturated solitons in a nonlinear dispersive system with instability and dissipation, Phys. Rev. Lett. 51 (1983) 381–383. [2] T. Kawahara, S. Toh, Pulse interaction in an unstable dissipative–dispersive nonlinear system, Phys. Fluids 31 (1988) 2103–2111. [3] D. Michelson, Steady solutions of the Kuramoto–Sivashinsky equation, Physica D 19 (1986) 89–111. [4] A.P. Hooper, R. Grimshaw, Travelling wave solutions of the Kuramoto–Shivashinsky equation, Wave Motion 10 (1988) 405–420. [5] Y.A. Demekhin, G.Yu. Tokarev, V.Ya. Shkadov, Hierarchy of bifurcations of space-periodic structures in a nonlinear model of active dissipative media, Physica D 52 (1991) 338– 361. [6] H.-C. Chang, E.A. Demekhin, D.I. Kopelevich, Laminarizing effects of dispersion in an active-dissipative nonlinear medium, Physica D 63 (1993) 299–320. [7] H.-C. Chang, Wave evolution on a falling film, Ann. Rev. Fluid Mech. 26 (1994) 103–136. [8] A. Nitzan, The critical behavior of nonequilibrium transitions in reacting diffusing systems, in: H. Haken (Ed.), Dynamics of Synergetic Systems, Springer, Berlin, 1980, pp. 119–131. [9] E. Ott, W.M. Manheimer, D.L. Book, J. Boris, Model equations for mode coupling saturation in unstable plasmas, Phys. Fluids 16 (1973) 855–862. [10] T. Kawahara, S. Toh, On some properties of solutions to a nonlinear equation including long-wavelength instability, in: A. Jeffrey (Ed.), Contributions to Nonlinear Wave Motion, Longman, New York, 1988, pp. 95–117. [11] S.B. Wineberg, J.F. McGrath, E.F. Gabl, L.R. Scott, C.E. Southwell, Implicit spectral methods for wave propagation problems, J. Comput. Phys. 97 (1991) 311–336. [12] J.P. Boyd, Spectral methods using rational basis functions on an infinite interval, J. Comput. Phys. 25 (1987) 113–143.