section0em2em subsection0em2em 11affiliationtext: University of Oxford, Andrew
Wiles Building Woodstock Rd, Oxford OX2 6GG, United Kingdom
<EMAIL_ADDRESS>
# Regularity for quasilinear vectorial elliptic systems through an iterative
scheme with numerical applications

Lukas Koch

###### Abstract

We consider an iterative procedure to solve quasilinear elliptic systems with
$p$-growth. The scheme was first considered by Koshelev in the quadratic case
$p=2$. We present numerical applications as well as applications to higher
regularity properties.

†† _2010 Mathematics Subject Classification:_ 35J57, 65N30
_Keywords._ Regularity, iterative Galerkin schemes, p-growth
_Acknowledgements._ L.K. was supported by the Engineering and Physical
Sciences Research Council [EP/L015811/1].

## 1 Introduction

Let $1<p<\infty$ and consider an open and bounded domain
$\Omega\subset\mathbb{R}^{n}$. This paper is concerned with problems of the
form: Seek $u\in W^{1,p}_{g}(\Omega)=g+W^{1,p}_{0}(\Omega,\mathbb{R}^{N})$
such that

$\displaystyle\textup{div}\,a(x,\textup{D}u)=$
$\displaystyle\,\textup{div}\,f\quad\text{ in }\Omega$ (1.1)

where ${a\colon\Omega\times\mathbb{R}^{N\times n}\to\mathbb{R}^{N\times n}}$
is a matrix-valued elliptic structure field satisfying controlled $p$-growth
conditions, $f\in L^{p^{\prime}}(\Omega,\mathbb{R}^{N\times n})$ and $g\in
W^{1,p}(\Omega,\mathbb{R}^{N})$. Here $p^{\prime}$ denotes the Hölder
conjugate of $p$ and (1.1) is to be understood in the sense of distributions.
To be precise, we make the following assumptions on $a$:

  1. (A1)[ass:A1]

For every $z\in\mathbb{R}^{N\times n}$, $a(\cdot,z)$ is measurable in $\Omega$
and, for almost every $x\in\Omega$, $a(x,\cdot)$ is continuously
differentiable in $\mathbb{R}^{N\times n}$.

  2. (A2)[ass:A2]

There are $\Lambda_{a}\geq\lambda_{a}>0$ and $\mu\geq 0$ such that, for almost
every $x\in\Omega$ and for all $z,\xi\in\mathbb{R}^{N\times n}$,

$\lambda_{a}(\mu^{2}+|z|^{2})^{\frac{p-2}{2}}|\xi|^{2}\leq\partial_{z}a(x,z)\xi\cdot\xi\leq\Lambda_{a}(\mu^{2}+|z|^{2})^{\frac{p-2}{2}}|\xi|^{2}.$

  3. (A3)[ass:A3]

There is $c>0$, such that, for all $x\in\Omega$, $z\in\mathbb{R}^{N\times n}$,
$|a(x,z)|\leq c\left(1+|z|^{p-1}\right)$.

For simplicity of presentation, in addition to
$\mathrm{\ref{ass:A1}-\ref{ass:A3}}$, we make the further assumption that

$\displaystyle\text{For almost every }x\in\Omega\text{ and every
}z\in\mathbb{R}^{N\times n},\,\partial_{z}a(x,z)\text{ is symmetric.}$ (A4)

######  Remark 1.1.

We have chosen to focus on systems of the form
$\textup{div}\,a(x,\textup{D}u)=\textup{div}\,f$, however with minor changes
to $\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ all the theory we develop applies
equally to systems of the form
$\textup{div}\,a(x,u,\textup{D}u)=\textup{div}\,f$ where
$a\colon\Omega\times\mathbb{R}^{N}\times\mathbb{R}^{N\times
n}\to\mathbb{R}^{N\times n}$. Further (A4) can be dropped if the results of
this paper are modified as indicated in Remark 3.3.

We will only concern ourselves with the vectorial case $N\geq 2$ as the theory
in the scalar case $N=1$ is different and much more can be said.

The study of elliptic systems of the form (1.1) is well established with a
very long list of important results. For an introduction we refer to [19] and
[21] as well as the classical [31]. It is well known that in general, without
stronger structural assumptions than ellipticity and controlled $p$-growth,
solutions to (1.1) will not be $C^{1,\alpha}$\- or even $C^{0,\alpha}$\-
regular. In fact, problems already arise in the simplest case when $p=2$. An
initial counterexample for systems with $a\equiv a(x,z)$ (with $n\geq 3$) was
developed in [12], where a quadratic system was given for which solutions are
not Hölder-continuous and, in fact, fail to be bounded. While the
$x$-dependence of the system in [12] is discontinuous, almost at the same time
an example was found of a higher-order elliptic operator with analytic
coefficients but discontinuous solutions [36]. We remark at this point that,
if $n=2$, solutions to (1.1) are Hölder continuous [39]. When $a\equiv
a(x,y,z)$ does indeed depend on $y$, it is possible to give examples where $a$
is analytic, but solutions to (1.1) are nowhere continuous, see [45, 27],
building on examples with analytic fields $a$ but discontinuous solutions in
[22]. In this set-up, an example of a system in optimal dimensions $n>2$,
$N=2$ with discontinuous solutions is given in [17]. However, the lack of
regularity occurs already without $x$-dependence of the coefficients. Even
when $a\equiv a(z)$ is quadratic and analytic, $C^{1,\alpha}$-regularity of
solutions need not hold, if $n>2$, as [40] shows.

Due to this lack of regularity, the general regularity theory of systems of
type (1.1) proceeds through notions of partial regularity, that is, regularity
outside of a context-dependent small relatively closed set. Since our focus
here lies on results holding everywhere in the domain, we refer the reader to
[19, 21, 38] for results and references in this direction.

Regularity results holding in the full domain are known only in a few special
cases. The radial Uhlenbeck structure $a\equiv a(\left|z\right|)z$, originally
developed in [48], is well known to guarantee full $C^{1,\alpha}$-regularity.
In [46], the result is generalised to fields of the form $g(z,z)z$, where
$g(\cdot,\cdot)$ is a symmetric, positive definite bilinear form. In
particular, this theory covers the model case of the $p$-Laplace operator
$a(x,y,z)=|z|^{p-2}z$. There now exists a vast literature concerning the
regularity theory of (1.1) if $a\equiv a(|z|)z$ and we refer the reader to
[20, 34, 35, 33, 15, 14] for further results and references.

A second direction of everywhere regularity results, in the case where
$a\equiv a(z)$, concerns the case where the modulus of continuity of
$\partial_{z}a(z)$ is not too large compared to the ellipticity constant of
$a(z)$. $C^{1,\alpha}_{\textup{loc}}$-regularity of solutions to (1.1) is
shown in this set-up in [9, 10, 11].

A further direction, closely related to the results of this paper, are results
of Cordes-Nirenberg type. For linear elliptic PDE, the classical Cordes-
Nirenberg results [41, 7] state that solutions of

$\displaystyle a_{ij}(x)\partial_{ij}u(x)=f(x)$ (1.2)

have interior $C^{1,\alpha}$ regularity, if $f\in L^{\infty}$ and
$|a_{ij}-\delta_{ij}|<\varepsilon$ for sufficiently small $\varepsilon$. That
is, $C^{1,\alpha}$-regularity is obtained when the coefficients $\\{a_{ij}\\}$
are sufficiently close to the identity matrix in an appropriate sense. Note
that the identity matrix is nothing but the coefficients corresponding to the
Laplacian. In [6] the result is extended to fully non-linear elliptic
equations with solutions in the sense of viscosity solutions. A result of a
similar spirit is [16], where local $BMO$-regularity for extremals of (a
suitable relaxed version of) $\textup{dist}_{K}^{2}(\textup{D}u)$ is shown.
Here $\textup{dist}_{K}(z)$ denotes the distance of $z\in\mathbb{R}^{N\times
n}$ to the compact set $K\subset\mathbb{R}^{N\times n}$.

We also comment on the case where $a(x,z)=a(x)z+b(x,z)z$ and $\|b(x,z)\|\leq
c(n)\lambda_{a}$. Here $a(x)$ is a positive definite matrix, with smallest
eigenvalue uniformly bounded below in $\Omega$ by $\lambda_{a}>0$ and $c(n)$
is a constant depending only on $n$. If $c(n)$ is sufficiently small, then
bounded solutions of (1.1) are $\alpha$-Hölder-continuous [49].

In this paper we obtain regularity results for fields $a$ that are close to a
suitable reference field $b$, see the start of Section 3 for a precise
definition of our notion of reference field. Suppose $a,b$ satisfy
LABEL:ass:A1 and LABEL:ass:A2, and let
$A(x,z)=\partial_{z}a(x,z),B(x,z)=\partial_{z}b(x,z)$. Denote by
$\lambda_{B^{-1}A}(x,z),\Lambda_{B^{-1}A}(x,z)$ the smallest and largest
eigenvalue of $B^{-1}(x,z)A(x,z)$, respectively. We note that then
$\lambda_{B^{-1}A}$ has a lower bound, and $\Lambda_{B^{-1}A}(x,z)$ has an
upper bound, that is independent of $x,z$. Denote these bounds by
$\lambda_{a,b}$ and $\Lambda_{a,b}$, respectively. We set

$\displaystyle
K_{a,b}=\frac{\Lambda_{a,b}-\lambda_{a,b}}{\lambda_{a,b}+\Lambda_{a,b}}$ (1.3)

We will operate in a regime where $K_{a,b}$ is sufficiently small and remark
that in all cases we are able to make the smallness assumption on $K_{a,b}$
explicit. In this regime we are able to transfer existence and regularity
results for solutions of the system
$\textup{div}\,b(x,\textup{D}u)=\textup{div}\,f$ to solutions of the system
$\textup{div}\,a(x,\textup{D}u)=\textup{div}\,f$.

The observation of [30] is that, if $a$ is a matrix-valued elliptic structure
field satisfying controlled quadratic growth, then considering the scheme

$\displaystyle\begin{cases}-\Delta u_{n+1}=-\Delta
u_{n}-\gamma\textup{div}\,a(x,\textup{D}u_{n})+\gamma\textup{div}\,f&\quad\text{
in }\Omega\\\ u_{n+1}=0&\quad\text{ on }\partial\Omega,\end{cases}$

with any choice of $u_{0}\in W^{1,2}_{g}(\Omega)$, gives a sequence
$\\{u_{n}\\}$ that converges in $W^{1,2}(\Omega)$ to a solution of (1.1). The
problem is to be understood in the sense of distributions. [30] further shows
convergence of the scheme in appropriate Morrey spaces, obtaining
$C^{0,\alpha}$\- and $C^{1,\alpha}$-regularity of solutions to (1.1) under
(sharp in the case of $C^{0,\alpha}$) smallness assumptions on $K$. We will
return to this point shortly.

We extend this result to matrix-valued elliptic structure fields with
$p$-growth.

######  Theorem 1.

Let $1<p<\infty$, $f\in L^{p^{\prime}}(\Omega,\mathbb{R}^{N\times n})$ and
$g\in W^{1,p}(\Omega,\mathbb{R}^{N})$. Suppose $a,b$ satisfy
$\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and (A4). Assume that the problem

$\textup{div}\,b(x,\textup{D}u)=\textup{div}\,F\text{ in }\Omega$

has a unique solution in $W^{1,p}_{g}(\Omega)$ for any choice of $F\in
L^{p^{\prime}}(\Omega)$. If

$\displaystyle\begin{cases}\dfrac{\Lambda_{b}}{\lambda_{b}}K_{a,b}\left(\dfrac{3+\sqrt{5}}{2}\right)^{\frac{p-2}{2p}}6^{p-2}(p-1)<1&\quad\text{
for }p\geq 2\\\\[20.0pt]
\dfrac{\Lambda_{b}}{\lambda_{b}}K_{a,b}\left(\dfrac{3-\sqrt{5}}{2}\right)^{\frac{(p-2)}{4}}\dfrac{2^{2-p}}{p-1}<1&\quad\text{
for }p\leq 2,\end{cases}$

then there is $\gamma>0$ such that, if ${u_{0}\in W^{1,p}_{g}(\Omega)}$, and
${u_{n+1}\in W^{1,p}_{g}(\Omega)}$ is inductively defined to be the weak
solution of

$\displaystyle\textup{div}\,b(x,\textup{D}u_{n+1})=\textup{div}\,(b(x,\textup{D}u_{n})-\gamma
a(x,\textup{D}u_{n})+\gamma f)\text{ in }\Omega$ (1.4)

in the class $W^{1,p}_{g}(\Omega)$, then $u_{n}\to u$ in $W^{1,p}(\Omega)$.
Here $u$ is the unique solution in $W^{1,p}_{g}(\Omega)$ of (1.1).

######  Remark 1.2.

We note that if $p=2$ and $b(x,z)=z$, we recover precisely the situation in
[30].

The iterative scheme (1.4) naturally lends itself as a numerical scheme to
solve (1.1) using finite element approximations. Take $g=0$ and suppose
$\Omega$ is a regular polytope. Consider a family of shape-regular
triangulations $\\{T_{h}\\}_{h\in(0,1]}$ of $\Omega$ of mesh-size $h$. Denote
by $X_{h}$ the space of continuous piecewise linear (subordinate to $T_{h}$)
functions. Take $u_{0}\in X_{h}$. Define $u_{n+1}^{h}\in X_{h}$ inductively as
a solution of the problem:

$\displaystyle\begin{cases}\int_{\Omega}b(x,\textup{D}u_{n+1}^{h})\cdot\textup{D}\phi\,\mathrm{d}x=\int_{\Omega}(b(x,\textup{D}u_{n}^{h})-\gamma
a(x,\textup{D}u_{n}^{h})+\gamma
f)\cdot\textup{D}\phi\,\mathrm{d}x\quad&\forall\phi\in X_{h}\\\
u_{n+1}^{h}=0&\text{ on }\partial\Omega\end{cases}$ (1.5)

Note that due to $\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and the theory of
monotone operators (1.5) is well-defined. We prove in Section 4 that under the
assumptions of Theorem 1 this scheme converges in $W^{1,p}(\Omega)$ to a
solution of the problem

$\displaystyle\begin{cases}\int_{\Omega}a(x,\textup{D}u_{h})\cdot\textup{D}\phi=\int_{\Omega}f\cdot\textup{D}\phi\,\mathrm{d}x\quad&\forall\phi\in
X_{h}\\\ u_{h}=0&\text{ on }\partial\Omega.\end{cases}$

Under the additional assumption that, if $F\in L^{p^{\prime}}(\Omega)$, then
there is $\alpha>0$ such that the equation
$\textup{div}\,b(x,\textup{D}u)=\textup{div}\,F$ has a solution $u$ in
$W^{1+\alpha,p}(\Omega)$, and, moreover, we have the following estimate for
some $c>0$:

$\displaystyle\|u\|_{W^{1+\alpha,p}(\Omega)}\leq
c\left(1+\|u\|_{W^{1,p}(\Omega)}+\|F\|_{L^{p^{\prime}}(\Omega)}^{1/(p-1)}\right),$
(1.6)

we show that $u_{h}\to u$ in $W^{1,p}(\Omega)$ as $h\to 0$ where $u\in
W^{1,p}_{g}(\Omega)$ is the solution of (1.1). We note that (1.6) is satisfied
by many fields of interest, e.g. the $p$-Laplacian $b(x,z)=|z|^{p-2}z$, [44].

If $p=2$ and $b(x,z)=z$, the proposed numerical scheme falls into the class of
iterative linearised Galerkin schemes studied in [26]. We refer to this
reference for an overview and further references regarding such schemes. An
advantage of the schemes studied in this paper, and of unified iterative
schemes of the form studied in [26] in general, is that it is possible to show
global convergence results at linear rate for a wide class of problems. This
should be compared with nonlinear Newton schemes for which convergence can
usually only be expected locally, albeit at quadratic rate. An example of a
scheme, fitting into the framework of unified iterative schemes, that has been
used to study some problems of $p$-Laplace type where $a(x,z)=a(x,|z|)z$ and
$a(x,|z|)$ has $p-2$-growth, is the so-called Kačanov iteration, originally
introduced in [28], and defined by solving

$\displaystyle\textup{div}\,a(x,|\textup{D}u_{n}|)\textup{D}u_{n+1}=\textup{div}\,f.$

Versions of this scheme have been studied in [50, 24, 18, 13, 25], but in
general have been restricted to the case $p\in[1,2]$. Note that the scheme
studied in this paper cannot be used to replace these schemes as it requires
the ability to solve (1.1) numerically for a reference field. However, given a
problem to which a known numerical scheme, such as the Kačanov iteration can
be applied, it enables to solve perturbations of the problem, to which the
known scheme cannot be applied directly. For example, in the case of the
Kačanov scheme, perturbations need not satisfy any structural assumption of
the form $a(x,z)=a(x,|z|)z$.

In [30] the iterative process is used to derive sharp conditions on the
dispersion $\Lambda_{a}/\lambda_{a}$ of elliptic systems with quadratic growth
that guarantee $C^{0,\alpha}$-regularity of solutions. An alternative proof
was suggested in [32]. Further, under additional assumptions, conditions for
$C^{1,\alpha}$-regularity of solutions are presented. In this spirit, we
obtain Cordes-Nirenberg type results with regards to Calderón-Zygmund
estimates for perturbations of fields of $p$-Laplace type. We recall a result
of [29]:

######  Theorem 2 (Theorem 1.6 in [29]).

Let $1<p<\infty$ and $\Omega\subset\mathbb{R}^{n}$ be a bounded domain with
$C^{1,\alpha}$-boundary for some $\alpha\in(0,1]$. Suppose that the
coefficients of $B$ belong to $VMO(\Omega)$ and that $F\in
L^{\frac{q}{p-1}}(\Omega)$ for some $q>p$. Then there exists a unique weak
solution $u\in W^{1,p}_{0}(\Omega)$ of
$\textup{div}\,b(x,\textup{D}u)=\textup{div}\,F$. Moreover, there is a
constant $C_{0}>0$ so that

$\int_{\Omega}|\textup{D}u|^{q}\,\mathrm{d}x\leq
C_{0}\int_{\Omega}|F|^{\frac{q}{p-1}}\,\mathrm{d}x.$

Here $C_{0}$ depends only on $n,p,q,\lambda,\Lambda$, the VMO data of $B$ and
$\partial\Omega$. In particular, this implies that $u\in W^{1,q}_{0}(\Omega)$.

We obtain the following extension:

######  Theorem 3.

Let $1<p<q<\infty$. Suppose $\Omega$ is a $C^{1,\alpha}$ domain for some
$0<\alpha\leq 1$. Take $f\in L^{q/(p-1)}(\Omega)$ and let $g=0$. Assume $b$
and $C_{0}$ are as in Theorem 2. Assume that $a$ satisfies assumptions
$\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and (A4) with $\mu=0$ and that the
assumptions of Theorem 1 are satisfied. If

$\dfrac{C_{0}^{(p-1)/q}\Lambda_{b}K_{a,b}}{(p-1)}<1,$

then the solution $u$ of (1.1) satisfies

$\displaystyle\|\textup{D}u\|_{L^{q}(\Omega)}\lesssim\left(1+\|f\|_{L^{q/(p-1)}}^{\frac{1}{p-1}}\right).$
(1.7)

Further we obtain a similar result with regards to weighted estimates for
perturbations of fields of $p$-Laplace type. We recall the following special
case of the main result of [42]:

######  Theorem 4 (c.f. Theorem 2.1 in [42]).

Let $1<p<q<\infty$ and let $w$ be an $A_{q/p}$ weight. Suppose $\Omega$ is a
$C^{1}$-domain. Suppose $B\in VMO(\Omega)$ and $B$ satisfies

$\displaystyle\lambda|\xi|^{2}\leq B(x)\xi\cdot\xi\leq\Lambda|\xi|^{2}.$

Then there exists a positive constant $C_{2}>0$ such that the following holds.
For a given vector field ${F\in L^{q(p-1)}_{w}(\mathbb{R}^{n})}$, there is a
unique weak solution ${u\in W^{1,p}_{0}(\Omega)}$ of
$\textup{div}\,b(x,\textup{D}u)=\textup{div}\,F$, satisfying $\textup{D}u\in
L^{q}_{w}(\Omega)$ with the estimate

$\displaystyle\|\textup{D}u\|_{L^{q}_{w}(\Omega)}\leq
C_{2}\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}.$

Here the constant $C_{2}$ depends only on $n,p,q,R,\lambda,\Lambda$ and
$[w]_{q/p}$.

We use the Koshelev iteration to obtain the following extension.

######  Theorem 5.

Let $1<p<q<\infty$ and $w\in A_{q/p}$. Assume $\Omega$ is a $C^{1}$ domain and
let $g=0$, $f\in L^{\frac{q}{p-1}}_{w}(\Omega)$. Suppose $b$ satisfies the
assumptions of Theorem 4. Let $C_{2}$ be the constant from Theorem 4. Suppose
$a,b$ satisfy the assumptions of Theorem 1. If

$\frac{C_{2}^{(p-1)/q}K_{a,b}\Lambda_{b}}{(p-1)}<1,$

then the solution $v$ of (1.1) satisfies the estimate

$\displaystyle\|\textup{D}v\|_{L^{q}_{w}(\Omega)}\leq
c(C,p,q)\left(1+\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}^{\frac{1}{p-1}}\right)$
(1.8)

######  Remark 1.3.

The result can easily be extended to Reifenberg-flat domains, where the
corresponding estimate for $b$, needed to replace Theorem 4 is due to [37].

The structure of this paper is as follows. In Section 2 we explain our
notation and present a number of preliminary results. In Section 3 we show
convergence in $W^{1,p}(\Omega)$ of the iterative process. In Section 4 we
present our numerical experiments regarding the iterative process. In Section
5 we use the iterative process to prove higher differentiability and weighted
estimates of perturbations to equations with $p$-Laplace structure.

## 2 Preliminaries and notation

$c$ will denote a a positive constant depending only on $\Omega,N,p$ that may
change from line to line. We write $a\lesssim b$ if there is $c>0$ depending
only on $\Omega,N,p$ such that $a\leq cb$.

Throughout $\Omega\subset\mathbb{R}^{n}$ will be a domain. Let $1\leq
p<\infty$. We denote by $L^{p}(\Omega)=L^{p}(\Omega,\mathbb{R}^{N})$ and
$W^{1,p}(\Omega)=W^{1,p}(\Omega,\mathbb{R}^{N})$ the usual Lebesgue and
Sobolev spaces respectively. $VMO(\Omega)$ denotes the space of maps with
vanishing mean oscillation. Further $p^{\prime}$ will denote the Hölder
conjugate of $p$. We also employ the standard fractional Sobolev spaces
$W^{k,p}(\Omega)$ where ${k\in(0,\infty)}$, whose theory can be found for
example in [47]. We recall in particular the following fact from [23]: if
$p\in(1,\infty)$, $\alpha>0$ and $\Omega$ is a regular polytope, then for
$u\in W^{1+\alpha,p}(\Omega)$ the best approximation $v$ to $u$ in the space
of continuous piecewise linear functions subordinate to a shape-regular
triangulation of $\Omega$ of mesh-size $h$ satisfies

$\displaystyle\|u-v\|_{W^{1,p}_{\infty}(\Omega)}\lesssim
h^{\alpha}\|u\|_{W^{1+\alpha,p}_{\infty}(\Omega)}.$ (2.1)

For $z\in\mathbb{R}^{n}$, $r>0$ we denote by $B_{r}(z)$ the open ball of
radius $r$ around $z$.

Given vectors $u,v\in\mathbb{R}^{n}$ we denote by $|v|$ the Euclidean norm and
we denote the inner product in $\mathbb{R}^{n}$ by $u\cdot v$.

Given a matrix $A\in\mathbb{R}^{n\times n}$ we denote by $\|A\|$ the operator
norm of $A$. If $A$ is positive definite, we denote its smallest eigenvalue by
$\lambda_{A}$ and its largest by $\Lambda_{A}$. In particular when $a\equiv
a(x,z)\colon\Omega\times\mathbb{R}^{N\times n}\to\mathbb{R}^{N\times n}$ is
such that $a(x,\cdot)$ is differentiable in $\mathbb{R}^{N\times n}$ for
almost every $x\in\Omega$, we view $\partial_{z}a(x,y,z)$ both as a matrix and
as a linear form. If $\partial_{z}a(x,y,z)$ is positive definite, we denote
$\lambda_{a}=\lambda_{A}$, $\Lambda_{a}=\Lambda_{A}$.

We recall the definition of Muckenhoupt weights. For a fixed $1<p<\infty$, we
say that a weight $w\colon\mathbb{R}^{n}\to[0,\infty)$ belongs to $A_{p}$ if
$w$ is locally integrable and there is a constant $C>0$ such that, for all
balls $B$ in $\mathbb{R}^{n}$, we have

$\fint_{B}w(x)\,\mathrm{d}x\left(\fint_{B}w^{-\frac{q}{p}}\,\mathrm{d}x\right)^{\frac{p}{q}}\leq
C<\infty.$

Given a weight $w$ on $\Omega$, we denote the weighted Lebesgue-spaces by
$L^{p}_{w}(\Omega)$.

Given $\mu\geq 0$, $p\in(1,\infty)$ and $v\in\mathbb{R}^{n}$ we write
$V_{\mu,p}(v)=(\mu^{2}+|v|^{2})^{\frac{p-2}{4}}v$. When the choice of $p$ is
clear from the context we suppress the index and write
$V_{\mu}(v)=V_{\mu,p}(v)$.

We recall some well-known tools for dealing with $p$-growth. The following
Lemma is standard, but the author has been unable to find a version in the
literature with explicit bounds that tend to $1$ as $\gamma\to 0$. Hence a
proof of the bounds of this type shown in the version below can be found in
the appendix.

######  Lemma 2.1.

Let $\xi,\eta\in\mathbb{R}^{m}$ for some $m>0$. For $\gamma\geq 0$ and
$\mu\geq 0$ we have

$\frac{1}{6^{\gamma}(2\gamma+1)}(\mu^{2}+|\eta|^{2}+|\eta-\xi|^{2})^{\gamma}\leq\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|^{2})^{\gamma}\,\mathrm{d}t\leq
2^{\gamma}(\mu^{2}+|\eta|^{2}+|\eta-\xi|^{2})^{\gamma}$

If $\gamma\in(-1/2,0]$ we have

$\displaystyle
2^{\gamma}(\mu^{2}+|\eta|^{2}+|\eta-\xi|^{2})^{\gamma}\leq\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|^{2})^{\gamma}\,\mathrm{d}t\leq\frac{1}{4^{\gamma}(\gamma+1)}(\mu^{2}+|\eta|^{2}+|\xi-\eta|^{2})^{\gamma}$

We also recall that the $V_{\mu}$-functional enjoys a Young-type inequality.

######  Lemma 2.2 (cf.[1], Lemma 2.3).

Let $x,y\in\mathbb{R}^{N\times n}$, $\mu\geq 0$ and $1\leq p$. Let
$\varepsilon>0$. Then with

$\displaystyle
C_{\varepsilon}=\max\left(\frac{1}{4\varepsilon},\frac{(p-1)^{p-1}}{p^{p}\varepsilon^{p-1}}\right),$

it holds that

$\displaystyle(\mu^{2}+|x|^{2})^{\frac{p-2}{2}}x\cdot
y\leq\varepsilon|V_{\mu}(x)|^{2}+C_{\varepsilon}|V_{\mu}(y)|^{2}.$

Further we note the following estimate:

######  Lemma 2.3.

For $a,b\in\mathbb{R}^{n}$ we have

$\displaystyle\frac{3-\sqrt{5}}{2}(|a|^{2}+|b|^{2})\leq|a|^{2}+|a-b|^{2}\leq\frac{3+\sqrt{5}}{2}(|a|^{2}+|b|^{2})$

We close this section by recalling a linear algebra result from [30]. Let
$A\in\mathbb{R}^{n\times n}$ be positive definite. Write $A=A^{+}A^{-}$ where
$A^{+},A^{-}$ are the symmetric and skew-symmetric part of $A$, respectively.
Set $C=A^{+}A^{-}-A^{-}A^{+}-(A^{-})^{2}$ and denote by $\sigma$ the largest
eigenvalue of $C$. Define

$\displaystyle K_{\gamma}\coloneqq\inf_{\gamma>0}\|I-\gamma A\|.$

Then we have

######  Lemma 2.4 (c.f. Lemma 1.1.2 in [30]).

The optimal constant $K_{\gamma}$ is achieved with the following choice:

$\displaystyle\begin{cases}K_{\gamma}^{2}=\dfrac{\sigma}{\sigma+\lambda_{A}^{2}},\qquad\gamma=\dfrac{\lambda_{A}}{\sigma+\lambda_{A}^{2}}\quad&\text{
if }\sigma\geq\dfrac{\lambda_{A}(\Lambda_{A}-\lambda_{A})}{2}\\\\[20.0pt]
K_{\gamma}^{2}=\dfrac{(\Lambda_{A}-\lambda_{A})^{2}+4\sigma}{(\Lambda_{A}+\lambda_{A})^{2}},\quad\gamma=\dfrac{2}{\Lambda_{A}+\lambda_{A}}\quad&\text{
if }\sigma\leq\dfrac{\lambda_{A}(\Lambda_{A}-\lambda_{A})}{2}\\\ \end{cases}$

######  Remark 2.5.

An inspection of the proof shows that the inequality $\|I-\gamma A\|\leq
K_{\gamma}$ still holds when $\lambda_{A},\Lambda_{A}$ and $\sigma$ are
replaced by any upper bound for the smallest and largest eigenvalue of $A$ and
the largest eigenvalue of $C$, respectively, in the definition of $K_{\gamma}$
and $\gamma$.

## 3 The iterative process

We call $b\colon\Omega\times\mathbb{R}^{N\times n}\to\mathbb{R}^{N\times n}$ a
reference field if the following holds: $b$ satisfies
$\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and (A4). In addition, whenever $F\in
L^{p^{\prime}}(\Omega,\mathbb{R}^{N\times n})$ and $g\in
W^{1,p}(\Omega,\mathbb{R}^{N})$, there exists a unique solution $u\in
W^{1,p}_{g}(\Omega)$ to the problem

$\displaystyle\textup{div}\,b(x,\textup{D}u)=\textup{div}\,F.$ (3.1)

(3.1) is to be understood in the sense of distributions.

######  Remark 3.1.

The $p$-Laplacian $|z|^{p-2}z$ is an example of such a field. We encourage the
reader to think of this case on a first reading.

Recall the definition of the iterative process. Let ${u_{0}\in
W^{1,p}_{g}(\Omega)}$ and $f\in L^{p^{\prime}}(\Omega)$. Take $\gamma>0$ to be
chosen later. Define $u_{n+1}\in W^{1,p}_{g}(\Omega)$ to be the weak solution
of the problem (3.1) with the choice

$\displaystyle F=b(x,\textup{D}u_{n})-\gamma a(x,\textup{D}u_{n})+\gamma f.$
(3.2)

Note that LABEL:ass:A3 ensures that $F\in
L^{p^{\prime}}(\Omega,\mathbb{R}^{N\times n})$, so that the sequence
$\\{u_{n}\\}$ is well-defined.

We want to show that $u_{n}\to u$ in $W^{1,p}(\Omega)$ where $u$ is a solution
of (1.1). The crucial observation to prove convergence is the following linear
algebra observation:

######  Lemma 3.2.

Suppose $A,B\in\mathbb{R}^{m\times m}$ are positive definite and symmetric.
Then, with the choice
$\gamma_{*}=\frac{2}{\Lambda_{B^{-1}A}+\lambda_{B^{-1}A}}$ and

$K=\frac{\Lambda_{B^{-1}A}-\lambda_{B^{-1}A}}{\Lambda_{B^{-1}A}+\lambda_{B^{-1}A}}<1,$

the estimate

$\|B-\gamma_{*}A\|\leq K\|B\|$

holds.

[Proof]By Lemma 2.4 there holds

$\displaystyle\|B-\gamma A\|\leq$ $\displaystyle\|B\|\|I-\gamma
B^{-1}A\|\leq\|B\|\frac{\Lambda_{B^{-1}A}-\lambda_{B^{-1}A}}{\Lambda_{B^{-1}A}+\lambda_{B^{-1}A}}=\|B\|K$

The moreover part follows immediately from Lemma 2.4.

######  Remark 3.3.

Using the remark after Lemma 2.4, whenever $a$ and $b$ are fields satisfying
$\mathrm{\ref{ass:A1}-\eqref{ass:A4}}$, we can use this Lemma to ensure the
existence of a $\gamma>0$ so that the inequality
$\|\partial_{z}b(x,z)-\gamma\partial_{z}a(x,z)\|\leq K_{a,b}\|B\|$ holds
uniformly in $x$ and $z$. We define

$\displaystyle K_{\gamma}\coloneqq\inf\\{c>0\colon\|B(x,z)-\gamma A(x,z)\|\leq
c\|B(x,z)\|\text{ for almost all }x\in\Omega,\text{ every
}z\in\mathbb{R}^{N\times n}\\}.$ (3.3)

In the case where $A,B$ are non-symmetric, the choice of $\gamma$ and $K$
needs to be modified with the obvious adaptions to the above proof coming from
Lemma 2.4. With these modifications the symmetry assumption (A4) may be
dropped in all results of this paper.

We are now able to prove Theorem 1, which we restate for the convenience of
the reader.

######  Theorem 6.

Suppose $b\colon\Omega\times\mathbb{R}^{N\times n}\to\mathbb{R}^{N}$ is a
reference field. Consider $a\colon\Omega\times\mathbb{R}^{N\times
n}\to\mathbb{R}^{N}$ satisfying $\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and (A4).
If

$\displaystyle\begin{cases}\dfrac{\Lambda_{b}}{\lambda_{b}}K_{a,b}\left(\dfrac{3+\sqrt{5}}{2}\right)^{\frac{p-2}{2p}}6^{p-2}(p-1)<1&\quad\text{
for }p\geq 2\\\\[20.0pt]
\dfrac{\Lambda_{b}}{\lambda_{b}}K_{a,b}\left(\dfrac{3-\sqrt{5}}{2}\right)^{\frac{(p-2)}{4}}\dfrac{2^{2-p}}{p-1}<1&\quad\text{
for }p\leq 2\end{cases}$

then (1.1) has a unique solution $u\in W^{1,p}_{g}(\Omega)$. Moreover, there
is $\gamma>0$ such that if the sequence ${u_{n}\in W^{1,p}_{g}(\Omega)}$ is
generated via (3.2), then $u_{n}\to u$ in $W^{1,p}(\Omega)$.

######  Remark 3.4.

From the proof it will be clear that, in fact, $\\{u_{n}\\}$ generated via
(3.2) converges to the weak solution of (1.1) whenever $\gamma>0$ is such that

$\displaystyle\begin{cases}R^{p-1}=\dfrac{\Lambda_{b}}{\lambda_{b}}K_{\gamma}\left(\dfrac{3+\sqrt{5}}{2}\right)^{\frac{p-2}{2p}}6^{p-2}(p-1)<1&\quad\text{
if }p\geq 2\\\
R^{p/2}=\dfrac{\Lambda_{b}}{\lambda_{b}}K_{\gamma}\left(\dfrac{3-\sqrt{5}}{2}\right)^{\frac{(p-2)}{4}}\dfrac{2^{2-p}}{p-1}<1&\quad\text{
if }p\leq 2.\end{cases}$

Here $K_{\gamma}$ is the constant from (3.3). Further, the convergence occurs
at linear rate $R$. To be precise, the following estimates hold:

$\displaystyle\|u_{n+1}-u_{n}\|_{W^{1,p}(\Omega)}\leq$ $\displaystyle
R^{n}\|u_{1}-u_{0}\|_{W^{1,p}(\Omega)}$
$\displaystyle\|u_{n+1}-u\|_{W^{1,p}(\Omega)}\leq$ $\displaystyle
R^{n}\|u_{1}-u\|_{W^{1,p}(\Omega)},$

where $u$ is the weak solution of (1.1).

Before presenting the proof we want to briefly outline the main idea. Consider
the sequence $(u_{n})$ generated via (3.2). Subtracting the equations defining
$u_{n+1}$ and $u_{n}$ we find

$\displaystyle\textup{div}\,\left(b(x,\textup{D}u_{n+1})-b(x,\textup{D}u_{n})\right)=\textup{div}\,\left(b(x,\textup{D}u_{n})-b(x,\textup{D}u_{n-1})-\gamma\left(a(x,\textup{D}u_{n})-a(x,\textup{D}u_{n-1})\right)\right)$

Using the mean-value theorem and denoting $A(x,z)=\partial_{z}a(x,z)$,
$B(x,z)=\partial_{z}b(x,z)$, we may rewrite this as

$\displaystyle\textup{div}\,\left(B(x,\xi_{n+1})\textup{D}(u_{n+1}-u_{n})\right)=\textup{div}\,\left((B-\gamma
A)(x,\tilde{\xi}_{n})\textup{D}(u_{n}-u_{n-1})\right).$

for some $\xi_{n+1}$ lying on the line segment between $\textup{D}u_{n+1}$ and
$\textup{D}u_{n}$ and $\tilde{\xi}_{n}$ lying on the line segment between
$\textup{D}u_{n}$ and $\textup{D}u_{n-1}$. If $\|B(x,z)-\gamma A(x,z)\|$ is
sufficiently small, uniformly in $x$ and $z$, then applying the ellipticity
assumption to bound the left-hand side from below, we are able to show that
for some $R<1$,

$\|\textup{D}(u_{n+1}-u_{n})\|_{L^{p}(\Omega)}\leq
R\|\textup{D}(u_{n}-u_{n-1})\|_{L^{p}(\Omega)}.$

We then conclude easily. [Proof]Let $\gamma$ be as in Lemma 3.2, so that,
using the notation of the lemma, $K_{\gamma}=K_{a,b}$. Throughout the proof we
write $A(x,z)=\partial_{z}a(x,z)$ and $B(x,z)=\partial_{z}b(x,z)$.

Let $n\geq 1$. Test the equations defining $u_{n+1}$ and $u_{n}$ against
$u_{n+1}-u_{n}$ to find:

$\displaystyle\int_{\Omega}\left(b(x,\textup{D}u_{n+1})-b(x,\textup{D}u_{n})\right)\textup{D}(u_{n+1}-u_{n})\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{\Omega}\left(b(x,\textup{D}u_{n})-b(x,\textup{D}u_{n-1})\right)\textup{D}(u_{n+1}-u_{n})-\gamma\left(a(x,\textup{D}u_{n})-a(x,\textup{D}u_{n-1})\right)\textup{D}(u_{n+1}-u_{n})\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{\Omega}\int_{0}^{1}\left(B(x,z_{n}(\theta))-\gamma
A(x,z_{n}(\theta))\right)\textup{D}(u_{n}-u_{n-1})\cdot\textup{D}(u_{n+1}-u_{n})\
d\theta\,\mathrm{d}x.$

Here $z_{n}(\theta)=(1-\theta)\textup{D}u_{n-1}+\theta\textup{D}u_{n}$. The
last line follows from an application of the mean value theorem. Now proceed
to estimate both sides of this equality.

We focus first on the case $p\geq 2$. Using again the mean-value theorem and
writing $z_{n+1}(\theta)=(1-\theta)\textup{D}u_{n+1}+\theta\textup{D}u_{n}$,
by assumption LABEL:ass:A2, the left-hand side gives:

$\displaystyle\int_{\Omega}\left(b(x,\textup{D}u_{n+1})-b(x,\textup{D}u_{n})\right)\textup{D}(u_{n+1}-u_{n})\
dx$ $\displaystyle=$
$\displaystyle\int_{\Omega}\int_{0}^{1}(B(x,z_{n+1}(\theta))\textup{D}(u_{n+1}-u_{n})\cdot\textup{D}(u_{n+1}-u_{n})\,\mathrm{d}\theta\,\mathrm{d}x$
$\displaystyle\geq$
$\displaystyle\lambda_{b}\int_{\Omega}\int_{0}^{1}(\mu^{2}+|z_{n+1}(\theta)|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n+1}-u_{n})|^{2}\,\mathrm{d}\theta\,\mathrm{d}x$
$\displaystyle\geq$
$\displaystyle\frac{\lambda_{b}}{6^{\frac{p-2}{2}}(p-1)}\int_{\Omega}(\mu^{2}+|\textup{D}u_{n+1}|^{2}+|\textup{D}u_{n+1}-\textup{D}u_{n}|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n+1}-u_{n})|^{2}\,\mathrm{d}x$

where the last line uses Lemma 2.1.

Note that $\|B(x,w)\|\leq\Lambda_{b}(\mu^{2}+|w|^{2})^{\frac{p-2}{2}}$ by
LABEL:ass:A2. Hence on the right-hand side by Lemma 3.2 and Lemma 2.3,

$\displaystyle\int_{\Omega}\int_{0}^{1}\left(B(x,z_{n}(\theta))-\gamma
A(x,z_{n}(\theta))\right)\textup{D}(u_{n}-u_{n-1})\cdot\textup{D}(u_{n+1}-u_{n})\,\mathrm{d}\theta\,\mathrm{d}x$
$\displaystyle\leq$
$\displaystyle\Lambda_{b}K\int_{\Omega}\int_{0}^{1}(\mu^{2}+|z_{n}(\theta)|^{2})^{\frac{p-2}{2}}|\textup{D}u_{n}-\textup{D}u_{n-1}||\textup{D}u_{n+1}-\textup{D}u_{n}|\,\mathrm{d}\theta\,\mathrm{d}x$
$\displaystyle\leq$
$\displaystyle\Lambda_{b}K\int_{\Omega}(\mu^{2}+|\textup{D}u_{n}|^{2}+|\textup{D}u_{n}-\textup{D}u_{n-1}|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n}-u_{n-1})||\textup{D}(u_{n+1}-u_{n})|\,\mathrm{d}x=I.$

We now apply Lemma 2.2 and Lemma 2.1 to find for any $\varepsilon>0$,

$\displaystyle I\leq$
$\displaystyle\Lambda_{b}K\left(C_{\varepsilon}\int_{\Omega}V_{(\mu^{2}+|\textup{D}u_{n}|^{2})^{1/2}}(\textup{D}(u_{n+1}-u_{n}))^{2}+\varepsilon\int_{\Omega}V_{(\mu^{2}+|\textup{D}u_{n}|^{2})^{1/2}}(\textup{D}(u_{n}-u_{n-1}))^{2}\right)$
$\displaystyle\leq$
$\displaystyle\Lambda_{b}K\Big{(}\varepsilon\int_{\Omega}V_{(\mu^{2}+|\textup{D}u_{n}|^{2})^{1/2}}(\textup{D}(u_{n}-u_{n-1}))^{2}$
$\displaystyle+C_{\varepsilon}\left(\frac{3+\sqrt{5}}{2}\right)^{\frac{p-2}{2}}\int_{\Omega}(\mu^{2}+|\textup{D}u_{n+1}|^{2}+|\textup{D}u_{n+1}-\textup{D}u_{n}|^{2})^{\frac{p-2}{2}}|\textup{D}u_{n+1}-\textup{D}u_{n}|^{2}\Big{)}.$

Combining these two estimates and re-arranging gives:

$\displaystyle\int_{\Omega}(\mu^{2}+|\textup{D}u_{n+1}|^{2}+|\textup{D}u_{n+1}-\textup{D}u_{n}|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n+1}-u_{n})|^{2}\,\mathrm{d}x$
$\displaystyle\leq$ $\displaystyle\frac{\varepsilon
C_{1}}{C_{0}-C_{\varepsilon}C_{2}}\int_{\Omega}(\mu^{2}+|\textup{D}u_{n}|^{2}+|\textup{D}(u_{n-1}-u_{n})|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n+1}-u_{n})|^{2}\,\mathrm{d}x.$

where $C_{0}=\frac{\lambda_{b}}{6^{\frac{p-2}{2}}(p-1)}$, $C_{1}=\Lambda_{b}K$
and $C_{2}=(\frac{3+\sqrt{5}}{2})^{\frac{p-2}{2}}\Lambda_{b}K$.

Optimising in $\varepsilon$ and using the hypothesis this shows by induction
that

$\int_{\Omega}(|\textup{D}u_{n}|^{2}+|\textup{D}(u_{n+1}-u_{n})|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n+1}-u_{n})|^{2}\,\mathrm{d}x\to
0$

at linear rate.

As $p\geq 2$, it immediately follows that $\\{u_{n}\\}$ converges in
$W^{1,p}(\Omega)$ at linear rate. Necessarily the limit $u$ lies in
$W^{1,p}_{g}(\Omega)$ and it is a solution of problem (1.1). Moreover
considering any solution $v\in W^{1,p}_{g}(\Omega)$ of (1.1) and considering
the above estimates with the starting point

$\displaystyle\int_{\Omega}(b(x,\textup{D}u_{n+1})-b(x,\textup{D}v))\textup{D}(u_{n+1}-v)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{\Omega}(b(x,\textup{D}u_{n})-b(x,\textup{D}v)-\gamma(a(x,\textup{D}u_{n})-a(x,\textup{D}v))+\gamma
f)\cdot\textup{D}(u_{n+1}-v)\,\mathrm{d}x$

we find that, for any $u_{0}\in W^{1,p}_{g}(\Omega)$, the corresponding
iterative process converges to $v$ in $W^{1,p}_{g}(\Omega)$. In particular,
the solution to (1.1) is unique.

We now consider $p<2$. Proceed as before to obtain

$\displaystyle\int_{\Omega}\left(b(x,\textup{D}u_{n+1})-b(x,\textup{D}u_{n})\right)\textup{D}(u_{n+1}-u_{n})\,\mathrm{d}x$
$\displaystyle\geq$
$\displaystyle\lambda_{b}\int_{\Omega}(\mu^{2}+|\textup{D}u_{n+1}|^{2}+|\textup{D}u_{n+1}-\textup{D}u_{n}|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n+1}-u_{n})|^{2}\,\mathrm{d}x$

and

$\displaystyle\int_{\Omega}\int_{0}^{1}\left(B(x,z(\theta))-\gamma
A(x,z(\theta))\right)\textup{D}(u_{n}-u_{n-1})\cdot\textup{D}(u_{n+1}-u_{n})\,\mathrm{d}\theta\,\mathrm{d}x$
$\displaystyle\leq$
$\displaystyle\frac{\Lambda_{b}K}{2^{p-2}(p-1)}\Bigg{(}\left(\frac{3-\sqrt{5}}{2}\right)^{\frac{p-2}{2}}C_{\varepsilon}\int_{\Omega}(\mu^{2}+|\textup{D}u_{n+1}|^{2}+|\textup{D}u_{n+1}-\textup{D}u_{n}|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n}-u_{n-1})|^{2}\,\mathrm{d}x$
$\displaystyle+\varepsilon\int_{\Omega}(\mu^{2}+|\textup{D}u_{n}|^{2}+|\textup{D}u_{n}-\textup{D}u_{n-1}|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n}-u_{n-1})|^{2}\,\mathrm{d}x\Bigg{)}.$

Combining these estimates and our assumptions, we again find that

$\int_{\Omega}(|\textup{D}u_{n}|^{2}+|\textup{D}u_{n+1}|^{2})^{\frac{p-2}{2}}|\textup{D}(u_{n+1}-u_{n})|^{2}\,\mathrm{d}x\to
0$

at linear rate.

Note that by Hölder there is $c>0$ such that

$\displaystyle\int_{\Omega}(|\textup{D}u_{n}|+|\textup{D}(u_{n+1}-u_{n})|)^{p-2}|\textup{D}(u_{n+1}-u_{n})|^{2}\,\mathrm{d}x$
$\displaystyle\geq$ $\displaystyle
c\frac{\|\textup{D}(u_{n+1}-u_{n})\|_{p}^{2}}{\|\textup{D}u_{n}\|_{p}^{2-p}+\|\textup{D}u_{n+1}\|_{p}^{2-p}}.$

If $\\{\textup{D}u_{n}\\}$ is a bounded sequence in $W^{1,p}$, convergence of
$\\{u_{n}\\}$ in $W^{1,p}(\Omega)$ and existence and uniqueness of solutions
to (1.1) follow by repeating the arguments of the case $p\geq 2$. For clarity
of presentation and as it essentially follows by repeating arguments similar
to those of this proof we postpone the proof of boundedness to Lemma 3.5.

######  Lemma 3.5.

Suppose the assumptions of Theorem 6 hold. Then $\\{u_{n}\\}$ is bounded in
$W^{1,p}(\Omega)$.

[Proof]

Let $\gamma$ be as in Lemma 3.2, so that, using the notation of the lemma,
$K_{\gamma}=K_{a,b}$. We write $A(x,z)=\partial_{z}a(x,z)$ and
$B(x,z)=\partial_{z}(x,z)$.

Let $n\geq 1$. Test the equation defining $u_{n+1}$ against $u_{n+1}$ to find:

$\displaystyle\int_{\Omega}\left(b(x,\textup{D}u_{n+1})-b(x,0)\right)\textup{D}u_{n+1}\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{\Omega}\left(b(x,\textup{D}u_{n})-b(x,0)\right)\textup{D}u_{n+1}$
$\displaystyle-\gamma\left(a(x,\textup{D}u_{n})-a(x,0)\right)\textup{D}u_{n+1}+\gamma
a(x,0)\textup{D}u_{n+1}\ dx$ $\displaystyle=$
$\displaystyle\int_{\Omega}\int_{0}^{1}\left(B(x,z_{n}(\theta))-\gamma
A(x,z_{n}(\theta))\right)\textup{D}u_{n}\cdot\textup{D}u_{n+1}+\gamma
a(x,0)\textup{D}u_{n+1}\,\mathrm{d}\theta\,\mathrm{d}x.$

Here $z_{n}(\theta)=\theta u_{n}$. The last line follows from an application
of the mean value theorem. Now proceed to estimate both sides of this equality
exactly as in Theorem 6.

We focus first on the case $p\geq 2$. Using again the mean-value theorem and
writing $z_{n+1}(\theta)=\theta u_{n+1}$, by assumption LABEL:ass:A2, the
left-hand side gives:

$\displaystyle\int_{\Omega}\left(b_{\alpha}(x,\textup{D}u_{n+1})-b_{\alpha}(x,0)\right)\textup{D}u_{n+1}\,\mathrm{d}x=$
$\displaystyle\int_{\Omega}\int_{0}^{1}(B(x,z_{n+1}(\theta))\textup{D}(u_{n+1})\cdot\textup{D}u_{n+1}\,\mathrm{d}\theta\,\mathrm{d}x$
$\displaystyle\geq$
$\displaystyle\lambda_{b}\int_{\Omega}\int_{0}^{1}(\mu^{2}+|z_{n_{1}}(\theta)|^{2})^{\frac{p-2}{2}}|\textup{D}u_{n+1}|^{2}\,\mathrm{d}\theta\,\mathrm{d}x$
$\displaystyle\geq$
$\displaystyle\frac{\lambda_{b}}{(p-1)6^{\frac{p-2}{2}}}\int_{\Omega}(\mu^{2}+|\textup{D}u_{n+1}|^{2})^{\frac{p-2}{2}}|\textup{D}u_{n+1}|^{2}\,\mathrm{d}x$

where the last line uses Lemma 2.1.

Note that $\|B(x,z)\|\leq\Lambda_{b}(\mu^{2}+|z|^{2})^{\frac{p-2}{2}}$ by
LABEL:ass:A2. Hence on the right-hand side by Lemma 3.2 and Lemma 2.1, as well
as LABEL:ass:A3

$\displaystyle\int_{\Omega}\int_{0}^{1}\left(B(x,z_{n}(\theta))-\gamma
A(x,z_{n}(\theta))\right)\textup{D}u_{n}\cdot\textup{D}u_{n+1}+\gamma
a(x,0)\cdot\textup{D}u_{n+1}\,\mathrm{d}\theta\,\mathrm{d}x$
$\displaystyle\leq$
$\displaystyle\Lambda_{b}K\int_{\Omega}\int_{0}^{1}(\mu^{2}+|z_{n}(\theta)|^{2})^{\frac{p-2}{2}}|\textup{D}u_{n}||\textup{D}u_{n+1}|+C\gamma|\textup{D}u_{n+1}|\mathrm{d\theta}\,\mathrm{d}x$
$\displaystyle\leq$
$\displaystyle\Lambda_{b}K\int_{\Omega}(\mu^{2}+|\textup{D}u_{n}|^{2})^{\frac{p-2}{2}}|\textup{D}u_{n}||\textup{D}(u_{n+1})|\,\mathrm{d}x+C\gamma\int_{\Omega}|\textup{D}u_{n+1}|\,\mathrm{d}x=I+II.$

We now apply Lemma 2.2, to find

$\displaystyle I\leq$
$\displaystyle\Lambda_{b}K\left(C_{\varepsilon}\int_{\Omega}V_{\mu}(\textup{D}u_{n+1})^{2}+\varepsilon\int_{\Omega}V_{\mu}(\textup{D}u_{n})^{2}\right).$

Young’s inequality also gives

$\displaystyle II\leq$ $\displaystyle
C\gamma\|\textup{D}u_{n+1}\|_{L^{p}(\Omega)}\leq\varepsilon_{1}\|\textup{D}u_{n+1}\|_{L^{p}(\Omega)}^{p}+C(\varepsilon_{1}).$

Combining these two estimates and re-arranging gives:

$\displaystyle\int_{\Omega}(\mu^{2}+|\textup{D}u_{n+1}|^{2})^{\frac{p-2}{2}}|\textup{D}u_{n+1}|^{2}\,\mathrm{d}x$
$\displaystyle\leq$ $\displaystyle\frac{\varepsilon
C_{1}}{C_{0}-C_{\varepsilon}C_{2}}\int_{\Omega}(\mu^{2}+|\textup{D}u_{n}|^{2})^{\frac{p-2}{2}}|\textup{D}u_{n})|^{2}\,\mathrm{d}x+\frac{1}{C_{0}-\varepsilon
C_{2}}(C(\varepsilon_{1})+\varepsilon_{1}\|\textup{D}u_{n+1}\|_{L^{p}(\Omega)}^{p}).$

where $C_{0}=\frac{\lambda_{b}}{6^{\frac{p-2}{2}}(p-1)}$ and
$C_{1}=C_{2}=\Lambda_{b}K$.

Optimising in $\varepsilon$, we can ensure that $\frac{\varepsilon
C_{1}}{C_{0}-C_{\varepsilon}C_{1}}<1$. Now choosing $\varepsilon_{1}$
sufficiently small, we conclude that

$\displaystyle\|V_{\mu}(\textup{D}u_{n+1})\|_{L^{2}(\Omega)}^{2}\leq\eta\|V_{\mu}(\textup{D}u_{n+1})\|_{L^{2}(\Omega)}^{2}+C$
(3.4)

for some constant $\eta<1$. The conclusion follows using induction.

For $p\leq 2$ we argue similarly to obtain (3.4) with the choices
$C_{0}=\lambda_{b}$, $C_{1}=\Lambda_{b}K$,
$C_{2}=\Lambda_{b}K(p-1)^{-1}4^{(2-p)2}$. Optimising in $\varepsilon$ we can
again ensure that

$\dfrac{\varepsilon C_{1}}{C_{0}-C_{\varepsilon}C_{1}}<1$

and use this to conclude as before.

## 4 Numerical analysis and experiments

Throughout this section we assume that $a(x,z),b(x,z)$ satisfy
$\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and (A4). We further assume that the
assumptions of Theorem 6 hold. In particular, this will imply that, with an
appropriate choice of $\gamma$, the iterative process (3.2) converges in
$W^{1,p}(\Omega)$ to a solution of (1.1) and so, in particular, the iterates
$\\{u_{n}\\}$ are uniformly bounded in $W^{1,p}(\Omega)$. We fix such a choice
of $\gamma$ from now on.

Recall that we study the following numerical scheme: Suppose $\Omega$ is a
regular polytope. Consider a sequence of shape-regular triangulations
$\\{T_{h}\\}_{h\in(0,1]}$ of $\Omega$ of mesh-size $h$. Denote by $X_{h}$ the
space of continuous piecewise linear (subordinate to $T_{h}$) functions.
Choose $u_{0}\in X_{h}$. Define $u_{n+1}^{h}\in X_{h}$ inductively as a
solution of the problem:

$\displaystyle\begin{cases}\int_{\Omega}b(x,\textup{D}u_{n+1}^{h})\cdot\textup{D}\phi\,\mathrm{d}x=\int_{\Omega}(b(x,\textup{D}u_{n}^{h})-\gamma
a(x,\textup{D}u_{n}^{h})+\gamma
f)\cdot\textup{D}\phi\,\mathrm{d}x&\quad\forall\phi\in X_{h}.\\\
u_{n+1}=0&\quad\text{ on }\partial\Omega\end{cases}$ (4.1)

Note that due to $\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and the theory of
monotone operators (4.1) is well-defined.

###  4.1 Analysis of the numerical scheme

In this section we present our main results regarding convergence, a-priori
and a-posteriori estimates for (4.1). We begin with the following convergence
result:

######  Theorem 7.

Let $u_{0}\in X_{h}$ and consider the sequence $\\{u_{n}^{h}\\}$ generated by
(4.1) starting from $u_{0}$. Then $u_{n}^{h}\to u_{h}$ in $W^{1,p}(\Omega)$
where $u_{h}\in X_{h}$ solves

$\displaystyle\begin{cases}\int_{\Omega}a(x,\textup{D}u_{h})\cdot\textup{D}\phi\,\mathrm{d}x=\int_{\Omega}f\cdot\textup{D}\phi\,\mathrm{d}x&\quad\forall\phi\in
X_{h}\\\ u_{h}=0&\quad\text{ on }\partial\Omega.\end{cases}$ (4.2)

Moreover there is $0<C_{1}<1$ such that

$\displaystyle\|u_{n+1}^{h}-u_{n}^{h}\|_{W^{1,p}(\Omega)}\leq
C_{1}^{n}\|u_{1}-u_{0}\|_{W^{1,p}(\Omega)}$

[Proof]This follows from carrying out line by line the proof of Theorem 6,
replacing $\\{u_{n}\\}$ with $\\{u_{n}^{h}\\}$ and $u$ with $u_{h}$
respectively.

We now wish to show that $u_{h}\to u$ in $W^{1,p}(\Omega)$ as $h\to 0$ where
$u\in W^{1,p}_{0}(\Omega)$ solves (1.1). In order to prove this we make the
following regularity assumption on $b$: If $F\in L^{p^{\prime}}(\Omega)$ then
there is $\alpha>0$ such that (3.1) has a solution $u$ in
$W^{1+\alpha,p}(\Omega)$ and moreover we have the following estimate for some
$c>0$:

$\displaystyle\|u\|_{W^{1+\alpha,p}(\Omega)}\leq
c\left(1+\|u\|_{W^{1,p}(\Omega)}+\|F\|_{L^{p^{\prime}}(\Omega)}^{1/(p-1)}\right)$
(4.3)

######  Remark 4.1.

(4.3) is satisfied for example when $b(x,z)=|z|^{p-2}z$ with any choice of
$\alpha>0$ such that ${\alpha<\min\left(1/(p-1)^{2},(p-1)^{2}\right)}$, see
[44].

We proceed to study the effect of decreasing the mesh-size $h$. We use the
notation of Theorem 7.

######  Theorem 8.

Set $g=0$. Assume the assumptions of Theorem 7 hold. Suppose moreover that
(4.3) is satisfied. Choose $u_{0}\in X_{h}$ and let $\\{u_{n}^{h}\\}$ be the
sequence generated by (4.1). Suppose $u\in W^{1,p}_{0}(\Omega)$ solves (1.1).
Then $u\in W^{1+\alpha,p}(\Omega)$ and $u^{h}\to u$ in $W^{1,p}(\Omega)$ as
$h\to 0$. Moreover we have the estimate

$\displaystyle\|u_{n+1}^{h}-u\|_{W^{1,p}(\Omega)}\leq
ch^{\frac{2\alpha}{\max(2,p)}}c\left(\|f\|_{L^{p^{\prime}}(\Omega)}\right)+C_{1}^{n}\|u_{1}-u_{0}\|_{W^{1,p}(\Omega)}.$

[Proof]Consider the iterative process $\\{u_{n}\\}$ started from $0$. By
Theorem 6, $u_{n}\to u$ in $W^{1,p}(\Omega)$. Using (4.3), we find

$\displaystyle\|u_{n}\|_{W^{1+\alpha,p}(\Omega)}\leq
c\left(1+\|u_{n}\|_{W^{1,p}(\Omega)}+\|b(x,\textup{D}u_{n})-\gamma
a(x,\textup{D}u_{n})+\gamma f\|_{L^{p^{\prime}}(\Omega)}^{1/(p-1)}\right)$

Repeating arguments from the proof of Lemma 3.5 we deduce using induction that
there is $c(\|f\|_{L^{p^{\prime}}(\Omega)})$ such that for all $n\geq 0$,

$\displaystyle\|u_{n}\|_{W^{1+\alpha,p}(\Omega)}\leq
c(\|f\|_{L^{p^{\prime}}(\Omega)}).$

Extracting a weakly convergent subsequence we conclude the same estimate holds
for $u$.

Let $v$ be the best approximation to $u$ in $X_{h}$. Using LABEL:ass:A2, the
fact that $u,u_{h}$ solve (1.1) and (4.2) respectively we find

$\displaystyle\lambda_{a}\int_{\Omega}(\mu^{2}+|\textup{D}u|^{2}+|\textup{D}u^{h}|^{2})^{\frac{p-2}{2}}|\textup{D}u-\textup{D}u^{h}|^{2}\,\mathrm{d}x$
$\displaystyle\leq$
$\displaystyle\int_{\Omega}(a(x,\textup{D}u)-a(x,\textup{D}u^{h}))\textup{D}(u-u^{h})\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{\Omega}(a(x,\textup{D}u)-a(x,\textup{D}u^{h}))\textup{D}(u-v)\,\mathrm{d}x$
$\displaystyle\lesssim$
$\displaystyle\Lambda_{a}\int_{\Omega}(\mu^{2}+|\textup{D}u|^{2}+|\textup{D}(u-u^{h})|^{2})^{\frac{p-2}{2}}|\textup{D}(u-u^{h})||\textup{D}(u-v)|\,\mathrm{d}x$
$\displaystyle\lesssim$ $\displaystyle
C_{\varepsilon}\Lambda_{a}\int_{\Omega}(\mu^{2}+|\textup{D}u|^{2}+|\textup{D}u_{h}|^{2})^{\frac{p-2}{2}}|\textup{D}(u-u_{h})|^{2}\,\mathrm{d}x$
$\displaystyle\quad+\varepsilon\Lambda_{a}\int_{\Omega}(\mu^{2}+|\textup{D}u|^{2}+|\textup{D}v|^{2})^{\frac{p-2}{2}}|\textup{D}(u-v)|^{2}\,\mathrm{d}x$

where to obtain the last line we have used Lemma 2.2. Choosing $\varepsilon$
sufficiently large we conclude using (2.1)

$\displaystyle I=$
$\displaystyle\int_{\Omega}(\mu^{2}+|\textup{D}u|^{2}+|\textup{D}u^{h}|^{2})^{\frac{p-2}{2}}|\textup{D}(u-u^{h})|^{2}$
$\displaystyle\lesssim$
$\displaystyle\int_{\Omega}(\mu^{2}+|\textup{D}u|^{2}+|\textup{D}v|^{2})^{\frac{p-2}{2}}|\textup{D}(u-v)|^{2}\,\mathrm{d}x$
$\displaystyle\leq$
$\displaystyle\|\textup{D}(u-v)\|_{L^{p}(\Omega)}^{2}\left(1+\|\textup{D}u\|_{W^{1,p}(\Omega)}^{p-1}+\|\textup{D}(u-v)\|_{W^{1,p}(\Omega)}^{p-1}\right)$
$\displaystyle\lesssim$ $\displaystyle
h^{2\alpha}\|u\|_{W^{1+\alpha,p}(\Omega)}\left(1+\|u\|_{W^{1,p}(\Omega)}^{p-1}+h^{\alpha(p-1)}\|u\|_{W^{1+\alpha,p}(\Omega)}^{p-1}\right)$
(4.4)

If $p\geq 2$, $I\geq\|\textup{D}(u-u_{h})\|_{L^{p}(\Omega)}^{p}$ whereas if
$p\leq 2$, we have by applying Hölder’s inequality
${I\geq\|\textup{D}(u-u^{h})\|_{L^{p}(\Omega)}^{2}(1+\|\textup{D}u\|_{L^{p}(\Omega)}+\|\textup{D}u^{h}\|_{L^{p}(\Omega)})}$.
Recalling the standard estimates
$\|\textup{D}u\|_{L^{p}(\Omega)}\leq\|f\|_{L^{p^{\prime}}(\Omega)}^{1/(p-1)}$
and
$\|\textup{D}u^{h}\|_{L^{p}(\Omega)}\leq\|f\|_{L^{p^{\prime}}(\Omega)}^{1/(p-1)}$,
we conclude by combining (4.1), Theorem 7 and the inequality

$\displaystyle\|u_{n+1}^{h}-u\|_{W^{1,p}(\Omega)}\leq\|u_{n+1}^{h}-u^{h}\|_{W^{1,p}(\Omega)}+\|u^{h}-u\|_{W^{1,p}(\Omega)}.$

We also have an a-posteriori error bound. The proof follows [26].

######  Proposition 4.2.

Assume the conditions of Theorem 7 hold. Then we have

$\displaystyle\|\textup{D}(u_{n}^{h}-u^{h})\|_{L^{p}(\Omega)}\leq
C(\|f\|_{L^{p^{\prime}}(\Omega)},\|\textup{D}u_{0}\|_{L^{p}(\Omega)})\|\textup{D}(u_{n}^{h}-u_{n-1}^{h})\|_{L^{p}(\Omega)}^{\frac{2}{\max(2,p)}}.$

[Proof]We compute using LABEL:ass:A2 and Lemma 2.2,

$\displaystyle\lambda\gamma\int_{\Omega}(\mu^{2}+|\textup{D}(u^{h}-u_{n-1}^{h})|^{2}+|\textup{D}u_{n-1}^{h}|^{2})^{\frac{p-2}{2}}|\textup{D}u^{h}-\textup{D}u_{n-1}^{h}|^{2}\,\mathrm{d}x$
$\displaystyle\lesssim$
$\displaystyle\gamma\int_{\Omega}(a(x,\textup{D}u^{h})-a(x,\textup{D}u_{n-1}^{h}))\cdot\textup{D}(u^{h}-u_{n-1}^{h})\,\mathrm{d}x$
$\displaystyle=$ $\displaystyle\int_{\Omega}\gamma
f\cdot\textup{D}(u^{h}-u_{n-1}^{h})+(b(x,\textup{D}(u_{n}^{h}))-b(x,\textup{D}u_{n-1}^{h}))-\gamma
f)\cdot\textup{D}(u^{h}-u_{n-1}^{h})\,\mathrm{d}x$ $\displaystyle\lesssim$
$\displaystyle\int_{\Omega}\left(\mu^{2}+|\textup{D}u_{n-1}^{h}|^{2}+|\textup{D}(u_{n}^{h}-u_{n-1}^{h})|^{2}\right)^{\frac{p-2}{2}}|\textup{D}(u_{n}^{h}-u_{n-1}^{h})||\textup{D}(u^{h}-u_{n-1}^{h})|\,\mathrm{d}x$
$\displaystyle\leq$
$\displaystyle\varepsilon\int_{\Omega}\left(\mu^{2}+|\textup{D}u_{n-1}^{h}|^{2}+|\textup{D}(u_{n}^{h}-u_{n-1}^{h})|^{2}\right)^{\frac{p-2}{2}}|\textup{D}(u_{n}^{h}-u_{n-1}^{h})|^{2}\,\mathrm{d}x$
$\displaystyle\quad+C_{\varepsilon}\int_{\Omega}(\mu^{2}+|\textup{D}(u^{h}-u_{n-1}^{h})|^{2}+|\textup{D}u_{n-1}^{h}|^{2})^{\frac{p-2}{2}}|\textup{D}u^{h}-\textup{D}u_{n-1}^{h}|^{2}\,\mathrm{d}x$

Thus choosing $\varepsilon$ sufficiently large, re-arranging, employing by now
standard arguments and recalling that $\\{u_{n}^{h}\\}$ is bounded uniformly
in $W^{1,p}(\Omega)$ we obtain

$\displaystyle\|\textup{D}(u^{h}-u_{n-1}^{h})\|_{L^{p}(\Omega)}\lesssim
c\left(\|\textup{D}u^{h}\|_{L^{p}(\Omega)},\|\textup{D}u_{0}\|_{L^{p}(\Omega)}\right)\|\textup{D}(u_{n}^{h}-u_{n-1}^{h})\|_{L^{p}(\Omega)}^{\frac{2}{\max(2,p)}}.$

By the triangle inequality, and using
$\|u^{h}\|_{W^{1,p}(\Omega)}\lesssim\|f\|_{L^{p^{\prime}}(\Omega)}^{1/(p-1)}$,
we conclude the desired estimate:

$\displaystyle\|\textup{D}(u_{n}^{h}-u^{h})\|_{L^{p}(\Omega)}\leq$
$\displaystyle\|\textup{D}(u_{n-1}^{h}-u^{h})\|_{L^{p}(\Omega)}+\|\textup{D}(u_{n}^{h}-u_{n-1}^{h})\|_{L^{p}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
c\left(\|f\|_{L^{p^{\prime}}(\Omega)},\|\textup{D}u_{0}\|_{L^{p}(\Omega)}\right)\|\textup{D}(u_{n}^{h}-u_{n-1}^{h})\|_{L^{2}(\Omega)}.$

We close this section by detailing a modification of the numerical scheme,
following the algorithm outlined in [26]. We strengthen our assumptions on
$\\{T_{h}\\}$ and assume $\\{T_{k}\\}_{k\in\mathbb{N}}$ is a sequence of
shape-regular triangulations with meshsize $h_{k}\to 0$ as $k\to\infty$.
Moreover, we assume that $\\{T_{k}\\}$ is obtained from $\\{T_{k-1}\\}$ by
refinement. Denote by $X^{k}$ the space of continuous piecewise linear
functions subordinate to $\\{T_{k}\\}$. We then consider the following
algorithm:

Set $u_{0}^{0}=0$. Set $k=0$. Choose a maximal mesh $k_{max}$ and a sequence
of tolerances $\textup{tol}_{k}\to 0$ as $k\to\infty$.

repeat

while $\|u_{n}^{k}-u_{n-1}^{k}\|>\textup{tol}_{k}$ do

solve $-\Delta u^{k}_{n+1}=-\Delta
u^{k}_{n}+\gamma\textup{div}\,(a(x,\textup{D}u^{k}_{n})-f)$ in $X^{k}$

$n\leftarrow n+1$

end while

$u_{0}^{k+1}\leftarrow u_{n}^{k}$

$k\leftarrow k+1$

until $k=k_{max}$.

Algorithm 1 A variant of the iteration scheme

From the results of this section it is clear that under the assumptions of
Theorem 8 Algorithm 1 converges to a solution of (1.1) as $k_{max}\to\infty$.

###  4.2 Numerical experiments

The results presented in this section are obtained using Firedrake [4, 5, 8,
43, 2, 3]. Throughout this section $\Omega=[0,1]^{3}$. We consider shape-
regular triangulations of $\Omega$ with uniformly spaced nodes at distance
$h=2^{-i}$. Note that in this set-up the mesh-size is $\sqrt{3}h$.

A linear example: We first consider a linear example where the exact solution
is known. We choose,

$f=\left(\begin{matrix}8\pi&8\pi&8\pi\\\ 10\pi&10\pi&10\pi\\\
2\pi&2\pi&2\pi\end{matrix}\right)(v,v,v)^{T}\quad\text{ where
}v=\left(\begin{matrix}cos(2\pi x)sin(2\pi y)sin(2\pi z)\\\ sin(2\pi
x)cos(2\pi y)sin(2\pi z)\\\ sin(2\pi x)sin(2\pi y)cos(2\pi
z)\end{matrix}\right)$

and consider the problem

$\displaystyle\begin{cases}-\textup{div}\,A\,\textup{D}u=\textup{div}\,f\quad&\text{
in }\Omega\\\ u=0&\text{ on }\partial\Omega\end{cases}\qquad\text{ where
}A=\left(\begin{matrix}1&1&2\\\ 0&2&3\\\ 0&0&1\end{matrix}\right).$

Note that with $\tilde{u}(x,y,z)=sin(2\pi x)sin(2\pi y)sin(2\pi z)$,

$u(x,y,z)=(\tilde{u}(x,y,z),\tilde{u}(x,y,z),\tilde{u}(x,y,z))^{T}$

is the exact solution of this problem.

For the iteration scheme we choose $u_{0}=0$, $b(x,z)=z$ and $\gamma=2/3$ or
$\gamma=1/2$. We solve each iteration step using GMRES with an incomplete LU
factorisation to precondition the problem. The iteration is terminated when
$\|u_{n+1}-u_{n}\|_{H^{1}(\Omega)}\leq 10^{-9}$.

We record the $H^{1}$-error of the numerical solution $v$ computed using the
iterative scheme (4.1) in Table 1. We also record the number of iterations
needed.

$h$ | $2^{-1}$ | $2^{-2}$ | $2^{-3}$ | $2^{-4}$ | $2^{-5}$
---|---|---|---|---|---
$\gamma=2/3$, error | $6.6924$ | $5.3767$ | $3.2456$ | $1.6786$ | $0.8436$
$\gamma=2/3$, iterations | $1$ | $22$ | $22$ | $22$ | $22$
$\gamma=1/2$, error | $6.6924$ | $5.3767$ | $3.2456$ | $1.6786$ | $0.8436$
$\gamma=1/2$, iterations | $1$ | $33$ | $33$ | $34$ | $34$
Table 1: $\|u-v\|_{H^{1}(\Omega)}$

We also compute the numerical solution $u^{h}$ directly using a
$LU$-factorisation. We record the $H^{1}$-distance between $u_{n}^{h}$ and
$u^{h}$ for two different choices of $\gamma$ in Fig.1.

(a) $\gamma=2/3$ (b) $\gamma=2/3$

(c) $\gamma=1/2$

Figure 1: $\|u_{n}^{h}-u_{h}\|_{H^{1}(\Omega)}$

A nonlinear example We also consider the nonlinear problem

$\displaystyle\begin{cases}-\textup{div}\,(1+|\textup{D}u|^{4})A\textup{D}u+|u|^{4}u=(y,x^{2},z^{2}+x^{2})^{T}\quad&\text{
in }\Omega\\\ u=0&\text{ on }\partial\Omega\end{cases}\qquad\text{with
}A=\left(\begin{matrix}1&3&5\\\ 0&2&4\\\ 0&0&1\end{matrix}\right).$

We set $b(x,z)=(1+|z|^{4})z$, $\gamma=0.65$ and $u_{0}=0$ and employ (4.1).
Each (non-linear) iteration step is solved using the standard ’solve’-method
in Firedrake. This utilises a nonlinear Newton linesearch scheme where the
linear step is computed using GMREs. Denote the solution obtained in this way
with $u_{dir}$. We record the $H^{1}$-distance
$\|u_{dir}-u_{n}^{h}\|_{H^{1}(\Omega)}$ in Fig.2.

Figure 2: $\|u_{dir}-u_{n}^{h}\|_{L^{2}(\Omega)}$

For this problem, we moreover compare computing times for (4.1) and Algorithm
(1) in Table 2. For (4.1) we use $\|u_{n+1}-u_{n}\|_{H^{1}(\Omega)}\leq
10^{-9}$ as our termination condition, while for Algorithm 1 we choose
$\textup{tol}_{k}=10^{-(i+4)}$ when $h=2^{-i}$.

$h$ | $2^{-1}$ | $2^{-2}$ | $2^{-3}$ | $2^{-4}$ | $2^{-5}$
---|---|---|---|---|---
iteration scheme | $0.201$ | $0.715$ | $3.931$ | $47.48$ | $821.4$
Algorithm 1 | $0.333$ | $1.000$ | $3.3365$ | $33.41$ | $576.9$
Table 2: runtime in seconds

######  Remark 4.3.

We do not apply a scheme optimised for the $p$-Laplacian in order to solve
each iteration step. Improving our computation in this way a further reduction
in runtime should be achieved.

## 5 Improved regularity results

In this section, we assume that

$b(x,z)=(B(x)z\cdot z)^{\frac{p-2}{2}}B(x)z,$

where $B$ is measurable and satisfies for $x\in\Omega$,
$\xi\in\mathbb{R}^{N\times n}$,

$\displaystyle\lambda|\xi|^{2}\leq B(x)\xi\cdot\xi\leq\Lambda|\xi|^{2}.$ (5.1)

We note that $b$ satisfies $\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ with $\mu=0$.
We will further assume that $B$ is symmetric almost everywhere in $\Omega$, so
that (A4) holds.

###  5.1 A first example: Calderón-Zygmund type estimates

We will use fields satisfying the assumptions of Theorem 2 as reference fields
and apply the iterative process to obtain the following result.

######  Theorem 9.

Let $1<p<\infty$. Suppose $f\in W^{1,q/(p-1)}(\Omega)$ and let $g=0$. Assume
$a$ satisfies $\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and (A4) with $\mu=0$.
Further suppose there is $B$ such that (5.1) and the assumptions of Theorem 2
are satisfied. Finally assume that with this choice of $b$ the assumptions of
Theorem 6 are satisfied. If

$C_{1}=\dfrac{C_{0}^{(p-1)/q}\Lambda_{b}K_{a,b}}{(p-1)}<1,$

where $K_{a,b}$ is given by (1.3), then the solution $u$ of (1.1) satisfies

$\displaystyle\|\textup{D}u\|_{L^{q}(\Omega)}\lesssim
1+\|f\|_{L^{q/(p-1)}}^{\frac{1}{p-1}}.$

Here $C_{0}$ is the constant arising in Theorem 2.

[Proof]

Let $\\{u_{n}\\}$ be the sequence generated by the iterative process with
$u_{0}=0$ and with the choice of $\gamma$ as in in Theorem 6. By Theorem 2,
Lemma 3.2 and LABEL:ass:A3, we have, for some $c>0$, the estimate

$\displaystyle\int_{\Omega}|\textup{D}u_{n+1}|^{q}\,\mathrm{d}x\leq$
$\displaystyle C_{0}\int_{\Omega}|b(x,\textup{D}u_{n})-\gamma
a(x,\textup{D}u_{n})+\gamma f|^{\frac{q}{p-1}}\,\mathrm{d}x$
$\displaystyle\leq$ $\displaystyle
C_{0}\int_{\Omega}|b(x,\textup{D}u_{n})-b(x,0)-\gamma(a(x,\textup{D}u_{n})+a(x,0))+\gamma
f|^{\frac{q}{p-1}}\,\mathrm{d}x+c$ $\displaystyle\leq$ $\displaystyle
C_{0}\Lambda_{b}^{\frac{q}{p-1}}K^{\frac{q}{p-1}}\int_{\Omega}\left(\int_{0}^{1}\theta^{p-2}\,\mathrm{d}\theta\right)^{\frac{q}{p-1}}|\textup{D}u_{n}|^{q}+C_{0}\int_{\Omega}|f|^{\frac{q}{p-1}}\,\mathrm{d}x+c$
$\displaystyle=$
$\displaystyle\frac{C_{0}\Lambda_{b}^{\frac{q}{p-1}}K^{\frac{q}{p-1}}}{(p-1)^{q/(p-1)}}\int_{\Omega}|\textup{D}u_{n}|^{q}+C_{0}\int_{\Omega}|f|^{\frac{q}{p-1}}\,\mathrm{d}x+c.$

Thus, as by assumption $C_{1}<1$, we find by induction

$\displaystyle\|\textup{D}u_{n+1}\|_{L^{q}(\Omega)}\lesssim\frac{C_{1}^{\frac{q}{p-1}}}{1-C_{1}^{\frac{q}{p-1}}}(1+\|f\|_{L^{q/(p-1)}}^{\frac{1}{p-1}}).$

Extracting a weakly convergent subsequence and noting that $u_{n}\to u$ in
$W^{1,p}(\Omega)$ by Theorem 6, where $u$ is the solution of (1.1), we find
the desired estimate holds.

###  5.2 A second example: Weighted estimates and Hölder continuity

We can use the Koshelev iteration to perturb Theorem 4 as follows:

######  Theorem 10.

Suppose $a$ is a field satisfying $\mathrm{\ref{ass:A1}-\ref{ass:A3}}$ and
(A4) with $\mu=0$ and moreover there is symmetric $B$ satisfying the
assumptions of Theorem 4 such that the assumptions of Theorem 6 apply with
this choice of $b(x,z)$. If

$\displaystyle
C_{3}=\dfrac{C_{2}^{\frac{p-1}{q}}K_{a,b}\Lambda_{b}}{(p-1)}<1,$

where $K_{a,b}$ is given by (1.3), then the solution $v$ of the boundary value
problem (1.1) satisfies the estimate

$\displaystyle\|\textup{D}v\|_{L^{q}_{w}(\Omega)}\lesssim
1+\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}^{\frac{1}{p-1}}.$

Here $C_{2}$ is the constant arising in Theorem 4.

[Proof]The proof is similar to the proof of Theorem 9. Let $\\{u_{n}\\}$ be
given by the iterative process with $u_{0}=0$ and $\gamma$ chosen as in
Theorem 6.

Note that $b(x,0)-\gamma a(x,0)+\gamma f\in L^{\frac{q}{p-1}}_{w}(\Omega)$.
Thus we find,

$\displaystyle\|\textup{D}u_{n+1}\|_{L^{q}_{w}(\Omega)}^{q}\leq$
$\displaystyle C_{2}\int_{\Omega}|b(x,\textup{D}u_{n})-\gamma
a(x,\textup{D}u_{n})+\gamma f|^{\frac{q}{p-1}}w(x)\,\mathrm{d}x$
$\displaystyle\leq$ $\displaystyle
C_{2}\int_{\Omega}|b(x,\textup{D}u_{n})-b(x,0)-\gamma(a(x,\textup{D}u_{n})-a(x,0))|^{\frac{q}{p-1}}w(x)\,\mathrm{d}x$
$\displaystyle+C_{2}|b(x,0)-\gamma a(x,0)+\gamma
f|^{\frac{q}{p-1}}w(x)\,\mathrm{d}x$ $\displaystyle\leq$ $\displaystyle
C_{2}K_{a,b}^{\frac{q}{p-1}}\Lambda_{b}^{\frac{q}{p-1}}2^{-\frac{q}{p-1}}\|\textup{D}u_{n}\|_{L^{q}_{w}(\Omega)}^{q}+c(C_{2},p,q)\left(1+\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}^{\frac{q}{p-1}}\right).$

Hence by induction, as by assumption $C_{3}<1$,

$\displaystyle\|\textup{D}u_{n+1}\|_{L^{q}_{w}(\Omega)}\lesssim\dfrac{C_{3}^{\frac{q}{p-1}}}{1-C_{3}^{\frac{q}{p-1}}}\left(1+\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}^{\frac{q}{p-1}}\right).$

Extracting a weakly convergent subsequence and noting that $u_{n}\to u$ in
$W^{1,p}(\Omega)$ by Theorem 6 where $u$ solves (1.1) we conclude the desired
estimate, after passing to the limit in the estimate.

We note that of particular interest is the choice $w(x)=|x|^{\alpha}$ which
can be used to obtain estimates in the the familiar Morrey spaces and hence
through the Morrey-Sobolev embedding allows to obtain continuity statements.
Recall the definition of the $L^{q,\theta}$-Morrey-norm:

$\displaystyle\|u\|_{L^{q,\theta}(\Omega)}=\sup_{0<r<\text{diam}(\Omega),z\in\Omega}r^{\frac{\theta-n}{q}}\|u\|_{L^{q}(B_{r}(z)\cap\Omega)},$

where $\theta\in(0,n)$. We assume that all the assumptions and the notation of
Theorem 10 hold and show how to deduce estimates in Morrey spaces.

Fix $z\in\Omega$, $r\in(0,\text{diam}(\Omega))$ and choose for
$\rho\in(0,\theta)$,

$\displaystyle
w(x)=\min\left(|x-z|^{-n+\theta-\rho},r^{-n+\theta-\rho}\right).$

Then by our previous work,

$\displaystyle\|\textup{D}u\|_{L^{q,\theta}(B_{r}(z)\cap\Omega)}^{q}\leq
r^{n-\theta+\rho}\|\textup{D}u\|_{L^{q}_{w}(B_{r}(z)\cap\Omega)}^{q}\leq
r^{n-\theta+\rho}c(C,p,q)\left(1+\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}^{\frac{q}{p-1}}\right).$

It remains to estimate $\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}$. For this we
proceed exactly as [37] but provide the argument here for the sake of
completeness. We will show that $\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}\leq
c\|f\|_{L^{q,\theta}(\Omega)}^{q}r^{-\rho}$, which will conclude the proof.

It is convenient to introduce $f^{\prime}$ where
$|f^{\prime}|^{p-2}f^{\prime}=f$. For $\alpha>0$ we denote the set
${E_{\alpha}=\\{x\in\Omega\colon|f^{\prime}|>\alpha\\}}$. Then we can write

$\displaystyle\|f\|_{L^{\frac{q}{p-1}}(\Omega)}^{q}=\|f^{\prime}\|_{L^{q}(\Omega)}^{q}=$
$\displaystyle
q\int_{0}^{\infty}\alpha^{q}\int_{E_{\alpha}}w(x)\,\mathrm{d}x\frac{d\alpha}{\alpha}$
$\displaystyle\leq$ $\displaystyle
q\int_{0}^{\infty}\alpha\int_{0}^{r^{-n+\theta-\rho}}\left|E_{\alpha}\cap
B_{\beta^{\frac{1}{-n+\theta-\rho}}}(z)\right|\,\mathrm{d}\beta\frac{d\alpha}{\alpha}.$

We now estimate the inner integral as follows:

$\displaystyle\int_{0}^{r^{n+\theta-\rho}}\left|E_{\alpha}\cap
B_{\beta^{\frac{1}{-n+\theta-\rho}}}(z)\right|d\beta\leq$
$\displaystyle\sum_{i=1}^{\infty}2^{-i}r^{-n+\theta-\beta}\left|E_{\alpha}\cap
B_{r2^{\frac{-i}{-n+\theta-\rho}}}(z)\right|$ $\displaystyle\leq$
$\displaystyle
2\int_{0}^{\frac{1}{2}r^{-n+\theta-\rho}}\beta\left|E_{\alpha}\cap
B_{\beta^{\frac{1}{-n+\beta-\rho}}}(z)\right|\frac{d\beta}{\beta}.$

Now returning to the original estimate and applying Fubini’s theorem we
conclude

$\displaystyle\|f\|_{L^{\frac{q}{p-1}}_{w}(\Omega)}^{q}\leq$ $\displaystyle
2q\int_{0}^{\infty}\alpha^{q}\int_{0}^{\frac{1}{2}r^{-n+\theta-\rho}}\beta\left|E_{\alpha}\cap
B_{\beta^{\frac{1}{-n+\theta-\rho}}}(z)\right|\frac{\,\mathrm{d}\beta}{\beta}\frac{\,\mathrm{d}\alpha}{\alpha}$
$\displaystyle\leq$ $\displaystyle
2q\|f\|_{L^{q,\theta}(\Omega)}^{q}\int_{0}^{\frac{1}{2}r^{-n+\theta-\rho}}\beta^{1+\frac{n-\theta}{-n+\theta-\rho}}\leq
c\|f\|_{L^{q,\theta}(\Omega)}^{q}r^{-\rho}.$

This gives the desired result.

##  Appendix A Proof of Lemma 2.1

We restate and prove Lemma 2.1 here.

######  Lemma A.1.

For $\gamma\geq 0$ and $\mu\geq 0$ we have

$\frac{1}{6^{\gamma}(2\gamma+1)}(\mu^{2}+|\eta|^{2}+|\eta-\xi|^{2})^{\gamma}\leq\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|^{2})^{\gamma}\,\mathrm{d}t\leq
2^{\gamma}(\mu^{2}+|\eta|^{2}+|\eta-\xi|^{2})^{\gamma}$

If $\gamma\in(-1/2,0]$ we have

$\displaystyle
2^{\gamma}(\mu^{2}+|\eta|^{2}+|\eta-\xi|^{2})^{\gamma}\leq\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|^{2})^{\gamma}\,\mathrm{d}t\leq\frac{1}{4^{\gamma}(\gamma+1)}(\mu^{2}+|\eta|^{2}+|\xi-\eta|^{2})^{\gamma}$

[Proof]We first consider $\gamma\geq 0$. Then by Young’s inequality,

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|^{2})^{\gamma}\,\mathrm{d}t\leq$
$\displaystyle
2^{\gamma}\int_{0}^{1}(\mu^{2}+|\eta|^{2}+t^{2}|\xi-\eta|^{2})^{\gamma}\,\mathrm{d}t$
$\displaystyle\leq$ $\displaystyle
2^{\gamma}(\mu^{2}+|\eta|^{2}+|\xi-\eta|^{2})^{\gamma}$

For the lower bound, note that by symmetry we may assume $|\eta|\geq|\xi|$.
Dividing both sides by $|\eta|^{\gamma}$ we see that moreover we may assume
$|\eta|=1$. Rotating coordinate axis if necessary we may even assume that
$\eta=e_{1}$, the unit vector in the first coordinate direction. Note that

$\displaystyle|\eta+t(\xi-\eta)|\geq|1-t|\xi-\eta||.$

Write $s=|\xi-\eta|$. Note that $s\in[0,2]$.

We first assume $s\leq 1$. Then as $\gamma\geq 0$,

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|^{\gamma}\,\mathrm{d}t\geq\int_{0}^{1}(\mu^{2}+(1-ts)^{2})^{\gamma}\,\mathrm{d}t\geq
2^{-\gamma}\int_{0}^{1}(\mu+1-ts)^{2\gamma}\,\mathrm{d}t$ $\displaystyle=$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)s}\left((\mu+1)^{2\gamma+1}-(\mu+1-s)^{2\gamma+1}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)s}(\mu+1)^{2\gamma}\left((\mu+1)\left(1-\left(\frac{\mu+1-s}{\mu+1}\right)^{2\gamma}\right)+s\left(\frac{\mu+1-s}{\mu+1}\right)^{2\gamma}\right)$
$\displaystyle\geq$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)}(\mu+1)^{2\gamma}\left((\mu+1)\left(1-\left(\frac{\mu+1-s}{\mu+1}\right)^{2\gamma}\right)+\left(\frac{\mu+1-s}{\mu+1}\right)^{2\gamma}\right)$
$\displaystyle\geq$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)}(\mu+1)^{2\gamma}$
$\displaystyle\geq$
$\displaystyle\frac{1}{4^{\gamma}(2\gamma+1)}(\mu^{2}+1+s^{2})^{\gamma}.$

Next assume $s\geq 1$. Then

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|^{\gamma}\,\mathrm{d}t\geq\int_{0}^{1}(\mu^{2}+(1-ts)^{2})^{\gamma}\,\mathrm{d}t$
$\displaystyle\geq$ $\displaystyle
2^{-\gamma}\int_{0}^{1/s}(\mu+1-ts)^{2\gamma}\,\mathrm{d}t+2^{-\gamma}\int_{1/s}^{1}(\mu+ts-1)^{2\gamma}$
$\displaystyle=$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)s}\left((\mu+1)^{2\gamma+1}-\mu^{2\gamma+1}+(\mu+s-1)^{2\gamma+1}-\mu^{2\gamma+1}\right)$
$\displaystyle\geq$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)s}(\mu+1+s)^{2\gamma}\left(\left(\frac{\mu+1}{\mu+1+s}\right)^{2\gamma}+(s-1)\left(\frac{\mu+s-1}{\mu+1+s}\right)^{2\gamma}\right)$

Define

$\displaystyle
f(s)=\frac{1}{s}\left(\left(\frac{\mu+1}{\mu+1+s}\right)^{2\gamma}+(s-1)\left(\frac{\mu+s-1}{\mu+1+s}\right)^{2\gamma}\right)$

Note that

$\displaystyle f^{\prime}(s)=$
$\displaystyle\frac{1}{s^{2}}\left((1-\frac{2}{1+\mu+s})^{2\gamma}+4\gamma(s-1)s\frac{(\mu+s-1)^{2\gamma-1}}{(\mu+s+1)^{2\gamma+1}}-\frac{(\mu+1)^{2\gamma}(1+\mu+s+2\gamma
s)}{(\mu+s+1)^{2\gamma+1}}\right)$ $\displaystyle\leq 0$

Hence we conclude

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|^{\gamma}\,\mathrm{d}t\geq$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)}(\mu+s)^{2\gamma}\left(\frac{\mu+1}{\mu+3}\right)^{2\gamma}$
$\displaystyle\geq$
$\displaystyle\frac{1}{6^{\gamma}(2\gamma+1)}(\mu+s)^{2\gamma}$

This completes the case $\gamma\geq 0$.

We now turn to $\gamma\in(-\frac{1}{2},0)$. We find

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\eta+(1-t)\xi|^{2})^{\gamma}\,\mathrm{d}t\geq$
$\displaystyle\int_{0}^{1}(\mu^{2}+2(|\eta|^{2}+|\xi-\eta|^{2})^{2})^{\gamma}\,\mathrm{d}t$
$\displaystyle\geq$ $\displaystyle
2^{\gamma}(\mu^{2}+|\eta|^{2}+|\eta-\xi|^{2})^{\gamma}.$

We turn to the upper bound. Again set $s=|\xi-\eta|$. First assume $s\leq
1$.Then

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|)^{\gamma}\,\mathrm{d}t\leq$
$\displaystyle\int_{0}^{1}(\mu^{2}+(1-ts)^{2})^{\gamma}\,\mathrm{d}t$
$\displaystyle\leq$ $\displaystyle
2^{-\gamma}\int_{0}^{1}(\mu+1-ts)^{2\gamma}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)s}\left((\mu+1)^{2\gamma+1}-(\mu+1-s)^{2\gamma+1}\right).$

Set

$\displaystyle
f(s)=\frac{1}{s}\left((\mu+1)^{2\gamma+1}-(\mu+1-s)^{2\gamma+1}\right)$

and note

$\displaystyle
f^{\prime}(s)=(1+\mu)((1+\mu-s)^{2\gamma}-(1+\mu)^{2\gamma})+2\gamma
s(1+\mu-s)^{2\gamma}\geq 0$

to conclude

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\eta+(1-t)\xi|^{2})^{\gamma}\,\mathrm{d}t\leq$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)}((\mu+1)^{2\gamma+1}-\mu^{2\gamma+1})$
$\displaystyle\leq$
$\displaystyle\frac{(\mu+1)^{2\gamma}}{2^{\gamma}(2\gamma+1)}$
$\displaystyle\leq$
$\displaystyle\frac{(\mu+1+s^{2})^{2\gamma}}{4^{\gamma}(2\gamma+1)}.$

If $s\geq 1$ we estimate

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|)^{\gamma}\,\mathrm{d}t\leq$
$\displaystyle\int_{0}^{1}(\mu^{2}+(1-ts)^{2})^{\gamma}\,\mathrm{d}t$
$\displaystyle\leq$ $\displaystyle
2^{-\gamma}\int_{0}^{1/s}(\mu+1-ts)^{2\gamma}+2^{-\gamma}\int_{1/s}^{1}(\mu+ts-1)^{2\gamma}$
$\displaystyle=$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)s}\left(\mu+1)^{2\gamma+1}-\mu^{2\gamma+1}+(\mu+s-1)^{2\gamma+1}-\mu^{2\gamma+1}\right)$
$\displaystyle\leq$
$\displaystyle\frac{1}{2^{\gamma}(2\gamma+1)s}\left((\mu+1)^{2\gamma}+(s-1)(\mu+s-1)^{2\gamma}\right)$

Define

$\displaystyle
g(s)=\frac{1}{s}\left((\frac{\mu+1}{\mu+s})^{2\gamma}+(s-1)\left(\frac{\mu+s-1}{\mu+s}\right)^{2\gamma}\right)$

and note that

$\displaystyle 0\leq
g^{\prime}(s)=\frac{1}{s^{2}(\mu+s-1)}(m+s)^{-1-2\gamma}(-(\mu+1)^{2\gamma}(\mu+s-1)(\mu+s+2\gamma
s)$
$\displaystyle\qquad+(m+s-1)^{2\gamma}(\mu^{2}+(1+2\gamma)(s-1)s+\mu(2s-1)))$
$\displaystyle\Leftrightarrow$
$\displaystyle(\frac{\mu+1}{\mu+s-1})^{2\gamma}>\frac{\mu^{2}+(1+2\gamma)s(s-1)+\mu(2s-1)}{(\mu+s-1)(\mu+s+2\gamma
s)}$

By Bernoulli’s inequality it suffices to check that

$\displaystyle\mu+s-1+2\gamma(s-2)\geq\frac{\mu^{2}+(1+2\gamma)s(s-1)+\mu(2s-1)}{\mu+s+2\gamma
s}$ $\displaystyle\Leftrightarrow$ $\displaystyle 2\gamma\mu(s-2)+2\gamma
s(s-2)+4\gamma^{2}s(s-2)\geq 0.$

It is straightforward to see that the last inequality holds. Thus $g(s)$ is
increasing and we can conclude

$\displaystyle\int_{0}^{1}(\mu^{2}+|t\xi+(1-t)\eta|)^{\gamma}\,\mathrm{d}t\leq$
$\displaystyle\frac{(\mu+s)^{2\gamma}}{2^{\gamma}(2\gamma+1)}$
$\displaystyle\leq$
$\displaystyle\frac{(\mu+s+1)^{2\gamma}}{4^{\gamma}(2\gamma+1)}.$

Acknowledgments: The author would like to thank Jan Kristensen for suggesting
the topic of this paper to him. He would also like to thank Ioannis
Papadopoulos and Pascal Heid for useful discussion and practical suggestions
regarding the numerical experiments.

## References

  * [1] Acerbi, E., and Mingione, G. Regularity results for stationary electro-rheological fluids.  Arch. Ration. Mech. Anal. 164 (2002), 213–259.
  * [2] Amestoy, P., Duff, I., L’Excellent, J.-Y., and Koster, J. A fully asynchronous multifrontal solver using distributed dynamic scheduling.  SIAM Journal on Matrix Analysis and Applications 23, 1 (2001), 15–41.
  * [3] Amestoy, P., Guermouche, A., L’Excellent, J.-Y., and Pralet, S. Hybrid scheduling for the parallel solution of linear systems.  Parallel Computing 32, 2 (2006), 136–156.
  * [4] Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W., Karpeyev, D., Kaushik, D., Knepley, M., May, D. A., McInnes, L. C., Mills, R. T., Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S., Zhang, H., and Zhang, H. PETSc Users Manual.  Tech. Rep. ANL-95/11 - Revision 3.11, Argonne National Laboratory, 2019\.
  * [5] Balay, S., Gropp, W., McInnes, L., and Smith, B. Efficient Management of Parallelism in Object Oriented Numerical Software Libraries.  In Modern Software Tools in Scientific Computing (1997), E. Arge, A. M. Bruaset, and H. P. Langtangen, Eds., Birkhäuser Press, pp. 163–202.
  * [6] Caffarelli, L. Interior a Priori Estimates for Solutions of Fully Non-Linear Equations.  Annals of Mathematics 130, 1 (1989), 189–213.
  * [7] Cordes, H. Über die erste Randwertaufgabe bei quasilinearen Differentialgleichungen in mehr als zwei Variablen.  Math. Ann. 131 (1956), 278–312.
  * [8] Dalcin, L., Paz, R., Kler, P., and Cosimo, A. Parallel distributed computing using Python.  Advances in Water Resources 34, 9 (2011), 1124–1139.
  * [9] Daněček, J., John, O., and Stará, J. Interior $C^{1,\gamma}$-regularity for weak solutions of nonlinear second order elliptic systems.  Math. Nachr. 276 (2004), 47–56.
  * [10] Daněček, J., John, O., and Stará, J. Structural conditions guaranteeing $C^{1,\gamma}$-regularity of weak solutions to nonlinear second-order elliptic systems.  Nonlinear Analysis 66 (2007), 288–300.
  * [11] Daněček, J., and Viszus, E. Regularity on the interior for some class of nonlinear second-order elliptic systems.  Electron. J. Diff. Eqns. 2013, 121 (2013), 1–17.
  * [12] De Giorgi, E. Un esempio di estremali discotinue per un problema variazionale di tipo ellitico.  Boll. Un. Mat. Ital 4, 1 (1968), 135–137.
  * [13] Diening, L., Fornasier, M., Tomasi, R., and Wank, M. A Relaxed Kačanov iteration for the $p$-poisson problem.  Numerische Mathematik 145 (2020), 1–34.
  * [14] Diening, L., Scharle, T., and Schwarzacher, S. Regularity for parabolic systems of Uhlenbeck type with Orlicz growth.  Journal of Mathematical Analysis and Applications 472 (2019), 46–60.
  * [15] Diening, L., Stroffolini, B., and Verde, A. Everywhere regularity for functionals with $\phi$-growth.  Manuscr. Math. 4 (2009), 449–481.
  * [16] Dolzmann, G., Kristensen, J., and Zhang, K. BMO and uniform estimates for multi-well problems.  Manuscr. Math. 140, 1-2 (2013), 83–114.
  * [17] Frehse, J. A discontinuous solution of a mildly nonlinear elliptic system.  Math. Z. 134 (1973), 229–230.
  * [18] Garau, E., Morin, P., and Zuppa, C. Convergence of an adaptive Kačanov FEM for quasi-linear problems.  Appl. Numer. Math. 61, 4 (2011), 512–529.
  * [19] Giaquinta, M., and Giusti, E. Differentiability of minima of nondifferentiable functionals.  Invent. Math. 72 (1983), 285–298.
  * [20] Giaquinta, M., and Modica, G. Remarks on the regularity of the minimizers of certain degenerate functionals.  Manuscr. Math. 57 (1986), 55–99.
  * [21] Giusti, E. Direct Methods in the Calculus of Variations.  World Scientific, 2003.
  * [22] Giusti, E., and Miranda, M. Un esempio di soluzioni discontinue per un problema di inimo relativo ad un integrale regolare del calcolo delle variazioni.  Boll. Un. Mat. Ital. 4, 1 (1968), 219–226.
  * [23] Guermond, J.-L., and Ern, A. Finite element quasi-interpolation and best approximation.  ESAIM: Mathematical Modelling and Numerical Analysis 51, 4 (2017), 1367–1385.
  * [24] Han, W., Jensen, S., and Shimansky, I. The Kačanov method for some nonlinear problems.  Appl. Numer. Meth. 24 (1997), 57–79.
  * [25] Heid, P., and Wihler, T. A modified Kačano iteration scheme with application to quasilinear diffusion models.  arXiv preprint arXiv:2101.10137 (2021).
  * [26] Heid, P., and Wihler, T. P. On the convergence of adaptive iterative linearized Galerkin methods.  Calcolo 57, 24 (2020).
  * [27] John, O., Malý, J., and Stará, J. Nowhere continuous solutions to elliptic systems.  Commentationes Mathematicae Universitatis Carolinae 30 (1989), 33–43.
  * [28] Kačanov, L. Variational methods of solution of plasticity problems.  J. Appl. Math. Mech. 23 (1959), 880–883.
  * [29] Kinnunen, J., and Zhou, S. A boundary estimate for nonlinear equations with discontinuous coefficients.  Differential and Integral Equations 14, 4 (2001), 475–492.
  * [30] Koshelev, A. Regularity Problem for Quasilinear Elliptic and Parabolic Systems.  Springer-Verlag Berlin Heidelberg, 1991.
  * [31] Ladyzhenskaya, O., and Ural’tseva, N. Linear and quasilinear elliptic equations.  Academic Press Inc., New York, 1968.
  * [32] Leonardi, S. Remarks on the regularity of solutions of elliptic systems.  In Nonlinear Applied Analysis, A. Sequeira, H. da Veiga, and J. H. Videman, Eds. Kluwer Academic/Plenum Publishers, New York, 1999, pp. 325–344.
  * [33] Lieberman, G. Hölder regularity for the gradients of solutions of degenerate parabolic systems.  Ukr. Math. Bull. 3 (2006), 352–373.
  * [34] Marcellini, P. Everywhere regularity for a class of elliptic systems without growth conditions.  Ann. Sc. Norm. Super. Pisa, Cl. Sci. 23, 1 (1996), 1–25.
  * [35] Marcellini, P., and Papi, G. Nonlinear elliptic systems with general growth.  J. Differ. Equations 221, 2 (2006), 412–443.
  * [36] Maz’ya, V. Examples of nonregular solutions of quasilinear elliptic equations with analytic coefficients.  Funkcional. Anal. i Priložen. 2, 3 (1968), 53–57.
  * [37] Mengesha, T., and Phuc, N. C. Global Estimates for Quasilinear Elliptic Equations on Reifenberg Flat Domains.  Arch. Ration. Mech. Anal. 203 (2012), 189–216.
  * [38] Mingione, G. Regularity of minima: an invitation to the dark side of the calculus of variations.  Appl. Math 51, 4 (2006), 355–426.
  * [39] Morrey, C. On the solutions of quasi-linear elliptic partial differential equations.  Transactions of the American Mathematical Society 43 (1938), 126–166.
  * [40] Nečas, J., John, O., and Stará, J. Counterexample to the regularity of weak solution of elliptic systems.  Commentationes Mathematicae Universitatis Carolinae 21 (1980), 145–154.
  * [41] Nirenberg, L. On a generalization of quasi-conformal mappings and its application to elliptic partial differential equations.  Annuals of Mathematics Studies 33 (1954), 95–100.
  * [42] Phuc, N. C. Weighted estimates for nonhomogeneous quasilinear equations with discontinuous coefficients.  Ann. Scuola Norm. Sup. Pisa Cl. Sci. 10, 5 (2011), 1–17.
  * [43] Rathgeber, F., Ham, D., Mitchell, L., Lange, M., Luporini, F., McRae, A., Bercea, G.-T., Markall, G., and Kelly, P. Firedrake: automating the finite element method by composing abstractions.  ACM Trans. Math. Softw. 43, 3 (2016), 24:1—-24:27.
  * [44] Simon, J. Regularite de la solution d’un probleme aux limites non lineaires.  Ann. fac. sci. Toulouse 3, 3-4 (1981), 247–274.
  * [45] Souček, J. Singular solutions to linear elliptic systems.  Commentationes Mathematicae Universitatis Carolinae 25 (1984), 273–281.
  * [46] Tolksdorf, P. Everywhere regularity for some quasi-linear systems with a lack of ellipticity.  Ann. Mat. Pura App. 134 (1983), 241–266.
  * [47] Triebel, H. Theory of Function Spaces, 1 ed.  Birkhäuser Basel, Basel, 1983.
  * [48] Uhlenbeck, K. Regularity for a class of non-linear elliptic systems.  Acta Math. 138 (1977), 219–240.
  * [49] Wiegner, M. Regularity theorems for nondiagonal elliptic systems.  Arkiv för Matematik 20, 1-2 (1982), 1–13.
  * [50] Zeidler, E. Nonlinear functional analysis and its applications.  Springer, New York, 1988.