On the numerical approximation of vectorial absolute minimisers in L∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^\infty $$\end{document}

Let Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} be an open set. We consider the supremal functional 1E∞(u,O):=‖Du‖L∞(O),O⊆Ωopen,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {E}_\infty (u,{\mathcal {O}})\, {:}{=}\, \Vert \text {D}u \Vert _{L^\infty ( {\mathcal {O}} )}, \ \ \ {\mathcal {O}} \subseteq \Omega \text { open}, \end{aligned}$$\end{document}applied to locally Lipschitz mappings u:Rn⊇Ω⟶RN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u : \mathbb {R}^n \supseteq \Omega \longrightarrow \mathbb {R}^N$$\end{document}, where n,N∈N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n,N\in \mathbb {N}$$\end{document}. This is the model functional of Calculus of Variations in L∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^\infty $$\end{document}. The area is developing rapidly, but the vectorial case of N≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\ge 2$$\end{document} is still poorly understood. Due to the non-local nature of (1), usual minimisers are not truly optimal. The concept of so-called absolute minimisers is the primary contender in the direction of variational concepts. However, these cannot be obtained by direct minimisation and the question of their existence under prescribed boundary data is open when n,N≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n,N\ge 2$$\end{document}. We present numerical experiments aimed at understanding the behaviour of minimisers through a new technique involving p-concentration measures.


Introduction
Let n, N ∈ N be fixed integers and Ω an open set in R n . In this paper we perform numerical experiments and make some theoretical observations regarding appropriately defined minimisers of the functional The area of Calculus of Variations in L ∞ has a relatively short history in Analysis, pioneered by G. Aronsson in the 1960s, see e.g. [4][5][6]. Nonetheless, since its inception it has attracted the interest of many mathematicians due to both the theoretical importance (see e.g. [2,10,11,16,22,28,38] and the expository texts [7,15,31]) as well as due to the relevance to various and diverse applications from electrical breakdown [26] to image processing [21] to polycrystals [13] and from conformal mappings [14,29] to game theory [9,37].
A major difficulty in the study of (1.1) is that, quite surprisingly, global minimisers are not truly optimal objects for such pointwise functionals, in the sense that "they can be minimised even further" and in general may not solve any kind of Euler-Lagrange equations. There are many simple examples even in one space dimension where a minimiser does not minimise with respect to its own boundary conditions on subdomains [31][Chapter 1]. The latter property is of course automatic for integral functionals. In essence, this owes to the "non-local" nature of (1.1): although we can see (1.1) as the limit case of the p-Dirichlet functional as p → ∞ when passing to the extreme case of p = ∞, the σ-additivity properties of the Lebesgue integral are lost. This was already realised by Aronsson, who introduced the next concept.
Here, and throughout this exposition, O acts as a fixed variable ranging through the open subsets of Ω.
Although the area is developing rapidly, the vectorial case of N ≥ 2 is still poorly understood. In particular, to date neither the existence nor the uniqueness of absolute minimisers with prescribed boundary data on ∂Ω are known, unless min{n, N } = 1, namely either when n = 1 (curves in R N ) or N = 1 (scalar functions on R n ), even in more general situations when we have a general Hamiltonian H(., u, Du) instead of |Du|. The difficulty lies in the fact that for the "localise" variational concept of Definition 1, the standard NoDEA Numerics for vectorial absolute minimisers Page 3 of 23 51 approach of direct minimisation ( [18,24,25]) of the functional (1.1) in the affine space for a fixed u 0 ∈ W 1,∞ (Ω; R N ) does not in general yield absolute minimisers, but merely global minimisers which minimise on Ω only. Hence, in the absence of alternative methods, the main vectorial tool is to resort to L p approximations and then let p → ∞. To this end, let (u p ) p be the family of p-minimisers in W 1,∞ u0 (Ω; R N ). The problem encountered in the attempt to pass to the limit in is that the weak* sequential compactness is not strong enough to obtain that the limit is an absolute minimiser. This is because that, although p-Harmonic maps are C 1,α , this is not uniform in p. It can be shown, by the Bhattacharya-DiBenedetto-Manfredi estimate [8], they are bounded in W 1,∞ (O) uniformly in p. Hence, up to a subsequence (p j ) ∞ 1 , there exists a u ∞ such that u p * − − u ∞ as p j → ∞. This obstacle can be bypassed when min{n, N } = 1 because of either the scalar nature of the competing functions and comparison methods, or because of the one-dimensionality of their domain of definition (see for instance [1,10,32]), although the mode of convergence itself cannot in general be strengthened. In particular, the functional (1.1) is only weakly* lower semicontinuous and not weakly* continuous.
Herein we present and perform numerical experiments based on a new method curently being investigated, where various existence and uniqueness results have been obtained. This approach is motivated by the paper [22] and relates to developments in mass transportation (see e.g. [19,20,40], underpinning also the approach followed in the recent paper [35], despite in different guises since the higher order case treated therein is quite special. We refrain from giving any details now, and instead expound on the main ideas of this method in some detail in Sect. 2.
The numerical method we employ is a finite element approximation, based on an earlier work of Barrett and Liu on numerical methods for elliptic systems [12]. Therein the authors prove that, for a fixed exponent p, the method converges to the respective p-harmonic mapping under certain regularity assumptions on the solution. We would like to stress that significant care must be taken with numerical computations using this approach because the underlying nonlinear system is ill-conditioned. This owes to the nonlinearity of the problem which grows exponentially with p. Work to overcome this issue includes, for example, the work of Huang, Li and Liu [27] where preconditioners based on gradient descent algorithms are designed and shown to work well for large p.
Our mechanism for designing approximations of absolute minimisers is based on the approximation of the L ∞ minimisation problem by a respective L p one. We then utilise the approach described in [33,34,39] where it was shown that by forming an appropriate limit we are able to select candidates for 51 Page 4 of 23 N. Katzourakis and T. Pryer NoDEA numerical approximation along a "good" sequence of solutions, the p-harmonic mappings.
The purpose of this work is to demonstrate some key properties of vectorial absolute minimisers using an analytically justifiable numerical scheme which currently is the only technique available to give insight into the limiting vector-valued problem. We note that our goal is not to construct an efficient approximation method for vectorial absolute minimisers; indeed this indirect approximation of the limiting problem is not computationally efficient.
We conclude by noting that, albeit the concept of absolute minimisers is the primary contender in the direction of variational concepts, is not the unique candidate in the vector case of N ≥ 2. This is connected to that fact that, when N ≥ 2, the associated Euler-Lagrange ∞-Laplace system admits smooth non-minimising solutions which are characterised by two distinct variational concepts, see e.g. [3,30,36].

Approximating absolute minimisers through p-concentration measures
In this section we present the recently proposed method to construct absolute minimisers. For the clarity of this discussion suppose that u ∈ C 1 (Ω; R N ). The central point of this approach is to bypass the difficulties caused by the lack of continuity of the essential-supremum with respect to weak convergence by attempting to write the supremum as an integral but for a different quirky measure, namely where σ is a Radon probability measure on the closure O, which in principle is O-dependent as well as u-dependent. Of course an infinity of different families of measures {σ O : O Ω} can be defined to perform this role: for instance, for any O Ω choose any point x O maximising |Du| on O, namely any point and then choose the Dirac mass σ O :=δ xO . However, this condition by itself does not suffice for our purposes and has to be coupled by another condition. Firstly, in order to utilise duality arguments with ease, let us modify (2.1) to and suppose in addition that the matrix-valued Radon measure Du σ is divergence-free on O, namely that 1 div Du σ = 0. The PDE (2.3) is to be understood in a sense stronger than the usual distributional sense, due to the emergence of the boundary. In fact, (2.3) has to be interpreted as namely in the space of those C 1 test maps φ : O −→ R n which vanish on ∂O and extend continuously together with their derivative on ∂O (but whose derivative may not vanish on ∂O). Although the pair of conditions (2.2)-(2.3) seem to be ad-hoc and unmotivated, the payoff is that it is quite easy to see that if we have a C 1 map 3), then in fact u is an absolute minimiser in the sense of Aronsson (Definition 1.1). Indeed, fix a test function φ ∈ C 1 0 (O; R N ) and for simplicity suppose that ∂O has zero Lebesgue n-measure. Then, we have where the last equality follows from the fact that σ O is a probability measure on O, that is σ O (O) = 1. Hence, u is indeed minimising on O, at least among smooth variations. The general case then needs some approximation arguments, appropriately adapted to L ∞ . What is less clear is how we can actually prove the existence of such objects in W 1,∞ u0 (Ω; R N ) with given boundary conditions. The idea is to use L p approximations in the following way: for each p > n, consider the minimisation problem By standard variational arguments (see e.g. [18,24]), this problem has a solution u p , which is a weak solution to the celebrated p- By juxtaposing (2.7) with (2.3), we see that upon devising appropriate analytic tools in order to pass in an appropriate weak* sense to some limit . This is particularly challenging, as on the one hand mass might be lost towards the boundary ∂O and measures on O rather than on O have to be considered, even though by definition σ O p (∂O) = 0. On the other hand, the main cause of additional ramifications is the possible "oscillations" of the pair of weakly* converging objects (Du p , σ O p ) which requires particular "compensated compactness" mechanisms to pass to the limit and show that the oscillations occur in such a way that cancel each other.
It remains to show that, at least formally, we obtain (2.2) in the limit as p → ∞. To see this, note that for any α ∈ (0, 1) and for , directly from the definition of σ O p we have the estimate Since we can check with simple arguments that (with semi-continuity methods similar to e.g. in [31][Ch. 8]), we see that, at least formally Numerics for vectorial absolute minimisers Page 7 of 23 51 Equality (2.8) implies that the limit probability measure σ O ∞ is supported on the arg-max set Hence, (2.2) indeed follows.

Numerical approximations of p-concentration measures
In this section we perform several numerical experiments in both the scalar and the vectorial case for various appropriately selected boundary conditions. These experiments demonstrate in a concrete fashion the plethora of possible behaviours of the limit measures σ O ∞ , depending on the shape of O and its position in relation to the level sets of the modulus of the gradient |Du ∞ |.
The method we use is a conforming finite element discretisation of the p-Laplace system analysed in [12] for fixed p. We will describe the discretisation and summarise extensive numerical experiments aimed at quantifying the behaviour of the limit measure σ O ∞ . We refer to [33,34,39] for analytic justification.
We let T be an admissible triangulation of Ω, namely, T is a finite collection of sets such that (1) K ∈ T implies K is an open triangle, (2) for any K, J ∈ T we have that K ∩ J is either ∅, a vertex, an edge, or the whole of K and J and (3) ∪ K∈T K = Ω.
The shape regularity constant of T is defined as the number where ρ K is the radius of the largest ball contained inside K and h K is the diameter of K. An indexed family of triangulations {T n } n is called shape regular if μ := inf n μ(T n ) > 0. In what follows we shall assume that all triangulations are shape-regular and quasiuniform. We let P 1 (T ) denote the space of piecewise linear polynomials over the triangulation T , i.e., P 1 (T ) = φ such that φ| K ∈ P 1 (K) (3.4) and introduce the finite element space to be the usual space of continuous piecewise linear polynomial functions.

Galerkin discretisation
We consider the Galerkin discretisation, to find This is a conforming finite element discretisation of the vectorial p-Laplacian system proposed in [12]. Existence and uniqueness of solution to (3.6) follows from examination of the p-functional Notice that (3.7) is strictly convex and coercive on W 1,p 0 (Ω, R N ) so we may apply standard arguments from the Calculus of Variations showing that the minimisation problem is well posed. Hence, there exists a u ∈ W 1,p u0 (Ω, R N ) such that E p (u, Ω) = min the same argument applies for the Galerkin approximation. Further, for fixed p, let {U p } be the finite element approximation generated by solving (3.6) (indexed by h) and u p , the weak solution of the p-Laplace system, then we have that (3.9)

Numerical experiments
Our numerical approximations are achieved using Galerkin approximations to the p-Laplacian for sufficiently high values of p. We focus on studying the behaviour solutions have as p increases which allow us to make various conjectures on the behaviour of their asymptotic limit as p → ∞.
The computation of p-harmonic mappings is an extremely challenging problem in its own right. The class of nonlinearity in the problem results in the algebraic system, which ultimately yields the finite element solution, being extremely badly conditioned. One method to tackle this class of problems is by using preconditioners based on descent algorithms [27]. For extremely large NoDEA Numerics for vectorial absolute minimisers Page 9 of 23 51 p, say p ≥ 10000, this may be required; however for our purposes we restrict our attention to p ∼ 100. This yields sufficient accuracy for the results we are illustrating. We emphasise that even the case p ∼ 100 is computationally difficult to handle. The numerical approximation we are using is based on a damped Newton solver. As it is well known, Newton solvers require a sufficiently close initial guess in order to converge. A reasonable initial guess for the p-Laplacian is given by numerically approximating with the q-Laplacian for q < p sufficiently close to p. This leads to an iterative process in the generation of the initial guess, i.e., we solve the 2-Laplacian as an initial guess to the 3-Laplacian which serves as an initial guess to the 4-Laplacian, and so on. In each of our experiments the number of nonlinear iterations required to achieve a relative tolerance of 10 −8 was achieved in less that 20 iterations.
Using the methodology of p-approximation which we advocate here, it has been numerically demonstrated in the scalar case the rates of convergence both in p and in h that we expect to achieve [39][Section 4]. It was noticed that these rates were dependant on the regularity of the underlying ∞-harmonic function and that In both cases as h is decreased, an increasing value of p is required to achieve optimal approximation (in h). This suggests a coupling p = Ch α is necessary to achieve convergence, where the α is determined by the regularity expected in u ∞ . We found experimentally that coupling p = h −1/2 worked well for the singular case (u ∞ ∈ C 1,1/3 ) and p = h −1 for the smooth case (u ∞ ∈ C ∞ ).
We mention that we do not have access to exact solutions in general. For the experiments in Sects. 3.3-3.5 we only provide boundary data and in principle the solutions to this problem are non-smooth and must be interpreted in an appropriate weak sense. For the scalar problem this is through the viscosity solution framework [17].

Scalar case: Aronsson solution
Here we examine the scalar problem with prescribed boundary data given by the celebrated Aronsson solution shape of the domain and its position in relation to the level sets of the gradient, but generally it is a sum of lower-dimensional concentration measures.

Scalar case: Eikonal solution
In this test we examine the scalar problem with prescribed boundary data given by the conic solution Notice this function is Eikonal. We study the measure σ O p (T ) with O taken such that 0 / ∈ O. Results are shown in Fig. 5 where we plot the domain O and σ O p (T ) for p = 2, 4, 10, 20, 100. Notice that, in contrast to the Aronsson case, the mass is not sent toward the boundary as p → ∞. In this case we conjecture that σ O is always absolutely continuous with respect to the Lebesgue measure L n O .

Vectorial case: Eikonal solution
In this test we examine the vectorial problem with prescribed boundary data given by u(x, y) = e ix − e iy , (3.14) where we use the notation e it = (cos t, sin t). Notice this function is Eikonal. We study the measure σ O p (T ) with O taken such that 0 / ∈ O. Results are shown in Fig. 6 where we plot the domain O and σ O p (T ) for p = 2, 4, 10, 20, 100. Notice that mass is not sent towards the boundary as p → ∞. We conjecture that again σ O is absolutely continuous with respect to the Lebesgue measure L n O .

Vectorial case: mixed boundary conditions
In this test we examine the vectorial problem with prescribed boundary data given by  Figure 2. As Fig. 1 Figure 6. A test to characterise σ O ∞ for the boundary data given by (3.14). We show the domain O and the measure σ O p (T ) for increasing values of p. In each of the subfigures b-f we plot 10 contours that are equally spaced between the minimum and maximum values of σ O p . Notice that as p increases these contours remain equally spaced indicating mass is not pushed to the boundary NoDEA Numerics for vectorial absolute minimisers Page 17 of 23 51 Figure 7. A test to characterise σ O ∞ for the boundary data given by (3.15) with λ = 1 2 . We show the domain O and the measure σ O p (T ) for increasing values of p. In each of the subfigures b-f we plot 10 contours that are equally spaced between the minimum and maximum values of σ O p . Notice that as p increases these contours become concentrated at the four points {x i } 3 i=0 . In this case we conjecture that σ O is given by the sum of four Dirac masses at the points {x i } 3 i=0 and otherwise it is absolutely continuous with respect to L n O , rescaled so that it is a probability on O