On the numerical approximation of $p$-Biharmonic and $\infty$-Biharmonic functions

In [KP16] (arXiv:1605.07880) the authors introduced a second-order variational problem in $L^{\infty}$. The associated equation, coined the $\infty$-Bilaplacian, is a \emph{third order} fully nonlinear PDE given by $\Delta^2_\infty u\, := (\Delta u)^3 | D (\Delta u) |^2 = 0.$ In this work we build a numerical method aimed at quantifying the nature of solutions to this problem which we call $\infty$-Biharmonic functions. For fixed $p$ we design a mixed finite element scheme for the pre-limiting equation, the $p$-Bilaplacian $\Delta^2_p u\, := \Delta(| \Delta u |^{p-2} \Delta u) = 0.$ We prove convergence of the numerical solution to the weak solution of $\Delta^2_p u = 0$ and show that we are able to pass to the limit $p\to\infty$. We perform various tests aimed at understanding the nature of solutions of $\Delta^2_\infty u$ and in 1-$d$ we prove convergence of our discretisation to an appropriate weak solution concept of this problem, that of $\mathcal D$-solutions.