The r-Hunter-Saxton equation, smooth and singular solutions and their approximation

In this work we introduce the r-Hunter-Saxton equation, a generalisation of the Hunter-Saxton equation arising as extremals of an action principle posed in L_r. We characterise solutions to the Cauchy problem, quantifying the blow-up time and studying various symmetry reductions. We construct piecewise linear functions and show that they are weak solutions to the r-Hunter-Saxton equation.


Introduction
The Hunter-Saxton (HS) equation is a 1+1 dimensional, variational, partial differential equation (PDE), originally introduced as a model for the propagation of waves in the director field of a nematic liquid crystal [HS91,HZ95a,HZ95b]. This problem arises in three different forms, related to each other formally through differentiation, where the subindices denote partial differentiation with respect to the corresponding independent variable. These equations can be considered over a domain Ω that is either the real line, the periodic interval or the interval [a, b] with zero boundary conditions u(a) = u(b) = 0.
The HS equation also has an important geometric interpretation as it describes the geodesic flow on the diffeomorphism group on the domain Ω, with the right-invariant metric defined through the H 1 -inner product Consider time-parameterised diffeomorphisms g ∈ Diff(Ω) with (5) g t (x, t) = u(g(x, t), t), then the HS equations can formally be interpreted as extreme points for the action principle and endpoint conditions δu(x, 0) = δu(x, T ) = 0. This action seeks a path of diffeomorphisms g(x, t) with g(x, 0) = g 1 (x) and g(x, T ) = g 2 (x) such that the distance functional l[u] is extremised. This problem is the one-dimensional version of problems that arise in computational anatomy [MTY06]. Then, (3) is obtained as the corresponding Euler-Poincaré equation for this action principle after observing that (5) implies that we must have constrained variations δu = w t + wu x − uw x for arbitrary w. The conserved energy then arises from Noether's theorem applied to time translation symmetry; we also obtain a Poisson structure for (3) through this route. The Euler-Poincaré derivation is a formal calculation that was made rigorous in [KM03]. The integrability of the HS equation was shown in [HZ94] by relating it the Camassa-Holm equation TP was partially supported through the EPSRC grant EP/P000835/1. CJC was partially supported through the EPSRC grant EP/R029423/1. This support is gratefully acknowledged. Indeed, it can be shown that (3) arises from a high frequency limit of (8) and, thus, has a bi-Hamiltonian structure [AH09]. Further connections can be made through the geometric interpretation as the Camassa-Holm equation describes geodesic flow with respect to the right-invariant metric [Kou99] (9) f, g = Ω f x g x + f g d x.
The HS equation has a conserved energy, which represents one of the two Hamiltonians of the problem, the other being The explicit control of these two functionals ensures some regularity of the solution. Solutions to this equation can even be written down explicitly. The HS equation given in the forms (1-2) has piecewise linear weak solutions that conserve this energy [HZ94]; these solutions can be described by a set of moving points plus the value of u at those points with linear interpolation in between. In [BC05] it was shown how to make sense of these piecewise solutions as weak solutions, and that weak solutions of HS are locally Lipschitz with respect to the initial conditions, making use of techniques from optimal transport. The locality is clear since it is possible to find piecewise linear solutions where two points with different u values collide in finite time, leading to a jump in the solution.
In this paper, we investigate the effect of modifying the distance functional l[u] such that corresponding to the W 1,r (Ω) distance, for r > 2 rather than for r = 2. To the author's knowledge, previous work on quantifying the nature of solutions to PDEs arising from energy functionals posed on L r is limited to the elliptic case. Indeed, a sequence of works initiated by [Aro65] focussed on the study of the Euler-Lagrange equations of (12) in the multi-dimensional setting as r → ∞, see [Kat15,Pry18] for an accessible overview. In [KP16] the authors obtained results for the vectorial analogy of (12) and in [KP16, KM17, KP18] a minimisation problem involving second derivatives was examined. In this work, for the first time, we examine the effect modifying the distance functional has on an evolution problem. For convenience we take r to be an even positive integer. One of the motivating reasons to modify the distance functional is that a certain amount of regularity is required for a solution to exist to (5) [DGM98], specifically u(x, t) is required to be in L 1 ([0, T ], W 1,∞ (Ω)). For the HS equations with r = 2, this is not satisfied in general and we see the loss of the diffeomorphism property when the piecewise-linear solutions blow up as above. With that in mind, we are interested in how the solutions to the resulting PDE behave as r → ∞. We call the resulting equation the r-Hunter-Saxton (r-HS) equation.
To extend the notion of piecewise linear solutions of the HS equation to the r-HS equation we use an optimal control formulation that arises when trying to optimise l[u] such that a set of points moving with u are transported from one configuration to another. We shall see that the optimal u is then piecewise linear. Following [CH09,GBR11], after eliminating u, we obtain a Hamiltonian system for the point locations q and their conjugate momenta p, with the conserved energy being equivalent to l [u]. For u in a finite dimensional state space, it can be shown that after eliminating p and q, u satisfies the Euler-Poincaré equation corresponding to u. However, in this particular case, the r-HS equation does not make sense because piecewise solutions do not have enough regularity. In this paper, we resolve this situation by showing that the piecewise linear solutions are weak solutions of another equation whose formal derivative gives the r-HS equation.
The rest of this paper is organised as follows. In Section 2 we derive the r-HS equation. In Section 3 we study some fundamental properties of the r-HS equation, specifically we characterise smooth solutions to the problem, give a blow-up criteria and give examples of special solutions arising from a symmetry reduction technique. In Section 4 we introduce piecewise-linear functions as solutions to an optimal control problem for points on an interval. In Section 5 we show that these piecewise-linear functions are weak solutions of an integrated form of the r-HS equation. In section 6 we calculate some numerical examples, and finally in 7 we provide a summary and outlook.

The r-Hunter-Saxton equation
In this section we introduce the r-HS equation, describe some of its properties.
Definition 2.1 (r-Hunter-Saxton equation). The r-Hunter-Saxton equation on the interval [a, b] in its most general form is given by with boundary conditions u(a) = u(b) = 0.
One should notice that, as with the (2)-HS equation, the general r-HS equation has various equivalent forms.
Proof. Firstly note that Similarly, Making use of these we see that Hence (14) is the formal derivative of (15). To see (16) note that Hence as required.
The equations can also be defined on the real line, or with periodic boundary conditions, but we concentrate on the boundary value problem in this paper for simplicity.
Proof. We have through an integration by parts making use of the boundary conditions u(a, t) = u(b, t) = 0. Following [HMR98], (5) implies that with w(x, 0) = w(x, T ) = 0. Substituting (27) into (24), we obtain which is satisfied by solutions of (13) for arbitrary w, as required.
Definition 2.4 (Weak integrated r-HS equation). Let W 1,r 0 (a, b) be the space of functions with for all test functions φ ∈ W 1,r 0 (a, b).

Symmetries and classical solutions
In this section we examine some symmetries and characterise classical solutions and their existence time for the r-HS equation. We take inspiration from the arguments in [HS91,§3] and show that, remarkably, many of the properties shown for the r = 2 case generalise for arbitrary r. We begin by examining a characteristic reduction. To that end, throughout this section we will assume u ∈ C 2 (Ω) be a classical solution of the Cauchy problem Also for this section and the remainder of this work, for clarity of exposition, we will continue under the assumption that r ≥ 2 is even.

Theorem 3.1 (Smooth solution characterisation).
Let H ∈ C 1 (R) be any function with H(0) = H (0) = 0 and, for k ∈ N, let G k ∈ C 2 (R) be defined by requiring Then, all classical solutions of (32) are given through the implicit relation Proof. Let ξ denote a characteristic curve with U (ξ, t) = u(X(ξ, t), t) and X satisfying the initial value problem and hence solves the second order initial value problem This can then be solved to show that This should be compared with the r = 2 case found in [HS91,§3]. The final result follows from solving the system for U as required.
Proof. The implicit function theorem guarantees that as long as X ξ = 0 there is a smooth solution of (32). Given .
Remark 3.3 (Relating to Burgers' equation). Notice the maximal smooth solution time of (32) increases linearly as r increases. In fact, it is exactly r-times the shock time for inviscid Burgers' equation. In particular, as r → ∞, the blow up time T → ∞. We will return to this point later in this work.

Symmetries. A Lie point symmetry of equation (32) is a flow
generated by a vector field such that u( x, y) is a solution of (32) whenever u(x, y) is a solution of (32). As usual, we denote by e X the Lie series ∞ k=0 k k! X k with X k = XX k−1 and X 0 = 1. To find the symmetries of the r-HS equation we are required to solve the infinitesimal invariance condition for the vector field (43). To do this we use the prolongation of X [BA08, Olv93,Ste89]. The infinitesimal symmetry condition decomposes to a large overdetermined system of linear PDEs for ξ 1 , ξ 2 and η known as determining equations. Note that this procedure is described in further detail in [PP19].
Proposition 3.4 (Infinitesimal invariance). The infinitesimal invariance condition is equivalent to the following system of 14 equations: Solutions of the overdetermined system of linear PDEs (44)-(47) will yield the algebra of the symmetry generators (43) of the r-HS equation.
Given (44)-(47) form an overdetermined system of linear partial differential equations it is possible that they only admit the trivial solution ξ 1 = ξ 2 = η = 0. This would imply that the only Lie symmetry of the r-HS equation is the identity transformation. In what follows we will see that this is not the case. We are able to obtain the Lie algebra for the symmetry generators the r-HS equation and thus, using the Lie series, derive the groups of Lie point symmetries.
Proposition 3.6 (Lie algebra generators). It follows that the solution (48) defines a six dimensional Lie algebra of generators where a basis is formed by the following vector fields Invariant solutions through symmetry reductions. We now state solutions that occur through symmetry reductions of (32) to ODEs by means of the algebra generators given in Proposition 3.6. We consider each generator seperately and examine some examples of solutions from each. X 1 . To begin notice that solutions of (32) that are invariant under the symmetry generated by X 1 are of the form u = f (t), which immediately yields u ≡ const as a trivial solution prescribed by the initial condition.
X 2 . Solutions invariant under the symmetry generated by X 2 are of the form u = f (x). The reduced equation is given by the following ODE: . This, in turn shows that either f ≡ const or f 2 x + rf f xx = 0 and must take the form for arbitrary constants c 1 , c 2 .
X 3 . Solutions invariant under the symmetry generated by X 3 are of the form u = xt −1 + f (t) with f prescribed by the initial condition.
X 4 . The quantities u and ξ = tx −1 are algebraic invariants of the Lie group generated by X 4 . Assume u = f (ξ) for non-constant f , then we obtain the reduced equation as the following ODE: The general solution of this ODE is not known.
X 5 . Solutions invariant under the symmetry generated by X 5 are of the form u = f (x)t −1 , for non-constant f . The reduced ODE is given by: whose solution is an inverse hypergeometric function, that is for constants c 1 , c 2 . X 6 . The quantities u and ξ = tx −1/r are algebraic invariants of the Lie group generated by X 6 . With u = f (ξ) we may derive the following reduced ODE: which has closed form solution (57) f (ξ) = c 2 exp 2r log(c 1 + ξ r ) +(−2r − 4) log(ξ) 2(r + 1) .

Piecewise linear solutions
In this section we consider an optimal control formulation that allows us to quantify piecewise linear solutions to the problem. We show that the computation of these solutions requires the solution of a nonlinear system and make various comparisons to the 2-HS equation. We leave the formal interpretation of these solutions to the next section.
Definition 4.1 (Optimal control problem). Let Q 1 (t), . . . , Q N (t) represent a moving set of points on the interval [a, b]. The W 1,r optimal control problem is to find Q 1 (t), . . . , Q N (t) and u ∈ W 1,r (a, b) such that for i = 1, . . . , N .
This problem has the following variational formulation.
This type of variational principle takes its name from variational principles of this form that can be used to derive equations of fluid dynamics. In [CH09,GBR11] these variational principles were considered in a general form where the Lagrange multipliers enforce dynamics (61) .
where Q is defined on a manifold M, u is defined on a Lie algebra X , and L u represents a Lie algebra action on M. In this case, it can be shown that P and Q can be eliminated, and u solves the corresponding Euler-Poincaré equation. In our case, we have (L u Q) i = u(Q i , t), which is a Lie algebra action of smooth vector fields on [a, b], However, W 1,r (a, b) is not closed under this bracket, and we see that we do not have enough regularity to complete this argument. However, as we shall find, it is still possible to make sense of these solutions.
Lemma 4.3. The optimal u to the variational problem (60) are given by piecewise linear functions, whose jumps occuring at x = Q i , i = 1, . . . , N , and jump condition Proof. The Euler-Lagrange equation corresponding to u is given by which is, in turn, the weak form of the r-Laplace equation where δ(x) is the Dirac measure. This means that u is piecewise linear, with jump conditions given by (62).
Notice that this means that the term (|u x | r−2 u x ) x u x in (13) does not make sense, since at x = Q i it is the product of a Dirac measure with a function that has a jump at the same location. Later we shall see that this problem is resolved since u solves (31) nevertheless.

Remark 4.4. Equation (64) defines an infinitesimally equivariant momentum map from
The equivariance explains why it is possible to eliminate P and Q and obtain an equation for u alone. This was explored in the context of peakon solutions of the Camassa-Holm equation in [HM05], where the infinitesimal equivariance is proved.
Lemma 4.5. The piecewise solution can be characterised through the points u i = u(Q i ), i = 1, . . . n with u 0 = 0, u n+1 = 0. Then, the set of coefficients { u i } n+1 i=0 solves the difference equation Proof. The piecewise linear solution has piecewise constant derivative, and substitution into (62) gives the result. Proof. Let The solution { u i } n+1 i=0 is then the minimiser of Writing ∆ u i = u i+1 − u i , i = 0, . . . , n, we see that {∆ u i } n+1 i=0 minimises the function subject to the constraint n i=1 ∆ u i = 0. The function is the sum of convex and linear functions, which is then also convex. Since the constraint defines a linear subspace, we are minimising a convex function over a finite dimensional vector space, which has a unique solution.
The solution does not have a closed form analytic expression but can be solved numerically using Newton's method. We denote the solution operator u = U (P, Q).
Theorem 4.7. The equations for P and Q satisfy . .
with the conserved energy .
as expected, after noting that the bracket vanishes since u satisfies (65). After noting that we see the equation for P is given by Note that when r = 2, the P equation specialises to which recovers an equation from the standard (r = 2) Hunter-Saxton case.
Remark 4.8. In the r = 2 case, this formula allows an ad hoc interpretation of mu x term applied the piecewise linear solutions, when m contains Dirac measures centred on points where u x jumps. This interpretation says that we can just take an average of u x from each side of the jump point, but this only coincidentally works for the r = 2 case, and is not true for r > 2.
Remark 4.9. The same (P, Q) dynamics is obtained from the following variational principle, This dynamics has a Q-dependent Lagrangian but trivial Lie algebra action of u on Q. The equivalence of these two formulations becomes clear after noticing that they result in the same Hamiltonian after Legendre transformation. However this second variational principle has more directly computable Euler-Lagrange equations, since the vanishing terms above simply do not appear in the first place.
Corollary 4.10. Consider initial conditions P i (0) = z r−1 i . Then, the r → ∞ limit of Equations (70-72) are given by This limit makes sense, since it keepsû i finite (for distinct Q i ).
Proof. First we examine the difference equation (65). After substituting z i into P i and taking the r − 1 root, we obtain which recovers (85) in the limit. Substituting z i into Equation (71) gives and hence .
which converges to zero as r → ∞, provided that ( The well-posedness of Equations (70-72) depends on the solveability of (72), which a linear system in the maxalgebra (the algebra where the addition operator is replaced by the max operator) that can be formulated as a linear program. The convergence of solutions of Equations (70-72) to solutions of Equations (84-85 is an open question, which we investigate through special cases and numerical solutions in Section 6.

Weak solutions of the r-HS equation
As previously discussed, it is natural to hope the results of [CH09,GBR11] allow us to show that the piecewise linear u derived in §4 satisfies the corresponding Euler-Poincaré equation, i.e., the r-Hunter-Saxton equation (13). However, since the corresponding Lie bracket is not closed in W 1,r , it is only closed in C ∞ , the results of those papers do not hold. Indeed, (13) is not well-defined for piecewise linear solutions, not even weakly, since the integral is not defined. However, remarkably, the piecewise linear solutions are solutions of (31). Proof. Let u(x, t) be a time-dependent piecewise linear function over the intervals [Q i (t), Q i+1 (t)], with u i = u(Q i (t)). Then, Now, since u is piecewise linear, u x is constant in every interval [Q i , Q i+1 ], hence and, upon splitting the integral into intervals and integrating by parts, Combining these together, we have where Taking since Φ i are arbitrary. Now we show that our singular solutions also satisfy this equation. First, define c(t) according to (65) and (71), we have

Now combine equation
which means that (100) by induction, and we have our required equation.

Numerical examples
In this section we compute a number of examples of piecewise linear solutions, obtained by integrating equations (70-72).
6.1. One point solutions. In the case of a single point Q 1 , we may use energy conservation to find an algebraic dependence between u 1 and Q 1 . Specialising to [a, b] = [0, 1] to simplify, the conserved energy is Since E is time-independent, we may compute it from the initial condition, and then assuming a positive root (it is a consequence of the one point P equation that u will stay positive if it is positive initially). We may then integrate the first-order equation . Q 1 = u(Q 1 ) (this must be done numerically in general). Since u is always positive, Q 1 is monotonically increasing, and so it makes sense to consider what happens when Q 1 → 1. Writing Q 1 = 1 − , we obtain the asymptotic formula (104) u ∼ E 1/r r 1/r (r−1)/r .
The equation indicating that the W 1,∞ -norm blows up in finite time.
Snapshots of a typical solution are shown in Figure 1. A comparison of the behaviour of single point solutions is made in Figure 2. We observe that the u − Q 1 relation approaches u = min(Q 1 , 1 − Q 1 ) as r → ∞. A plot of the derivative u x in [q, 1] is given in Figure 3. We observe that whilst the solution always blows up in finite time, the blow up time is later for higher r, and is roughly proportional to r, as predicted by the asymptotic approximation. This is also compatible with the intuition that the r → ∞ limit should preserve the W 1,∞ norm.
For the single point case, the r → ∞ limit suggested in Corollary 4.10 results in the equations We also deduce that . z = 0 from taking the limit in the p-equation. For initial conditions Q 1 (0) > 0.5, we obtain i.e., the point does not attain its supremum in finite time, and hence the solution does not blow up in finite time.
The infinite limit solutions are compared with the finite r solutions in Figure (4).
6.2. Two point solutions. For two or more points, unless a symmetric solution is sought, there is no closed form for the u − q relation and equation (72) must be solved iteratively using Newton's method. The result can then be used in an explicit time-integrator such as the standard 4th-order Runge-Kutta scheme that we used for these examples. In this section, we consider three cases, as follows.
The numerical solutions are shown in Figure 5. We note that finite time singularities occur in all three cases, even the chasing collision (in contrast to the same case for the Camassa-Holm equation, where the solitons transfer momentum at a distance and the solution does not blow up). Finally, we investigate the approximation of smooth solutions by piecewise-linear solutions and their subsequent evolution in time. We take 101 equispaced points on the interval [0, 1] and interpolate the function u(x) = sin(2πx) to produce an initial condition for {Q i } 100 i=0 , {P i } 100 i=0 . We observe a strong jump in the derivative emerging for large r. The results are shown in Figure 6. |u + x | r = 2 r = 4 r = 6 r = 8 r = 10 r = 12 r = 14 r = 16 r = 18 r = 20

Summary and outlook
In this paper we introduced an extension of the Hunter-Saxton equation, which we call the r-Hunter-Saxton equation. These equations are the Euler-Poincaré equations with reduced Lagrangian given by the W 1,r norm, and are associated with the evolution of geodesics in the diffeomorphism group with metric defined by the W 1,r norm. This replaces the linear Laplacian relating momentum and velocity by a nonlinear r-Laplacian. We introduced an optimal control variational principle for parameterising finite dimensional subspaces of the solutions of PDEs, and derived a Hamiltonian system for a finite dimensional set of points Q 1 , . . . , Q n plus their conjugate momenta. This corresponds to piecewise linear functions u with jumps in the derivative at each of these points. Although these piecewise functions are not smooth enough to satisfy the r-Hunter-Saxton equation, we showed that remarkably they are still weak solutions of an integrated form of the r-Hunter-Saxton equation.
There are plenty of interesting open questions about the r-Hunter-Saxton equations and the piecewise linear solutions in the limit as r → ∞. We have presented some evidence that solutions of Equations (70-72) converge to solutions of Equations (84-85) (before blow-up). We have also presented evidence that solutions of (84-85) exist for all times, whilst solutions of (70-72) blow up in finite time.
It is very interesting to link the solutions to (84-85) back to the r → ∞ limit of the optimal control problem in Definition 4.1. This is dangerous, since it involves exchanging the limits r → ∞ and the limits defining variational derivatives. However, it is tantalising that the solutions to (84-85) preserve their W 1,∞ norm.
Another interesting question is whether solutions of (70-72) can be used to approximate smooth solutions of the r-Hunter-Saxton equation. The control provided over the W 1 p norm of the numerical solution should provide a useful tool for this.
A final direction of enquiry is, what is the correct r → ∞ limit of the r-Hunter-Saxton equation? Can the solutions of the limiting equation be approximated by (84-85)? Does the limiting solution preserve the W 1 ∞ norm? If so this would provide a way to establish existence for arbitrary time intervals. Do the limiting solutions generate geodesics for diffeomorphisms in 1D under the W 1 ∞ -norm? All of these intriguing questions will be the subject of future work.