Weak periodic solutions and numerical case studies of the Fornberg-Whitham equation

Spatially periodic solutions of the Fornberg-Whitham equation are studied to illustrate the mechanism of wave breaking and the formation of shocks for a large class of initial data. We show that these solutions can be considered to be weak solutions satisfying the entropy condition. By numerical experiments, we show that the breaking waves become shock-wave type in the time evolution.


Introduction
We denote by T = R/Z the one-dimensional torus group and identify functions on T with 1periodic functions on R. The Fornberg-Whitham equation for the wave height u : T×[0, ∞[→ R as a function of a spatial variable x ∈ T and time t ≥ 0 reads u t + uu x + 1 − ∂ 2 x −1 u x = 0 (x ∈ T, t > 0) (1) and is supplied with the initial condition u(x, 0) = u 0 (x) (x ∈ T).
Well-posedness results, local in time, with strong solutions in Sobolev and Besov spaces have been obtained for both cases, periodic and non-periodic, in [12,13]. We cannot always expect globally defined strong solutions. This was predicted with a sketch of proof by [23] and later proved rigorously in [3,10,14]. These results proved the existence of wave-breaking, i.e., that the solution u remains bounded but u x becomes singular, if the initial data displays sufficient asymmetry in its slope. However two questions seem to have remained unanswered: (i) What can we say about the solution after a singularity has emerged? (ii) What is the nature of the singularity?
Our first goal is to prove existence of globally defined weak solutions which extend beyond the time of wave-breaking and singularity formation. The second goal of the present paper is to show numerically that the singularity developing in such a solution looks very similar to a shock-wave solution of the inviscid Burgers equation. We emphasize that for the case of the real line in place of the torus, similar investigations have been carried out in [7], where the Fornberg-Whitham equation was formulated as and named Burgers-Poisson system (hence we have been unaware of that paper until almost completion of our current paper).

Weak entropy solutions
A well-known machinery exists for weak solution concepts for nonlinear partial differential equations in the form of hyperbolic conservation laws, but our Equation (1) involves a nonlocal term and the question is whether or not this can be harmful to global existence. As we will prove, this is not a big hurdle.
In the literature on nonlinear conservation laws, one of the fundamental references is Kružkov's classic paper [20], where he considered equations of the form u t + (ϕ(u)) x + ψ(u) = 0 (x ∈ R, t > 0).
He proved, under mild assumptions on the functions ϕ, ψ and on the initial data u 0 , that a weak solution exists globally in time and that it is unique in an appropriate class of functions. His setting is different from ours in two respects: First, we have x belonging to the one-dimensional torus, while Kružkov described the case with x on the real line. Second, in Kružkov's paper φ and ψ are ordinary functions, while our equation involves a non-local dependence on u in ψ. The first difference causes no problem. The second issue indeed forces us to alter some of the technicalities along the way, but a close examination of the very lucid presentation of the proofs in [20] reveals that the essence of most arguments can be applied to our equation with only slight modifications, which we will indicate in the sequel. Note first that we may write where the convolution is in the x-variable only and with kernel function given by K : T → R is given by K(x) = (e x + e 1−x )/(2(e − 1)) = √ e e−1 cosh(x − 1 2 ) for 0 ≤ x ≤ 1 (see, e.g., [14,Section 3]). Note that K is continuous but is not C 1 on T. Note also that the derivative K ′ is not continuous but is bounded.
We will relate to Kružkov's notation in [20] as closely as possible, but will switch the function arguments from (t, x) to (x, t) and occasionally write u(t) to denote the function x → u(x, t) with a frozen t. Compared with the main terms of the equation as labeled by Kružkov we set t))(x). The term involving ϕ conforms perfectly with the specifications from [20] and requires no extra consideration at all. For the linear, but non-local, term given by the operator ψ we will describe below suitable adaptations in the proof of uniqueness and make note of an alternative a priori estimate in the proof of existence. Overall in our situation, spatial periodicity, i.e., the compactness of T, simplifies several estimates along the way in following the various proofs of key results in [20]. Moreover, we have L ∞ (T) ֒→ L 1 (T).
We use the convolution representation K ′ * u in place of Kružkov's term ψ(u) in the following definition of the weak solution concept (where we also incorporate the initial condition into the basic inequality).
is called a weak entropy solution of (1)-(2) if the following (3) holds true: for any λ ∈ R and for any nonnegative C 1 -function of compact support in T × R.
As is well-known, upon putting λ = ± sup |u(x, t)| we may deduce that a weak entropy solution is also a weak (distributional) solution in the sense that in the integro-differential equation (1) the term uu x may be interpreted as ∂ x (u 2 )/2, since u ∈ L ∞ . Thus, we have

Uniqueness
We first show that the solution of (3) is unique. Existence will be proved later.
To show uniqueness of weak entropy solution we need an adaptation of [20, Theorem 1]noting that for large R we simply have K = [0, T 0 ]×T and S τ = T in that statement-and of its proof to our situation with the non-local linear term ψu = K ′ * u. This requires an appropriate replacement of the constant γ in [20,Equation (3.1)] and an alternative argument in the course of the proof of [20,Theorem 1] for the transition from [20,Equation (3.12)] to the inequalities at the bottom of page 228 and top of page 229 in [20], since in our case we cannot directly have a pointwise estimate calling on the mean value theorem for a classical differentiable function ψ. Instead, with two weak solutions u and v with initial values u 0 and v 0 , respectively, we obviously have which implies the following replacement of the next to last inequality on page 228 in [20] (with mollifier δ h and cut-off χ ε as chosen by Kružkov) Therefore, we replace the inequality stretching in [20] from the bottom of page 228 to the top of page 229 by Then, sending ρ → 0, we arrive at (Here, the L 1 -continuity of the weak solution is used.) Now Gronwall's inequality implies the following uniqueness result.

Existence
We will establish the global existence here under the condition that u 0 ∈ L ∞ (T).
The key idea in [20] is to consider a parabolic regularization of (1) in the form with a small parameter ε > 0. Our aim is to show that, at least for a subsequence of ε → 0, the corresponding solutions (all with the same initial value u 0 ) converge strongly in L 1 and are uniformly bounded in L ∞ . For any ε > 0, we will show that a strong solution u to (4) with initial value exists uniquely and globally in time. Due to the nonlocal term, we cannot directly call on the same references that Kružkov uses in [20] on page 231, but the desired result can be shown by a careful iteration of a standard contraction argument, which we describe in Appendix. The solution u of (4) and (5) belongs to C([0, T ]; L 2 (T)) ∩ C(]0, T ]; H 1 (T)) and, in fact, is smooth Let T > 0 be arbitrary and consider the unique solution of (4) and (5). For every ε > 0 we denote now the solution of (4)-(5) by u ε . We now have to find an alternative way to obtain Kružkov's basic estimate (4.6) in [20], which he got from a direct application of the maximum principle. To arrive at our analogue of the basic estimate in (8) below, we proceed as follows: Multiplying (4) by u ε and integrating over the spatial variable x yields , and integrating by parts on the right-hand side, yields (thanks to periodicity) In particular, we have for every t ∈ [0, T ] the a priori estimate which is independent of ε. From (7) and since ( for some constant c ′ independent of ε and of t. We now interpret the original equation in the form and apply the theorem on page 230 in [18] to obtain as a coefficient to the first order derivative and also note that in our case the zero order coefficient vanishes. Although the text in [18] assumes more regularity of the coefficients and data, continuity is sufficient.) We then let s → 0 to obtain This inequality is our analogue of Kružkov's basic estimate (4.6) in [20] and we may now return to follow along his lines again more closely (note that we are closest to what he classifies as 'Case A' in his paper on page 230 in next to the last paragraph) to establish a uniform modulus of L 1 -continuity. In fact, [20, Equation (4.7)] for w(x, t) := u ε (x + z, t) − u ε (x, t) holds in our case with e i = 0 and replacing the term cw by K ′ * w, preserving the Lipschitz continuity properties noted on top of page 232 in [20] (in our case even globally on the compact torus).
Lemma 5 in [20] is applicable in essentially the same way in establishing the key modulus of continuity estimate [20,Equation (4.15)] (with an appropriate proof variant taking into account the convolution term and showing (4.13) directly from (1.3) there). Thus, we have at least sketched how everything is in place and sufficient to apply Kružkov's method in proving existence from compactness in L 1 (via boundedness and equicontinuity of the L 1 -norm, [1,Theorem 4.26]) and uniform L ∞ -bounds. In combination with the above uniqueness result, we have the following statement.

Numerical experiments by finite differences
We now study numerical solutions and will observe that many of these exhibit singularities of shock-wave type. The solutions have been computed numerically under periodic boundary conditions and employing Godunov's finite difference method. We employ the following finite difference method with the uniform grid sizes h = ∆x, τ = ∆t. The nonlinear term is discretized by where the numerical flux is defined by Here u n k is the approximation for u(kh, nτ ), and f (u) = u 2 /2. This is the Godunov method as explicated in [4].
The nonlocal term is discretized by the following central difference scheme: We set N = 1000, h = 1/N, and τ = 0.4h/q, where q is the typical size of the initial data. We first test the following two initial data: The corresponding profiles of u from data1 are shown in Figure 1. As the profile moves to the right, the formation of a shock is clearly visible. After the emergence of the shock, the jump height at the discontinuity decreases as t increases, namely, if ξ(t) denotes the position of the discontinuity at time t, the difference of the one-sided limits u(ξ(t) − 0) − u(ξ(t) + 0) is a decreasing function of t. The solution with data2 is shown in Figure 2. In this case, even multiple shocks are formed and merging of some of the shocks over time can be observed. These and other experiments suggest the development of wave-breaking singularities into shock discontinuities. In fact, as far as we have computed, only shock-type singularities in wave solutions have been found. Moreover, the computations also suggest that global solutions exist, if the initial data are small-the recent result in [16, Theorem 1.5] on Fornberg-Whitham-type equations with nonlinear term ∂ x (u p /p), p ≥ 5, seems to support our observation.

Remarks on shock conditions and asymptotic properties
In our experiments all initial functions of the following form u 0 = 1≤k≤3 [a k cos(2kπx) + b k sin(2kπx)] which are periodic in x, but not necessarily symmetric, produce a shock spontaneously, if its L ∞ -norm is large. If the initial data is small, the numerical solution exists for quite a long time. For instance, if u(x, 0) = q cos(2πx), a shock wave was observed for q > 0.015. For q < 0.01 the solution seems to exist for ever. For 0.01 < q < 0.015, our experiments are indecisive. It may exists for ever, or it may have a shock after quite a long, numerically undetectable time has passed. Note that the guaranteed time of existence according to [12,Theorem 1] is inverse proportional to u 0 H 2 . The paper mentioned above, [16] discussing a whole class of Fornberg-Whitham-type equations, contains also detailed information on the existence time.
The propagation speed of the shock is the same as in the case of the inviscid Burgers equation Indeed, for any w in L 2 , the function (1 − ∂ 2 x ) −1 w x belongs to H 2 (T), hence is continuous, and therefore does not contribute to the Rankine-Hugoniot jump relation Here, c is the propagation speed of the shock and u + and u − denote the limit values of u from the right and the left at the discontinuity, respectively. Although the formation of a shock is quite similar in both equations, the asymptotic behavior may be different. In the case of the Burgers equation, the jump discontinuity vanishes in the limit and the solution converges to a constant function as t → ∞ (for a proof see [21] or [5, Theorems 6.4.9 and 11.12.5]). On the other hand, the Fornberg-Whitham equation possesses many traveling wave solutions that obviously do not decay and it is not yet conclusively shown whether or not for the discontinuous solutions displayed above the jump heights definitely decay to zero. We are not aware of a decisive result whether discontinuous periodic traveling waves exist, although we discuss below a negative result in the case of a single shock. (Existence of non-periodic discontinuous wave solutions has been shown in [7,15].) Another feature of smooth solutions u to the inviscid Burgers equation is that it preserves the values of maxima and minima of u(t, ·) as t progresses. The Fornberg-Whitham equation does not keep those values constant due to the presence of the nonlocal term. However, our experiments suggest that, although max x u(x, t) and min x u(x, t) is not a constant function of t, they are nearly constant, or, at least they seem to stay bounded as t → ∞.
The number of peaks of any solution to the inviscid Burgers equation is non-increasing, while we observe that in case of our equation they may increase (and decrease) as t varies.
Although several of the minor discrepancies exist as indicated, the solutions of (1) have some similarity with solutions of the inviscid Burgers equation in any time interval until and little after the formation of a shock discontinuity.

Traveling waves
In this subsection we first investigate the stability of traveling wave solutions and then the non-existence of such with a single shock.
Consider a continuously differentiable traveling wave u(x, t) = U(x − ct) with speed c > 0. Inserting this into (1) and integrating once we deduce that the profile function U satisfies the equation (The constant of integration may be set to zero without losing generality, if U and c are suitably normalized.) If we define V by where we choose the negative sign for the root, since only this connects U(x) = 0 with V (x) = 0. We obtain the differential equation which possesses V ≡ 0 as trivial solution for all c. Linearization of (11) at V ≡ 0 yields In terms of the Fourier coefficients (V (m)) m∈Z for the 1-periodic function V , this means 1 + (2πm) 2 − 1 c ·V (m) = 0 for every m ∈ Z. We obtain for every n ∈ N a nontrivial solution to the linearized equation in the form V n (x) := cos(2πnx), if the speed attains the bifurcation value c = c n := 1 1 + (2πn) 2 , otherwise we are left with V ≡ 0 as the only solution.
In the case of n = 1, nontrivial bifurcating solutions to the nonlinear differential equation exist in the interval c 1 ≈ 0.0247 < c < 0.02695. As c tends to the upper limit, min(c 2 + 2V ) tends to zero, and the profile U = c − √ c 2 + 2V tends to a function with a corner singularity. This Lipshitz continuous traveling wave, called peakon, has already been known to exist for a long time (cf. [8,25]) and occurs also as part of the analysis in [26] . Here we have computed these traveling waves in the context above, i.e., as the solutions of the boundary value problem (11) and illustrate some of them in Figure 3. In order to investigate stability, we picked the solution corresponding to c = 0.0255, input it as initial data, and computed the time dependent solution shown in Figure 4, which illustrates that the fixed wave profile travels at constant speed. In the time interval used by us, 0 < t < 30, the effect of numerical viscosity is invisible, but for long time intervals, it is expected to influence the numerical solution. We consider a perturbation of the initial data in the form U(x) + δd(x), where δ is 5% of the amplitude of U and d(x) is given as follows: Note that the latter disturbance is asymmetric. It turns out that in the time interval 0 ≤ t ≤ 300 the solution stays in a small neighborhood of the original solution U and only a small oscillation could be observed as is shown in Figure 5, where the points (u n 300 , u n 600 ) with n corresponding to 0 ≤ t ≤ 300 were plotted. If the initial perturbation is null, these describe a closed curve; once the initial disturbance is added, the corresponding curve departs from the closed orbit, but remains in a certain neighborhood. These and similar experiments may be interpreted as support for the conjecture of stability of the traveling wave with speed c = 0.0255. More experiments carried out for the case c = 0.0265 gave similar results.  (12). The points (u n 300 , u n 600 ) with n corresponding to 0 ≤ t ≤ 300 are plotted.

Nonexistence of traveling waves with a single shock
In view of the solutions shown earlier that created shocks after wave breaking, with jump height decreasing as time progresses, it is natural to ask whether periodic traveling wave solutions with jump discontinuities exist. For the non-periodic case it has been shown that discontinuous traveling waves with single jumps exist (see [7,15]). Contrary to this, there is no periodic traveling wave with a single jump, as we will argue in the following. Suppose U is the profile function for a discontinuous traveling wave solution that is piecewise C 1 on the torus T, i.e., C 1 except for a single point x 0 ∈ T where U as well as U ′ possess one-sided limits. Since the Fornberg-Whitham equation is invariant under translations in the x-variable, we may restrict to the case x 0 = 0, thus we may think of the profile function as a C 1 -function U : [0, 1] → R with U(0) = U(1). The Rankine-Hugoniot shock condition requires We come back to Equation (10) for the profile function U, but this time keep track of the constant of integration, for later convenience in the form which also shows that Z = (Y 2 /2) ′ = −K ′ * Y is continuous as a function on the torus and piecewise C 1 , i.e. can be thought of as C 1 -function Z : [0, 1] → R with We may therefore apply (1 − ∂ 2 x ) and arrive at We collect the equations for Y and Z in the following first-order system for the C 1 -functions Y and Z on [0, 1]. The requirements on U, including the Rankine-Hugoniot condition, now read while for Z we have the periodic continuity condition (13). We split the further analysis into two cases depending on the value of β: Case 2β + 1 ≤ 0: Equation (15) implies that Z ′ ≥ 0. By (13), this leaves only the option that Z ′ = 0, hence 0 ≤ (Y + 1) 2 = 2β + 1 ≤ 0, so that Y would have to be constant (equal to −1), which is a contradiction to (16).
Since 2β + 1 > 0, the vector field defining the right-hand sides of (14) and (15) has equilibirum points in (−1 ± √ 2β + 1, 0) and the qualitative analysis in [7, Section 3] may be applied to show that no trajectory satisfying all the above specifications can exist (the quantities u, E, d used there correspond here to Y , Z, β, respectively).
In summary, we have shown that there is no periodic traveling wave with a single shock. The question remains whether there exist traveling wave solutions with two or more shock discontinuities (and being piecwise C 1 ). Reviewing the above line of arguments, the reasoning in the first case, 2β + 1 ≤ 0, seems to essentially stay valid (with monotonicity arguments on subintervals instead), whereas in the second case, 2β + 1 > 0, the crucial consequence that Y (and hence Z) has to vanish somewhere is lost, if Y is allowed to have an additional jump discontinuity. We have to leave this issue open for potential future analysis and note that, even in that case, Z still has to be a continuous function on all of T thanks to the relation

Concluding remarks
The Godunov method is of first order and in further numerical studies one might want to employ a method of higher order such as the ENO (Essentially Non-Oscillatory) scheme for better accuracy. Furthermore, the computation of (1 − ∂ 2 x ) −2 ∂ x u from u in our approach may contain a significant truncation error. Therefore, there admittedly is a lot of room for improvement in the numerical experiments. Nevertheless, we do expect our numerical solutions to correctly show several main qualitative features. For instance, the almost generic emergence of shocks in the wave solution u appears to be undoubted and our experiments strongly indicate that many of the traveling wave solutions are stable. On the other hand, the large time behavior of solutions in general cannot be assessed quantitatively by our method and we clearly lack theoretical insight.
In summary, we are left with (at least) the following questions which we were unable to answer in terms of a rigorous analysis so far: (i) Existence of periodic traveling waves with jump discontinuities (i.e, non-decreasing shocks), although the case of a single jump could be ruled out.
(ii) Boundedness of the spatial minimum and maximum of u(t) as a function of t and independent of the existence time T.
(iii) Wave breaking in more generic cases than those covered by the asymmetry condition in the wave breaking theorems.
(iv) Global (in time) existence of strong solutions for (generic) initial data with small Sobolev norm.
We will show that this equation has a unique solution for every u 0 ∈ L 2 (T). The proof is carried out by a successive approximation as follows: u (0) (t) = e −tA u 0 , u (n+1) (t) = F (u (n) )(t) for n = 0, 1, . . . and the solution is then obtained by showing that F is a contraction mapping with respect to a suitable metric. In the sequel we will denote by c or c ′ various constants that may depend on ǫ but not on t. Furthermore, let γ 0 be a constant such that where here and hereafter denotes the L 2 -norm. Let R > 0 and suppose that we are given a continuous map t → w(t), [0, T ] → L 2 (T), which takes its values for t > 0 in H 1 (T) and satisfies We then have w(t) L ∞ ≤ c w(t) 1/2 ∂ x w(t) 1/2 ≤ cRt −1/4 and it is not difficult to deduce Similarly, where γ 0 is as above and B(·, ·) denotes Euler's beta function. Now we may take any R > max{ u 0 , γ 0 u 0 }. Then we may choose T small enough such that any w satisfying (18) implies that (18) holds also with F (w) in place of w.
If then the corresponding closed ball B R of radius R around 0 in Y T is mapped into itself (upon noting that t → ∂ x F (u)(t) is continuous ]0, T ] → L 2 for every u ∈ Y T ). Our next task is to show that the map F is a contraction with respect to the norm * (for some T > 0 and on bounded subsets). The proof is similar to the above estimates. Indeed, c w * w − z * s −3/4 + c z * w − z * s −3/4 + c w − z * ds ≤ c ′ t 1/4 w * + z * w − z * + ct w − z * .
and similarly for ∂ x F (u) − ∂ x F (z) . Thus, for T sufficiently small, F becomes a contractive mapping on B R and we obtain a unique solution u ∈ Y T of (17).
We emphasize that here T depends on R and ǫ, but not on any individual u 0 as far as u 0 satisfies max{ u 0 , γ 0 u 0 } < R. Since the spatial L 2 -norm in the solution is conserved over time, the solution may be continued any number of times, thus we obtain a solution globally in time. Therefore, we conclude the global unique existence of a solution to the parabolic equation (4). The following additional features of this solution are used in the proof of the weak solution: Due to the properties of the heat kernel, u ∈ C ∞ (T× ]0, T ]) and moreover, if u 0 ∈ L ∞ , although t → u(t) may be discontinuous at t = 0 as a map into L ∞ , the real-valued map t → u(t) L ∞ is continuous on the closed interval [0, T ].