\ensuremath{\beta}\xspaceegin{document}

\maketitle
\ensuremath{\beta}\xspaceegin{abstract}
  We present a continuous finite element method for some examples of
  fully nonlinear elliptic equation.  A key tool is the discretisation
  proposed in Lakkis \& Pryer (2011) allowing us to work directly on
  the strong form of a linear PDE. An added benefit to making use of
  this discretisation method is that a \emph{recovered (finite
    element) Hessian} is a biproduct of the solution process. We build
  on the linear basis and ultimately construct two different
  methodologies for the solution of second order fully nonlinear PDEs.
  Benchmark numerical results illustrate the convergence properties of
  the scheme for some test problems as well as the Monge--Ampère
  equation and the Pucci equation.
\end{abstract}

      \section{Introduction}
      \label{sec:intro}
      
      Fully nonlinear PDEs arise in many areas, including differential
      geometry (the \MA equation), mass transportation (the
      Monge--Kantorovich problem), dynamic programming (the Bellman
      equation) and fluid dynamics (the geostrophic equations).  The
      computer approximation of the solutions of such equations is thus an
      important scientific task.  There are at least three main difficulties
      apparent to someone attempting to derive numerical methods for fully
      nonlinear equations: first, the strong nonlinearity on the highest
      order derivative which generally precludes a variational formulation,
      second, a fully nonlinear equation does not always admit a classical
      solution, even if the problem data is smooth, and the solution has to
      sought in a generalised sense (e.g., viscosity solutions), which is
      bound to slow down convergence rates, and third, a common problem in
      nonlinear solvers, the exact solution may not be unique and
      constraints, such as convexity requirements must be included in the
      constraints to ensure uniqueness.
      
      Regardless of the problems, the \emph{numerical approximation of fully
        nonlinear second order elliptic equations}, as described in
      \cite{CaffarelliCabre:1995}, have been the object of considerable recent
      research, particularly for the case of \MA of which
      \cite{OlikerPrussner:1988, LoeperRapetti:2005, DeanGlowinski:2006,
        FengNeilan:09:vanishing, Oberman:2008, Awanou:10, DavydovSaeed:12:techreport,
        BrennerGudiNeilanSung:11, Froese:11} are selected examples.
      
      For more general classes of fully nonlinear equations some methods
      have been presented, most notably, at least from a theoretical view
      point, in \cite{Bohmer:2008} where the author presents a $\cont{1}$
      finite element method shows stability and consistency (hence
      convergence) of the scheme, following a classical ``finite
      difference'' approach outlined by \cite{Stetter:1973} which requires a
      high degree of smoothness on the exact solution. From a practical
      point of view this approach presents difficulties, in that the
      $\cont{1}$ finite elements are hard to design and complicated to
      implement, in \cite{DavydovSaeed:10} a useful overview of
      Bézier-Bernestein splines in two spatial dimensions is provided and a
      full implementation in \cite{DavydovSaeed:12:techreport}. Similar
      difficulties are encountered in finite difference methods and the
      concept of \emph{wide-stencil} appears to be useful, for example by
      \cite{KuoTrudinger:92,KuoTrudinger:05,Oberman:2008,Froese:11}.
      
      In \cite{FengNeilan:09:mixed,FengNeilan:09:vanishing,Awanou:10}
      the authors give a method in which they approximate the general second
      order fully nonlinear PDE by a sequence of fourth order quasilinear
      PDEs. These are quasilinear biharmonic equations which are discretised
      via mixed finite elements, or using high-regularity elements such as
      splines.  In fact for the \MA equation, which admits two solutions, of
      which one is convex and another concave, this method allows for the
      approximation of both solutions via the correct choice of a
      parameter. On the other hand although computationally less expensive
      than $\cont{1}$ finite elements (an alternative to mixed methods for
      solving the biharmonic problem), the mixed formulation still results
      in an extremely large algebraic system and the lack of maximum
      principle for general fourth order equations makes it hard to apply
      vanishing viscosity arguments to prove convergence.  A somewhat
      different approach, based on $C^0$-penalty, has been recently proposed
      by \cite{BrennerGudiNeilanSung:11}, as well as ``pseudo time'' one by
      \cite{Awanou:11}.
      
      It is worth citing also a \emph{least square} approach described by
      \cite{DeanGlowinski:2006}.  This method consists in minimising the
      mean-square of the residual, using a Lagrange multiplier method. Also
      here a fourth order elliptic term appears in the energy.
      
      In this paper, we depart from the above proposed methods and explore a
      more ``direct'' approach by applying the \emph{nonvariational finite
        element method}, introduced in \cite{LakkisPryer:2011}, as a solver
      for the Newton iteration directly derived from the PDE.  To be more
      specific, consider the following model problem
      \ensuremath{\beta}\xspaceegin{equation}
          \label{eq:modelproblem}
          \nlop{u} := \nlfunk(\ensuremath{\vec{H}}ess u) - f = 0
      \end{equation}
      with homogeneous Dirichlet boundary conditions where $\funk f\W\ensuremath{\rho}\xspaceeals$
      is prescribed function and $\funk F\symm\ensuremath{\rho}\xspaceeals$ is a real-valued
      algebraic function of symmetric matrixes, which provides an elliptic
      operator in the sense of \cite{CaffarelliCabre:1995}, as explained
      below in Definition \ensuremath{\rho}\xspaceef{def:ellipticity}.  The method we propose,
      consists in applying a Newton's method, given below by equation
      \eqref{eq:linearisedequation} of the PDE \eqref{eq:modelproblem},
      which results in a sequence of linear nonvariational elliptic PDEs
      that fall the framework of the nonvariational finite element method
      (NVFEM) proposed in \cite{LakkisPryer:2011}.  The results in this
      paper are computational, so despite not having a complete proof of
      convergence, we test our algorithm various problems that are
      specifically constructed to be well posed.  In particular, we test our
      method on the \MA problem, which is the de-facto benchmark for
      numerical methods of fully nonlinear elliptic equations.  This is in
      spite of \MA having an extra complication, which is conditional
      ellipiticity (the operator is elliptic only if the function is convex
      or concave.  A crucial, empirically observed feature of our method is
      that the convexity (or concavity) is automatically preserved if one
      uses $\poly2$ elements or higher.  For $\poly1$ elements this is not
      true and the scheme must be stabilized by reenforcing convexity (or
      concavity) at each timestep.  This was achieved in \cite{Pryer:2010}
      using a semidefinite programming method.  In a different spirit, but
      somewhat reminiscent, a stabilization procedure was obtained in
      \cite{BrennerGudiNeilanSung:11} by adding a penalty term.
      
      The rest of this paper is set out as follows. In
      \S\ensuremath{\rho}\xspaceef{sec:notation_and_discretisation} we introduce some notation,
      the model problem, discuss its ellipticity and Newton's method, which
      yields a sequences of nonvariational linearised PDE's. In
      \S\ensuremath{\rho}\xspaceef{sec:ndfem} we review of the nonvariational finite element
      method proposed in \cite{LakkisPryer:2011} and apply it to discretise
      the nonvariational linearised PDE's in Newton's method. In
      \S\ensuremath{\rho}\xspaceef{sec:unconstrained-fnl-pde} we numerically demonstrate the
      performance of our discretisation on a class of fully nonlinear PDE,
      those that are elliptic and well posed without constraining our
      solution to a certain class of functions.  In \S\ensuremath{\rho}\xspaceef{sec:ma} we turn
      to conditionally elliptic problems by dealing with the prime example
      of such problems, i.e., \MA.  We apply the discretisation to the \MA
      equation making use of the work \cite{Aguilera:2008} to check
      \emph{finite element convexity} is preserved at each iteration.
      Finally in \S\ensuremath{\rho}\xspaceef{sec:Pucci} we address the approximation of Pucci's
      equation, which is another important example of fully nonlinear
      elliptic equation.
      
      All the numerical experiments for this research, were carried out
      using the DOLFIN\xspace interface for \fenics \cite{LoggWells:2010} and
      making use of Gnuplot\xspace and \paraview for the graphics.

      \section{Notation}
      \label{sec:notation_and_discretisation}
      \subsection{Functional set-up}
      Let $\W\subset \R{d}$ be an open and bounded Lipschitz domain. We
      denote $\leb{2}(\W)$ to be the space of square (Lebesgue) integrable
      functions on $\W$ together with {its} inner product
      $\ltwop{v}{w} := \int_\W v w$ and norm $\Norm{v} :=
      \Norm{v}_{\leb2(\W)} = \ltwop{v}{v}^{1/2}$. 
      We denote by
      $\ensuremath{\,\mathrm{d}}uality{v}{w}$ the action of a distribution $v$ on the function
      $w$. 

      
      We use the convention that the derivative $\D u$ of a function
      $u:\W\to\ensuremath{\rho}\xspaceeals$ is a row vector, while the gradient of $u$, $\nabla u$
      is the derivatives transpose (an element of $\ensuremath{\rho}\xspaceeals^d$, representing
      $\D u$ in the canonical basis). Hence
      \ensuremath{\beta}\xspaceegin{equation}
        \nabla u 
        =
        \Transpose{\left(\D u\ensuremath{\rho}\xspaceight)}.
      \end{equation}
      For second derivatives, we follow the common innocuous abuse of
      notation whereby the Hessian of $u$ is denoted as $\ensuremath{\vec{H}}ess u$ (instead
      of the more consistent $\D\nabla  u$) and is represented by a $d\times d$
      matrix.
      
      The standard Sobolev spaces are ~\cite{Ciarlet:1978,Evans:1998}
      \ensuremath{\beta}\xspaceegin{gather}
        \sobh{k}(\W)
        := 
        \sob{k}{2}(\W) 
        = 
        \ensemble{\phi\in\leb2(\W)}
                 {\sum_{\norm{\vec\ensuremath{\alpha}\xspacelpha} \leq k}
                   \D^{\vec\ensuremath{\alpha}\xspacelpha}\phi\in\leb2(\W)},
                 \\
        \hoz := \text{closure of }\cont{\infty}_0(\W) \text{ in } \sobh{1}(\W)
      \end{gather}
      where $\vec\ensuremath{\alpha}\xspacelpha = \{ \ensuremath{\alpha}\xspacelpha_1,...,\ensuremath{\alpha}\xspacelpha_d\}$ is a
      multi-index, $\norm{\vec\ensuremath{\alpha}\xspacelpha} = \sum_{i=1}^d\ensuremath{\alpha}\xspacelpha_i$ and
      derivatives $\D^{\vec\ensuremath{\alpha}\xspacelpha}$ are understood in a weak sense. 
      
      We consider the case when the model problem \eqref{eq:modelproblem} is
      uniformly elliptic in the following sense.
      
      \ensuremath{\beta}\xspaceegin{Defn}[{ellipticity \cite{CaffarelliCabre:1995}}]
        \label{def:ellipticity}
        The operator $\nlop\cdot$ in Problem (\ensuremath{\rho}\xspaceef{eq:modelproblem}) is 
        called \emph{elliptic on $\cC\subseteq\symm$} if and only if
        for each $\geomat M\in\cC$ there exist $\Lambda\geq\lambda>0$, that
        may depend on $\geomat M$ such that
        \ensuremath{\beta}\xspaceegin{equation}
          \label{eqn:conditional-ellipticity-of-nl}
          \lambda \sup_{\norm{\geomat \xi} = 1}\norm{\geomat N\geomat \xi}
          \leq 
          \nlfunk(\geomat M+\geomat N) - \nlfunk(\geomat M) 
          \leq
          \La \sup_{\norm{\geomat \xi} = 1}\norm{\geomat N\geomat \xi}
          \Foreach \geomat N\in\symm.
        \end{equation}
        If the largest possible set $\cC$ for which
        (\ensuremath{\rho}\xspaceef{eqn:conditional-ellipticity-of-nl}) is satisfied is a proper
        subset of $\symm$ we say that $\nlop\cdot$ is \emph{conditionally elliptic}.
        
        The operator $\nlop\cdot$ in Problem (\ensuremath{\rho}\xspaceef{eq:modelproblem}) is called
        to be \emph{uniformly elliptic} if and only if for some
        $\lambda,\La>0$, called \emph{ellipticity constants}, we have
        \ensuremath{\beta}\xspaceegin{equation}
          \label{eq:ellipticity-of-nl}
          \lambda \sup_{\norm{\geomat \xi} = 1}\norm{\geomat N\geomat \xi}
          \leq 
          \nlfunk(\geomat M+\geomat N) - \nlfunk(\geomat M) 
          \leq
          \La \sup_{\norm{\geomat \xi} = 1}\norm{\geomat N\geomat \xi}
          \Foreach \geomat N,\geomat M\in\symm.
        \end{equation}
      \end{Defn}

      
      If $\nlfunk$ is differentiable (\ensuremath{\rho}\xspaceef{eq:ellipticity-of-nl}) can be
      obtained from conditions on the derivative of $\nlfunk$. A generic
      $\geomat M\in\symm$ is written as
      \ensuremath{\beta}\xspaceegin{equation}
        \geomat M =
        \squmatmd md
        ;
      \end{equation}
      so the derivative of $\nlfunk$ in the direction $\geomat N$ is given by
      \ensuremath{\beta}\xspaceegin{equation}
        \D\nlfunk(\geomat M)\geomat N
        =
        \frob{\nlfunk'(\geomat M)}{\geomat N}
      \end{equation}
      where the \emph{derivative matrix} $\nlfunk'(\geomat M)$ is defined by
      \ensuremath{\beta}\xspaceegin{equation}
        \nlfunk'(\geomat M)
        :=
        \squmatmd{\partial\nlfunk(\geomat M)/\partial m}{d}

        .
      \end{equation}

      

            Suppose $\nlfunk$ is differentiable. Then
            (\ensuremath{\rho}\xspaceef{eqn:conditional-ellipticity-of-nl}) is satisfied if and only if
            for each $\geomat M\in\cC$ there exists $\mu>0$ such that
            \ensuremath{\beta}\xspaceegin{equation} 
              \label{eq:ellipticity-for-derivative}
              \Transpose{\geovec\xi}\nlfunk'(\geomat M)\geovec\xi 
              \geq 
              \mu
              \norm{\geovec \xi}^2 
              \Foreach 

              \geovec\xi\in\R d.  
            \end{equation}
            Furthermore $\cC=\symm$ and $\mu$ is independent of $\geomat M$
            if and only if (\ensuremath{\rho}\xspaceef{eq:ellipticity-of-nl}) is satisfied.

      \ensuremath{\beta}\xspaceegin{Hyp}[smooth elliptic operator]
        \label{hyp:smooth-elliptic-operator}
        In the remainder of this paper we shall assume that
        $\nlop\cdot$ is conditionally elliptic on $\cC$ and
        \ensuremath{\beta}\xspaceegin{equation}
          \nlfunk \in \cont{1}(\cC).
        \end{equation}
        Unless otherwise stated we will also assume that $\cC=\symm$.
      \end{Hyp}

      \subsection{Newton's method}
      \label{sec:newtons-method}
      The smoothness assumption \ensuremath{\rho}\xspaceef{hyp:smooth-elliptic-operator} allows
      to apply Newton's method to solve Problem (\ensuremath{\rho}\xspaceef{eq:modelproblem}).
      
      Given the initial guess $u^0\in\cont2(\W)$, with $\ensuremath{\vec{H}}ess u^0\in\cC$, for
      each $n\in\NO$, find $u^{n+1}\in\cont2(\W)$ with $\ensuremath{\vec{H}}ess u^{n+1}\in\cC$ such that
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:newtonsmethod}
        \Dnlop{u^n}\left(u^{n+1} - u^n\ensuremath{\rho}\xspaceight) = -\nlop{u^n},
      \end{equation}
      where $\Dnlop u$ indicates the (Fréchet) derivative, which is formally
      given by
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:newtonsmethod2}
        \ensuremath{\beta}\xspaceegin{split}
          {\Dnlop u}{v}
          &=
          \lim_{\epsilon\ensuremath{\rho}\xspaceightarrow 0} 
          \frac{\nlop{u + \epsilon v} - \nlop{u}}{\epsilon}
          \\ 
          &= 
          \lim_{\epsilon\ensuremath{\rho}\xspaceightarrow 0}
          \frac{\nlfunk(\ensuremath{\vec{H}}ess u + \epsilon \ensuremath{\vec{H}}ess v) - \nlfunk(\ensuremath{\vec{H}}ess u)}{\epsilon}
          \\ 
          &=
          \nlfunk'(\ensuremath{\vec{H}}ess u) : \ensuremath{\vec{H}}ess v,
        \end{split}
      \end{equation}
      for each $v\in\cont2(\W)$.
      Combining (\ensuremath{\rho}\xspaceef{eq:newtonsmethod}) and (\ensuremath{\rho}\xspaceef{eq:newtonsmethod2}) then
      results in the following nonvariational sequence of linear PDEs. Given
      $u^0$ for each $n\in\naturals_0$ find $u^{n+1}$ such that
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:system-of-nonvar-pdes}
        \nlfunk'(\ensuremath{\vec{H}}ess u^n) : \ensuremath{\vec{H}}ess\left(u^{n+1}-u^n\ensuremath{\rho}\xspaceight) = f - \nlfunk(\ensuremath{\vec{H}}ess u^n).
      \end{equation}

      
      The PDE (\ensuremath{\rho}\xspaceef{eq:system-of-nonvar-pdes}) comes naturally in a nonvariational
      form. If we attempted to rewrite into a variational form, in order, say, to
      apply a ``standard'' Galerkin method, we would introduce an advection term
      which would depend on derivatives of $\nlfunk'$, \ie for generic $v,w$
      \ensuremath{\beta}\xspaceegin{equation}
        \frob{\nlfunk'(\ensuremath{\vec{H}}ess v)}{\ensuremath{\vec{H}}ess w}
        =
        \ensuremath{\,\mathrm{d}}iv\qb{\nlfunk'(\ensuremath{\vec{H}}ess v) \nabla w}- \ensuremath{\,\mathrm{d}}iv\qb{ \nlfunk'(\ensuremath{\vec{H}}ess v) }\nabla w.
      \end{equation}
      where the matrix-divergence is taken row-wise:
      \ensuremath{\beta}\xspaceegin{equation}
        \ensuremath{\,\mathrm{d}}iv\qb{\nlfunk'(\ensuremath{\vec{H}}ess v(\vec x))} 
        :=
        \qp{
          \sum_{i=1}^d \frac{\partial}{\partial x_i} 
          \qb{\matentry{F'}i1(\ensuremath{\vec{H}}ess v(\vec x))}
          , 
          \ensuremath{\,\mathrm{d}}otsc
          ,
          \sum_{i=1}^d \frac{\partial}{\partial x_i} 
          \qb{\matentry{F'}id(\ensuremath{\vec{H}}ess v(\vec x))}

        }
      \end{equation}
      and the chain rule provides us, for each $\ensuremath{\rho}\xspaceangefromto j1d$, with
      \ensuremath{\beta}\xspaceegin{equation}
          \sum_{i=1}^d \frac{\partial}{\partial x_i} 
          \qb{\matentry{F'}ij(\ensuremath{\vec{H}}ess v(\vec x))}
          =
          \sumifromto{k,l}1d
          \partial_{k,l}\matentry{F'}ij
          (\ensuremath{\vec{H}}ess v(\vec x))
          \partial_{ikl}v(\vec x).
      \end{equation}
      This procedure is undesirable for many reasons. Firstly it requires
      $F$ to be twice differentiable and it involves a third order
      derivative of the functions $u^{n+1}$ and $u^n$ appearing in
      (\ensuremath{\rho}\xspaceef{eq:newtonsmethod}).  Moreover, the ``variational'' reformulation
      could very well result in the problem becoming advection dominated and
      unstable for conforming FEM, as was manifested in numerical examples
      for the linear equation \cite[\S 4.2]{LakkisPryer:2011}.  In order to
      avoid these problems, we here propose the use of the nonvariational
      finite element method described next.

      \section[The NVFEM]{The nonvariational finite element method}
      \label{sec:ndfem}
      

      In \cite{LakkisPryer:2011} we proposed the \emph{nonvariational finite
        element method} (NVFEM) to approximate the solution of problems of
      the form (\ensuremath{\rho}\xspaceef{eq:system-of-nonvar-pdes}).  We review here the NVFEM
      and explain how to use it in combination with the Newton method to
      derive a practical Galerkin method for the numerical approximation of
      Problem (\ensuremath{\rho}\xspaceef{eq:modelproblem})'s solution.

      \subsection{Distributional form of (\ensuremath{\rho}\xspaceef{eq:system-of-nonvar-pdes}) and generalised Hessian}
      Let $\A \in \Le{\infty}(\W)^{d\times d}$ and for each $\geovec
      x\in\W$, let $\A(\geovec x)\in\symm$, the space of bounded, symmetric,
      positive definite, $d \times d$ matrices and $\funk f\W\ensuremath{\rho}\xspaceeals$. The
      \emph{Dirichlet linear nonvariational elliptic problem} associated
      with $\A$ and $f$ is
      \ensuremath{\beta}\xspaceegin{equation}
        \frob{\A}\ensuremath{\vec{H}}ess u=f
        \AND
        \ensuremath{\rho}\xspaceestriction u{\ensuremath{\beta}\xspaceoundary\W}=0.
      \end{equation}
      Testing this equation, and assuming $u\in\sobh2(\W)\cap\hoz$ such
      that $\ensuremath{\rho}\xspaceestriction{\nabla u}{\ensuremath{\beta}\xspaceoundary\W}\in\leb2(\ensuremath{\beta}\xspaceoundary\W)$,
      we may write it as
      \ensuremath{\beta}\xspaceegin{equation}
        \ensuremath{\beta}\xspaceegin{split}
          \label{Problem}
          \ltwop{\frob{\A}{\ensuremath{\vec{H}}ess u}}{\phi} 
          &= \ltwop{f}{\phi}
          \qquad \Foreach
          \phi\in\cont\infty_0(\W).
        \end{split}
      \end{equation}
      To allow a Galerkin type discretisation of (\ensuremath{\rho}\xspaceef{Problem}), we need to
      restrict the test functions $\phi$ to finite element function spaces
      that are generally \emph{not subspaces} of $\sobh2(\W)$.  So before
      restricting, we need to extend and we use a traditional
      distribution-theory (or generalised-functions) approach.

      Given a function $v\in\sobh2(\W)$ and let $\geovec
      n:\partial\W\to\ensuremath{\rho}\xspaceeals^d$ be the outward pointing normal of $\W$ then the
      Hessian of $v$, $\ensuremath{\vec{H}}ess v$ satisfies the following identity:
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:generalised-hessian}
        \ltwop{\ensuremath{\vec{H}}ess v}{\phi}
        = 
        -
        \int_\W \nabla v \otimes \nabla \phi
        +
        \int_{\partial\W} \nabla v \otimes \geovec n\,\phi 
        \Foreach \phi\in\sobh1(\W),
      \end{equation}
      where $\geovec a\otimes\geovec b:=\geovec a\transposevec b$ for
      $\vec a,\vec b$ column vectors in $\R d$.  If $v\in\sobh{1}(\W)$
      with $\ensuremath{\rho}\xspaceestriction{\nabla v}{\ensuremath{\beta}\xspaceoundary\W}\in\sobh{-1/2}(\partial\W)$ the
      right-hand side of (\ensuremath{\rho}\xspaceef{eq:generalised-hessian}) still makes sense
      and defines $\ensuremath{\vec{H}}ess v$ as an element in the dual of $\sobh1(\W)$ via

      \ensuremath{\beta}\xspaceegin{equation}

        \ensuremath{\,\mathrm{d}}uality{\ensuremath{\vec{H}}ess v}{\phi}
        :=
        -
        \int_\W \nabla v \otimes \nabla \phi
        +
        \int_{\partial\W} \nabla v \otimes \geovec n\,\phi
        \Foreach \phi\in\sobh1(\W)
        ,
      \end{equation}
      where $\ensuremath{\,\mathrm{d}}uality\cdot\cdot$ denotes the duality action on $\sobh1(\W)$
      from its dual.  We call $\ensuremath{\vec{H}}ess v$ the \emph{generalised Hessian} of
      $v$, and assuming that the coefficient tensor $\A$ is in
      $\cont0(\W)^{d\times d}$, for the product with a distribution to make
      sense, we now seek $u\in\sobhz1(\W)$ such that $\ensuremath{\rho}\xspaceestriction{\nabla
        u}{\ensuremath{\beta}\xspaceoundary\W}\in\sobh{-1/2}(\W)$ and whose generalised Hessian
      satisfies
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:linear-model}
        \ensuremath{\,\mathrm{d}}uality{\frob\A\ensuremath{\vec{H}}ess v}{\phi}
        =
        \ltwop f\phi
        \Foreach
        \phi\in\sobh1(\W).
      \end{equation}

      \ensuremath{\beta}\xspaceegin{comment}

        The mixed formulation of the model problem (\ensuremath{\rho}\xspaceef{Problem}) we consider is to seek
        the pair $(u, \ensuremath{\vec{H}}[u]) \in \sobh{1}(\W)\cap\sobh{1}(\partial\W)
        \times \leb{2}(\W)^{d\times d}$ such that
        \ensuremath{\beta}\xspaceegin{gather}
          \label{eq:mixed-1}
          \ltwop{\ensuremath{\vec{H}}[u]}{\phi} 
          +
          \int_\W{\nabla u}\otimes{\nabla \phi} 
          -
          \int_{\partial \W} {\nabla u}\otimes{\geovec{n} \ \phi}
          =
          \geomat 0
          \\
          \label{eq:mixed-2}
          \ltwop{\frob{\A}{\ensuremath{\vec{H}}[u]}}{\psi} 
          =
          \ltwop{f}{\psi}
          \Foreach (\phi,\psi)\in\sobh1(\W)\times \leb{2}(\W).
        \end{gather}

      \end{comment}

      \subsection{Finite element discretisation and finite element Hessian}
      We discretise (\ensuremath{\rho}\xspaceef{eq:linear-model}) for simplicity
      with a standard piecewise polynomial approximation for test and trial
      spaces for both problem variable, $U$, and auxiliary (mixed-type) variable,
      $\ensuremath{\vec{H}}[U]$. Let $\T{}$ be a conforming, shape regular
      triangulation of $\W$, namely, $\T{}$ is a finite family of sets such
      that
      \ensuremath{\beta}\xspaceegin{enumerate}
      \item $K\in\T{}$ implies $K$ is an open simplex (segment for $d=1$,
        triangle for $d=2$, tetrahedron for $d=3$),
      \item for any $K,J\in\T{}$ we have that $\closure K\meet\closure J$ is
        a full subsimplex (i.e., it is either $\emptyset$, a vertex, an
        edge, a face, or the whole of $\closure K$ and $\closure J$) of both
        $\closure K$ and $\closure J$ and
      \item $\union{K\in\T{}}\closure K=\closure\W$.
      \end{enumerate}
      We use the convention where $\funk h\W\ensuremath{\rho}\xspaceeals$ denotes the
      \emph{meshsize function} of $\T{}$, i.e.,
      \ensuremath{\beta}\xspaceegin{equation}
        h(\vec{x}):=\max_{\closure K\ni \vec x}h_K.
      \end{equation}
      We introduce the \emph{finite element spaces}
      \ensuremath{\beta}\xspaceegin{gather}
        \label{eqn:def:finite-element-space}
        \fes
        :=\ensemble{\Phi \in \sobh1(\W)}{\Phi\vert_{K} \in \poly p\Forall
          K\in\T{} \AND \Phi\in\cont{0}(\W)},
        \\
        \feszero
        :=\fes \cap \hoz,
      \end{gather}
      where $\poly k$ denotes the linear space of polynomials in $d$
      variables of degree no higher than a positive integer $k$. We
      consider $p\geq 1$ to be fixed and denote by $\Nzero :=
      \ensuremath{\,\mathrm{d}}im{\feszero}$ and $N := \ensuremath{\,\mathrm{d}}im{\fes}$. 

      The discretisation of problem then reads: Find $(U, \ensuremath{\vec{H}}[U]) \in
      \feszero \times \fes^{d\times d}$
      such that
      \ensuremath{\beta}\xspaceegin{equation}
        \ensuremath{\beta}\xspaceegin{gathered}
          \label{eq:discrete-linear}
          \ltwop{\ensuremath{\vec{H}}[U]}{\Phi} 
          = 
          -
          \int_\W{\nabla U}\otimes{\nabla \Phi} 
          +
          \int_{\partial \W} {\nabla U}\otimes{\geovec{n} \ \Phi}
          \Foreach \Phi\in\fes,
          \\
          \ltwop{\frob{\A}{\ensuremath{\vec{H}}[U]}}{\Psi} 
          =
          \ltwop{f}{\Psi}
          \Foreach \Psi\in\feszero.
        \end{gathered}
      \end{equation}
      For an algebraic formulation of (\ensuremath{\rho}\xspaceef{eq:discrete-linear}) we refer
      the reader to \cite[\S 2]{LakkisPryer:2011}.  Note that this
      discretisation can be interpreted as a mixed method whereby the
      first (matrix) equation defines the \emph{finite element Hessian} and the
      second (scalar) equation approximates the original PDE (\ensuremath{\rho}\xspaceef{Problem}).

      \subsection{Two discretisation stategies of \eqref{eq:modelproblem}}
        \label{sec:two-nonlinear-finite-element-methods}
        The finite element Hessian allows us two discretisation strategies.
        The first strategy, detailed in \S\ensuremath{\rho}\xspaceef{sec:unconstrained-fnl-pde},
        consists in applying Newton first to set-up
        (\ensuremath{\rho}\xspaceef{eq:system-of-nonvar-pdes}) and then using the NVFEM
        (\ensuremath{\rho}\xspaceef{eq:discrete-linear}) to solve each step.  A second strategy
        becomes possible, upon noting that given $U\in\fes$ the finite
        element Hessian $\ensuremath{\vec{H}}[U]$ is a regular function,\footnote{A
          generalised function $v$ is a \emph{regular function}, or just
          \emph{regular}, if it can be represented by a Lebesgue measurable
          function $f\in\lebloc1$ such that $\ensuremath{\,\mathrm{d}}uality v\phi=\int_\W f\phi$
          for all $\phi\in\cont\infty_0(\W)$.  We follow the customary and
          harmless abuse in identifying $v$ with $f$.}  which the
        generalised Hessian $\ensuremath{\vec{H}}ess U$ might fail to be.  This allows to
        apply nonlinear functions such as $F$ to $\ensuremath{\vec{H}}[U]$ and consider the
        following \emph{fully nonlinear finite element method} (FNFEM)
        \ensuremath{\beta}\xspaceegin{equation}
        \ensuremath{\beta}\xspaceegin{gathered}
          \label{eq:fully-nonlinear-FEM}
          \ltwop{\ensuremath{\vec{H}}[U]}{\Phi} 
          = 
          -
          \int_\W{\nabla U}\otimes{\nabla \Phi} 
          +
          \int_{\partial \W} {\nabla U}\otimes{\geovec{n} \ \Phi}
          \Foreach \Phi\in\fes,
          \\
          \ltwop{F(\ensuremath{\vec{H}}[U])}\Psi
          =
          \ltwop f\Psi
          \Foreach \Psi\in\fez.
        \end{gathered}
        \end{equation}
        Of course, in order to solve the second equation, a
        finite-dimensional Newton method may be necessary (but this
        strategy leaves the door open for other nonlinear solvers, e.g.,
        fixed point iterations).  A finite element code based on this idea
        will be tested in \S\ensuremath{\rho}\xspaceef{sec:Pucci} to solve the Pucci equation.
      
        In summary the finite element Hessian allows both paths in the
        following diagram:
        \ensuremath{\beta}\xspaceegin{equation}
          \small
          \label{eqn:dia:Newton-discretizations}
          \ensuremath{\beta}\xspaceegin{tikzpicture}
            \path (0,3) node(X){fully nonlinear PDE \eqref{eq:modelproblem}};
            \path (6,3) node(N){nonvariational linear PDE's \eqref{eq:system-of-nonvar-pdes}};
            \path (0,0) node(F){fully nonlinear FE discretization \eqref{eq:fully-nonlinear-FEM}};
            \path (6,.5) node(D){\color adiscrete linear 1};
            \path (6,0) node(L){\color bdiscrete linear 2};
            \ensuremath{\,\mathrm{d}}raw[-stealth,color=a] (X)--(N) node[pos=0.5,above]{Newton};
            \ensuremath{\,\mathrm{d}}raw[-stealth,color=a] (N)--(D) node[pos=0.5,above,sloped]{NVFEM};
            \ensuremath{\,\mathrm{d}}raw[-stealth,color=b] (X)--(F) node[pos=0.5,above,sloped]{FNFEM};
            \ensuremath{\,\mathrm{d}}raw[-stealth,color=b] (F)--(L) node[pos=0.5,above]{Newton};
          \end{tikzpicture}
        \end{equation}
        Although the diagram in (\ensuremath{\rho}\xspaceef{eqn:dia:Newton-discretizations}) does
        not generally commute, if the function $F$ is algebraically
        accessible, then it is commutative.  By ``algebraically accessible''
        we mean a function that can be computed in a finite number of
        algebraic operations or inverses thereof.  In this paper, we use
        only algebraically accessible nonlinearities, but, in principle
        assuming derivatives are available, our methods could be extended to
        algebraically inaccessible nonlinearities, such as Bellman's (or
        Isaacs's) operators involving optimums over infinite families, e.g.,
        \ensuremath{\beta}\xspaceegin{equation}
          F(\geomat M):=\inf_{\ensuremath{\alpha}\xspacelpha\in\cA}\frob{\geomat L_{\ensuremath{\alpha}\xspacelpha}}{\geomat M}
          \quad
          \qpbig{\OR
          F(\geomat M):=\inf_{\ensuremath{\alpha}\xspacelpha\in\cA}\sup_{\ensuremath{\beta}\xspaceeta\in\cB}\frob{\geomat L_{\ensuremath{\alpha}\xspacelpha,\ensuremath{\beta}\xspaceeta}}{\geomat M}}
          ,
        \end{equation}
        where $\ensemble{\geomat L_{\ensuremath{\alpha}\xspacelpha}}{\ensuremath{\alpha}\xspacelpha\in\cA}$ 
        (or $\ensemble{\geomat L_{\ensuremath{\alpha}\xspacelpha,\ensuremath{\beta}\xspaceeta}}{(\ensuremath{\alpha}\xspacelpha,\ensuremath{\beta}\xspaceeta)\in\cA\times\cB}$)
        is a family of elliptic operators.

      \section{The discretisation of unconstrained fully nonlinear PDEs}
      \label{sec:unconstrained-fnl-pde}
      
      In this section we detail the application of the method reviewed in
      \S\ensuremath{\rho}\xspaceef{sec:ndfem} to the fully nonlinear model problem
      \eqref{eq:modelproblem}.  Many fully nonlinear elliptic PDEs must be
      constrained in order to admit a unique solution. For example the \MAD
      is elliptic and admits a unique solution in the cone of convex (or
      concave) functions when $f>0$ (or $f<0$, respectively).  Before we
      turn our attention to the more complicated constrained PDE's in
      \S\ensuremath{\rho}\xspaceef{sec:ma} and we illustrate the Newton--NVFEM method in the simplest light.
      In this section we study fully nonlinear PDEs which have no such
      constraint.

      \ensuremath{\beta}\xspaceegin{Hyp}[unconditionally elliptic linearisation]
        \label{ass:well-posed}
        We assume, in this section, that the Newton-step linearisation
        (\ensuremath{\rho}\xspaceef{eq:system-of-nonvar-pdes}) is elliptic.  For this assumption
        to hold, it is sufficient to assume uniform ellipticity, i.e.,
        (\ensuremath{\rho}\xspaceef{eq:ellipticity-for-derivative}) with $\cC=\symm$ and $\mu>0$
        independent of $\geomat M$.
      \end{Hyp}

      \subsection{The Newton-NVFEM method}
        \label{The:nlfem}
        Suppose we are given a BVP of the form, finding
        $u\in\sobh2(\W)\cap\hoz$ such that
        \ensuremath{\beta}\xspaceegin{equation}
          \label{eq:nonlinearproblem}
          \ensuremath{\beta}\xspaceegin{split}
            \nlop{u} = \nlfunk(\ensuremath{\vec{H}}ess u) - f &= 0 \qquad \text{ in } \W,     
          \end{split}
        \end{equation}
        which satisfies Assumption \ensuremath{\rho}\xspaceef{ass:well-posed}.
      
        Upon applying Newton's method to approximate the solution of problem
        (\ensuremath{\rho}\xspaceef{eq:nonlinearproblem}) we obtain a sequence of functions
        $(u^n)_{n\in\naturals_0}$ solving the following linear equations in
        nonvariational form,
        \ensuremath{\beta}\xspaceegin{equation}
          \label{eq:linearisedequation}
          \frob{\N(\ensuremath{\vec{H}}ess u^n)}{\ensuremath{\vec{H}}ess u^{n+1}} = g(\ensuremath{\vec{H}}ess u^n)
        \end{equation}
        where
        \ensuremath{\rho}\xspaceenewcommand{\ensuremath{\theta}\xspaceefootnote}{\fnsymbol{footnote}} 	
        \ensuremath{\beta}\xspaceegin{gather}
          \label{eq:newton_prob_data:N}
          \N(\geomat X) := F'(\geomat X),\\
          \label{eq:newton_prob_data:g}
          g(\geomat X) 
          :=
          f - \nlfunk(\geomat X) + \frob{F'(\geomat X)}{\geomat X}.
        \end{gather}
      
        The nonlinear finite element method to approximate
        (\ensuremath{\rho}\xspaceef{eq:linearisedequation}) is: given an initial guess $U^0 :=
        \Pi_0 u^0$ for each $n\in\naturals_0$ find $(U^{n+1},
        \ensuremath{\vec{H}}[U^{n+1}])\in \feszero\times\fes^{d\times d}$ such that
        \ensuremath{\beta}\xspaceegin{equation}
          \ensuremath{\beta}\xspaceegin{gathered}
            \label{eq:discretenewtonsmethod}
            \ltwop{\ensuremath{\vec{H}}[U^{n+1}]}{\Phi} 
            +
            \int_\W{\nabla U^{n+1}}\otimes{\nabla \Phi} 
            -
            \int_{\partial \W} {\nabla U^{n+1}}\otimes{\geovec{n}}\,\Phi
            =
            \geomat 0
            \Foreach \Phi\in\fes
            \\
            \AND
            \ltwop{\frob{\N(\ensuremath{\vec{H}}[U^{n}])}{\ensuremath{\vec{H}}[U^{n+1}]}}{\Psi} 
            =
            \ltwop{g(\ensuremath{\vec{H}}[U^n])}{\Psi}
            \Foreach \Psi\in\feszero.
          \end{gathered}
        \end{equation}

        \ensuremath{\beta}\xspaceegin{Obs}[initial conditions]
          
        \end{Obs}

      \ensuremath{\beta}\xspaceegin{comment}
      We now give an algorithm for the general method.
      \subsection{$\Algoname{The NVFEM for a class of fully nonlinear problems}$}
      \label{alg:fully-nonlinear}
      \ensuremath{\beta}\xspaceegin{algorithmic}
        \operatorname{re}quire $(\T{}_0, u^0, p, N, K_{\max}, \tol)$ 
        \Ensure $(U)$ the
        NVFE solution of (\ensuremath{\rho}\xspaceef{eq:discretenewtonsmethod})
        \State $k = 0$
        \While{$k \leq K_{\max}$} 
        \State $\fes_k = \Algoname{FE Space}(\T{}_k, p)$ 
        \If{$k = 0$}
        \State $U^0 = \laginterpol{\fes_0} u^0$
        \State $\ensuremath{\vec{H}}[U^0] = \Algoname{Hessian Recovery}(U^0,\fes_0)$
        \EndIf
        \State $n = 0$
        \While{$n \leq N$}
        \State $[U^{n+1}, \ensuremath{\vec{H}}[U^{n+1}]] = \Algoname{NVFEM}(\fes_k, \geomat N(\ensuremath{\vec{H}}[U^n]), g(\ensuremath{\vec{H}}[U^n]))$
        \If{$\Norm{U^{n+1}-U^n} \leq \tol$}
        \State break
        \EndIf
        \State $n = n+1$
        \EndWhile
        \State $\T{}_{k+1} = \Algoname{Global Refine}(\T{}_k)$
        \State $k = k+1$
        \EndWhile
      \end{algorithmic}
      \end{comment}

      

      

      \subsection{Numerical experiments: a simple example}

      In this section we detail numerical experiments aimed at demonstrating
      the application of \eqref{eq:discretenewtonsmethod} to a simple model
      problem.

      \ensuremath{\beta}\xspaceegin{Example}[a simple fully nonlinear PDE]
        \label{ex:fullynl-abs}
        The first example we consider is a fully nonlinear PDE with a very
        smooth nonlinearity. The problem is
        \ensuremath{\beta}\xspaceegin{equation}
            \ensuremath{\beta}\xspaceegin{split}
            \nlop{u} := \sin{\Delta u} + 2\Delta u -f &= 0 \text{ in } \W,
            \\
            u &= 0 \text{ on } \pd{}\W.
            \end{split}

        \end{equation}
        which is specifically constructed to be uniformly elliptic. Indeed
        \ensuremath{\beta}\xspaceegin{equation}
          \nlfunk'(\ensuremath{\vec{H}}ess u)
          = 
          \left(\cos{\Delta u} + 2\ensuremath{\rho}\xspaceight)\eye.
        \end{equation}
        which is uniformly positive definite.
        The Newton linearisation of the problem is then: Given $u^0$, for
        $n\in\naturals_0$ find $u^{n+1}$ such that
        \ensuremath{\beta}\xspaceegin{equation}
          \frob{\left(\cos{\Delta u^n} + 2\ensuremath{\rho}\xspaceight)\eye}
               {\ensuremath{\vec{H}}ess (u^{n+1} - u^n)}
               =
               f - \sin{\Delta u^n} - 2\Delta u^n.
        \end{equation}
        and our approximation scheme is nothing but
        \ensuremath{\rho}\xspaceef{eq:discretenewtonsmethod} with
        \ensuremath{\beta}\xspaceegin{gather}
          \N(\geomat X) = \left(\cos{\trace{\geomat X}} + 2\ensuremath{\rho}\xspaceight)\geomat I
          \\
          g(\geomat X) = f - \sin{\trace{\geomat X}} - 2\trace{\geomat X}.
        \end{gather}
        Figure \ensuremath{\rho}\xspaceef{Fig:nonlin-regular-abs-lap} details a numerical
        experiment on this problem when $d=2$ and when $\W =[-1,1]^2$ is a
        square which is triangulated using a criss-cross mesh.
      \end{Example}

      \ensuremath{\beta}\xspaceegin{Rem}[simplification of Example \ensuremath{\rho}\xspaceef{ex:fullynl-abs}]
        \label{rem:variational-simplification}
        Example \ensuremath{\rho}\xspaceef{ex:fullynl-abs} can be simplified considerably by noticing that
        \ensuremath{\beta}\xspaceegin{equation}
          \int_\W \tr{\ensuremath{\vec{H}}[U]} \Phi = \int_\W \Transpose{\qp{\nabla U}}\nabla \Phi
          \Foreach \Phi\in\feszero.
        \end{equation}
        This coincides with the definition of the \emph{discrete Laplacian}
        and makes the NVFEM coincide with the standard conforming FEM.  This
        observation applies to all fully nonlinear equations, with
        nonlinearity of the form \eqref{eq:modelproblem} with
        \ensuremath{\beta}\xspaceegin{equation}
          F(\geomat M):=a(\tr{\geomat M}),
        \end{equation}
        for some given $a$.  This class of problems, can be solved using a
        variational finite element method and can be used for comparison
        with our method. Note that in
        \cite{JensenSmears:12:techreport:Finite} the authors use this
        together with a localisation argument in order to prove convergence
        of a finite element method for a specific class of
        Hamilton--Jacobi--Bellman equation.  Their method coincides with ours,
        for an appropriate choice of quadrature.
      \end{Rem}
      

      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption[]{\label{Fig:nonlin-regular-abs-lap} Numerical experiments
          for Example \ensuremath{\rho}\xspaceef{ex:fullynl-abs}. Choosing $f$ appropriately such
          that $u(\geovec x) = \exp\qp{-10\norm{\geovec x}^2}$. We use an
          initial guess $u^0 = 0$ and run the iterative procedure until
          $\Norm{U^{n+1} - U^n} \leq 10^{-8}$, setting $U := U^M$
          the final Newton iterate of the sequence. Here we are plotting
          log--log error plots together with experimental convergence rates
          for $\leb{2}(\W), \sobh1(\W)$ error functionals for the problem
          variable, $U$, and an $\leb{2}(\W)$ error functional for the
          auxiliary variable, $\ensuremath{\vec{H}}[U]$. Notice that there is a
          ``superconvergence'' of the auxiliary variable for both
          approximations.}

        \subfigure[][{Taking $\fes$ to be the space of
            piecewise linear functions on $\W$ ($p = 1$). Notice that $\Norm{u -
              U^M} = \Oh(h^2)$, $\norm{u - U^M}_1 = \Oh(h)$ and $\Norm{\ensuremath{\vec{H}}ess u
              - \ensuremath{\vec{H}}[U^M]} = \Oh(h^{0.5})$. }]{
          \label{Fig:sin-lap-p1}
          \includegraphics[scale=\figscale,width=0.47\figwidth]{Figures/fullynonlinear-abs-sin-linears}
        } 
 \subfigure[][{Taking $\fes$ to be the space of piecewise
            quadratic functions on $\W$ ($p = 2$). Notice that $\Norm{u -
              U^M} = \Oh(h^3)$, $\norm{u - U^M}_1 = \Oh(h^2)$ and
            $\Norm{\ensuremath{\vec{H}}ess u - \ensuremath{\vec{H}}[U^M]} = \Oh(h^{1.5})$}]{
          \label{Fig:sin-lap-p2}
          \includegraphics[scale=\figscale,width=0.47\figwidth]{Figures/fullynonlinear-abs-sin-quadratics}
        }

      \end{figure}

      \ensuremath{\beta}\xspaceegin{Example}[nonvariational example]
        \label{eq:fullynl-unsolvable}
        This is a simple example where the \emph{variational trick}
        mentioned in Remark \ensuremath{\rho}\xspaceef{rem:variational-simplification} cannot be applied. We fix
        $d=2$ and consider the problem
        \ensuremath{\beta}\xspaceegin{equation}
          \ensuremath{\beta}\xspaceegin{split}
            \nlop{u}:=
            \qp{\partial_{11} u}^3
            +
            \qp{\partial_{22} u}^3
            +
            \partial_{11} u
            +
            \partial_{22} u
            -
            f
            &=
            0 \text{ in } \W
            \\
            u &= 0 \text{ on } \partial\W.
          \end{split}
        \end{equation}
        The approximation scheme is then \eqref{eq:discretenewtonsmethod} with
        \ensuremath{\beta}\xspaceegin{gather}
          \N(\geomat X) :=
          \ensuremath{\beta}\xspaceegin{bmatrix}
            3\geomat X_{11}^2 + 1 & 0\\
            0 & 3\geomat X_{22}^2 +1
          \end{bmatrix}
          \\
          g(\geomat X)
          := 
          f + 2\qp{\geomat X_{11}^3 + \geomat X_{22}^3}.
        \end{gather}
        Figure \ensuremath{\rho}\xspaceef{Fig:nonlin-krylov} details a numerical experiment on
        this problem in the case $d=2$ and $\W = [-1,1]^2$ triangulated with
        a criss-cross mesh. A similar example is also studied in
        \cite[Ex 5.2]{DavydovSaeed:12:techreport} using B\"ohmers method.
      

      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption[]{\label{Fig:nonlin-krylov} Numerical experiments
          for Example \ensuremath{\rho}\xspaceef{eq:fullynl-unsolvable}. Choosing $f$ appropriately such
          that $u(\geovec x) = \exp\qp{-10\norm{\geovec x}^2}$. We use an
          initial guess $u^0 = 0$ and run the iterative procedure until
          $\Norm{U^{n+1} - U^n} \leq 10^{-8}$, setting $U := U^M$
          the final Newton iterate of the sequence. Here we are plotting
          log--log error plots together with experimental convergence rates
          for $\leb{2}(\W), \sobh1(\W)$ error functionals for the problem
          variable, $U$, and an $\leb{2}(\W)$ error functional for the
          auxiliary variable, $\ensuremath{\vec{H}}[U]$. Notice that there is a
          ``superconvergence'' of the auxiliary variable for both
          approximations.}

        \subfigure[][{Taking $\fes$ to be the space of
            piecewise linear functions on $\W$ ($p = 1$). Notice that $\Norm{u -
              U^M} = \Oh(h^2)$, $\norm{u - U^M}_1 = \Oh(h)$ and $\Norm{\ensuremath{\vec{H}}ess u
              - \ensuremath{\vec{H}}[U^M]} \ensuremath{\alpha}\xspacepprox \Oh(h^{1.5})$. }]{
          \label{Fig:sin-lap-p1}
          \includegraphics[scale=\figscale,width=0.47\figwidth]{Figures/fullynonlinear-krylov-exp-linears}
        } 
 \subfigure[][{Taking $\fes$ to be the space of piecewise
            quadratic functions on $\W$ ($p = 2$). Notice that $\Norm{u -
              U^M} = \Oh(h^3)$, $\norm{u - U^M}_1 = \Oh(h^2)$ and
            $\Norm{\ensuremath{\vec{H}}ess u - \ensuremath{\vec{H}}[U^M]} = \Oh(h^{1.3})$}]{
          \label{Fig:sin-lap-p2}
          \includegraphics[scale=\figscale,width=0.47\figwidth]{Figures/fullynonlinear-krylov-exp-quadratics}
        }

      \end{figure}

      
      \end{Example}
      

      

      \section{The \MAD problem}
      \label{sec:ma}
      
      In this section we propose a numerical method for the \MAD (MAD) problem
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:ma}
        \ensuremath{\beta}\xspaceegin{split}
          \ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}ess u}
          &=
          f \text{ in }\W
          \\
          u 
          &= 
          g \text{ on }\partial \W.
        \end{split}
      \end{equation}
      Our numerical experiments exhibit robustness of our method when
      computing (smooth) classical solutions of the MAD equation.  Most
      importantly we noted the following facts:
      \ensuremath{\beta}\xspaceegin{enumerate}[(i)\ ]
      \item 
        the use of $\poly p$ elements with $p\geq2$ is essential
        as $\poly1$ do not work,
      \item 
        the convexity of the Newton iterates is
        conserved throughout the computation, in a similar way to the
        observations in \cite{loeper-rapetti:05}, where the authors prove this
        convexity-conservation property.
      \end{enumerate}
      Our observations are purely
      empirical from computations, which leaves an interesting open problem
      of proving this property.

      

      
      \ensuremath{\beta}\xspaceegin{Rem}[the MAD problem fails to satisfy Assumption \ensuremath{\rho}\xspaceef{ass:well-posed}]
        To clarify Assumption \ensuremath{\rho}\xspaceef{ass:well-posed} for the MAD problem
        (\ensuremath{\rho}\xspaceef{eq:ma}), in view of the characteristic expansion of
        determinant if $\geomat X,\geomat Y \in \symm$
        \ensuremath{\beta}\xspaceegin{equation}
          \ensuremath{\,\mathrm{d}}et\qp{\geomat X + \epsilon \geomat Y} 
          =
          \ensuremath{\,\mathrm{d}}et{\geomat X} + \epsilon \frob{\cof\qp{\geomat X}}{\geomat Y} + \Oh(\epsilon^2),
        \end{equation}
        where $\cof{\geovec X}$ is the
        matrix of cofactors of $\geovec X$. Hence
        \ensuremath{\beta}\xspaceegin{equation}
          \label{eq:der-of-det-is-cof}
          \nlfunk'(\geomat X) = \cof{\geomat X}
          .
        \end{equation}
        This implies that the linearisation of MAD is only well posed if we
        restrict the class of functions we consider to those $u$ that
        satisfy
        \ensuremath{\beta}\xspaceegin{equation}
          \label{eqn:MAD-ellipticity-via-derivative}
          \Transpose{\geovec \xi} \cof \ensuremath{\vec{H}}ess u\, \geovec \xi 
          \geq
          \lambda \norm{\geovec\xi}^2 \Foreach \geovec\xi\in\ensuremath{\rho}\xspaceeals^d
        \end{equation}
        for some ($u$-dependent) $\lambda>0$.
        Note that (\ensuremath{\rho}\xspaceef{eqn:MAD-ellipticity-via-derivative}) is equivalent to the
        following two conditions as well
        \ensuremath{\beta}\xspaceegin{equation}
          \ensuremath{\beta}\xspaceegin{gathered}
            \Transpose{\geovec \xi} \ensuremath{\vec{H}}ess u\, \geovec \xi 
            \geq
            \lambda \norm{\geovec\xi}^2 \Foreach \geovec\xi\in\ensuremath{\rho}\xspaceeals^d
            \\
            $u$ \text{ is strictly convex. }
          \end{gathered}
        \end{equation}
      
        \cite{LoeperRapetti:2005} have shown that for the \emph{continuous}
        (infinite dimensional) Newton method described in
        \ensuremath{\rho}\xspaceef{sec:newtons-method}, given an strictly convex initial guess
        $u^0$, each iterate $u^n$ will be convex.  It is crucial that this
        property is preserved at the discrete level, as it guarantees the
        solvability of each iteration in the \emph{discretised} Newton
        method.  For this it the right notion of convexity turns out to be
        the \emph{finite element convexity} as developed in
        \cite{Aguilera:2008}.  In \cite{Pryer:2010}, an intricate method
        based on semidefinite programming provided a way to constrain the
        solution in the case of $\poly1$ elements.  Here we observe that
        the finite element convexity is automatically preserved, provided we
        use $\poly2$ or higher conforming elements.
      \end{Rem}

      \subsection{Newton's method applied to \MA}
      
      In view of (\ensuremath{\rho}\xspaceef{eq:der-of-det-is-cof}) it is clear that
      \ensuremath{\beta}\xspaceegin{equation}
        \Dnlop{u}v = \frob{\cof{\ensuremath{\vec{H}}ess u}}{\ensuremath{\vec{H}}ess v}.
      \end{equation} 
      Applying the methodology set out in \S\ensuremath{\rho}\xspaceef{sec:unconstrained-fnl-pde}
      we set
      \ensuremath{\beta}\xspaceegin{gather}
        \label{eq:maN}
        \N(\ensuremath{\vec{H}}ess u^n) = \cof{\ensuremath{\vec{H}}ess u^n},
        \\
        \label{eq:mag}
        g(\ensuremath{\vec{H}}ess u^n) = f - \ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}ess u^n} + \frob{\cof{\ensuremath{\vec{H}}ess
            u^n}}{\ensuremath{\vec{H}}ess u^n},
      \end{gather}
      
      \ensuremath{\beta}\xspaceegin{Rem}[{relating cofactors to determinants}]
        \label{rem:cof-to-det}
        For a generic (twice differentiable) function $v$ it holds
        that
        \ensuremath{\beta}\xspaceegin{equation}
          d \ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}ess v} = \frob{\cof{\ensuremath{\vec{H}}ess v}}{\ensuremath{\vec{H}}ess v}.
        \end{equation}
        Using this formulation we could construct a simple fixed point
        method for the \MA equation.
      \end{Rem}
      In view of Remark \ensuremath{\rho}\xspaceef{rem:cof-to-det} $g$ can be further
      simplified
      \ensuremath{\beta}\xspaceegin{equation}
        \ensuremath{\beta}\xspaceegin{split}
          g(\ensuremath{\vec{H}}ess u^n) &= f - \ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}ess u^n} + \frob{\cof{\ensuremath{\vec{H}}ess
              u^n}}{\ensuremath{\vec{H}}ess u^n}
          \\
          &= f + (d-1) \ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}ess u^n}.
        \end{split}
      \end{equation}
      Newton's method reads: Given $u^0$ for each $n\in\naturals_0$ find
      $u^{n+1}$ such that
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:linearised-ma}
        \frob{\N(\ensuremath{\vec{H}}ess u^n)}{\ensuremath{\vec{H}}ess u^{n+1}} 
        =
        g(\ensuremath{\vec{H}}ess u^n).
      \end{equation}
      
      \subsection{Numerical experiments}
      
      In this section we study the numerical behaviour of the scheme
      presented in Definition \ensuremath{\rho}\xspaceef{The:nlfem} applied to the MAD problem.
      
      We present a set of benchmark problems constructed from the problem
      data such that the solution to the \MA equation is known. We fix $\W$
      to be the square $S = [-1,1]^2$ or $[0,1]^2$ (specified in the
      problem) and test convergence rates of the discrete solution to the
      exact solution.
      
      Figures \ensuremath{\rho}\xspaceef{Fig:EOC-ma1}--\ensuremath{\rho}\xspaceef{Fig:EOC-ma2} details the various
      experiments and shows numerical convergence results for each of the
      problems studied as well as solution plots, it is worthy of note that
      each of the solutions seems to be convex, however this is not
      necessarily the case. They are all though \emph{finite element convex}
      \cite{Aguilera:2008}. In each of these cases the Dirichlet boundary
      values are not zero. The implementation of nontrivial boundary
      conditions is described in \cite[\S 3.6]{LakkisPryer:2011} or in more
      detail in \cite[\S 4.4]{Pryer:2010}.
      
      \ensuremath{\beta}\xspaceegin{Rem}[choosing the ``right'' initial guess]
        As with any Newton method we require a starting guess, not just for
        $U^0$ but also of $\ensuremath{\vec{H}}[U^0]$. Due to the mild nonlinearity with the
        previous example an initial guess of $U^0\equiv 0$ and $\ensuremath{\vec{H}}[U^0]
        \equiv \geomat 0$ was sufficient. The initial guess to the MAD
        problem must be more carefully sought. 
        
        Since we restrict our solution to the space of convex functions, it
        is prudent for the initial guess to also be convex. Moreover we must
        rule out constant and linear functions over $\W$, since the Hessian
        of these objects would be identically zero, destroying ellipticity
        on the initial Newton step. Hence we specify that the initial guess
        to (\ensuremath{\rho}\xspaceef{eq:linearised-ma}) must be strictly convex.  Rather than
        postprocessing the finite element Hessian from a initial project
        (although this is an option) to initialise the algorithm we solve a
        linear problem using the nonvariational finite element
        method. Following a trick, described in \cite{DeanGlowinski:2003},
        we chose $U^0$ to be the standard $\fes$-finite element
        approximation of $u^0$ such that
        \ensuremath{\beta}\xspaceegin{gather}
          \label{eq:laplace-initial-guess}
          \Delta u^0 = 2\sqrt{f} \text{ in } \W
          \\
          u^0 = g \text{ on } \partial\W.
        \end{gather}

      \end{Rem}
      
      \ensuremath{\beta}\xspaceegin{Rem}[degree of the FE space]
        \label{rem:degree-of-fe-space}
        In the previous example the lowest order convergent scheme was found
        by taking $\fes$ to be the space of piecewise linear functions
        ($p=1$). For the MAD problem we require a higher approximation
        power, hence we take $\fes$ to be the space of piecewise
        quadratic functions, \ie $p=2$. 
      
        Although the choice of $p=1$ gives a stable scheme, convergence is
        not achieved. This can be characterised by \cite[Thm
          3.6]{Aguilera:2008} that roughly says you require more
        approximation power than what piecewise linear functions provide to
        be able to approximate all convex functions. Compare with Figure \ensuremath{\rho}\xspaceef{Fig:EOC-ma4}.
      \end{Rem}
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption{ \label{Fig:EOC-ma1} Numerical results for the MAD problem
          on the square $S = [-1,1]^2$. We choose the problem data $f$ and
          $g$ appropriately such that the solution is the radially symmetric
          function $u(\geovec x) = \exp\qp{{\norm{\geovec x}^2}/{2}}$. We
          plot the finite element solution together with a log--log error
          plot for various error functionals as in Figure
          \ensuremath{\rho}\xspaceef{Fig:nonlin-regular-abs-lap}. Note for $p=2$ the $\leb{2}(\W)$
          error rate of convergence is suboptimal, this is in agreement with
          the numerical examples produced in \cite{BrennerGudiNeilanSung:2011}}
        \ensuremath{\beta}\xspaceegin{center}
              \subfigure[][{The FE approximation to the function $u(\geovec x) = \exp\qp{\frac{\norm{\geovec
                x}^2}{2}}$.}]{

        \includegraphics[scale=\figscale,width=0.47\figwidth]{Figures/ma-solution-exp}
      }  

          \subfigure[][{Log--log error plot for $\poly{2}$ Lagrange FEs.}]{
            \label{Fig:ma-radial}
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/EOC-MA-radial}
          } 
        \end{center}
      \end{figure}
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption{ \label{Fig:EOC-ma3} A FE--convexity test for the numerical example given in \ensuremath{\rho}\xspaceef{Fig:EOC-ma1}. We plot $\ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}[U]}$ together with the principle minor of $\ensuremath{\vec{H}}[U]$.}
        \ensuremath{\beta}\xspaceegin{center}
              \subfigure[][{The principal minor of $\ensuremath{\vec{H}}[U]$, an approximation to the Hessian of the function $u(\geovec x) = \exp\qp{\frac{\norm{\geovec x}^2}{2}}$.}]{
                \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/principle-minor-p2-exp.jpg}
      }  

          \subfigure[][{The determinant of $\ensuremath{\vec{H}}[U]$.}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/determinant-p2-exp.jpg}
          } 
        \end{center}
      \end{figure}
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption{ \label{Fig:EOC-ma4} Numerical results for the MAD problem
          on the square $S = [-1,1]^2$. We choose the problem data $f$ and
          $g$ appropriately such that the solution is the radially symmetric
          function $u(\geovec x) = \exp\qp{{\norm{\geovec x}^2}/{2}}$.}
        \ensuremath{\beta}\xspaceegin{center}
              \subfigure[][{The $\poly{1}$ FE approximation to the function $u(\geovec x) = \exp\qp{\frac{\norm{\geovec
                x}^2}{2}}$.}]{
                \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/p1-test-solution.jpg}
      }  

          \subfigure[][{The error $u-U$ plotted as a function over $\W$. Note the 
      FE approximation does not converge in this case, see Remark \ensuremath{\rho}\xspaceef{rem:degree-of-fe-space}}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/p1-test-error.jpg}
          } 
        \end{center}
      \end{figure}
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption{\label{Fig:EOC-ma2} Numerical results for the MAD problem
          on the square $S=[-1,1]^2$. Choosing $f$ and $g$ appropriately
          such that the solution is $u(\geovec x) = -\qp{2 - x_1^2 -
            x_2^2}^{1/2}$. Note the function has singular derivatives on the
          corners of $S$. We plot the finite element solution together with
          a log--log error plot for various error functionals as in Figure
          \ensuremath{\rho}\xspaceef{Fig:nonlin-regular-abs-lap}.}
        \ensuremath{\beta}\xspaceegin{center}
          \subfigure[][{The FE approximation to the function $u(\geovec x)
              = -\qp{2 - x_1^2 - x_2^2}^{1/2}$.}]{

        \includegraphics[scale=\figscale,width=0.47\figwidth]{Figures/ma-solution-ballie}
      } 
 
         \subfigure[][{Log--log error plot for $\poly{2}$ Lagrange FEs.}]{
            \label{Fig:ma-ballie}
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/EOC-MA-ballie}
          } 
      \end{center}
      \end{figure}
      
      \subsection{Nonclassical solutions}
      \label{sec:nonclassical-monge}
      The numerical examples given in Figures \ensuremath{\rho}\xspaceef{Fig:EOC-ma1}--\ensuremath{\rho}\xspaceef{Fig:EOC-ma2} both describe the numerical approximation of classical solutions to the MAD problem. In the case of Figure \ensuremath{\rho}\xspaceef{Fig:EOC-ma1} $u\in\cont{\infty}(\closure\W)$ whereas in Figure \ensuremath{\rho}\xspaceef{Fig:EOC-ma2} $u\in\cont{\infty}(\W)\cap\cont{0}(\closure{\W})$. We now take a moment to study less regular solutions, \ie viscocity solutions which are not classical. In this test we the solution
      \ensuremath{\beta}\xspaceegin{equation}
        u(\geovec x) = \norm{\geovec x}^{2\ensuremath{\alpha}\xspacelpha}
      \end{equation}
      for $\ensuremath{\alpha}\xspacelpha\in (1/2, 3/4)$. The solution $u(\geovec x) \notin \sobh{2}(\W)$. In Figures \ensuremath{\rho}\xspaceef{Fig:EOC-alpha-0-55}--\ensuremath{\rho}\xspaceef{Fig:EOC-alpha-0-7} we vary the value of $\ensuremath{\alpha}\xspacelpha$ and study the convergence properties of the method. We note that the method fails to find a solution for $\ensuremath{\alpha}\xspacelpha \leq 1/2$. Finally in Figure \ensuremath{\rho}\xspaceef{Fig:monge-adaptive} we conduct an adaptive experiment based on a gradient recovery aposteriori estimator. The recovery estimator we make use of is the Zienkiewicz--Zhu patch recovery technique see
      \cite{publication1}, \cite[\S 2.4]{Pryer:2010} or \cite[\S
        4]{Ainsworth:2000} for further details.
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption{\label{Fig:EOC-alpha-0-55} Numerical results for the MAD problem
          on the square $S=[-1,1]^2$. Choosing $f$ and $g$ appropriately
          such that the solution is $u(\geovec x) = \norm{\geovec x}^{2\ensuremath{\alpha}\xspacelpha}$, with $\ensuremath{\alpha}\xspacelpha=0.55$. Note the function is singular at the origin. We plot the finite element solution together with a log--log error plot for various error functionals as in Figure
          \ensuremath{\rho}\xspaceef{Fig:nonlin-regular-abs-lap}.}
        \ensuremath{\beta}\xspaceegin{center}
          \subfigure[][{The FE approximation to the function $u(\geovec x) =
              \norm{\geovec x}^{2\ensuremath{\alpha}\xspacelpha}$, with $\ensuremath{\alpha}\xspacelpha=0.55$.}]{

        \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/monge-singular-alpha-0-55.jpg}
      } 
 
         \subfigure[][{Log--log error plot for $\poly{2}$ Lagrange FEs.}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/monge-singular-alpha-0-55-convergence}
          } 
      \end{center}
      \end{figure}
      
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption{\label{Fig:EOC-alpha-0-6} Numerical results for the MAD problem
          on the square $S=[-1,1]^2$. Choosing $f$ and $g$ appropriately
          such that the solution is $u(\geovec x) = \norm{\geovec x}^{2\ensuremath{\alpha}\xspacelpha}$, with $\ensuremath{\alpha}\xspacelpha=0.6$. Note the function is singular at the origin. We plot the finite element solution together with a log--log error plot for various error functionals as in Figure
          \ensuremath{\rho}\xspaceef{Fig:nonlin-regular-abs-lap}.}
        \ensuremath{\beta}\xspaceegin{center}
          \subfigure[][{The FE approximation to the function $u(\geovec x) =
              \norm{\geovec x}^{2\ensuremath{\alpha}\xspacelpha}$, with $\ensuremath{\alpha}\xspacelpha=0.6$.}]{

        \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/monge-singular-alpha-0-6.jpg}
      } 
 
         \subfigure[][{Log--log error plot for $\poly{2}$ Lagrange FEs.}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/monge-singular-alpha-0-6-convergence}
          } 
      \end{center}
      \end{figure}
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption{\label{Fig:EOC-alpha-0-7} Numerical results for the MAD problem
          on the square $S=[-1,1]^2$. Choosing $f$ and $g$ appropriately
          such that the solution is $u(\geovec x) = \norm{\geovec x}^{2\ensuremath{\alpha}\xspacelpha}$, with $\ensuremath{\alpha}\xspacelpha=0.7$. Note the function is singular at the origin. We plot the finite element solution together with a log--log error plot for various error functionals as in Figure
          \ensuremath{\rho}\xspaceef{Fig:nonlin-regular-abs-lap}.}
        \ensuremath{\beta}\xspaceegin{center}
          \subfigure[][{The FE approximation to the function $u(\geovec x) =
              \norm{\geovec x}^{2\ensuremath{\alpha}\xspacelpha}$, with $\ensuremath{\alpha}\xspacelpha=0.7$.}]{

        \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/monge-singular-alpha-0-7.jpg}
      } 
 
         \subfigure[][{Log--log error plot for $\poly{2}$ Lagrange FEs.}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/monge-singular-alpha-0-7-convergence}
          } 
      \end{center}
      \end{figure}
      

      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption[]{Numerical results for a solution to the \MAD equation with
          $f$ and $g$ appropriately
          such that the solution is $u(\geovec x) = \norm{\geovec x}^{2\ensuremath{\alpha}\xspacelpha}$, with $\ensuremath{\alpha}\xspacelpha=0.55$.
          We choose $p=2$, and use an adaptive
          scheme based on Z--Z gradient recovery. The mesh is refined
          correctly about the origin. Note that when $\ensuremath{\,\mathrm{d}}im{\fes} = 20,420$ the adaptive solution achieves $\Norm{u - U^M} \ensuremath{\alpha}\xspacepprox 0.0078$, the uniform solution given in Figure \ensuremath{\rho}\xspaceef{Fig:EOC-alpha-0-55} satisfies $\Norm{u-U^M} \ensuremath{\alpha}\xspacepprox 0.18$ using the same number of degrees of freedom. Using the adaptive strategy both $\Norm{u-U^M}$ and $\norm{u - U^M}_1$ converge like $\Oh(N^{-1})$. {\label{Fig:monge-adaptive} }}
        \ensuremath{\beta}\xspaceegin{center}
          \subfigure[][{
              Adaptive mesh
          }]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/monge-singular-alpha-0-55-adaptive-mesh.jpg}
          }
          

          \subfigure[][{Log--log error plot for $\poly{2}$ Lagrange FEs.}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/monge-singular-alpha-0-55-adaptive-convergence.pdf}
          }
        \end{center}
      \end{figure}
      
      
      

      \section{Pucci's equation}
      \label{sec:Pucci}
      In this section we look to discretise the nonlinear problem, in this
      case Pucci's equation as a system of nonlinear equations.  Pucci's
      equation arises as a linear combination of Pucci's extremal
      operators. It can nevertheless be written in an algebraically
      accessible form, without the need to compute the eigenvalues.

      \ensuremath{\beta}\xspaceegin{Defn}[Pucci's extremal operators \cite{CaffarelliCabre:1995}]
        Let $\geomat N \in \symm$ and $\eig{\geomat N}$ be it's spectrum,
        then the extremal operators are
        \ensuremath{\beta}\xspaceegin{gather}
          \cM(\geomat N) := \sum_{\lambda_i \in \sigma(\geomat N)} \ensuremath{\alpha}\xspacelpha_i \lambda_i = 0,
        \end{gather}
        with $\ensuremath{\alpha}\xspacelpha_i\in\ensuremath{\rho}\xspaceeals$. The maximal (minimal) operator, commonly
        denoted $\cM^+$ ($\cM^-$), has coefficients that satisfy
        \ensuremath{\beta}\xspaceegin{equation}
          0 < \ensuremath{\alpha}\xspacelpha_1 \leq \ensuremath{\,\mathrm{d}}ots \leq \ensuremath{\alpha}\xspacelpha_n 
          \qquad
          \qp{ \ensuremath{\alpha}\xspacelpha_1 \geq \ensuremath{\,\mathrm{d}}ots \geq \ensuremath{\alpha}\xspacelpha_n > 0}
        \end{equation}
        respectively.
      \end{Defn}

      \subsection{The planar case and uniform ellipticity}
      In the case $d=2$ the normalised Pucci's equation reduces to finding $u$ such that
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:Pucci}
        \ensuremath{\alpha}\xspacelpha \lambda_2 + \lambda_1 = 0
      \end{equation}
      where $\geomat N := \ensuremath{\vec{H}}ess u$. Note that if $\ensuremath{\alpha}\xspacelpha = 1$
      (\ensuremath{\rho}\xspaceef{eq:Pucci}) reduces to the Poisson--Dirichlet problem. This can
      be easily seen when reformulating the problem as a second order PDE
      \cite{DeanGlowinski:2005}. Making use of the characteristic
      polynomial, we see
      \ensuremath{\beta}\xspaceegin{equation}
        \ensuremath{\beta}\xspaceegin{split}
          \lambda_i = \frac{\Delta u \pm \qp{\qp{\Delta u}^2 - 4\ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}ess u}}^{1/2}}{2} \qquad i = 1,2.
        \end{split}
      \end{equation}
      Thus Pucci's equation can be written as
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:pucci-pde}
        0 = \qp{\ensuremath{\alpha}\xspacelpha + 1}\Delta u + \qp{\ensuremath{\alpha}\xspacelpha - 1}\qp{\qp{\Delta u}^2 - 4\ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}ess u}}^{1/2},
      \end{equation}
      which is a nonlinear combination of \MA and Poisson problems.  However
      owing to the Laplacian terms, and unlike the \MAD problem, Pucci's
      equation is (unconditionally) uniformly elliptic for

        \ensuremath{\beta}\xspaceegin{equation}
          \qp{\tr{\geomat X}}^2 - 4\ensuremath{\,\mathrm{d}}et{\geomat X} 
          \geq
          0
          \Foreach
          \geomat X \in \ensuremath{\rho}\xspaceeals^{2\times 2}.
        \end{equation}

      
      
      The discrete problem we use is a direct approximation of
      (\ensuremath{\rho}\xspaceef{eq:pucci-pde}), we seek $\qp{U, \ensuremath{\vec{H}}[U]}$ such that 
      \ensuremath{\beta}\xspaceegin{gather}
          \int_\W \ensuremath{\beta}\xspaceigg(\qp{\ensuremath{\alpha}\xspacelpha + 1}\tr{\ensuremath{\vec{H}}[U]} + \qp{\ensuremath{\alpha}\xspacelpha - 1}\qp{\qp{\tr{\ensuremath{\vec{H}}[U]}}^2 - 4\ensuremath{\,\mathrm{d}}et{\ensuremath{\vec{H}}[U]}}^{1/2}\ensuremath{\beta}\xspaceigg)\Phi = 0
          \\
          \ltwop{\ensuremath{\vec{H}}[U]}{\Psi} 
          = 
          -
          \int_\W{\nabla U}\otimes{\nabla \Psi} 
          +
          \int_{\partial \W} {\nabla U}\otimes{\geovec{n} \ \Psi}
          \Foreach (\Phi,\Psi)\in\feszero \times \fes.
      \end{gather}
      
      The result is a nonlinear system of equations which was solved using a
      algebraic Newton method.
      
      \subsection{Numerical experiments}
      
      We conduct numerical experiments to be compared with those of
      \cite{DeanGlowinski:2005}. The first problem we consider is a
      classical solution of Pucci's equation (\ensuremath{\rho}\xspaceef{eq:Pucci}). Let $\geovec
      x = \Transpose{\qp{x, y}}$, then the function
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:Pucci-solution}
        u(\geovec x) = -\qp{\qp{\qp{x + 1}^2 + \qp{y+1}^2}^{\qp{1-\ensuremath{\alpha}\xspacelpha}/2}}
      \end{equation}
      solves Pucci's equation almost everywhere away from $(x, y) = (-1,-1)$
      with $g := u|_{\partial\W}$. Let $\T{}$ be an irregular triangulation
      of $\W = [-0.95, 1]^2$. In Figure \ensuremath{\rho}\xspaceef{Fig:pucci-converence} we detail
      a numerical experiement considering the case $\ensuremath{\alpha}\xspacelpha \in [2,5]$.
      
      \ensuremath{\beta}\xspaceegin{figure}
        \caption[]{Numerical results for a classical solution to Pucci's
          equation (\ensuremath{\rho}\xspaceef{eq:Pucci-solution}). As with the case of the MAD
          problem we choose $p=2$. We use a Newton method to solve the
          algebraic system until the residual of the problem (see
          \cite[c.f.]{Kelley:1995}) is less than $10^{-10}$ (which is
          overkill to minimise Newton error effects). We plot log--log error
          plots with experimental orders of convergence, for various norms
          and values of $\ensuremath{\alpha}\xspacelpha$.{\label{Fig:pucci-converence} }}
        \ensuremath{\beta}\xspaceegin{center}
          \subfigure[][{
          $\ensuremath{\alpha}\xspacelpha = 2$
          }]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-classical-alpha-2.pdf}
          }
          

          \subfigure[][{$\ensuremath{\alpha}\xspacelpha =3$}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-classical-alpha-3.pdf}
          }
          

          \subfigure[][{
              $\ensuremath{\alpha}\xspacelpha = 4$
          }]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-classical-alpha-4.pdf}
          }
          

          \subfigure[][{$\ensuremath{\alpha}\xspacelpha = 5$}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-classical-alpha-5.pdf}
          }
        \end{center}
      \end{figure}
      
      We also conduct a numerical experiment to be compared with
      \cite{Oberman:2008}. In this problem we consider a solution of Pucci's
      equation with a piecewise defined boundary. Let $\W = [-1,1]^2$ and
      the boundary data be given as
      \ensuremath{\beta}\xspaceegin{equation}
        \label{eq:pucci-pw-boundary}
        g(\geovec x):= 
        \ensuremath{\beta}\xspaceegin{cases}
          1 \qquad \text{ when } \ensuremath{\alpha}\xspacebs{x} \geq \frac{1}{2} \AND \ensuremath{\alpha}\xspacebs{y} \geq \frac{1}{2}
          \\
          0 \qquad \text{ otherwise.}
        \end{cases}
      \end{equation}
      Figure \ensuremath{\rho}\xspaceef{Fig:pucci-piecewise} details the numerical experiment on
      this problem with various values of $\ensuremath{\alpha}\xspacelpha$.
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption[]{Numerical results for a solution to Pucci's equation with
          a piecewise defined boundary condition
          (\ensuremath{\rho}\xspaceef{eq:pucci-pw-boundary}). We choose $p=2$, use a Newton method
          to solve the algebraic system until the residual of the problem is
          less than $10^{-10}$. We plot the solution for various values of
          $\ensuremath{\alpha}\xspacelpha$ as well as a cross section through the coordinate
          axis. Notice that the solution becomes extremely badly behaved as
          $\ensuremath{\alpha}\xspacelpha$ increases. {\label{Fig:pucci-piecewise} }}
        \ensuremath{\beta}\xspaceegin{center}
          \subfigure[][{
          $\ensuremath{\alpha}\xspacelpha = 2$
          }]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-pw-alpha-2.jpg}
          }
          

          \subfigure[][{$\ensuremath{\alpha}\xspacelpha =3$}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-pw-alpha-3.jpg}
          }
          

          \subfigure[][{
              $\ensuremath{\alpha}\xspacelpha = 4$
          }]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-pw-alpha-4.jpg}
          }
          

          \subfigure[][{$\ensuremath{\alpha}\xspacelpha = 5$}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-pw-alpha-5.jpg}
          }
          

          \subfigure[][{Cross section about $x = 0$ for each value of $\ensuremath{\alpha}\xspacelpha$.}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-pw-cross-section.jpg}
          }
        \end{center}
      \end{figure}
      
      Since the solution to the Pucci's equation with piecewise boundary
      (\ensuremath{\rho}\xspaceef{eq:pucci-pw-boundary}) is clearly singular near the
      discontinuities we have also conducted an adaptive experiment based on
      a gradient recovery aposteriori estimator (as in \S\ensuremath{\rho}\xspaceef{sec:nonclassical-monge}). As can be seen from Figure
      \ensuremath{\rho}\xspaceef{Fig:pucci-adaptive} we regain qualitively similar results using
      far fewer degrees of freedom.
      
      \ensuremath{\beta}\xspaceegin{figure}[h]
        \caption[]{Numerical results for a solution to Pucci's equation with
          a piecewise defined boundary condition
          (\ensuremath{\rho}\xspaceef{eq:pucci-pw-boundary}). We choose $p=2$, and use an adaptive
          scheme based on Z--Z gradient recovery. The mesh is refined
          correctly about the jumps on the
          boundary. {\label{Fig:pucci-adaptive} }}
        \ensuremath{\beta}\xspaceegin{center}
          \subfigure[][{
          Finite element solution for $\ensuremath{\alpha}\xspacelpha = 3$
          }]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-pw-adaptive-alpha-3.jpg}
          }
          

          \subfigure[][{Adaptive mesh}]{
            \includegraphics[scale=\figscale,width=0.47\figwidth]{fullynonlinear/Figures/pucci-pw-adaptive-mesh-alpha-3.jpg}
          }
        \end{center}
      \end{figure}
      
      
      

      \section{Conclusions and outlook}
      
      In this work we have proposed a novel numerical scheme for fully
      nonlinear and generic quasilinear PDEs. The scheme was based on a
      previous work for nonvariational PDEs (those given in nondivergence
      form) \cite{LakkisPryer:2011}.
      
      We have illustrated the application of the method for a simple, non
      physically motivated example, moving on to more interesting problems,
      that of the \MA equation and Pucci's equation. 
      
      For Pucci's equation we numerically showed convergence and conducted
      experiments which may be compared with previous numerical studies.
      
      We demonstrated for classical solutions to the \MA equation the method
      is robust again showing numerical convergence. For less regular
      viscocity solutions we have found that the method must be augmented
      with a penalty term in a similar light to
      \cite{BrennerGudiNeilanSung:2011}.
      
      We postulate that the method is better suited to a discontinuous
      Galerkin framework which is the subject of ongoing research.
      

\ensuremath{\beta}\xspaceibliographystyle{abbrvnat}
              \ensuremath{\beta}\xspaceegin{thebibliography}{30}
              \providecommand{\natexlab}[1]{#1}
              \providecommand{\url}[1]{\texttt{#1}}
              \expandafter\ifx\csname urlstyle\endcsname\ensuremath{\rho}\xspaceelax
                \providecommand{\ensuremath{\,\mathrm{d}}oi}[1]{doi: #1}\else
                \providecommand{\ensuremath{\,\mathrm{d}}oi}{doi: \ensuremath{\beta}\xspaceegingroup \urlstyle{rm}\Url}\fi
              
              \ensuremath{\beta}\xspaceibitem[Aguilera and Morin(2009)]{Aguilera:2008}
              N.~E. Aguilera and P.~Morin.
              \newblock On convex functions and the finite element method.
              \newblock \emph{SIAM J. Numer. Anal.}, 47\penalty0 (4):\penalty0 3139--3157,
                2009.
              \newblock ISSN 0036-1429.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1137/080720917}.
              \newblock URL \url{http://dx.doi.org/10.1137/080720917}.
              
              \ensuremath{\beta}\xspaceibitem[Ainsworth and Oden(2000)]{Ainsworth:2000}
              M.~Ainsworth and J.~T. Oden.
              \newblock \emph{A posteriori error estimation in finite element analysis}.
              \newblock Pure and Applied Mathematics (New York). Wiley-Interscience [John
                Wiley \& Sons], New York, 2000.
              \newblock ISBN 0-471-29411-X.
              
              \ensuremath{\beta}\xspaceibitem[Awanou(2010)]{Awanou:10}
              G.~Awanou.
              \newblock Spline element method for the {Monge--Amp{\`e}re} equation.
              \newblock online preprint, Northern Illinois University, 2010.
              \newblock URL \url{http://www.math.niu.edu/~awanou/SplineMonge7.pdf}.
              
              \ensuremath{\beta}\xspaceibitem[Awanou(2011)]{Awanou:11}
              G.~Awanou.
              \newblock Pseudo time continuation and time marching methods for
                {Monge--Amp{\`e}re} type equations.
              \newblock online preprint, January 2011.
              \newblock URL \url{http://www.math.niu.edu/~awanou/MongePseudo05.pdf}.
              
              \ensuremath{\beta}\xspaceibitem[B{\"o}hmer(2008)]{Bohmer:2008}
              K.~B{\"o}hmer.
              \newblock On finite element methods for fully nonlinear elliptic equations of
                second order.
              \newblock \emph{SIAM J. Numer. Anal.}, 46\penalty0 (3):\penalty0 1212--1249,
                2008.
              \newblock ISSN 0036-1429.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1137/040621740}.
              \newblock URL \url{http://dx.doi.org/10.1137/040621740}.
              
              \ensuremath{\beta}\xspaceibitem[Brenner et~al.(2011{\natexlab{a}})Brenner, Gudi, Neilan, and
                Sung]{BrennerGudiNeilanSung:11}
              S.~C. Brenner, T.~Gudi, M.~Neilan, and L.-y. Sung.
              \newblock {$C^1$} penalty methods for the fully nonlinear monge-amp{{\`e}}re
                equation.
              \newblock \emph{Math. Comp.}, 80:\penalty0 1979--1995, March
                2011{\natexlab{a}}.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1090/S0025-5718-2011-02487-7}.
              
              \ensuremath{\beta}\xspaceibitem[Brenner et~al.(2011{\natexlab{b}})Brenner, Gudi, Neilan, and
                Sung]{BrennerGudiNeilanSung:2011}
              S.~C. Brenner, T.~Gudi, M.~Neilan, and L.-y. Sung.
              \newblock {$C^0$} penalty methods for the fully nonlinear {M}onge-{A}mp\`ere
                equation.
              \newblock \emph{Math. Comp.}, 80\penalty0 (276):\penalty0 1979--1995,
                2011{\natexlab{b}}.
              \newblock ISSN 0025-5718.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1090/S0025-5718-2011-02487-7}.
              \newblock URL \url{http://dx.doi.org/10.1090/S0025-5718-2011-02487-7}.
              
              \ensuremath{\beta}\xspaceibitem[Caffarelli and Cabr{\'e}(1995)]{CaffarelliCabre:1995}
              L.~A. Caffarelli and X.~Cabr{\'e}.
              \newblock \emph{Fully nonlinear elliptic equations}, volume~43 of
                \emph{American Mathematical Society Colloquium Publications}.
              \newblock American Mathematical Society, Providence, RI, 1995.
              \newblock ISBN 0-8218-0437-5.
              
              \ensuremath{\beta}\xspaceibitem[Ciarlet(1978)]{Ciarlet:1978}
              P.~G. Ciarlet.
              \newblock \emph{The finite element method for elliptic problems}.
              \newblock North-Holland Publishing Co., Amsterdam, 1978.
              \newblock ISBN 0-444-85028-7.
              \newblock Studies in Mathematics and its Applications, Vol. 4.
              
              \ensuremath{\beta}\xspaceibitem[Davydov and Saeed(2010)]{DavydovSaeed:10}
              O.~Davydov and A.~Saeed.
              \newblock Stable splitting of bivariate splines spaces by bernstein-b{\'e}zier
                methods.
              \newblock online preprint, University of Strathclyde, Department of Mathematics
                and Statistics, University of Strathclyde, Glasgow, Scotland GB, 11 2010.
              \newblock URL
                \url{http://personal.strath.ac.uk/oleg.davydov/stable_splBB.html}.
              \newblock to appear on LNCS.
              
              \ensuremath{\beta}\xspaceibitem[Davydov and Saeed(2012)]{DavydovSaeed:12:techreport}
              O.~Davydov and A.~Saeed.
              \newblock Numerical solution of fully nonlinear elliptic equations by
                {B{\"o}hmer}'s method.
              \newblock Technical report, University of Strathclyde, 26 Richmond Street,
                Glasgow GB-G1 1XH, Scotland, UK, January 2012.
              \newblock URL
                \url{http://personal.strath.ac.uk/oleg.davydov/fully_nonlin.html}.
              
              \ensuremath{\beta}\xspaceibitem[Dean and Glowinski(2003)]{DeanGlowinski:2003}
              E.~J. Dean and R.~Glowinski.
              \newblock Numerical solution of the two-dimensional elliptic {M}onge-{A}mp\`ere
                equation with {D}irichlet boundary conditions: an augmented {L}agrangian
                approach.
              \newblock \emph{C. R. Math. Acad. Sci. Paris}, 336\penalty0 (9):\penalty0
                779--784, 2003.
              \newblock ISSN 1631-073X.
              
              \ensuremath{\beta}\xspaceibitem[Dean and Glowinski(2005)]{DeanGlowinski:2005}
              E.~J. Dean and R.~Glowinski.
              \newblock On the numerical solution of a two-dimensional {P}ucci's equation
                with {D}irichlet boundary conditions: a least-squares approach.
              \newblock \emph{C. R. Math. Acad. Sci. Paris}, 341\penalty0 (6):\penalty0
                375--380, 2005.
              \newblock ISSN 1631-073X.
              
              \ensuremath{\beta}\xspaceibitem[Dean and Glowinski(2006)]{DeanGlowinski:2006}
              E.~J. Dean and R.~Glowinski.
              \newblock Numerical methods for fully nonlinear elliptic equations of the
                {M}onge-{A}mp\`ere type.
              \newblock \emph{Comput. Methods Appl. Mech. Engrg.}, 195\penalty0
                (13-16):\penalty0 1344--1386, 2006.
              \newblock ISSN 0045-7825.
              
              \ensuremath{\beta}\xspaceibitem[Evans(1998)]{Evans:1998}
              L.~C. Evans.
              \newblock \emph{Partial differential equations}, volume~19 of \emph{Graduate
                Studies in Mathematics}.
              \newblock American Mathematical Society, Providence, RI, 1998.
              \newblock ISBN 0-8218-0772-2.
              
              \ensuremath{\beta}\xspaceibitem[Feng and Neilan(2009{\natexlab{a}})]{FengNeilan:09:mixed}
              X.~Feng and M.~Neilan.
              \newblock Mixed finite element methods for the fully nonlinear
                {M}onge-{A}mp{\`e}re equation based on the vanishing moment method.
              \newblock \emph{SIAM J. Numer. Anal.}, 47\penalty0 (2):\penalty0 1226--1250,
                2009{\natexlab{a}}.
              \newblock ISSN 0036-1429.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1137/070710378}.
              \newblock URL \url{http://dx.doi.org/10.1137/070710378}.
              
              \ensuremath{\beta}\xspaceibitem[Feng and Neilan(2009{\natexlab{b}})]{FengNeilan:09:vanishing}
              X.~Feng and M.~Neilan.
              \newblock Vanishing moment method and moment solutions for fully nonlinear
                second order partial differential equations.
              \newblock \emph{J. Sci. Comput.}, 38\penalty0 (1):\penalty0 74--98,
                2009{\natexlab{b}}.
              \newblock ISSN 0885-7474.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1007/s10915-008-9221-9}.
              \newblock URL \url{http://dx.doi.org/10.1007/s10915-008-9221-9}.
              
              \ensuremath{\beta}\xspaceibitem[Froese(2011)]{Froese:11}
              B.~D. Froese.
              \newblock A numerical method for the elliptic monge-amp{\`e}re equation with
                transport boundary conditions.
              \newblock Technical report, 01 2011.
              \newblock URL \url{http://arxiv.org/abs/1101.4981v1}.
              
              \ensuremath{\beta}\xspaceibitem[Jensen and Smears(2012)]{JensenSmears:12:techreport:Finite}
              M.~Jensen and I.~Smears.
              \newblock Finite element methods with artificial diffusion for
                hamilton-jacobi-bellman equations.
              \newblock Technical report, 01 2012.
              \newblock URL \url{http://arxiv.org/abs/1201.3581v2}.
              
              \ensuremath{\beta}\xspaceibitem[Kelley(1995)]{Kelley:1995}
              C.~T. Kelley.
              \newblock \emph{Iterative methods for linear and nonlinear equations},
                volume~16 of \emph{Frontiers in Applied Mathematics}.
              \newblock Society for Industrial and Applied Mathematics (SIAM), Philadelphia,
                PA, 1995.
              \newblock ISBN 0-89871-352-8.
              \newblock With separately available software.
              
              \ensuremath{\beta}\xspaceibitem[Kuo and Trudinger(1992)]{KuoTrudinger:92}
              H.~J. Kuo and N.~S. Trudinger.
              \newblock Discrete methods for fully nonlinear elliptic equations.
              \newblock \emph{SIAM J. Numer. Anal.}, 29\penalty0 (1):\penalty0 123--135,
                1992.
              \newblock ISSN 0036-1429.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1137/0729008}.
              \newblock URL \url{http://dx.doi.org/10.1137/0729008}.
              
              \ensuremath{\beta}\xspaceibitem[Kuo and Trudinger(2005)]{KuoTrudinger:05}
              H.-J. Kuo and N.~S. Trudinger.
              \newblock Estimates for solutions of fully nonlinear discrete schemes.
              \newblock In \emph{Trends in partial differential equations of mathematical
                physics}, volume~61 of \emph{Progr. Nonlinear Differential Equations Appl.},
                pages 275--282. Birkh{\"a}user, Basel, 2005.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1007/3-7643-7317-2_20}.
              \newblock URL \url{http://dx.doi.org/10.1007/3-7643-7317-2_20}.
              
              \ensuremath{\beta}\xspaceibitem[Lakkis and Pryer(2011)]{LakkisPryer:2011}
              O.~Lakkis and T.~Pryer.
              \newblock A finite element method for second order nonvariational elliptic
                problems.
              \newblock \emph{SIAM J. Sci. Comput.}, 33\penalty0 (2):\penalty0 786--801,
                2011.
              \newblock ISSN 1064-8275.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1137/100787672}.
              \newblock URL \url{http://dx.doi.org/10.1137/100787672}.
              
              \ensuremath{\beta}\xspaceibitem[Loeper and Rapetti(2005{\natexlab{a}})]{LoeperRapetti:2005}
              G.~Loeper and F.~Rapetti.
              \newblock Numerical solution of the {M}onge-{A}mp\`ere equation by a {N}ewton's
                algorithm.
              \newblock \emph{C. R. Math. Acad. Sci. Paris}, 340\penalty0 (4):\penalty0
                319--324, 2005{\natexlab{a}}.
              \newblock ISSN 1631-073X.
              
              \ensuremath{\beta}\xspaceibitem[Loeper and Rapetti(2005{\natexlab{b}})]{loeper-rapetti:05}
              G.~Loeper and F.~Rapetti.
              \newblock Numerical solution of the {M}onge-{A}mp\`ere equation by a {N}ewton's
                algorithm.
              \newblock \emph{C. R. Math. Acad. Sci. Paris}, 340\penalty0 (4):\penalty0
                319--324, 2005{\natexlab{b}}.
              \newblock ISSN 1631-073X.
              
              \ensuremath{\beta}\xspaceibitem[Logg and Wells(2010)]{LoggWells:2010}
              A.~Logg and G.~N. Wells.
              \newblock D{OLFIN}: automated finite element computing.
              \newblock \emph{ACM Trans. Math. Software}, 37\penalty0 (2):\penalty0 Art. 20,
                28, 2010.
              \newblock ISSN 0098-3500.
              \newblock \ensuremath{\,\mathrm{d}}oi{10.1145/1731022.1731030}.
              \newblock URL \url{http://dx.doi.org/10.1145/1731022.1731030}.
              
              \ensuremath{\beta}\xspaceibitem[Oberman(2008)]{Oberman:2008}
              A.~M. Oberman.
              \newblock Wide stencil finite difference schemes for the elliptic
                {M}onge-{A}mp\`ere equation and functions of the eigenvalues of the
                {H}essian.
              \newblock \emph{Discrete Contin. Dyn. Syst. Ser. B}, 10\penalty0 (1):\penalty0
                221--238, 2008.
              \newblock ISSN 1531-3492.
              
              \ensuremath{\beta}\xspaceibitem[Oliker and Prussner(1988)]{OlikerPrussner:1988}
              V.~I. Oliker and L.~D. Prussner.
              \newblock On the numerical solution of the equation {$(\partial\sp 2z/\partial
                x\sp 2)(\partial\sp 2z/\partial y\sp 2)-((\partial\sp 2z/\partial x\partial
                y))\sp 2=f$} and its discretizations. {I}.
              \newblock \emph{Numer. Math.}, 54\penalty0 (3):\penalty0 271--293, 1988.
              \newblock ISSN 0029-599X.
              
              \ensuremath{\beta}\xspaceibitem[Pryer(2010)]{Pryer:2010}
              T.~Pryer.
              \newblock Recovery methods for evolution and nonlinear problems.
              \newblock \emph{DPhil Thesis, University of Sussex}, 2010.
              
              \ensuremath{\beta}\xspaceibitem[Stetter(1973)]{Stetter:1973}
              H.~J. Stetter.
              \newblock \emph{Analysis of discretization methods for ordinary differential
                equations}.
              \newblock Springer-Verlag, New York, 1973.
              \newblock Springer Tracts in Natural Philosophy, Vol. 23.
              
              \end{thebibliography}

\end{document}