On the Numerical Approximation of $\infty$-Harmonic Mappings

Given a map $u : \Omega \subseteq \mathbb{R}^n \longrightarrow \mathbb{R}^N$, the $\infty$-Laplacian is the system \[ \label{1} \Delta_\infty u \, :=\, \Big(\text{D}u \otimes \text{D}u + |\text{D}u|^2 [\text{D}u]^\bot \! \otimes I \Big) : \text{D}^2 u\, = \, 0. \tag{1} \] \eqref{1} is the model system of vectorial Calculus of Variations in $L^\infty$ and arises as the"Euler-Lagrange"equation in relation to the supremal functional \[ \label{2} E_\infty(u,\Omega)\, :=\, \| \text{D}u \|_{L^\infty(\Omega)}. \tag{2} \] The scalar case of \eqref{1} has been introduced by Aronsson in the 1960s and by now is relatively classical and well understood. The general system \eqref{1} has been discovered and studied by the first author in a series of recent papers. Supremal functionals are fundamental for applications because they provide more realistic models as opposed to conventional integral models. Herein we provide numerical approximations of solutions to the Dirichlet problem when $n=2$ and $N=2,3$ for certain carefully selected boundary data on the unit square. Our experiments demonstrate interesting and unexpected phenomena occurring in $L^\infty$ and provide insights on the structure of general solutions and the natural separation to phases they present.


Introduction
Let n, N ∈ N and Ω an open set in R n . Given a smooth map u : Ω ⊆ R n −→ R N with components (u 1 , ..., u N ) , the gradient matrix map Du : Ω ⊆ R n −→ R N n is denoted by Du = (D i ..n and the hessian tensor D 2 u : Ω ⊆ R n −→ R N n 2 s is denoted by D 2 u = (D 2 ij u α ) α=1...N i,j=1...n . We use R N n and R N n 2 s to denote respectively the matrix space and the symmetric tensor space in which Du and D 2 u are valued. In this paper we are interested in the numerical approximation of solutions to the ∞-Laplacian which is defined on smooth maps as the following PDE system (1.1) ∆ ∞ u := Du ⊗ Du + |Du| 2 [Du] ⊥ ⊗ I : D 2 u = 0.
In the above, "⊗" denotes the usual tensor product, "|Du|" denotes the Euclidean norm of the gradient in the matrix space R N n , that is "I" is the identity matrix and "[Du] ⊥ " denotes the orthogonal projection on the subspace of R N which is the orthogonal complement of the range of the linear mapping Du(x) : R n −→ R N for fixed x ∈ Ω: (1.2) [Du] ⊥ := Proj (R(Du)) ⊥ .
The notation ":" symbolises a contraction with respect to 3 indices as above which extends the usual Euclidean inner product in R N n . In index form with respect to the standard bases, (1.1) reads N β=1 n i,j=1 D i u α D j u β + |Du| 2 [Du] ⊥ αβ δ ij D 2 ij u β = 0, α = 1, ..., N.
The so-called supremal functional (1.3) and the PDE system (1.1) are the archetypal model objects of vectorial Calculus of Variations in the space L ∞ . In the scalar case N = 1, the system (1.1) simplifies to the following single equation (1.4) Du ⊗ Du : D 2 u = n i,j=1 since the second term involving the orthogonal projection (1.2) vanishes identically. The field of Calculus of Variations in the space L ∞ has been pioneered by Aronsson in the 1960s [A1]- [A6] who was the first to derive and study (1.4) is relation to variational problems arising from (1.3). Since then, the field has undergone an explosion of interest within the PDE community, we refer e.g. to [B, C, K] and references therein. Until the early 2010s, all considerations where restricted to the scalar case of N = 1. The general vectorial case of (1.1) as well as the study of the associated PDE systems arising from more general first order functionals (1.5) E ∞ (u, Ω) := H(·, u, Du) L ∞ (Ω) has been pioneered by the first author in a series of recent papers [K1]- [K9]. Except for the intrinsic mathematical interest of the field, L ∞ functionals are very important for applications as well since we obtain more realistic variational models as opposed to their integral counterparts: for example in mechanics it is preferable to know to maximum tolerance of a specific material to distractive forces rather than the average tolerance.
A basic difficulty arising already in the scalar case is that (1.4) is degenerate elliptic and in non-divergence form and classical approaches to define and study weak solutions via integration-by-parts fail. In general the Dirichlet problem can not be solved in the class of smooth functions since Aronsson himself demonstrated the existence of singular solutions in [A6, A7] which however are minimising for the functional. The theory of viscosity solutions of Crandall-Ishii-Lions (see [CIL, C] and for a pedagogical introduction we refer to [K]) proved to be the appropriate framework to study scalar L ∞ variational problems.
In the vectorial case, N ≥ 2, singular solutions of (1.1) still appear (see [K1, K3]). A further difficulty associated to (1.1) which is not present when N = 1 is that the projection [Du] ⊥ may be discontinuous even for C ∞ maps u, whence the nonlinear operator ∆ ∞ of (1.2) may have discontinuous coefficients even when applied to C ∞ maps. This happens because the range of the gradient matrix Du(x) ∈ R N n may not be constant throughout the domain Ω. A simple example of C ∞ smooth solution to (1.1) from the unit rectangle of R 2 into R 2 is u(x, y) = e ix − e iy ( [K1]). The rank of Du is 2 away from the diagonal {x = y} because the partial derivative D x u, D y u are linearly independent but on the diagonal the rank jumps to 1 and (1.2) becomes discontinuous thereon. This is a general phenomenon and the solutions of (1.1) always come in phases of different rank values of the gradient with interfaces separating them (see [K2, K3]). The appropriate duality-free notion of generalised solution for (1.1) has very recently been proposed by the first author in [K7] and is based on the probabilistic representation of the hessian which does not exist classically via Young measures valued the compactification of the space R N n 2 s . In this setting, among other far-reaching results, we have been able to prove existence to the Dirichlet problem for (1.1) and for the equations arising from (1.5) for n = 1. Subsequent developments in the context of this theory appear in [K8]- [K11].
In the scalar case some numerical schemes have been proposed for the direct approximation of viscosity solutions of (1.4). In [O,EO,O1] Oberman uses techniques from Barles and Souganidis [BS] for the approximation of fully nonlinear PDEs to construct difference schemes for the ∞-Laplacian. See also [LP1] where the second author constructs an h-adaptive finite element scheme based on a residual error indicator and the method derived in [LP]. Herein we are concerned with numerical experiments which provide further insights on the understanding of ∞-Harmonic mappings. To the best of our knowledge, the experiments we perform have not been attempted in the vectorial case. We consider a set of five different boundary conditions on the unit rectangle Ω = (−1, 1) 2 ⊆ R 2 with target either R 2 or R 3 (Subsections 4.2-4.5). The method we follow is based on the approximation of the L ∞ system (1.1) by the respective L p Euler-Lagrange system, that is the p-Laplacian (1.6) ∆ p u := div |Du| p−2 Du = 0, for large p ∈ N. Our numerical scheme of choice is the finite element method. We utilise the approach described in [P] for the scalar version of (1.1) where it was shown that by forming an appropriate limit we are able to select candidates for numerical approximation along a "good" sequence of solutions, the p-Harmonic mappings. This approach has been analytically justified by the first author in the scalar case of (1.4) in [P]. Herein we justify its application to the full vectorial case of (1.1). 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 [BL]. 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 resulting nonlinear system is extremely badly 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 [HLL] where preconditioners based on gradient descent algorithms are designed and shown to work well for p up to 1000. We circumvent the need for such preconditioners by choosing our boundary data carefully such that |Du| ≈ 1 over the domain.
A word of caution: to be perfectly clear about this paper, the purpose of our work is to demonstrate some key properties of ∞-Harmonic mappings by using an analytically justifiable scheme which currently is the only technique available to give insight into the limiting vectorial problem. We note that our goal is not to construct an efficient approximation method for the ∞-Laplace system; indeed this indirect approximation of the ∞-Laplacian system by variational problems is not computationally efficient.
Our results exhibit very interesting phenomena arising in the p-limiting case. More specifically, as p increases the image of the solutions tends to "flatten" and they behave like minimal surfaces. If the boundary condition includes two components which have ranks equal to 1, 2 (Subsections 4.2, 4.3), then the solutions tend to achieve the maximum possible rank throughout the domain. Moreover, as p increases the angle between the 2 partial derivative vectors appears to approach a constant value throughout without any interfaces inside the domain. If we prescribe boundary data which have rank equal to 1 (Subsection 4.4), then as p increases the solutions tend to "break" and become of rank equal to 1 for p = ∞ without interfaces, while for all finite p there is a region whereon the rank of the gradient is 2 and nontrivial interfaces appear. However, in general interfaces may be formed and they may not be either smooth or with locally Euclidean topology: in subsections 4.6 and 4.5 we use as boundary data the explicit ∞-Harmonic maps constructed in [K2] whose interfaces are either rectangular or with triple junctions. Our results show that the p-Harmonic maps approach the explicit solutions as p increases, forming interfaces with these shapes.
We conclude this introduction by noting that in the paper [SS] the authors derived a different more singular multi-valued version of "∞-Laplacian" which describes optimal Lipschitz extensions. In our setting this amounts to changing in (1.3) from the Euclidean norm "| · |" we are using on R N n to the nonsmooth operator norm A = max |a|=1 |Aa|.

Basics on Vectorial Calculus of Variations in L ∞ and its Fundamental Equations
2.1. L p approximations as p → ∞ of the L ∞ equations. The nomenclature ∞-Laplacian of (1.1) owes to its very derivation as the limit of the p-Laplacian (1.6) as p → ∞. In addition, the respective functionals also approximate the L ∞ functional, if rescaled appropriately, in that for any W 1,∞ (Ω) function we have In the vectorial case N ≥ 2 the derivation of the full system (1.1) from (1.6) was first performed by the first author in [K1]. We recall here the formal derivation which is the scaffolding we employ for the numerical approximations to our solutions. Suppose u : Ω ⊆ R n −→ R N is a smooth map. We rewrite the p-Laplacian (1.6) is index form as By distributing derivatives, we have N β=1 n i,j=1 and we may normalise and contract the derivatives in the first summand to find In vector notation this means 3), we lose information and we formally obtain only the system Du ⊗ Du : D 2 u = 0 which is one of the components of (1.1). In the scalar case, however, this idea is correct and we obtain the full equation (1.4). In the general vectorial case, we have the information that the two summands of (2.3) are opposite and in particular |Du| 2 ∆u is tangential to the image u(Ω) ⊆ R N . In order to retain this information, for any fixed x ∈ Ω we split R N as the direct orthogonal sum of the range of Du(x) : R n −→ R N and of its complement and by recalling (1.2), we also set By utilising (2.4) and (1.2), we split the system (2.3) as follows Note now that the term in the first bracket is tangential since for each x ∈ Ω it is valued in R(Du(x)), while the second term is orthogonal to the first and is valued in R(Du(x)) ⊥ . Hence, the two summands are linearly independent. We choose to renormalise (2.5) by multiplying the second summand by p − 2. Then, after this normalisation we get By letting p → ∞ in (2.6) we obtain the full system (1.1). A byproduct of this derivation is that (1.1) actually can be decoupled to a pair of systems (tangential & normal) which are independent of one another: The above arguments require too much regularity and too strong estimates in order to make sense rigorously for classical C 2 solutions, but they are very instructive of the approach we follow. In the scalar case N = 1 the above method of studying the L ∞ equations by utilising the asymptotic limits of the L p equations as p → ∞ has proved to be very popular and successful. This idea which dates back to Aronsson has been effectively put into action by applying the theory of Viscosity Solutions to the ∞-Laplacian (see e.g. [C, K] and references therein) which almost by definition is very stable under limits. Further, in view of the uniqueness theorems available in the scalar case, all subsequential limits as p → ∞ give rise to a viscosity solution of the limit equation.
In the vectorial case N ≥ 2 the situation is much more complicated since there is no effective counterpart of Viscosity Solutions stable under limits which would allow to prove existence with elementary estimates. Motivated partly by the equations arising in L ∞ , the first author has recently proposed in [K7] a new duality-free theory of generalised solutions which applies to general fully nonlinear systems of any order. In particular, it allows to make sense of (1.1) in the appropriate regularity class of W 1,∞ (Ω, R N ) mappings. Among other results, in [K7] we proved existence of solution to the Dirichlet problem for (1.1) and in [K9] we proved variational characterisations of ∞-Harmonic maps in terms of the functional (1.3). The idea behind this new notion of so-called D-solutions is briefly explained in the next subsection.
However, as it has been proved in [K5] due to vectorial obstructions it is impossible to obtain uniqueness for the Dirichlet problem to (1.1) even in the class of C ∞ solutions and for n = N = 2. Evidence provided in [K5] and herein suggest that the method of L p approximations as p → ∞ provides a selection principle of "good" solution to (1.1) which has been conjectured in [K7] that it is unique.
We now state an existence result for asymptotic limits of p-Harmonic maps as p → ∞ needed later. The proof can be found e.g. in [P] and is a minor vectorial extension of standard results on limits of p-Harmonic functions as p → ∞. Since the proof utilises only arguments involving norms, the proof in the vectorial case is essentially the same as in the scalar case given in the book [K]. Let us also remind here the standard notion of weak solutions to the p-Laplacian: a map u ∈ W 1,p g (Ω, R N ) is weakly p-Harmonic when 2.2. Theorem (Existence of limits of p-Harmonic maps as p → ∞). Let {u p } ∞ p=1 denote a sequence of weak solutions to the p-Laplace system (1.6) with u p ∈ W 1,p g (Ω, R N ). Then, there exists a subsequence such that as p → ∞ that sequence converges uniformly to a mapping u ∞ ∈ W 1,∞ (Ω, R N ). Namely, We emphasise that the map u ∞ is a candidate ∞-Harmonic mapping, that is a generalised solution to (1.1). This is true in the scalar case N = 1 in the sense of Viscosity Solutions of Crandall-Ishii-Lions. In the vectorial case, it is true as well at least in the case of the ODE system arising from (1.5) when n = 1 in the sense of D-solutions of the first author (see [K8]). We conjecture to be true in the case of (1.1) as well, but this is not a consequence of the current results of [K7] since the method of the existence proof was based on an ad-hoc method (an analytic counterpart of Gromov's "Convex Integration" for a differential inclusion) rather than on p-Harmonic approximations. A complete proof of this conjecture, at least to date, eludes us, but recently we have made significant progress in this regard.
2.3. Generalised solutions of the L ∞ equations. For the sake of completeness we briefly motivate here the definition of generalised solutions to (1.1). Since we do not utilise it in an essential manner in this paper, we refrain from giving all the details which can be found in [K7] and [K8]- [K11]. The idea applies to general fully nonlinear systems of any order and allows for merely measurable solutions. It is based on the following observation: a map u : The equation (2.10) is a mere restatement of the system (1.1) where we just change the viewpoint and instead of considering the hessian as a classical map D 2 u : Ω ⊆ R n −→ R N n 2 s , we instead view it as a probability valued map given by the Dirac mass at the hessian: This change of viewpoint turns out to be very fruitful since by attaching one point and compactifying R N n 2 s , we have that for any perhaps nonsmooth map, there always exist limits of the difference quotients in the appropriate space of probability-valued maps which may not be the concentration measures δ D 2 u but instead more general probability valued maps D 2 u : Ω ⊆ R n −→ P(R N n 2 s ∪ {∞}) called diffuse hessians. More precisely, if D 1,h denotes the first difference quotient operator on R n , our generalised hessians are the subsequential limits of the form The respective notion which generalises (2.10) is called D-solutions to the ∞-Laplacian (1.1) and is the primary notion of generalised solution in this context for the vectorial case.
(b) The graph of the parametrisation function K of (II) used in 2.12.
(c) Phases and interfaces of the smooth ∞-Harmonic map given by the explicit formula 2.12 when K is as in (I).
(d) Phases and interfaces of the smooth ∞-Harmonic map given by the explicit formula 2.12 when K is as in (II).
The above examples show that in the vectorial case very complicated phenomena can arise even for smooth solutions. We will further examine these phenomena in Section 4 numerically by studying p-Harmonic mappings for increasing values of p using boundary data provided by 2.12 with choices of K as in (I), (II).
3. Numerical Approximations of ∞-Harmonic mappings for n = N = 2 and n = 2 < N = 3 In this section we describe the technique we use to approximate ∞-Harmonic mappings. The method we use is a conforming finite element discretisation of the p-Laplacian analysed in [BL] for fixed p. We will describe the discretisation and justify its application to the problem at hand by studying the behaviour as the meshsize parameter tends to zero and as p gets large. To that end we will also extend the results given in [P] to the vectorial case.
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 Further, we define h : Ω → R to be the piecewise constant meshsize function of T given by A mesh is called quasiuniform when there exists a positive constant C such that max x∈Ω h ≤ C min x∈Ω h.
In what follows we shall assume that all triangulations are shape-regular and quasiuniform. We let P k (T ) denote the space of piecewise polynomials of degree k over the triangulation T ,i.e., (3.4) P k (T ) = φ such that φ| K ∈ P k (K) and introduce the finite element space to be the usual space of continuous piecewise polynomial functions of degree k.
3.1. Definition (L 2 (Ω) projection operator). The L 2 (Ω) projection operator, P h : It is well known that this operator satisfies the following approximation properties for v ∈ W 1,p (Ω) 3.2. Galerkin discretisation. We consider the Galerkin discretisation of (2.2), to find U ∈ [V] N with U | ∂Ω = P h g such that (3.9) This is a conforming finite element discretisation of the vectorial p-Laplacian system proposed in [BL].
3.3. Proposition (existence and uniqueness of solution to (3.9)). There exists a unique solution of both the weak formulation (2.8) and the Galerkin approximation (3.9).
Proof . Existence and uniqueness of this problem follows from examination of the p-functional Notice that (3.10) 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 u ∈ W 1,p g (Ω, R N ) such that (3.11) E p (u, Ω) = min v∈W 1,p 0 (Ω,R N ) Noticing that the weak problem (2.8) is the weak Euler-Lagrange equation for (3.11) we have equivalence of (2.8) and (3.11) as such the weak formulation is also well posed. To see this for the Galerkin approximation notice that [V] N ⊂ W 1,p (Ω, R N ). As such the minimisation problem over [V] N is equivalent to (3.9) and the same argument applies as in the continuous case.
3.4. Theorem (convergence of the discrete scheme to weak solutions). For fixed p let {U p } be the finite element approximation generated by solving (3.9) (indexed by h) and u p , the weak solution of (2.8), then we have that Proof We begin by noting the discrete weak formulation (3.9) is equivalent to the minimisation problem: Using this, we immediately have In view of the stability of the L 2 projection in W 1,p (Ω) [CT] we have uniformly in h. Hence by weak compactness there exists a (weak) limit to the finite element sequence, which we will call u * . Due to the weak semicontinuity of E p (·, Ω) we have In addition, in view of the approximation properties of P h given in Definition 3.1 we have for any v ∈ C ∞ (Ω, R N ) that Using the fact that U is a discrete minimiser of (3.13) we have Now, as v was generic we may use density arguments and that u p was the unique minimser to conclude u * = u p , concluding the proof.
3.5. Theorem (convergence). Let U p be the Galerkin solution of (3.9) and u ∞ the candidate ∞-Harmonic mapping. Then, along a subsequence we have Proof The proof consists of combining Theorems 2.2 and 3.4 and noticing that

Numerical experiments
In this section we summarise extensive numerical experiments which focus on quantifying the structure of solutions to the ∞-Laplacian PDE system (1.1). This is 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 → ∞.

4.1.
Remark (practical computation of (3.9) for large 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 [HLL]. For extremely large p, say p ≥ 10000, this may be required; however for our purposes we restrict our attention to p ∼ 1000. This yields sufficient accuracy for the results we want to illustrate.
We emphasise that even the case p ∼ 1000 is computationally tough to handle. The numerical approximation we are using is based on a 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.
To tie into the explicit examples given in Figure 1 we are particularly interested in the rank of the solution. Except for the intrinsic interest, this relates directly to a deeper understanding of the solutions to ∞-Laplace system since the coefficients of (1.1) are discontinuous on the interfaces of the solution. We compute this by calculating det(DU ) and representing the areas Ω 2 , Ω 1 and S of (2.11) by plotting contours of the function det(DU ).

4.2.
Solutions (−1, +1) 2 ⊆ R 2 −→ R 2 with mixed rank-two and rank-one boundary data. In this test we construct approximations of solutions of (1.1) with mixed rank-two and rank-one boundary data. We take This gives us rank-one data in the quadrant x < 0 and y < 0 and rank-two data elsewhere. The results are illustrated in Figure 2. 4.3. Solutions (−1, +1) 2 ⊆ R 2 −→ R 3 with mixed rank-two and rank-one boundary data. In this experiment we examine an extension to the example given in Section 4.2 to the case N = 3. We study the image of the problem with boundary data given as 4.4. Solutions (−1, +1) 2 ⊆ R 2 −→ R 2 with rank-one boundary data. In this test we construct mappings with rank-one boundary data. We take where e x = (1, 0) and e y = (0, 1) denote the unit vectors in the x and y directions. Notice that g ∈ C 0 (∂Ω, R 2 ). We examine the image and rank of the numerical approximation for various values of p in Figures 4 and 5. 4.5. Solutions (−1, +1) 2 ⊆ R 2 −→ R 2 with boundary data an explicit ∞-Harmonic map with triple junction interface. In this test we construct boundary data which give rise to the example of an explicit smooth ∞-Harmonic mapping with a triple junction interface as illustrated in Figure 1c. We take    We also take    The numerical experiment is given in Figure 6. 4.6. Solutions (−1, +1) 2 ⊆ R 2 −→ R 2 with boundary data an explicit ∞-Harmonic map with rectangular interface. In this test we constuct boundary data which give rise to the example of an explicit The numerical experiment is given in Figure 7.       . An illustration of the rank of the solution to the vectorial p-Laplacian with the boundary conditions given in 2.12 for the rectangular interface illustrated in Figure 1c for various increasing values of p. Here we plot det(DU ) and associated contour lines. These are plotted at increments of 0.05. Notice as p increases, the structure illustrated in Figure  1d becomes more pronounced.