Numerical conservative solutions of the Hunter--Saxton equation

In the article a convergent numerical method for conservative solutions of the Hunter--Saxton equation is derived. The method is based on piecewise linear projections, followed by evolution along characteristics where the time step is chosen in order to prevent wave breaking. Convergence is obtained when the time step is proportional to the square root of the spatial step size, which is a milder restriction than the common CFL condition for conservation laws.


Introduction
The Hunter-Saxton (HS) equation is given by It was derived, in differentiated form, from the nonlinear variational wave equation ψ tt − c(ψ)(c(ψ)ψ x ) x = 0 as an asymptotic model of the director field of a nematic liquid crystal [13]. Furthermore, the Hunter-Saxton equation is the high frequency limit of the Camassa-Holm equation [6]. It is completely integrable [14] and can be interpreted as a geodesic flow [17].
Another main property is that weak solutions are not unique, see e.g. [15,16]. The main reason being the following: Solutions of (1.1) may experience wave breaking in finite time, i.e., u x → −∞ pointwise while the energy u x (t, ·) 2 remains uniformly bounded and the solution u stays continuous. Furthermore, a finite amount of energy is concentrated on a set of measure zero.
We illustrate wave breaking with an example by considering a peakon solution -a soliton-like solution that is continuous and piecewise linear in space. It is not a classical solution. Indeed the function is not differentiable at the break points between the linear segments. with 0 ≤ t < 2. Note that for t < 2,

Example 1.1 (Wave breaking for peakons). A particular peakon solution that illustrates wave breaking is given by
shrinks to a single point (see Figure  1). One can check that the function u remains uniformly bounded and uniformly Hölder continuous with exponent 1 2 on [0, 2] × R. It is possible to extend weak solutions past wave breaking in various ways, see [1,2,3,9,19]. One could ignore the part of the solution that blows up. That amounts to continuing the solution in Example 1.1 as u(t, x) = 0 for all t ≥ 2. Such solutions are called (energy) dissipative and are unique [5,7]. A different approach is to "reinsert" the energy after wave breaking to get (energy) conservative solutions. To extend the solution in Example 1.1 as a conservative solution we let the formula defining u hold for t ≥ 2 as well. Uniqueness of conservative solutions is only known in several special cases [25,26]. The different solution concepts mimic the ones for some closely related equations: the Camassa-Holm equation [4], the nonlinear variational wave equation [21], and various generalizations of these equations.
From now on we focus on weak solutions that preserve the energy, that is conservative solutions. It has been shown in [2] that there exists a Lipschitz continuous semigroup of weak conservative solutions to (1.1). Existence of solutions is proved using Lagrangian coordinates and characteristics. Note that the curves describing the position of the break points in Example 1.1 are examples of characteristic curves.
To prolong the solution past wave breaking and to attempt to overcome the nonuniqueness of weak solutions past wave breaking, we include the cumulative energy F as part of the solution. The HS equation is then reformulated as with appropriate initial conditions and the conditions F (t, x) = µ(t, (−∞, x)) for some positive, finite Radon measure µ(t, ·), x (t, x) dx = µ ac (t, (a, b)), (1.3c) where µ ac (t, ·) is the absolutely continuous part of µ(t, ·). A closer look at the imposed conditions reveals that one challenge is to find a numerical method that respects condition (1.3c). The key to overcome this difficulty is to consider (1.2) with the slightly more general conditions (1.3a), (1.3b), and This new system is a reformulation of the so-called two-component Hunter-Saxton (2HS) system, which not only generalizes the HS equation, but can also be studied using the same methods and ideas, see [9] and [19].
then this property will be preserved along characteristics and no wave breaking takes place. Furthermore note that applying a piecewise linear projection operator to pairs (u, F ) satisfying (1.3c) yields pairs (ũ,F ) satisfying (1.4). Thus using the method of characteristics and piecewise linear projection operators as building blocks for a numerical method seems to be a good choice.
1.1. Numerical methods for the Hunter-Saxton equation. Despite receiving a considerable amount of attention theoretically, relatively little numerical work has been done on the Hunter-Saxton equation. In [11] a finite difference method was constructed and proved to converge to dissipative solutions. In [22] and [23] discontinuous Galerkin methods were introduced, followed by a convergence proof in the dissipative case but not in the conservative case. More recently, a geometric numerical integrator, based on the complete integrability of (1.1), was introduced and studied in [18]. The method seems to converge to the conservative solutions, but no proof was presented. In [20] a difference method that converges for smooth solutions of a modified Hunter-Saxton equation in the periodic setting was introduced. The analysis in [20] does not apply to our setting since the method relies crucially on the modification of the equation, and even for (1.1) the periodic case and the real line case are essentially different [24].
In this paper, we contribute to this line of research by introducing a convergent numerical method for the conservative solutions of the Hunter-Saxton equation (1.2). The method, introduced at the beginning of Section 2, is inspired by Godunovtype methods for conservation laws and is based on piecewise linear projections, followed by evolution along characteristics forward in time. As for finite difference (and volume) schemes for conservation laws, where one limits the time step ∆t to prevent shocks from occurring, we limit the size of ∆t to prevent wave breaking (see (2.3)). In contrast to the situation for conservation laws, we get the improved bound ∆t ≤ C √ ∆x for some C that depends on the initial data.
After establishing some a priori bounds of the numerical solutions in Section 2.1, we show in Section 2.2 that the numerical approximation converges with a rate of O( √ ∆x) to the unique solution of (1.2) whenever the solution is Lipschitz continuous. We also prove the existence of a convergent subsequence of the proposed numerical method in the general case, which converges to a weak solution preserving F . Unfortunately, the present lack of a satisfactory uniqueness theory for conservative solutions of (1.1) prevents us from guaranteeing that the sequence as a whole converges to the unique conservative solution. However, we perform numerical experiments towards the end of the paper, see Section 3, showing that the numerical approximations seem to converge towards the desired solutions also in the case of non-Lipschitz solutions.

Numerical conservative solutions of the Hunter-Saxton equation
For the (to be defined) numerical solutions to approximate conservative solutions of the HS equation, we will require that they mimic certain aspects of these solutions. In particular, we will design a method such that the numerical approximations are pairs (u, F ) in a suitable function space D, which resembles the one used for the 2HS system in [19]: Remark 2.2. Given a pair (u, F ) ∈ D, there exists a positive finite Radon measure µ, such that F (x) = µ((−∞, x)).
Let T t be the conservative solution operator associated to (1.2), as defined in [19], mapping every initial data (u, F ) to the corresponding solution at time t. For continuous and piecewise linear initial data (u, F ), the conservative solution of (1.2) takes a particularly simple form as long as no wave breaking takes place: The solution is again continuous and piecewise linear and the breakpoints x j (t) travel along characteristics, i.e. along the curves x j (t) given by we get x j (t)) = F (0, x j (0)), (2.2b) with linear interpolation between the breakpoints. Thus the equations (2.2) implicitly define the solution operator T t in the case of continuous and piecewise linear initial data (u, F ).
Turning our attention once more towards Example 1.1, we see that the two curves describe the position of the breakpoints. Furthermore, at the breaking time t = 2 we have x 1 (t) = x 2 (t). In the general case of a continuous and piecewise linear initial data (u, F ), wave breaking occurs at times t where at least two break points coincide, i.e., x j (t) = x k (t) for some j = k.
Using the above observations, we will now derive the numerical scheme. The idea is to use piecewise linear projection operators P ∆x to project the solution at each time step, and T ∆t to evolve the solution one time step ∆t ahead. To improve the readability, we define points in space and time t n = n∆t, n ∈ N, Definition 2.3. Define the projection operator P ∆x : D → D so that (ū,F ) = P ∆x (u, F ) is given byū with linear interpolation in between grid points ∆xZ.
Remark 2.4. The operator P ∆x is well defined since it is assumed that F is (left) continuous, and thus one can evaluate F at any point.
Assume now that the time step ∆t is so small that no wave breaking occurs as the piecewise linear approximation is evolved from one time step to the next. Then the scheme is defined by (U 0 , F 0 ) = P ∆x (u 0 , F 0 ) and (U n+1 , F n+1 ) = P ∆x T ∆t (U n , F n ) for n ≥ 0.
We will need to interpret the numerical solution as a function from [0, ∞) × R to R × R + .
That is, we follow the solution along lines x = x j from one time step to the next, and interpolate linearly in between.
After each evolution ∆t forward in time, the solution is projected onto the space of continuous piecewise linear functions. As multiple peakons can be glued together to form multipeakons, which solve (1.2), we can continue computing the solution forward in time after each projection.
Remark 2.6. Note that as the numerical approximation consists of linear interpolations between grid points and solving exactly between time steps, F ∆x satisfies We introduce a CFL-like condition that ensures that characteristic curves x j (t) do not collide as long as we evolve the equations less than ∆t. In particular, the condition prevents wave breaking, which occurs when x j (t) = x j+1 (t) for some j ∈ Z and t > 0. We arrive at the following bound on ∆t in terms of the initial data and the grid length ∆x. The condition is not a true CFL condition in the sense that characteristics may travel past several cells [x j , x j+1 ] during one time step.
Definition 2.7 (CFL-like condition). We require that ∆t satisfies Note that (2.3) is less restrictive than the CFL conditions used for conservation laws, which reads ∆t < C∆x for some C depending on the initial data and the particular flux function.
Remark 2.8. In the upcoming proofs we will use Note that we could have chosen any fixed α ∈ (0, 1] to take the step from (2.3) to (2.4) (with 1 replaced by α). As a consequence the least distance between characteristics x j (t) and x j+1 (t), starting from neighboring grid points, would be given by β(α)∆x and could be computed using (2.1).
Similarly to the forward characteristics governed by (2.1), there are characteristics backwards in time. In particular, we can associate to any grid point (x j , τ ) with t n ≤ τ ≤ t n+1 , the unique point (t n , ξ n j (τ )) given by ). Remark 2.10. The numerical scheme can be written in the more familiar form where the backward characteristic from x i at t n+1 satisfies ξ n i (∆t) ∈ [x j , x j+1 ], see (2.5).

2.1.
A priori bounds of the numerical solutions. In this section, we prove certain a priori bounds of the proposed method, which are needed to prove convergence. We begin with some preliminary results on the projection operator P ∆x .
Proof. For any grid point which proves the first inequality. Next, we have From the L 1 -estimate one can obtain the L 2 -estimate, To prove that the numerical approximation converges, we wish to employ the Arzelà-Ascoli theorem to ensure convergence of a subsequence of u ∆x , and subsequently a version of the Kolmogorov compactness theorem to get convergence of a subsequence of F ∆x . To invoke the Arzelà-Ascoli theorem, we need u ∆x to be uniformly equicontinuous and equibounded. For the Kolmogorov compactness theorem we need that F ∆x is of uniformly bounded total variation, that F ∆x (t, ·) is continuous in t in the L 1 (R)-norm, and that F ∆x (t, ·) does not escape to infinity as ∆x tends to zero. First we establish some immediate properties of the solutions (u ∆x , F ∆x ).
The bounds on u ∆x (t, x) and F ∆x (t, x) follow from (2.2) and Definition 2.5. Since both (2.2) and the projection operator preserve the monotonicity of F , we have that F ∆x is monotone increasing. Continuity follows from the fact that characteristics emanating from different grid points are at least 1 2 ∆x apart as long as the time step is controlled by (2.4). We To begin with let t = 0. Since u ∆x (0, ·) and F ∆x (0, ·) are both piecewise linear and continuous it suffices to show the result for x j ≤ a ≤ b ≤ x j+1 . By assumption one has that b a u 2 Now, let t = t n + τ , and denote by τ → (ũ(τ ),F (τ )) the conservative solution with initial data (u ∆x (t n ), F ∆x (t n )). Furthermore, assume that (u ∆x (t n ), F ∆x (t n )) satisfies (2.7c). Then we have for each spatial grid point x j thatũ(τ, since this property is preserved along characteristics. Applying the projection operator, we can follow the same lines as in the case t = 0, to obtain that (2.7c) holds for all t ∈ [t n , t n+1 ]. By . Let x j− be the closest gridpoint to a from below, and let x j+ be the closest gridpoint to b from above. Then F ∆x, Next we show that also F ∆x,x (∆t, ·) is compactly supported. By (2.1), we have Here it is essential that Since ∆x = 4F ∞ ∆t 2 , we have that (k + 1)∆x = ∆x + 4F ∞ (k∆t)∆t. From the interpolation between temporal grid points we get The total variation estimate follows from the fact that it holds for conservative solutions, and that the projection operator can only reduce the total variation.
Remark 2.13 (Spatial Hölder continuity). An immediately derivable property of the numerical solution from (2.7c) is spatial Hölder continuity of u ∆x : In order to obtain temporal Hölder continuity for u ∆x we will need to compare a numerical solution with itself several time steps ahead. Lemma 2.14. For each i, n, k there are non-negative constants β ink j such that Proof. We prove the lemma by induction on k. First note that the statement is trivially true for k = 0. Then assume that it holds for k = l. We show that it must then hold for k = l + 1 as well. We have that SinceC ∆x is greater than the C ∆x in the inductive assumption, we get The computation for U n+k i is analogous. Indeed, we have Next is an important corollary which provides a discrete Hölder continuity estimate for the numerical solution u ∆x .
Proof. Using Lemma 2.14, we compute and thus, remembering Remark 2.13, (2.8c), and (2.4), We are now ready to prove that for each T > 0 the solutions u ∆x are uniformly Hölder continuous on [0, T ] × R. Uniform Hölder continuity implies equicontinuity, which is necessary for the Arzelà-Ascoli theorem. Lemma 2.16 (Hölder continuity). Let 0 ≤ t, s ≤ T and x, y ∈ R, then Proof. Assume first that t n ≤ s < t ≤ t n+1 and x j ≤ x ≤ x j+1 . We start by adding and subtracting u ∆x (s, x) and obtain Then, we have by definition, Note that at the spatial grid points x l the solution u ∆x (t, x l ) equals the conservative solution given by (2.2) with initial data (u ∆x (t n , ·), F ∆x (t n , ·)) evolved t − t n < ∆t forward in time. For conservative solutions given by (2.2) we do have Hölder continuity with the constant C depending on F ∞ , u 0 ∞ , and T only. To be more specific it has been shown in the proof of [9, Theorem 3.14] that for all j ∈ Z. Hence, we have We look at the general case t n−1 ≤ s ≤ t n ≤ t n+k ≤ t ≤ t n+k+1 , and x j−1 ≤ y ≤ x j ≤ x j+l ≤ x ≤ x j+l+1 . Then, by Corollary 2.15, we have,

NUMERICAL CONSERVATIVE SOLUTIONS OF THE HUNTER-SAXTON EQUATION 13
To use the Kolmogorov compactness theorem we need uniform regularity of F ∆x in t.
Otherwise, we have that there exist j − and j + , possibly equal, such that The number of terms in the above sum is bounded from above by |ξ n j (t) − ξ n j (s)|. Direct calculations yield ξ n j (t) − ξ n j (s) ≤ u ∆x (t n , ξ n j (s))(s − t n ) − u ∆x (t n , ξ n j (t))(t − t n ) and therefore Due to condition (2.4) characteristics (forward as well as backward) from neighbouring grid points have a minimum distance of 1 2 ∆x. Hence for each j ′ , the maximal number of backward characteristics ξ n j (∆t) ending up in [x j ′ , x j ′ +1 ] equals two. Hence we have the bound The general case 0 ≤ s < t ≤ T can now be found by iteration over time steps and using condition (2.4),

2.2.
Convergence of the numerical solutions. In this section we prove that there exists a convergent subsequence of (u ∆x , F ∆x ), and that the limit is a conservative weak solution of (1.2), which satisfies condition (1.4). First we rigorously define, as in [2,9,19], conservative weak solutions.
, for all T ≥ 0, (u(t), F (t)) ∈ D for all t ≥ 0, where µ(t) is the finite positive Radon measure with F (t, ·) as its distribution function, see Definition 2.1.
We prove the existence of a convergent subsequence of (u ∆x , F ∆x ).
Theorem 2.19. To any initial data (u 0 , F 0 ) ∈ D such that µ 0 has compact support, there exists a convergent subsequence of (u ∆x , F ∆x ). The convergence is in C [0, T ], L 1 (R) , pointwise a.e. in x for F ∆x , and uniform on [0, T ] × R for u ∆x . Moreover, the limit (u, F ) satisfies Here the relation between the positive Radon measure µ 0 and F 0 is given by F 0 (x) = µ 0 ((−∞, x)).
Proof. We have from Lemma 2.12 that the family u ∆x is uniformly bounded on [0, T ]×R and that u ∆x,x (t, ·) has compact support for all t ∈ [0, T ]. Furthermore, by Lemma 2.16 u ∆x is uniformly equicontinuous. Hence the conditions for the Arzelà-Ascoli theorem are satisfied and there exists a convergent subsequence (u ∆x,i , F ∆x,i ) of (u ∆x , F ∆x ) such that u ∆x,i converges to some u ∈ L ∞ ([0, T ] × R) for each T > 0. The limit of u ∆x,i is bounded and Hölder continuous with the same constants as the individual u ∆x,i . Next we show that the limit u satisfies u x (t, ·) ∈ L 2 (R) for all t ∈ [0, T ]. By construction we have that u ∆x,x (t, ·) L 2 ≤ √ F ∞ for all t ∈ [0, T ]. Thus there exists a subsequence (u ∆x,i k (t, ·), F ∆x,i k (t, ·)) of (u ∆x,i (t, ·), F ∆x,i (t, ·)), so that u ∆x,x,i k (t, ·) converges weakly to v in L 2 (R). Thus for any φ ∈ C ∞ c (R) we have and v(·) = u x (t, ·). Thus we have u ∆x,x,i k (t, ·) ⇀ u x (t, ·) in L 2 (R). A closer look reveals that the above argument shows that every weakly convergent subsequence has the same limit and therefore u ∆x,x,i (t, ·) ⇀ u x (t, ·) in L 2 (R) for all t ∈ [0, T ]. Combining Lemma 2.12 and [10, Theorem 12], we obtain, that for each t ∈ [0, T ], there exists a subsequence (u ∆x,ij (t, ·), F ∆x,ij (t, ·)) of (u ∆x,i (t, ·), F ∆x,i (t, ·)) such that F ∆x,ij (t, ·) converges pointwise everywhere and in the L 1 -norm to a function of bounded variation. Following the lines of the proof of [12,Theorem A.11] and taking into account Lemma 2.17, it then follows that there exists a subsequence of (u ∆x,i , F ∆x,i ), for which the F ∆x converge in C [0, T ], L 1 (R) . Furthermore, denoting the limit by F , we even have pointwise almost everywhere convergence of a further subsequence to F .
Last but not least, we have a look at the connection between u x and F . Denote by (ũ ∆x ,F ∆x ) the very last subsequence of (u ∆x , F ∆x ), then we have that b a u 2 Since, we can find two sequences a j ↓ a and b j ↑ b such that lim ∆x→0 F ∆x (t, a j ) = F (t, a j ) and We still need to prove that the limit of the convergent subsequence is a conservative weak solution in the sense of Definition 2.18. Proof. It remains to show that the integrals (2.9) and (2.10) hold. We compute the integrals for (u ∆x , F ∆x ) as follows Here (ũ n (τ, x),F n (τ, x)) denotes the conservative solution given by (2.1) and (2.2) with initial data (u ∆x (t n , x), F ∆x (t n , x)). Since the conservative solution is a weak solution we get Recalling that u ∆x (t n + τ ) = P ∆xũn (τ ) yields where we applied Proposition 2.11 in the last step as follows Here T φ = inf{t ≥ 0 | φ(s, ·) 2 + φ t (s, ·) 2 = 0 for all s ≥ t}. Note that T φ is finite since φ has compact support. We now turn our attention to the first sum. Recall that u ∆x (t n +τ ) = P ∆xũn (τ ), F ∆x (t n + τ ) = P ∆xFn (τ ), and keep Proposition 2.11 and Lemma 2.12 in mind.

Direct calculations then yield
Since φ has compact support we end up with In particular, we have By letting ∆x → 0, we end up with (2.9) as the subsequence of (u ∆x , F ∆x ), constructed in Theorem 2.19, converges to (u, F ) uniformly in [0, T ] × R for u ∆x and in C [0, T ], L 1 (R) for F ∆x .
In a similar fashion we demonstrate that the second integral equation (2.10) must be satisfied as ∆x → 0 as well. We have

NUMERICAL CONSERVATIVE SOLUTIONS OF THE HUNTER-SAXTON EQUATION 19
Since conservative solutions are weak solutions we get Summation over n then yields From the estimates in Proposition 2.11 we get The second term can be estimated as follows It remains to prove that the first term tends to zero as ∆x → 0. We compute Finally, we end up with (2.12) Since for every t we have that F ∆x (t, ·) → F (t, ·) almost everywhere, the convergence of (2.12) to (2.10) as ∆x → 0 follows from the proof of Proposition 7.19 in [8] and the fact that u ∆x → u in C([0, T ] × R).
A satisfactory uniqueness theory for conservative weak solutions of the Hunter-Saxton equation would have ensured that all limits in Theorem 2.19 are equal, and thus that the sequence as a whole converges. Unfortunately, uniqueness of conservative solutions is unknown at the present time. On the other hand it is known that if the initial data (u 0 , F 0 ) is Lipschitz continuous, also the solution (u(t, ·), F (t, ·)) will be Lipschitz continuous for all t ∈ [0, − 2 min(u0,x) ), at least. In particular, as Example 1.1 and 3.2 indicate, when wave breaking occurs u(t, ·) may be Hölder, but not Lipschitz continuous and F may even be discontinuous.
In the next theorem, we consider the special case of weak conservative solutions (u, F ) such that (u(t, ·), F (t, ·)) are Lipschitz continuous for all t ∈ [0, T ].
Proof. For any conservative solution (u, F ) ∈ [W 1,∞ ([0, T ] × R)] 2 with initial data (u 0 , F 0 ), the characteristic equation is uniquely solvable. Furthermore, the classical method of characteristics implies that the solution is unique and given by Introduce (ũ 0 ,F 0 )(t) = T t P ∆x (u 0 , F 0 ) and recall that x j (t) denotes the characteristic starting at the grid point x j , then (ũ 0 (t, x j (t)),F 0 (t, x j (t))) = (u(t, x j (t)), F (t, x j (t))) for all j ∈ N.

Numerical examples
In this section we perform two experiments to see whether the scheme in Definition 2.5 converges to the desired conservative solution. We compare the numerical solutions with known, exact solutions in two cases, namely a peakon and a cusp. These two have been selected since the exact solutions represent two distinct challenges for the numerical solver. The peakon solution experiences wave breaking at t = 2 and all energy is concentrated in a single point. Thus F becomes discontinuous while u becomes a constant function. The cusp solution, on the other hand, experiences wave breaking at each time t with t ∈ [0, 3], but only an infinitesimal amount of energy concentrates at any given time.
In our examples we assume that the initial data u 0 is constant outside some finite interval [a, b]. By (2.2) and Lemma 2.12, we then obtain that at each time t, both u(t, ·) and u ∆x (t, ·) will be constant outside some finite interval [a(t), b(t)]. Thus for any T > 0, by choosing the computational domain accordingly, we can ensure that u(t, ·) and u ∆x (t, ·) are constant outside the computational domain for all t ∈ [0, T ]. We look at the peakon example first.

Example 3.1 (Peakon solution). We have initial data
with the exact solution Figure 2 the numerical solution (u ∆x , F ∆x ) is computed and compared with the exact solution at t = 0, t = 2, and t = 4. Figure 3 shows the error, when compared to the exact solution, and we see that the method captures the wave breaking phenomena in this example. with the exact solution  3 , 1 + t + 1 3 t 2 < x. In Figure 4 the numerical solution (u ∆x , F ∆x ) is computed and compared with the exact solution at t = 0, t = 2, and t = 4. Figure 5 shows the error when compared to the exact solution.
From Figure 3 and 5, we see that the best we can hope for in terms of convergence rates in L ∞ for u and L 1 for F in the general case is O √ ∆x . As uniqueness of conservative solutions is still an open problem, proving any form of convergence rate of the numerical method seems to be extremely challenging. Remark 3.3. Clearly we cannot expect a better convergence order than one half in the L ∞ -norm for u, since there exists u 0 ∈ D such that u 0 −u 0∆x ∞ = √ F ∞ √ ∆x, see Proposition 2.11.  Figure 4. The functions u ∆x (left) and F ∆x , (right) in the case of cusp initial data, plotted at t = 0, t = 1.93, and t = 4. Here ∆x = 1/4 and ∆t ≈ 0.148. Note the slight discrepancy between the numerical solution and the exact solution in the variable F already at t = 0 due to the projection operator being applied to the numerical initial data.