# High-order QMC nonconforming FEMs for nearly incompressible planar
stochastic elasticity equations††thanks: School of Mathematics and Statistics,
University of New South Wales, Sydney, Australia††thanks: This work was
supported by the Australian Research Council grant DP220101811.

J. Dick, T. Le Gia, W. McLean, K. Mustapha, T. Tran

###### Abstract

In a recent work (Dick et al, QMC sparse grid Galerkin finite element methods
for linear elasticity equations with uncertainties,
https:arxiv.org/pdf/2310.06187.pdf), we considered a linear stochastic
elasticity equation with random Lamé parameters which are parameterized by a
countably infinite number of terms in separate expansions. We estimated the
expected values over the infinite dimensional parametric space of linear
functionals ${\mathcal{L}}$ acting on the continuous solution ${\bf u}$ of the
elasticity equation. This was achieved by truncating the expansions of the
random parameters, then using a high-order quasi-Monte Carlo (QMC) method to
approximate the high dimensional integral combined with the conforming
Galerkin finite element method (FEM) to approximate the displacement over the
physical domain $\Omega.$ In this work, as a further development of
aforementioned article, we focus on the case of a nearly incompressible linear
stochastic elasticity equation. To serve this purpose, in the presence of
stochastic inhomogeneous (variable Lamé parameters) nearly compressible
material, we develop a new locking-free symmetric nonconforming Galerkin FEM
that handles the inhomogeneity. In the case of nearly incompressible material,
one known important advantage of nonconforming approximations is that they
yield optimal order convergence rates that are uniform in the Poisson
coefficient. Proving the convergence of the nonconforming FEM leads to another
challenge that is summed up in showing the needed regularity properties of
${\bf u}$. For the error estimates from the high-order QMC method, which is
needed to estimate the expected value over the infinite dimensional parametric
space of ${\mathcal{L}}{\bf u},$ we are required here to show certain
regularity properties of ${\bf u}$ with respect to the random coefficients.
Some numerical results are delivered at the end.

##  1 Introduction

In this work we explore and analyze the application of the quasi-Monte Carlo
(QMC) method combined with a new nonconforming Galerkin finite element method
(FEM) to solve a linear elastic model with uncertainties [14, 20, 21, 22, 29].
More specifically, we consider the case where the properties of the elastic
inhomogeneous material (that is, the Lamé parameters) are varying spatially in
an uncertain way (see (1.2)). The equation governing small elastic
deformations of a body $\Omega$ in $\mathbb{R}^{2}$ with polygonal boundary
can be written as

$-\nabla\cdot\boldsymbol{\sigma}({\boldsymbol{r}};{\bf
u}({\boldsymbol{x}},{\boldsymbol{r}}))={\bf f}({\boldsymbol{x}})\quad\text{for
}{\boldsymbol{x}}\in\Omega,~{}~{}{\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}}),$
(1.1)

subject to homogeneous Dirichlet boundary conditions ${\bf
u}(\cdot,{\boldsymbol{r}})={\bf 0}$ on $\Gamma:=\partial\Omega$ with
${\boldsymbol{y}}$ and ${\boldsymbol{z}}$ being parameter vectors describing
randomness. The parametric Cauchy stress tensor
$\boldsymbol{\sigma}({\boldsymbol{r}};\cdot)\in[L^{2}(\Omega)]^{2\times 2}$ is
defined as

$\boldsymbol{\sigma}({\boldsymbol{r}};{\bf
u}(\cdot,{\boldsymbol{r}}))=\lambda(\cdot,{\boldsymbol{z}})\Big{(}\nabla\cdot{\bf
u}(\cdot,{\boldsymbol{r}})\Big{)}I+2\mu(\cdot,{\boldsymbol{y}}){\boldsymbol{\varepsilon}}({\bf
u}(\cdot,{\boldsymbol{r}}))\quad{\rm on}\quad\Omega$

with ${\bf u}(\cdot,{\boldsymbol{r}})$ being the displacement vector field and
the symmetric strain tensor ${\boldsymbol{\varepsilon}}({\bf
u}):=\frac{1}{2}(\nabla{\bf u}+(\nabla{\bf u})^{T}).$ Here, ${\bf f}$ is the
body force per unit volume, and $I$ is the identity tensor. The gradient
($\nabla$) and the divergence ($\nabla\cdot$) are understood to be with
respect to the physical variable ${\boldsymbol{x}}\in\Omega$. The random
parameters $\mu$ and $\lambda$ are expressed in the following separate
expansions:

$\mu({\boldsymbol{x}},{\boldsymbol{y}})=\mu_{0}({\boldsymbol{x}})+\sum_{j=1}^{\infty}y_{j}\psi_{j}({\boldsymbol{x}})~{}~{}{\rm
and}~{}~{}\lambda({\boldsymbol{x}},{\boldsymbol{z}})=\Lambda\widehat{\lambda}({\boldsymbol{x}},{\boldsymbol{z}})=\Lambda\Big{(}\widehat{\lambda}_{0}({\boldsymbol{x}})+\sum_{j=1}^{\infty}z_{j}\phi_{j}({\boldsymbol{x}})\Big{)},$
(1.2)

where ${\boldsymbol{y}}=(y_{j})_{j\geq 1}$ and
${\boldsymbol{z}}=(z_{k})_{k\geq 1}$ belong to
$U:=(-\frac{1}{2},\frac{1}{2})^{\mathbb{N}}$ consisting of a countable number
of parameters $y_{j}$ and $z_{k}$, respectively, which are assumed to be
i.i.d. uniformly distributed. Here, the constant $\Lambda$ satisfies
$1\leq\Lambda<\infty$, and $\\{\psi_{j}\\}$ and $\\{\phi_{j}\\}$ are
orthogonal basis functions for $L^{2}(\Omega)$. For the well-posedness of the
elastic problem (1.1), we assume that $\mu_{\min}\leq\mu\leq\mu_{\max}$ and
$\widehat{\lambda}_{\min}\leq\widehat{\lambda}\leq\widehat{\lambda}_{\max}$ on
$\Omega\times U,$ for some positive constants $\mu_{\min}$, $\mu_{\max}$,
$\widehat{\lambda}_{\min},$ and $\widehat{\lambda}_{\max}.$

To ensure that $\mu$ and $\widehat{\lambda}$ are well-defined for all
parameters ${\boldsymbol{y}}$, ${\boldsymbol{z}}\in U$, we impose that

$\mu_{0},\,\widehat{\lambda}_{0}\in
L^{\infty}(\Omega),\quad\sum_{j=1}^{\infty}\|\psi_{j}\|_{L^{\infty}(\Omega)}<\infty~{}~{}{\rm
and}~{}~{}\sum_{j=1}^{\infty}\|\phi_{j}\|_{L^{\infty}(\Omega)}<\infty\,.$ (A1)

As in [13], in this work, we are interested in efficient approximation of the
expected value

$\mathcal{I}(\mathcal{L}({\bf u})):=\int_{U}\int_{U}\mathcal{L}\big{(}{\bf
u}(\cdot,{\boldsymbol{r}})\big{)}\,d{\boldsymbol{r}}$ (1.3)

of $\mathcal{L}({\bf u}(\cdot,{\boldsymbol{r}}))$ for a certain linear
functional $\mathcal{L}:{\bf L}^{2}(\Omega)\to\mathbb{R}$, with respect to the
random vector variables
${\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}}),$
$d{\boldsymbol{r}}=d{\boldsymbol{y}}\,d{\boldsymbol{z}}$, $d{\boldsymbol{y}}$
and $d{\boldsymbol{z}}$ are the uniform probability measures on $U$, and ${\bf
u}(\cdot,{\boldsymbol{r}})$ (displacement) is the random vector function that
solves problem (1.1). The main focus is the elastic material that is permitted
to be nearly incompressible, that is, the Poisson ratio of the elastic
material $\nu$ approaches $1/2$, which is not the case in [13]. So, in the
current analysis, we seek estimates that are independent of the (possibly very
large) constant factor $\Lambda$ in (1.2). Note that, as a special case, if
$\lambda$ is a constant multiple of $\mu$, that is, the randomnesses of the
Lamé parameters are due to the ones in the Young modulus $E$, then
$\nu=\frac{\lambda}{2(\lambda+\mu)}$ is constant, see [21, 22].

To estimate $\mathcal{I}(\mathcal{L}({\bf u}))$, we have to deal with three
sources of errors: a dimension truncation error from truncating the infinite
expansions in (1.2), a QMC quadrature error from approximating the high-
dimensional integral, and a nonconforming Galerkin error from approximating
the continuous solution over the physical domain $\Omega$. We explore these
steps and state the main results of this work (Theorem 2.3) in the next
section. Prior to this, we discuss the estimation of the expected value of a
given high dimensional function via a two different high-order QMC rules. The
first one is based on a tensor product structure of two separate high-order
rules for the random coefficients arising from the expansions for $\lambda$
and $\mu$, respectively, see Theorem 2.1. The second high-order QMC rule is
based on simulating both $\lambda$ and $\mu$ simultaneously using one family
of QMC points. Under some additional regularity assumptions, the latter QMC
rule leads to a better rate of convergence, see Theorem 2.2. We refer to [13]
for the proofs of Theorems 2.1 and 2.2.

The achieved regularity properties in [13] are not applicable here due to the
$\Lambda$ dependence. So, in Section 3, we show in detail the required
regularity properties of the continuous solution ${\bf u}$ (that are $\Lambda$
independent) with respect to both the random parameters ${\boldsymbol{y}}$ and
${\boldsymbol{z}}$, and the spatial variable ${\boldsymbol{x}}.$ These results
are needed to guarantee the convergence of the errors from both the QMC
integration and the nonconforming finite element discretization. For the
convergence of the latter, we require the parameters
$\mu(\cdot,{\boldsymbol{y}}),\,\widehat{\lambda}(\cdot,{\boldsymbol{z}})\in
W^{1,\infty}(\Omega)$ for every ${\boldsymbol{y}},\,{\boldsymbol{z}}\in U$. To
have this, we additionally assume that

$\mu_{0},\,\widehat{\lambda}_{0}\in
W^{1,\infty}(\Omega),~{}~{}\sum_{j=1}^{\infty}\|\nabla\psi_{j}\|_{L^{\infty}(\Omega)}<\infty~{}~{}{\rm
and}~{}~{}\sum_{j=1}^{\infty}\|\nabla\phi_{j}\|_{L^{\infty}(\Omega)}<\infty.$
(A2)

The error from truncating the infinite series expansion in (1.2) is
investigated in Section 4 (Theorem 4.1). More precisely, if ${\bf u}_{{\bf
s}}$, ${\bf s}:=(s_{1},s_{2})\in\mathbb{N}^{2}$, is the solution of (1.1)
obtained by truncating the infinite expansions in (1.2) at $j=s_{1}$ and
$k=s_{2}$, then the task is to estimate the error between ${\bf u}$ and ${\bf
u}_{{\bf s}}$. In other words, ${\bf u}_{{\bf s}}$ is the solution of (1.1)
corresponding to the truncated (finite) Lamé parameters
$\mu_{s_{1}}({\boldsymbol{x}},{\boldsymbol{y}}^{s_{1}})=\mu({\boldsymbol{x}},{\boldsymbol{y}}^{s_{1}})$
and
$\lambda_{s_{2}}({\boldsymbol{x}},{\boldsymbol{z}}^{s_{2}})=\lambda({\boldsymbol{x}},{\boldsymbol{z}}^{s_{2}}),$
where ${\boldsymbol{y}}^{s_{1}}=(y_{1},y_{2},\cdots,y_{s_{1}},0,0,\cdots)$ and
${\boldsymbol{z}}^{s_{2}}=(z_{1},z_{2},\cdots,z_{s_{2}},0,0,\cdots)$. To
minimize the errors (see Theorem 4.1), we assume that the functions $\psi_{j}$
and $\phi_{j}$ are ordered so that $\|\psi_{j}\|_{L^{\infty}(\Omega)}$ and
$\|\phi_{j}\|_{L^{\infty}(\Omega)}$ are nonincreasing, that is,

$\|\psi_{j}\|_{L^{\infty}(\Omega)}\geq\|\psi_{j+1}\|_{L^{\infty}(\Omega)}~{}~{}{\rm
and}~{}~{}\|\phi_{j}\|_{L^{\infty}(\Omega)}\geq\|\phi_{j+1}\|_{L^{\infty}(\Omega)},\quad{\rm
for}~{}~{}j\geq 1.$ (A3)

We may also require that (this assumption is needed to show a high-order QMC
error as well)

$\sum_{j=1}^{\infty}\|\psi_{j}\|^{p}_{L^{\infty}(\Omega)}<\infty~{}~{}{\rm
and}~{}~{}\sum_{j=1}^{\infty}\|\phi_{j}\|^{q}_{L^{\infty}(\Omega)}<\infty,\quad{\rm
for~{}some}~{}~{}0<p,\,q<1\,.$ (A4)

To avoid the locking and to deal with the variability of the Lamé parameters,
we develop a new symmetric low-order (piecewise linear) nonconforming Galerkin
FEM for approximating the parametric solution ${\bf
u}(\cdot,{\boldsymbol{r}})$ of (1.1) over the physical domain $\Omega$ in
Section 5, and prove error estimates. It is known that conforming FEMs
(studied in [13]) suffer from so called locking (or non-robustness) [1] in the
nearly incompressible case. The locking phenomenon was avoided by using
different approaches; these include nonconforming and mixed Galerkin FEMs [5,
6, 7, 8, 9, 15] for the case of constant Lamé parameters. The nonconforming
FEM in the existing literature is not applicable for the case of variable Lamé
parameters. Section 6 is devoted to numerical simulations. In four different
examples, we numerically illustrate the achieved theoretical convergence
results.

###### Remark 1.1

Apart from the regularity results in Theorem 3.2, the achieved results can be
extended to the case of homogeneous mixed boundary conditions of the form
${\bf u}(\cdot,{\boldsymbol{r}})=0$ on $\Gamma_{D}$ and ${\bf
n}\cdot\boldsymbol{\sigma}({\boldsymbol{r}};{\bf
u}(\cdot,{\boldsymbol{r}}))=0$ on $\Gamma_{N}$, where $\Gamma_{D}$ and
$\Gamma_{N}$ are disjoint, with $\Gamma_{D}\cup\Gamma_{N}=\Gamma$ and measure
of $\Gamma_{D}>0$, for every ${\boldsymbol{r}}\in U\times U$. Here, ${\bf n}$
is the outward unit normal column vector to $\Gamma_{N}.$

##  2 Estimations of the expected value

To estimate the expected value $\mathcal{I}(\mathcal{L}({\bf u}))$ in (1.3),
we initially approximate ${\bf u}$ by ${\bf u}_{\bf s}$ (${\bf
s}:=(s_{1},s_{2})\in\mathbb{N}^{2}$) which is the solution of (1.1)
corresponding to the truncated $\mu_{s_{1}}$ and $\lambda_{s_{2}}$. Then, with
$U_{i}=[0,1]^{s_{i}}$ being of fixed dimensions $s_{i}$ for $i=1,2,$ we intend
to approximate $\mathcal{I}(\mathcal{L}({\bf u}))$ by

$\mathcal{I}_{\bf s}(\mathcal{L}({\bf u}_{{\bf
s}})):=\int_{U_{2}}\int_{U_{1}}F({\boldsymbol{r}})d{\boldsymbol{r}},\quad{\rm
where}~{}~{}F({\boldsymbol{r}}):=\mathcal{L}\Big{(}\bf{\bf u}_{\bf
s}\big{(}\cdot,{\boldsymbol{r}}-{\bf\frac{1}{2}}\big{)}\Big{)}\,.$ (2.1)

The shifting of the coordinates by ${\bf\frac{1}{2}}$ translates $U_{i}$ to
$\big{[}-\frac{1}{2},\frac{1}{2}\big{]}^{s_{i}}$ for $i=1,2$. The idea behind
using this shift is due to the second step where we approximate the high
dimensional integrals in (2.1) using _deterministic, interlaced high-order
polynomial lattice QMC rules_ of the form

$\mathcal{I}_{\bf s}(\mathcal{L}({\bf u}_{\bf s}))\approx\mathcal{I}_{{\bf
s},{\bf N}}(\mathcal{L}({\bf u}_{\bf
s})):=\frac{1}{N_{1}\,N_{2}}\sum_{k=0}^{N_{2}-1}\sum_{j=0}^{N_{1}-1}F({\boldsymbol{r}}_{j,k})~{}~{}{\rm
with}~{}~{}{\bf N}=(N_{1},N_{2})\,,$ (2.2)

where ${\boldsymbol{r}}_{j,k}=({\boldsymbol{y}}_{j},{\boldsymbol{z}}_{k})$,
and the QMC points
$\\{{\boldsymbol{y}}_{0},\ldots,{\boldsymbol{y}}_{N_{1}-1}\\}\in U_{1}$ and
$\\{{\boldsymbol{z}}_{0},\ldots,{\boldsymbol{z}}_{N_{2}-1}\\}\in U_{2}$ (the
QMC points are always constructed over the unit cube $[0,1]^{s}$ making the
shift by $1/2$ necessary). With $i=1,2,$ we generate the $N_{i}:=b^{m_{i}}$
($m_{i}\geq 1$ is integer and $b$ is prime) QMC points using a _generating
vector_ of polynomials with degree less than $m_{i}$ with its coefficients
taken from a finite field ${\mathbb{Z}}_{b}$; see [13, Section 6] for full
details.

######  Theorem 2.1 ([13, Theorem 6.1])

Let $\boldsymbol{\chi}=(\chi_{j})_{j\geq 1}$ and
$\boldsymbol{\varphi}=(\varphi_{j})_{j\geq 1}$ be two sequences of positive
numbers such that $\sum_{j=1}^{\infty}\chi_{j}^{p}$ and
$\sum_{j=1}^{\infty}\varphi_{j}^{q}$ being finite for some $0<p,\,q<1.$ Let
$\boldsymbol{\chi}_{s_{1}}=(\chi_{j})_{1\leq j\leq s_{1}}$,
$\boldsymbol{\varphi}_{s_{2}}=(\varphi_{j})_{1\leq j\leq s_{2}}$,
$\alpha\,:=\,\lfloor 1/p\rfloor+1$, and $\beta\,:=\,\lfloor 1/q\rfloor+1$.
Assume that there exists a positive constant $c$ such that

$|\partial^{\boldsymbol{\alpha}}_{\boldsymbol{y}}F({\boldsymbol{r}})|\leq
c|\boldsymbol{\alpha}|!\boldsymbol{\chi}_{s_{1}}^{\boldsymbol{\alpha}}\quad{\rm
and}\quad|\partial^{\boldsymbol{\beta}}_{\boldsymbol{z}}F({\boldsymbol{r}})|\leq
c|\boldsymbol{\beta}|!\boldsymbol{\varphi}_{s_{2}}^{\boldsymbol{\beta}},$
(2.3)

for all ${\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}})\in U_{1}\times
U_{2}$, $\boldsymbol{\alpha}\in\\{0,1,\ldots,\alpha\\}^{s_{1}}$, and
$\boldsymbol{\beta}\in\\{0,1,\ldots,\beta\\}^{s_{2}}$. Then one can construct
two interlaced high-order polynomial lattice rules of order $\alpha$ with
$N_{1}=b^{m_{1}}$ points, and of order $\beta$ with $N_{2}=b^{m_{2}}$ points,
with $|m_{1}q-m_{2}p|<1$, using a fast component-by-component (CBC) algorithm,
such that the quadrature error $|\mathcal{I}_{\bf s}(\mathcal{L}({\bf u}_{\bf
s}))-\mathcal{I}_{{\bf s},{\bf N}}(\mathcal{L}({\bf u}_{\bf s}))|\leq
C\,N^{-\frac{1}{p+q}}$, where $N=N_{1}N_{2}$ is the total number of QMC
quadrature points. The constant $C$ depends on $p,q,$ and $b$, but is
independent of $s_{1}$, $s_{2}$, $m_{1}$ and $m_{2}$.

In the next theorem, we state the error from approximating the integral in
(2.1) by a QMC rule in dimension $s_{1}+s_{2}$ directly. In this approach we
combine the different weights arising from simulating $\mu$ and $\lambda$. The
proof follows directly from [12, Theorem 3.1].

###### Theorem 2.2

Let $\boldsymbol{\chi}$ and $\boldsymbol{\varphi}$ be the two sequences
introduced in Theorem 2.1, and let $\gamma\,:=\min(\lfloor 1/p\rfloor,\lfloor
1/q\rfloor)+1.$ For any
$\boldsymbol{\gamma}\in\\{0,1,\ldots,\gamma\\}^{s_{1}+s_{2}}$, assume that $F$
satisfies

$|\partial^{\boldsymbol{\gamma}}_{{\boldsymbol{r}}}{F}({\boldsymbol{r}})|=|\partial^{\boldsymbol{\gamma}_{1},\boldsymbol{\gamma}_{2}}_{{\boldsymbol{y}},{\boldsymbol{z}}}{F}({\boldsymbol{y}},{\boldsymbol{z}})|\leq
c|\boldsymbol{\gamma}|!\boldsymbol{\chi}_{s_{1}}^{\boldsymbol{\gamma}_{1}}\boldsymbol{\varphi}_{s_{2}}^{\boldsymbol{\gamma}_{2}},~{}~{}{\rm
for~{}any}~{}~{}{\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}})\in
U_{1}\times U_{2}=[0,1]^{s_{1}+s_{2}},$ (2.4)

where the vectors $\boldsymbol{\gamma}_{1}$ and $\boldsymbol{\gamma}_{2}$ are
formed from the first $s_{1}$ and last $s_{2}$ components of
$\boldsymbol{\gamma},$ respectively, and the constant $c$ is independent of
${\boldsymbol{r}}$, $s_{1}$, $s_{2}$, and of $p$ and $q.$ Then one can
construct an interlaced polynomial lattice rule of order $\gamma$ with
$N=b^{m}$ points using a fast CBC algorithm so that

$|\mathcal{I}_{\bf s}(\mathcal{L}({\bf u}_{\bf s}))-\mathcal{I}_{{\bf s},{\bf
N}}(\mathcal{L}({\bf u}_{\bf s}))|=\Big{|}\mathcal{I}_{\bf s}(\mathcal{L}({\bf
u}_{\bf s}))-\frac{1}{N}\sum_{n=0}^{N-1}\mathcal{L}({\bf u}_{\bf
s}(\boldsymbol{r}_{n}-\frac{1}{2}))\Big{|}\leq C\,N^{-\min(1/p,1/q)}\,,$

where the generic constant $C$ depends on $b,p$ and $q$, but is independent of
${\bf s}$ and $m.$

In the third step of our error analysis, we approximate the unknown truncated
solution ${\bf u}_{\bf s}$ over the physical domain $\Omega$ by ${\bf
u}_{{{\bf s}},h}\in{\bf V}_{h}$ (see Section 5 for the definition of the
finite element space ${\bf V}_{h}$). However, before summarising these steps
in the upcoming theorem, we introduce the following spaces and notations. Let
$\|{\bf w}\|$ denote the norm of a vector field ${\bf w}$ in ${\bf
L}^{2}(\Omega)=[L^{2}(\Omega)]^{2}$, and let $\|{\bf w}\|_{\bf V}$ and $\|{\bf
w}\|_{\bf H}$ denote the norms of ${\bf w}$ in ${\bf
V}=[H^{1}_{0}(\Omega)]^{2}$ and ${\bf H}=[H^{2}(\Omega)]^{2}$, respectively.
Here, $H^{1}_{0}(\Omega)$ and $H^{2}(\Omega)$ are the usual Sobolev spaces of
scalar-valued functions on $\Omega$. Finally, ${\bf V}^{*}$ is the dual space
of ${\bf V}$ with respect to the inner product in ${\bf L}^{2}(\Omega)$.

In the next theorem, the generic constant $C$ is independent of the finite
element mesh size $h$, the dimensions $s_{1}$ and $s_{2},$ the number of QMC
point $N$, and the Lamé parameter factor $\Lambda$. However, $C$ may depend on
$\Omega$ and other parameters including $p,q,b,$ the upper/lower bounds of
$\mu$ and $\widehat{\lambda}$, and the upper bounds of
$\|\nabla\mu\|_{L^{\infty}(\Omega)}$ and
$\|\nabla\widehat{\lambda}\|_{L^{\infty}(\Omega)}.$ Write
${\boldsymbol{r}}^{\bf s}=({\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})$
for the truncated version of
${\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}})$.

###### Theorem 2.3

Assume that $\Omega$ is convex and (A1)–(A4) are satisfied. For every
${\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}})\in U\times U,$ let ${\bf
u}\in{\bf H}$ be the solution of problem (1.1) and ${\bf u}_{{\bf s},h}$ be
the nonconforming finite element solution, defined as in (5.4) with
${\boldsymbol{r}}^{\bf s}$ in place of ${\boldsymbol{r}}$ for $s_{1},s_{2}\geq
1.$ Let $\alpha\,:=\,\lfloor 1/p\rfloor+1$ and $\beta\,:=\,\lfloor
1/q\rfloor+1$, where $p$ and $q$ are those in (A4). Then one can construct two
interlaced high-order polynomial lattice rules of order $\alpha$ with $N_{1}$
points, and of order $\beta$ with $N_{2}$ points, such that

$|\mathcal{I}(\mathcal{L}({\bf u}))-\mathcal{I}_{{\bf s},{\bf
N}}(\mathcal{L}({\bf u}_{{\bf s},h}))|\leq
C\,\Big{(}s_{1}^{1-1/p}+s_{2}^{1-1/q}+N^{-1/(p+q)}\Big{)}\|{\bf f}\|_{{\bf
V}^{*}}\|\mathcal{L}\|_{{\bf V}^{*}}+C\,h^{2}\|{\bf f}\|\|\mathcal{L}\|\,.$

Further, if we use a QMC rule in dimension $s_{1}+s_{2}$ directly with $N$
points as in Theorem 2.2, then the above estimate remains valid with
$N^{-\min(1/p,1/q)}$ in place of $N^{-1/(p+q)}$.

Proof We start our proof by noticing that the error
$|\mathcal{I}(\mathcal{L}({\bf u}))-\mathcal{I}_{{\bf s},{\bf
N}}(\mathcal{L}({\bf u}_{{\bf s},h}))|$ is bounded by

$|\mathcal{I}(\mathcal{L}({\bf u}))-\mathcal{I}_{{\bf s}}(\mathcal{L}({\bf
u}_{\bf s}))|+|\mathcal{I}_{{\bf s}}(\mathcal{L}({\bf u}_{\bf
s}))-\mathcal{I}_{{\bf s},{\bf N}}(\mathcal{L}({\bf u}_{\bf
s}))|+|\mathcal{I}_{{\bf s},{\bf N}}(\mathcal{L}({\bf u}_{\bf s}-{\bf u}_{{\bf
s},h}))|.$ (2.5)

Since $d{\boldsymbol{y}}$ and $d{\boldsymbol{z}}$ are the uniform probability
measures with i.i.d. uniformly distributed parameters on $U$, by using Theorem
4.1 we obtain

$\mathcal{I}(\mathcal{L}({\bf u}))-\mathcal{I}_{{\bf s}}(\mathcal{L}({\bf
u}_{\bf s}))=\int_{U}\int_{U}\mathcal{L}\Big{(}{\bf
u}(\cdot,{\boldsymbol{r}})-{\bf u}_{{\bf s}}(\cdot,{\boldsymbol{r}}^{\bf
s})\Big{)}\,d{\boldsymbol{r}}\leq
C\,\Big{(}s_{1}^{1-1/p}+s_{2}^{1-1/q}\Big{)}\|{\bf f}\|_{{\bf
V}^{*}}\|\mathcal{L}\|_{{\bf V}^{*}}\,.$

To bound the second term in (2.5), which is the QMC quadrature error, we apply
Theorem 2.1 with $F({\boldsymbol{r}}):=\mathcal{L}\Big{(}\bf{\bf u}_{\bf
s}\big{(}\cdot,{\boldsymbol{r}}-{\bf\frac{1}{2}}\big{)}\Big{)}.$ Since
$|\partial_{{\boldsymbol{r}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}F({\boldsymbol{r}})|\leq\|\mathcal{L}\|_{{\bf
V}*}\,\|\partial_{{\boldsymbol{r}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}_{\bf s}\big{(}\cdot,{\boldsymbol{r}}-{\bf\frac{1}{2}}\big{)}\|_{{\bf V}}$,
the needed regularity conditions in (2.3) are satisfied owing to Theorem 3.3.
Therefore,

$|\mathcal{I}_{{\bf s}}(\mathcal{L}({\bf u}_{\bf s}))-\mathcal{I}_{{\bf
s},{\bf N}}(\mathcal{L}({\bf u}_{\bf s}))|\leq C\,N^{-\frac{1}{p+q}}\|{\bf
f}\|_{{\bf V}^{*}}\|\mathcal{L}\|_{{\bf V}^{*}}\,.$

To estimate the third term in (2.5), we refer to Theorem 5.2 by taking the
vectors ${\boldsymbol{y}}$ and ${\boldsymbol{z}}$ with $y_{j}=z_{k}=0$ for
$j>s_{1}$ and $k>s_{2}$. Hence, the proof of this theorem is complete. $\Box$

##  3 Regularity properties

This section is devoted to showing some useful regularity properties of the
weak solution with respect to both the physical variable ${\boldsymbol{x}}$
and parametric variables ${\boldsymbol{y}}$ and ${\boldsymbol{z}}$. Following
[13], the weak formulation of parameter-dependent equation (1.1) is: Find
${\bf u}(\cdot,{\boldsymbol{r}})\in{\bf V}$ so that

$\mathcal{B}({\boldsymbol{r}};{\bf u}(\cdot,{\boldsymbol{r}}),{\bf
v})=\ell({\bf v}),\quad\text{for all~{}~{}${\bf v}\in{\bf V}$ and for every
${\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}})\in U\times U$,}$ (3.1)

where the linear functional $\ell$ and the bilinear operator $\mathcal{B}$ are
defined by

$\ell({\bf v})=\int_{\Omega}{\bf f}\cdot{\bf
v}\,d{\boldsymbol{x}}\quad\text{and}\quad\mathcal{B}({\boldsymbol{r}};{\bf
u}(\cdot,{\boldsymbol{r}}),{\bf
v})=\int_{\Omega}[2\mu\,{\boldsymbol{\varepsilon}}({\bf
u}(\cdot,{\boldsymbol{r}})):{\boldsymbol{\varepsilon}}({\bf
v})+\lambda\nabla\cdot{\bf u}(\cdot,{\boldsymbol{r}})\nabla\cdot{\bf
v}]\,d{\boldsymbol{x}}.$

The colon operator is the inner product between tensors. Since
$\max\\{\frac{\|\nabla\cdot{\bf
w}\|}{\sqrt{2}},\|{\boldsymbol{\varepsilon}}({\bf w})\|\\}\leq\|\nabla{\bf
w}\|$,

$\mathcal{B}({\boldsymbol{r}};{\bf v},{\bf w})\leq C_{\lambda,\mu}\|\nabla{\bf
v}\|\,\|\nabla{\bf w}\|\leq C_{\lambda,\mu}\|{\bf v}\|_{{\bf V}}\,\|{\bf
w}\|_{{\bf V}},~{}~{}{\rm for}~{}~{}{\bf v},{\bf w}\in{\bf V}.$ (3.2)

So, $\mathcal{B}({\boldsymbol{r}};\cdot,\cdot)$ is bounded over ${\bf
V}\times{\bf V}$. The coercivity property of
$\mathcal{B}({\boldsymbol{r}};\cdot,\cdot)$ on ${\bf V}$ follows from the
nonnegativity assumptions on $\mu$ and $\lambda$ and Korn’s inequality as

$\mathcal{B}({\boldsymbol{r}};{\bf v},{\bf v})\geq
2\mu_{\min}\|{\boldsymbol{\varepsilon}}({\bf v})\|^{2}\geq C\mu_{\min}\|{\bf
v}\|_{{\bf V}}^{2},\quad{\bf v}\in{\bf V}.$ (3.3)

Owing to these two properties, and since $|\ell({\bf v})|\leq\|{\bf f}\|_{{\bf
V}^{*}}\|{\bf v}\|_{{\bf V}}$, an application of the Lax-Milgram theorem
completes the proof of the next theorem.

###### Theorem 3.1

For every $f\in{\bf V}^{*}$, and for every ${\boldsymbol{r}}\in U\times U$,
problem (3.1) has a unique solution.

###  3.1 Regularity in ${\boldsymbol{x}}$

For the convergence of the nonconforming finite element solution, we derive
next some regularity estimates of the parametric solution of (3.1). For the
planar linear elasticity problem with constant Lamé parameters $\lambda$ and
$\mu$, these estimates were proved in [28, Theorem A.1] for domains with
smooth boundary. Later on, some ideas were borrowed from [28] to show these
results for the case of convex polygonal domains, see [7, Lemma 2.2]. Our
results extend the ones in the existing literature to the case of variable
random parameters $\lambda$ and $\mu$. Some ideas from [7, Lemma 2.2] are used
in our proof.

###### Theorem 3.2

Assume that (A2) is satisfied. Then, for every $f\in{\bf V}^{*}$, and for
every ${\boldsymbol{r}}\in U\times U$, the parametric weak solution ${\bf
u}={\bf u}(\cdot,{\boldsymbol{r}})$ of problem (3.1) satisfies the estimate

$\|{\bf u}\|_{{\bf V}}+\Lambda\|\nabla\cdot{\bf u}\|\leq C\,\|{\bf f}\|_{{\bf
V}^{*}}\,.$ (3.4)

Furthermore, if $\Omega$ is a convex domain (or if $\Omega$ has a smooth
boundary) and ${\bf f}\in{\bf L}^{2}(\Omega)$, then ${\bf
u}(\cdot,{\boldsymbol{r}})\in{\bf V}\cap{\bf H}$, which is a strong solution
of (1.1), and the following estimate holds true

$\|{\bf u}\|_{\bf H}+\|(\lambda+\mu)\nabla\cdot{\bf u}\|_{{\bf V}}\leq
C\,\|{\bf f}\|,\quad{\rm for~{}every}~{}~{}{\boldsymbol{r}}\in U\times U\,.$
(3.5)

The constant $C$ in (3.4) depends on $\Omega$, $\mu$ and $\widehat{\lambda}$,
while in (3.5), $C$ depends (linearly) on
$\|\nabla\mu(\cdot,{\boldsymbol{y}})\|_{L^{\infty}(\Omega)}$ and
$\|\nabla\widehat{\lambda}(\cdot,{\boldsymbol{z}})\|_{L^{\infty}(\Omega)}$ as
well. In both (3.4) and (3.5), $C$ is independent of $\Lambda$.

Proof From the coercivity property in (3.3) and the weak formulation in (3.1),
we have

$C\mu_{\min}\|{\bf u}\|_{{\bf V}}^{2}\leq
2\mu_{\min}\|{\boldsymbol{\varepsilon}}({\bf
u})\|^{2}\leq\mathcal{B}({\boldsymbol{r}};{\bf u},{\bf u})=\ell({\bf
u})\leq\|{\bf f}\|_{{\bf V}^{*}}\|{\bf u}\|_{{\bf V}},$ (3.6)

for every ${\boldsymbol{r}}\in U\times U$. Since ${\bf
u}(\cdot,{\boldsymbol{r}})\in{\bf V},$ by [3, Theorem 3.1 and (3.2)], there
exists ${\bf u}^{*}\in{\bf V}$ such that $\nabla\cdot{\bf u}=\nabla\cdot{\bf
u}^{*}$ and $\|{\bf u}^{*}\|_{{\bf V}}\leq C\|\nabla\cdot{\bf u}\|$. Using
this in the weak formulation (3.1), after some calculations, we reach

$\Lambda\widehat{\lambda}_{\min}\|\nabla\cdot{\bf
u}\|^{2}\leq\int_{\Omega}\lambda\nabla\cdot{\bf u}\nabla\cdot{\bf
u}^{*}\,\,d{\boldsymbol{x}}=\ell({\bf
u}^{*})-\int_{\Omega}2\mu{\boldsymbol{\varepsilon}}({\bf
u}):{\boldsymbol{\varepsilon}}({\bf u}^{*})\,\,d{\boldsymbol{x}}\\\ \leq
C\|{\bf f}\|_{{\bf V}^{*}}\|{\bf u}^{*}\|_{{\bf
V}}+2\mu_{\max}\|{\boldsymbol{\varepsilon}}({\bf
u})\|\,\|{\boldsymbol{\varepsilon}}({\bf u}^{*})\|\leq C\Big{(}\|{\bf
f}\|_{{\bf V}^{*}}+2\mu_{\max}\|{\bf u}\|_{{\bf V}}\Big{)}\|\nabla\cdot{\bf
u}\|\,.$ (3.7)

Thus, $\Lambda\widehat{\lambda}_{\min}\|\nabla\cdot{\bf u}\|\leq
C\Big{(}\|{\bf f}\|_{{\bf V}^{*}}+2\mu_{\max}\|{\bf u}\|_{{\bf V}}\Big{)}$.
Combining this estimate with the inequality $\|{\bf u}\|_{{\bf
V}}\leq\frac{C}{\mu_{\min}}\|{\bf f}\|_{{\bf V}^{*}}$, which is due to (3.6),
completes the proof of the regularity estimate in (3.4).

We turn now to the case of a convex domain $\Omega$. Since ${\bf f}\in{\bf
L}^{2}(\Omega)$ and since the bilinear operator
$\mathcal{B}({\boldsymbol{r}};\cdot,\cdot)$ is ${\bf V}-$elliptic for every
${\boldsymbol{r}}\in U\times U$, the solution ${\bf
u}(\cdot,{\boldsymbol{r}})$ of (3.1) is in the space ${\boldsymbol{H}}$ owing
to the elliptic regularity theory, and consequently, ${\bf
u}(\cdot,{\boldsymbol{r}})$ is a strong solution of (1.1). For more details
about strongly elliptic systems and regularity, see [24, Theorem 4.18].

To show the regularity estimate in (3.5), we use first the identity
$\nabla\cdot(g{\bf A})=g\nabla\cdot{\bf A}+\nabla g\cdot{\bf A}$ (with $g\in
W^{1,\infty}(\Omega)$, ${\bf A}=[a_{ij}]$ is an $m$-by-$m$ matrix-valued
function, and so $\nabla g\cdot{\bf A}=[\nabla g\cdot{\bf A}_{i}]$ is an
$m$-by-$1$ column vector where ${\bf A}_{i}$ is the $i$th row of ${\bf A}$),
and notice that

$\nabla\cdot(\boldsymbol{\sigma}({\bf
u},\lambda,\mu))=\lambda\nabla\cdot\Big{(}\nabla\cdot{\bf
u}\,I\Big{)}+2\mu\nabla\cdot{\boldsymbol{\varepsilon}}({\bf
u})+\nabla\lambda\cdot\Big{(}\nabla\cdot{\bf
u}\,I\Big{)}+2\nabla\mu\cdot{\boldsymbol{\varepsilon}}({\bf u})\,.$

However $\nabla\cdot(\nabla\cdot{\bf u}I)=\nabla(\nabla\cdot{\bf u})$ and
$2\nabla\cdot{\boldsymbol{\varepsilon}}({\bf u})=\Delta{\bf
u}+\nabla(\nabla\cdot{\bf u})$, so problem (1.1) can be reformulated as: for
${\boldsymbol{r}}\in U\times U,$

$-\mu\Delta{\bf u}-(\lambda+\mu)\nabla(\nabla\cdot{\bf u})={\bf
f}+\nabla\lambda\,\Big{(}\nabla\cdot{\bf
u}\Big{)}+2\nabla\mu\cdot{\boldsymbol{\varepsilon}}({\bf u})=:\mu\,{\bf G}.$
(3.8)

For later use, it is not hard to observe that

$\|\nabla(\nabla\cdot{\bf u})\|\leq\|\eta^{-1}\Delta{\bf u}\|+\|\eta^{-1}{\bf
G}\|~{}~{}{\rm
with}~{}~{}\eta=\frac{\lambda+\mu}{\mu}=\frac{\lambda}{\mu}+1\,.$ (3.9)

By dividing now both sides of (3.8) by $\mu$ and letting
$p=-\eta\nabla\cdot{\bf u}\in H^{1}(\Omega)$, after some simplifications, we
obtain

$-\Delta{\bf u}+\nabla p={\bf\widetilde{G}}:={\bf
f}/\mu+2\nabla\mu\cdot{\boldsymbol{\varepsilon}}({\bf u})/\mu+\nabla\cdot{\bf
u}\,\Big{(}\lambda/\mu^{2}\Big{)}\nabla\mu\,.$ (3.10)

We borrow now the idea used in [7, Lemma 2.2], specifically (2.19) and (2.20)
therein,

$\exists~{}{\bf u}^{*}\in{\bf H}\cap{\bf V}~{}{\rm
such~{}that}~{}~{}\nabla\cdot{\bf u}=\nabla\cdot{\bf u}^{*}~{}{\rm
and}~{}\|{\bf u}^{*}\|_{\bf H}\leq C_{0}\|\nabla\cdot{\bf u}\|_{{\bf V}},$
(3.11)

where the generic constant $C_{0}$ depends on $\Omega$ only. (This can be
guaranteed by using [3, Theorem 3.1 and (3.2)].) Hence, by setting ${\bf
w}={\bf u}-{\bf u}^{*}$, problem (3.10) can be reformulated in (homogeneous)
Stokes form as: find ${\bf w}\in{\bf H}\cap{\bf V}$ and $p\in H^{1}(\Omega)$
such that

$-\Delta{\bf w}+\nabla p={\bf\widetilde{G}}-\Delta{\bf u}^{*},\quad{\rm
with}~{}~{}\nabla\cdot{\bf w}=0\,.$ (3.12)

Thanks to the regularity estimates for the Stokes problem in [23, Theorem 2],
we have

$\|{\bf w}\|_{\bf H}+\|\nabla p\|\leq C_{0}\|{\bf\widetilde{G}}-\Delta{\bf
u}^{*}\|,~{}~{}\text{assuming that $\Omega$ is convex.}$

If $\Omega$ is a ${\mathcal{C}}^{1,1}$ bounded domain in $\mathbb{R}^{d}$
(with $d\in\\{2,3\\}$), a similar result was proved in [2, Theorem 3].
Furthermore, for a ${\mathcal{C}}^{2}$ class domain, we refer the reader to
[16, Theorem IV.6.1] for similar regularity results. Recall that ${\bf u}={\bf
w}+{\bf u}^{*}$ and $p=-\eta\nabla\cdot{\bf u}$, thence

$\|{\bf u}\|_{\bf H}+\|\nabla(\eta\nabla\cdot{\bf u})\|\leq
C_{0}\Big{(}\|{\bf\widetilde{G}}\|+\|{\bf u}^{*}\|_{\bf H}\Big{)}\,.$

However, from (3.11), the regularity estimate in (3.4), and the pre-estimate
in (3.9), we have

$\|{\bf u}^{*}\|_{\bf H}\leq C_{0}\|\nabla\cdot{\bf u}\|_{{\bf V}}\leq
C_{0}\|\nabla(\nabla\cdot{\bf u})\|+C\Lambda^{-1}\|{\bf f}\|\leq
C_{0}(\|\eta^{-1}\Delta{\bf u}\|+\|\eta^{-1}{\bf G}\|)+C\Lambda^{-1}\|{\bf
f}\|,$

and consequently,

$\|{\bf u}\|_{\bf H}+\|\nabla(\eta\nabla\cdot{\bf u})\|\leq
C_{0}(\|\eta^{-1}\Delta{\bf u}\|+\|\eta^{-1}{\bf G}\|)+C\Lambda^{-1}\|{\bf
f}\|+C\|{\bf\widetilde{G}}\|\,.$

From the estimate in (3.4) together with the embedding of ${\bf
L}^{2}(\Omega)$ into ${\bf V}^{*}$, we have

$\|\eta^{-1}{\bf G}\|\leq
C\Big{(}1+\|\nabla\widehat{\lambda}\|_{L^{\infty}(\Omega)}/\widehat{\lambda}_{\min}\Big{)}\|{\bf
f}\|~{}~{}{\rm and}~{}~{}\|{\bf\widetilde{G}}\|\leq
C\Big{(}1+\widehat{\lambda}_{\max}/\widehat{\lambda}_{\min}\Big{)}\|{\bf
f}\|\,.$

Combining the above two equations, we conclude that

$\|{\bf u}\|_{\bf H}+\|\nabla(\eta\nabla\cdot{\bf u})\|\leq
C_{0}\|\eta^{-1}\Delta{\bf u}\|+C\|{\bf f}\|\leq
C_{0}\frac{\mu_{\max}}{\widehat{\lambda}_{\min}\Lambda}\|\Delta{\bf
u}\|+C\|{\bf f}\|\leq\frac{1}{2}\|{\bf u}\|_{\bf H}+C\|{\bf f}\|,$

for a sufficiently large $\Lambda$. Canceling the common terms yields

$\|{\bf u}\|_{\bf H}+2\|\nabla(\eta\nabla\cdot{\bf u})\|\leq C\|{\bf f}\|\,.$
(3.13)

To finalize the proof of (3.5), it is sufficient to show that

$\|\nabla((\lambda+\mu)\nabla\cdot{\bf u})\|\leq C\|{\bf f}\|\,,$ (3.14)

and then, merging this with (3.13) and the regularity estimate in (3.4). In
fact, the bound in (3.14) is a consequence of the inequality

$\|\nabla((\lambda+\mu)\nabla\cdot{\bf u})\|=\|\nabla(\mu\eta\nabla\cdot{\bf
u})\|\leq C\|\nabla(\eta\nabla\cdot{\bf u})\|+C\|\eta\nabla\cdot{\bf u}\|\leq
C\|\nabla(\eta\nabla\cdot{\bf u})\|+C\Lambda\|\nabla\cdot{\bf u}\|\,,$

and the achieved estimates in (3.13) and (3.4). $\Box$

###  3.2 Regularity in random coefficients

For the QMC error bounds, we prove in the next theorem a mixed partial
derivatives estimate of the parametric displacement ${\bf u}$ with respect to
the random variables $y_{j}$ and $z_{k}$. Here, ${\mathcal{S}}$ denotes the
set of (multi-index) infinite vectors $\boldsymbol{\alpha}=(\alpha_{j})_{j\geq
1}$ with nonnegative integer entries such that
$|\boldsymbol{\alpha}|:=\sum_{j\geq 1}\alpha_{j}<\infty.$ For
$\boldsymbol{\alpha}=(\alpha_{j})_{j\geq 1}$ and
$\boldsymbol{\beta}=(\beta_{j})_{j\geq 1}$ in $\mathcal{S},$ the mixed partial
derivative with respect to ${\boldsymbol{y}}$ and ${\boldsymbol{z}}$ is
denoted by

$\partial_{{\boldsymbol{r}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}:=\partial_{{\boldsymbol{z}}}^{\boldsymbol{\beta}}\partial_{{\boldsymbol{y}}}^{\boldsymbol{\alpha}}=\frac{\partial^{|\boldsymbol{\beta}|}}{\partial_{z_{1}}^{\beta_{1}}\partial_{z_{2}}^{\beta_{2}}\cdots}\frac{\partial^{|\boldsymbol{\alpha}|}}{\partial_{y_{1}}^{\alpha_{1}}\partial_{y_{2}}^{\alpha_{2}}\cdots},\quad{\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}})\,.$

###### Theorem 3.3

Under assumptions (A1) and (A3), for every ${\bf f}\in{\bf V}^{*}$, every
${\boldsymbol{r}}\in U\times U$ and every
$\boldsymbol{\alpha},\boldsymbol{\beta}\in\mathcal{S}$, the solution ${\bf
u}(\cdot,{\boldsymbol{r}})$ of the parametric weak problem (3.1) satisfies

$\big{\|}{\boldsymbol{\varepsilon}}\big{(}\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\big{)}\big{\|}+\Lambda\big{\|}{\boldsymbol{\varepsilon}}\big{(}\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\big{)}\big{\|}\leq(|\boldsymbol{\alpha}|+|\boldsymbol{\beta}|)!\
\widetilde{\bf b}^{\boldsymbol{\alpha}}\widehat{\bf
b}^{\boldsymbol{\beta}}\Big{(}\|{\boldsymbol{\varepsilon}}({\bf
u})\|+\Lambda\|\nabla\cdot{\bf u}\|\Big{)},$ (3.15)

where

$\widetilde{\bf b}^{\boldsymbol{\alpha}}=\prod_{i\geq
1}(\widetilde{b}_{i})^{\alpha_{i}},\quad\widetilde{b}_{j}=\Big{(}\frac{2C_{1}}{\widehat{\lambda}_{\min}}\Big{(}1+\frac{\mu_{\max}}{\mu_{\min}}\Big{)}+\frac{1}{\mu_{\min}}\Big{)}\|\psi_{j}\|_{L^{\infty}(\Omega)},$

and

$\widehat{\bf b}^{\boldsymbol{\beta}}=\prod_{i\geq
1}(\widehat{b}_{i})^{\beta_{i}},\quad\widehat{b}_{j}=\Big{(}\frac{1}{\widehat{\lambda}_{\min}}\Big{(}1+\sqrt{2}C_{1}\frac{\mu_{\max}}{\mu_{\min}}\Big{)}+\frac{1}{\sqrt{2}\mu_{\min}}\Big{)}\|\phi_{j}\|_{L^{\infty}(\Omega)}\,.$

The constant $C_{1}$ occurring in the proof depends on $\Omega$ only.
Consequently,

$\|\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf V}}\leq
C(|\boldsymbol{\alpha}|+|\boldsymbol{\beta}|)!\ \widetilde{\bf
b}^{\boldsymbol{\alpha}}\widehat{\bf b}^{\boldsymbol{\beta}}\|{\bf f}\|_{{\bf
V}^{*}},$ (3.16)

where the constant $C$ depends on $\Omega$, $\mu$ and $\widehat{\lambda}$, but
is independent of $\Lambda$.

Proof Differentiating both sides of (3.1) with respect to the variables
$y_{j}$ and $z_{k}$, we find the following recurrence after some calculations

$\mathcal{B}({\boldsymbol{r}};\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u},{\bf
v})=-\int_{\Omega}\Big{[}2\sum_{\boldsymbol{\alpha}}\alpha_{j}\psi_{j}\,{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha}-{\bf
e}_{j},\boldsymbol{\beta}}{\bf u}):{\boldsymbol{\varepsilon}}({\bf
v})+\Lambda\sum_{\boldsymbol{\beta}}\beta_{k}\phi_{k}\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}-{\bf
e}_{k}}{\bf u}\nabla\cdot{\bf v}\Big{]}\,d{\boldsymbol{x}},$ (3.17)

where $\sum_{\boldsymbol{\alpha}}=\sum_{j,\alpha_{j}\neq 0}$ (that is, the sum
over the nonzero indices of $\boldsymbol{\alpha}$), and ${\bf
e}_{i}\in\mathcal{S}$ denotes the multi-index with entry $1$ in position $i$
and zeros elsewhere. Choose ${\bf
v}=\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}$ in (3.17), then from the coercivity property in (3.3) we have

$2\mu_{\min}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})\|^{2}\leq\mathcal{B}({\boldsymbol{r}};\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u},\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})\leq
2\sum_{\boldsymbol{\alpha}}\alpha_{j}\|\psi_{j}\|_{L^{\infty}(\Omega)}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha}-{\bf
e}_{j},\boldsymbol{\beta}}{\bf
u})\|\,\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})\|\\\
+\Lambda\sum_{\boldsymbol{\beta}}\beta_{k}\,\|\phi_{k}\|_{L^{\infty}(\Omega)}\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}-{\bf
e}_{k}}{\bf
u}\|\,\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\|\,,$

and thus, using the inequality $\|\nabla\cdot{\bf
v}\|\leq\sqrt{2}\,\|{\boldsymbol{\varepsilon}}({\bf v})\|$ ($2$ is the
physical dimension) gives

$\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})\|\leq\frac{1}{\mu_{\min}}\Big{(}\sum_{\boldsymbol{\alpha}}\alpha_{j}\|\psi_{j}\|_{L^{\infty}(\Omega)}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha}-{\bf
e}_{j},\boldsymbol{\beta}}{\bf
u})\|+\Lambda\sum_{\boldsymbol{\beta}}\beta_{k}\,\frac{\|\phi_{k}\|_{L^{\infty}(\Omega)}}{\sqrt{2}}\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}-{\bf
e}_{k}}{\bf u}\|\Big{)}\,,$ (3.18)

for every ${\boldsymbol{r}}\in U\times U$. To proceed in our proof, we argue
as in (3.7) and deduce that there exists
$(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*}\in{\bf V}$ such that
$\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}=\nabla\cdot(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*}$ and
$\|(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*}\|_{{\bf V}}\leq
C_{1}\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\|$. Using this and the inequality
$\|{\boldsymbol{\varepsilon}}((\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*})\|\leq\|(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*}\|_{{\bf V}}$ in the weak formulation (3.17) with ${\bf
v}=(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*}$, we obtain

$\Lambda\widehat{\lambda}_{\min}\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\|^{2}\leq\int_{\Omega}\lambda\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\nabla\cdot(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*}\,\,d{\boldsymbol{x}}\leq
2\sum_{\boldsymbol{\alpha}}\alpha_{j}\|\psi_{j}\|_{L^{\infty}(\Omega)}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha}-{\bf
e}_{j},\boldsymbol{\beta}}{\bf
u})\|\,\|{\boldsymbol{\varepsilon}}((\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*})\|\\\
+\Lambda\sum_{\boldsymbol{\beta}}\beta_{k}\,\|\phi_{k}\|_{L^{\infty}(\Omega)}\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}-{\bf
e}_{k}}{\bf
u}\|\,\|\nabla\cdot(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*}\|+2\mu_{\max}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})\|\,\|{\boldsymbol{\varepsilon}}((\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})^{*})\|\\\
\leq\Big{(}2C_{1}\sum_{\boldsymbol{\alpha}}\alpha_{j}\|\psi_{j}\|_{L^{\infty}(\Omega)}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha}-{\bf
e}_{j},\boldsymbol{\beta}}{\bf
u})\|+\Lambda\sum_{\boldsymbol{\beta}}\beta_{k}\,\|\phi_{k}\|_{L^{\infty}(\Omega)}\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}-{\bf
e}_{k}}{\bf u}\|\\\
+2C_{1}\mu_{\max}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})\|\Big{)}\,\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\|\,.$

Dividing both sides by
$\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\|$ and using (3.18) yield

$\Lambda\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\|\leq\frac{2C_{1}}{\widehat{\lambda}_{\min}}\Big{(}1+\frac{\mu_{\max}}{\mu_{\min}}\Big{)}\sum_{\boldsymbol{\alpha}}\alpha_{j}\|\psi_{j}\|_{L^{\infty}(\Omega)}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha}-{\bf
e}_{j},\boldsymbol{\beta}}{\bf u})\|\\\
+\frac{\Lambda}{\widehat{\lambda}_{\min}}\Big{(}1+\sqrt{2}C_{1}\frac{\mu_{\max}}{\mu_{\min}}\Big{)}\sum_{\boldsymbol{\beta}}\beta_{k}\,\|\phi_{k}\|_{L^{\infty}(\Omega)}\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}-{\bf
e}_{k}}{\bf u}\|\,.$

Combining this estimate with (3.18) leads to

$\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})\|+\Lambda\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\|\leq\sum_{\boldsymbol{\alpha}}\alpha_{j}\widetilde{b}_{j}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha}-{\bf
e}_{j},\boldsymbol{\beta}}{\bf
u})\|+\Lambda\sum_{\boldsymbol{\beta}}\beta_{k}\,\widehat{b}_{k}\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}-{\bf
e}_{k}}{\bf u}\|\,.$ (3.19)

To proceed in our proof, we use the induction hypothesis on
$n:=|\boldsymbol{\alpha}+\boldsymbol{\beta}|$. From (3.19), it is clear that
(3.15) holds true when $|\boldsymbol{\alpha}+\boldsymbol{\beta}|=1$. Now,
assume that (3.15) is true for $|\boldsymbol{\alpha}+\boldsymbol{\beta}|=n$,
and the task is to show it for $|\boldsymbol{\alpha}+\boldsymbol{\beta}|=n+1.$

From (3.19) and the induction hypothesis, we have

$\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u})\|+\Lambda\|\nabla\cdot\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}\|\leq
n!\Big{(}\sum_{\boldsymbol{\alpha}}\alpha_{j}\widetilde{b}_{j}\widetilde{\bf
b}^{\boldsymbol{\alpha}-{\bf e}_{j}}\widehat{\bf
b}^{\boldsymbol{\beta}}\|{\boldsymbol{\varepsilon}}({\bf
u})\|+\Lambda\sum_{\boldsymbol{\beta}}\beta_{k}\,\widehat{b}_{k}\widetilde{\bf
b}^{\boldsymbol{\alpha}}\widehat{\bf b}^{\boldsymbol{\beta}-{\bf
e}_{k}}\|\nabla\cdot{\bf u}\|\Big{)}\\\ \leq n!\widetilde{\bf
b}^{\boldsymbol{\alpha}}\widehat{\bf
b}^{\boldsymbol{\beta}}\Big{(}\sum_{\boldsymbol{\alpha}}\alpha_{j}+\sum_{\boldsymbol{\beta}}\beta_{k}\Big{)}(\|{\boldsymbol{\varepsilon}}({\bf
u})\|+\Lambda\|\nabla\cdot{\bf u}\|)\,.$

Since
$\sum_{\boldsymbol{\alpha}}\alpha_{j}+\sum_{\boldsymbol{\beta}}\beta_{k}=n+1$,
the proof of (3.15) is completed.

Finally, since
$\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))\|\geq
C\,\|\partial_{{\boldsymbol{y}},{\boldsymbol{z}}}^{\boldsymbol{\alpha},\boldsymbol{\beta}}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf V}}$ (by Korn’s
inequality) and since $\|{\boldsymbol{\varepsilon}}({\bf
u})\|+\Lambda\|\nabla\cdot{\bf u}\|\leq\|\nabla{\bf
u}\|+\Lambda\|\nabla\cdot{\bf u}\|\leq C\|{\bf f}\|_{{\bf V}^{*}}$ (by (3.4)),
we obtain (3.16) from (3.15). $\Box$

##  4 A truncated problem and error estimates

In this section, we investigate the error from truncating the first and second
sums in (1.2) at $s_{1}$ and $s_{2}$ terms, respectively, for some
$s_{1},s_{2}\in\mathbb{N}.$ First, let us define the truncated weak
formulation problem: for every ${\boldsymbol{r}}^{\bf s}\in U\times U,$ find
${\bf u}_{\bf s}(\cdot,{\boldsymbol{r}}^{\bf s})\in{\bf V}$, such that

$\mathcal{B}({\boldsymbol{r}}^{\bf s};{\bf u}_{\bf
s}(\cdot,{\boldsymbol{r}}^{\bf s}),{\bf v})=\ell({\bf v}),\qquad\forall~{}{\bf
v}\in{\bf V}.$ (4.1)

By Theorem 3.1, the variational problem (4.1) has a unique solution. In the
next theorem, we estimate $\mathcal{L}({\bf u}-{\bf u}_{\bf s})$. This is
needed for measuring the QMC finite element error in (2.5).

###### Theorem 4.1

Let ${\bf u}$ and ${\bf u}_{\bf s}$ be the solutions of the parametric weak
problems (3.1) and (4.1), respectively. Assume that $|\mathcal{L}({\bf
w})|\leq\|\mathcal{L}\|_{{\bf V}^{*}}\|{\bf w}\|_{\bf V}$ for any ${\bf
w}\in{\bf V}$. Under Assumption (A1), for every ${\bf f}\in{\bf V}^{*}$, every
${\boldsymbol{r}}\in U\times U$, and every ${\bf s}\in\mathbb{N}^{2}$, we have

$|\mathcal{L}({\bf u}(\cdot,{\boldsymbol{r}}))-\mathcal{L}({\bf u}_{\bf
s}(\cdot,{\boldsymbol{r}}^{\bf s}))|\leq C\,\widetilde{C}\,\|{\bf f}\|_{{\bf
V}^{*}}\|\mathcal{L}\|_{{\bf V}^{*}},$ (4.2)

where $C$ depends on $\Omega,$ $\mu$, and $\widehat{\lambda}$, and
$\widetilde{C}=\sum_{j\geq
s_{1}+1}\|\psi_{j}\|_{L^{\infty}(\Omega)}+\sum_{j\geq
s_{2}+1}\|\phi_{j}\|_{L^{\infty}(\Omega)}\,.$ Moreover, if $\sum_{j\geq
s_{1}+1}\|\psi_{j}\|_{L^{\infty}(\Omega)}\leq C_{p}s_{1}^{1-1/p}$ and
$\sum_{j\geq s_{2}+1}\|\phi_{j}\|_{L^{\infty}(\Omega)}\leq C_{q}s_{2}^{1-1/q}$
for some $0<p,\,q<1,$ then $\widetilde{C}\leq
C_{p,q}\Big{(}s_{1}^{1-\frac{1}{p}}+s_{2}^{1-\frac{1}{q}}\Big{)}.$ Note that,
by using the Stechkin inequality, a similar bound can be observed if (A3) and
(A4) are satisfied instead.

Proof From the variational formulations in (3.1) and (4.1), we have

$\mathcal{B}({\boldsymbol{r}}^{\bf s};{\bf u}_{\bf
s}(\cdot,{\boldsymbol{r}}^{\bf s})-{\bf u}(\cdot,{\boldsymbol{r}}),{\bf
v})=\mathcal{B}({\boldsymbol{r}}-{\boldsymbol{r}}^{\bf s};{\bf
u}(\cdot,{\boldsymbol{r}}),{\bf v})\,.$ (4.3)

By following the steps in (3.2) and using the achieved estimate in (3.4), we
have

$|\mathcal{B}({\boldsymbol{r}}-{\boldsymbol{r}}^{\bf s};{\bf
u}(\cdot,{\boldsymbol{r}}),{\bf v})|\leq C\Big{(}\sum_{j\geq
s_{1}+1}\|\psi_{j}\|_{L^{\infty}(\Omega)}\,\|{\bf
u}(\cdot,{\boldsymbol{r}})\|_{{\bf V}}+\Lambda\sum_{j\geq
s_{2}+1}\|\phi_{j}\|_{L^{\infty}(\Omega)}\,\|\nabla\cdot{\bf
u}(\cdot,{\boldsymbol{r}})\|\Big{)}\,\|{\bf v}\|_{{\bf V}}\\\ \leq
C\,\widetilde{C}\Big{(}\|{\bf u}(\cdot,{\boldsymbol{r}})\|_{{\bf
V}}+\Lambda\,\|\nabla\cdot{\bf u}(\cdot,{\boldsymbol{r}})\|\Big{)}\,\|{\bf
v}\|_{{\bf V}}\leq C\,\widetilde{C}\|{\bf f}\|_{{\bf V}^{*}}\|{\bf v}\|_{{\bf
V}}\,.$

Now, by mimicking the proof of (3.4) with (4.3) instead of the weak
formulation in (3.1), and the above estimate instead of $|\ell({\bf
v})|\leq\|{\bf f}\|_{{\bf V}^{*}}\|{\bf v}\|_{{\bf V}},$ we obtain

$\|{\bf u}(\cdot,{\boldsymbol{r}})-{\bf u}_{\bf s}(\cdot,{\boldsymbol{r}}^{\bf
s})\|_{{\bf V}}\leq C\,\widetilde{C}\|{\bf f}\|_{{\bf V}^{*}}\,.$ (4.4)

From the linearity and boundedness properties of $\mathcal{L}$,

$|\mathcal{L}({\bf u}(\cdot,{\boldsymbol{r}})-{\bf u}_{\bf
s}(\cdot,{\boldsymbol{r}}^{\bf s}))|\leq C\|{\bf
u}(\cdot,{\boldsymbol{r}})-{\bf u}_{\bf s}(\cdot,{\boldsymbol{r}}^{\bf
s})\|_{{\bf V}}\|\mathcal{L}\|_{{\bf V}^{*}}\,,$

and hence, to complete the proof, we insert (4.4) in this equation. $\Box$

##  5 Nonconforming FEM

In this section, we develop and analyze a low-order nonconforming Galerkin FEM
for approximating the solution of the model problem (3.1) (and consequently,
of problem (1.1)) over the physical domain $\Omega$, and also investigate the
numerical errors. We introduce a family of regular quasi-uniform
triangulations $\mathcal{T}_{h}$ of the domain $\overline{\Omega}$ and let
$h=\max_{K\in\mathcal{T}_{h}}(h_{K})$, where $h_{K}$ denotes the diameter of
the element $K$. For convenience, we introduce the following notations:
${\mathcal{E}}=\\{E_{1},E_{2},\ldots,E_{q_{h}}\\}$ is the set of all edges of
the interelement boundaries of the elements in ${\mathcal{T}}_{h}$, while
${\mathcal{F}}=\\{F_{1},F_{2},\ldots,F_{b_{h}}\\}$ is the set of all remaining
edges that are on the boundary of $\Omega.$ The corresponding midpoints of the
edges in ${\mathcal{E}}$ and ${\mathcal{F}}$ are denoted by
$\xi_{1},\xi_{2},\ldots,\xi_{q_{h}}$ and
$\eta_{1},\eta_{2},\ldots,\eta_{b_{h}}$, respectively. The space of linear
polynomials is denoted by $P_{1}.$ The nonconforming finite element space is
then defined by

${\bf V}_{h}=\\{{\bf w}_{h}\in{\bf L}^{2}(\Omega):~{}{\bf
w}_{h}|_{K}\in[P_{1}]^{2},~{}\forall~{}K\in\mathcal{T}_{h},\\\ ~{}{\bf
w}_{h}{\rm~{}is~{}continuous~{}at}~{}\xi_{i}~{}{\rm and}~{}{\bf
w}_{h}(\eta_{j})=0,~{}{\rm for}~{}~{}1\leq i\leq q_{h},~{}1\leq j\leq
b_{h}\\}.$

Developing a well-posed nonconforming FEM that maintains the symmetric
structure of the variational formulation of problem (1.1) (which is important
at both a theoretical and practical level) becomes more technically
challenging when the Lamé parameter $\mu$ is allowed to depend on the spatial
variable ${\boldsymbol{x}}$, so an innovative touch is needed here. To
elaborate, from the definitions of $\boldsymbol{\sigma}$ and
${\boldsymbol{\varepsilon}},$ and by using the identity
$\nabla\cdot(\mu\nabla{\bf u})^{T}=\nabla\mu\cdot(\nabla{\bf
u})^{T}+\mu\nabla(\nabla\cdot{\bf u})$ (which holds true provided that the
mixed first derivatives of the entries of ${\bf u}$ are commutative), we
notice that

$\nabla\cdot\boldsymbol{\sigma}({\bf u})=\nabla(\lambda\nabla\cdot{\bf
u})+\nabla\cdot(\mu\nabla{\bf u})+\nabla\cdot(\mu\nabla{\bf u})^{T}\\\
=\nabla(\lambda\nabla\cdot{\bf u})+\nabla\cdot(\mu\nabla{\bf
u})+\nabla\mu\cdot(\nabla{\bf u})^{T}+\mu\nabla(\nabla\cdot{\bf u})\,.$

Therefore, since $\mu\nabla(\nabla\cdot{\bf u})=\nabla(\mu\nabla\cdot{\bf
u})-\nabla\mu\nabla\cdot{\bf u}$, we may reformulate (1.1) by rewriting the
principle part in divergence form as: for every ${\boldsymbol{r}}\in U\times
U$ and every ${\boldsymbol{x}}\in\Omega,$

$-\nabla\cdot(\mu\nabla{\bf u})-\nabla((\lambda+\mu)\nabla\cdot{\bf
u})+\nabla\mu\,\nabla\cdot{\bf u}-\nabla\mu\cdot(\nabla{\bf u})^{T}={\bf
f}\,.$ (5.1)

When the Lamé parameters $\mu$ and $\lambda$ are constants, (5.1) reduces to
the following strong form: $-\mu\Delta{\bf
u}-(\lambda+\mu)\nabla(\nabla\cdot{\bf u})={\bf f}$. The weak formulation and
the nonconforming FEM can then be defined accordingly; see for example [7].
For the weak formulation of (5.1) when $\mu$ and $\lambda$ vary with
${\boldsymbol{x}}$, we take the ${\bf L}^{2}(\Omega)$ inner product with a
test function ${\bf v}\in{\bf V}$, then use the divergence theorem and the
homogeneous Dirichlet boundary condition on ${\bf v}$. The variational form of
problem (5.1) is then defined as: find ${\bf u}(\cdot,{\boldsymbol{r}})\in{\bf
V}$ such that

$\widehat{\mathcal{B}}({\boldsymbol{r}};{\bf u},{\bf v})=\ell({\bf
v}),\quad\text{for all}\quad{\bf v}\in{\bf V},$ (5.2)

where

$\widehat{\mathcal{B}}({\boldsymbol{r}};{\bf u},{\bf
v}):=\int_{\Omega}[\mu\,\nabla{\bf u}:\nabla{\bf
v}+(\mu+\lambda)\nabla\cdot{\bf u}\nabla\cdot{\bf
v}+(\nabla\mu\,\nabla\cdot{\bf u}-\nabla\mu\cdot(\nabla{\bf u})^{T})\cdot{\bf
v}]\,d{\boldsymbol{x}}\,.$

Here, the last term $\int_{\Omega}(\nabla\mu\,\nabla\cdot{\bf
u}-\nabla\mu\cdot(\nabla{\bf u})^{T})\cdot{\bf v}\,d{\boldsymbol{x}}$ of the
integrand equals

$\int_{\Omega}\nabla\mu\cdot({\bf K}_{u_{1}}v_{2}+{\bf
K}_{v_{1}}u_{2})\,d{\boldsymbol{x}},~{}~{}{\rm with}~{}~{}{\bf
u}=\begin{bmatrix}u_{1}\\\ u_{2}\end{bmatrix},~{}~{}{\bf
v}=\begin{bmatrix}v_{1}\\\
v_{2}\end{bmatrix},~{}~{}~{}{\boldsymbol{x}}=(x_{1},x_{2}),~{}~{}{\rm
and}~{}~{}{\bf K}_{g}=\begin{bmatrix}-g_{x_{2}}\\\ g_{x_{1}}\end{bmatrix},$

where the sub index $x_{i}$ denotes the partial derivative with respect to
$x_{i}$ for $i=1,2.$ Then, by applying again the divergence theorem and using
the homogeneous Dirichlet boundary conditions $u_{2}=0$ on $\Gamma$, assuming
that the mixed first order derivatives of $\mu$ exist and commute, we
reformulate the operator $\widehat{\mathcal{B}}({\boldsymbol{r}};\cdot,\cdot)$
in a symmetric form as:

$\widehat{\mathcal{B}}({\boldsymbol{r}};{\bf u},{\bf
v}):=\int_{\Omega}[\mu\,\nabla{\bf u}:\nabla{\bf
v}+(\mu+\lambda)\nabla\cdot{\bf u}\nabla\cdot{\bf v}+\nabla\mu\cdot({\bf
K}_{u_{1}}v_{2}+{\bf K}_{v_{1}}u_{2})]\,d{\boldsymbol{x}}\,.$ (5.3)

Alternatively, $\widehat{\mathcal{B}}({\boldsymbol{r}};{\bf u},{\bf v})$ can
be extracted from $\mathcal{B}({\boldsymbol{r}};{\bf u},{\bf v})$ by
reformulating the latter as

$\mathcal{B}({\boldsymbol{r}};{\bf u},{\bf v})=\int_{\Omega}[\mu\,\nabla{\bf
u}:\nabla{\bf v}+(\mu+\lambda)\nabla\cdot{\bf u}\nabla\cdot{\bf v}-\mu\,{\bf
K}_{u_{1}}\cdot\nabla v_{2}-\mu\,{\bf K}_{v_{1}}\cdot\nabla
u_{2})]\,d{\boldsymbol{x}}\,,$

then applying the divergence theorem to the last two terms to obtain
$\widehat{\mathcal{B}}({\boldsymbol{r}};{\bf u},{\bf v})$, assuming that the
mixed derivatives of $u_{1}$ and of $v_{1}$ exist and commute.

Motivated by (5.2), we define the nonconforming FEM as: find ${\bf
u}_{h}(\cdot,{\boldsymbol{r}})\in{\bf V}_{h}$ such that

$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf u}_{h},{\bf v}_{h})=\ell({\bf
v}_{h}),\quad\text{for all ${\bf v}_{h}\in{\bf V}_{h}$ and for every
${\boldsymbol{r}}\in U\times U,$}$ (5.4)

where the broken bilinear form
$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};\cdot,\cdot)$ is defined as
$\widehat{\mathcal{B}}({\boldsymbol{r}};\cdot,\cdot)$ in (5.3), but with
$\nabla_{h}$ and $\nabla_{h}\cdot$ in place of $\nabla$ and $\nabla\cdot$,
respectively. The discrete gradient operator $\nabla_{h}$ is defined as:
$\nabla_{h}{\bf v}|_{K}=\nabla({\bf v}|_{K})$ for any $K\in\mathcal{T}_{h}$,
and the discrete divergence operator $\nabla_{h}\cdot$ is defined similarly.

The existence and uniqueness of the solution ${\bf
u}_{h}(\cdot,{\boldsymbol{r}})$ of (5.4) is guaranteed if the operator
$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};\cdot,\cdot)$ is positive-definite
on ${\bf V}_{h}$. To show this, we introduce the following broken norms on the
finite dimensional space ${\bf V}_{h}$: for ${\bf v}\in V_{h},$

$\|{\bf v}\|_{h}^{2}=\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla{\bf
v}:\nabla{\bf v}\,d{\boldsymbol{x}}~{}~{}{\rm and}~{}~{}\|{\bf
v}\|_{\mathcal{B}_{h}}^{2}=\sum_{K\in\mathcal{T}_{h}}\int_{K}[\mu\,\nabla{\bf
v}:\nabla{\bf v}+(\mu+\lambda)(\nabla\cdot{\bf v})^{2}]\,d{\boldsymbol{x}}\,.$

There exists a positive constant $B$ (depending on $\Omega$) such that

$\|{\bf v}\|\leq B\,\|{\bf v}\|_{h}\leq B/\sqrt{\mu_{\min}}\|{\bf
v}\|_{\mathcal{B}_{h}},\quad{\rm for~{}all}~{}~{}{\bf v}\in{\bf V}_{h};$ (5.5)

see [4] for the proof of the first inequality. Noting that $\|{\bf v}\|_{h}=0$
means that ${\bf v}$ is a piecewise constant vector on each
$K\in\mathcal{T}_{h}$, and hence, by using the zero boundary conditions
together with continuity at midpoints, it follows that ${\bf v}\equiv{\bf 0}.$

Assume that $\max\\{|\mu_{x_{1}}|,|\mu_{x_{2}}|\\}\leq c_{0}$ for some
positive constant $c_{0}$. By using (5.5), we get

$2\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla\mu\cdot{\bf
K}_{v_{1}}v_{2}|\,d{\boldsymbol{x}}\leq
2c_{0}\sum_{K\in\mathcal{T}_{h}}\int_{K}(|v_{1_{x_{1}}}|+|v_{1_{x_{2}}}|)|v_{2}|\,d{\boldsymbol{x}}\\\
\leq
c_{0}\sum_{K\in\mathcal{T}_{h}}\int_{K}(|v_{1_{x_{1}}}|^{2}+|v_{1_{x_{2}}}|^{2}+2|v_{2}|^{2})\,d{\boldsymbol{x}}\\\
\leq c_{0}\Big{(}\|{\bf v}\|_{h}^{2}+2\|{\bf v}\|^{2}\Big{)}\leq
c_{0}(1+2B^{2})/\mu_{\min}\,\|{\bf v}\|_{\mathcal{B}_{h}}^{2}\,.$

Consequently, for ${\bf v}\in{\bf V}_{h}$,

$(1-c_{0}(1+2B^{2})/\mu_{\min})\,\|{\bf
v}\|_{\mathcal{B}_{h}}^{2}\leq\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf
v},{\bf v})\leq(1+c_{0}(1+2B^{2})/\mu_{\min})\,\|{\bf
v}\|_{\mathcal{B}_{h}}^{2}\,.$ (5.6)

Therefore, $\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};\cdot,\cdot)$ is
positive-definite over ${\bf V}_{h}$ whenever $c_{0}<\mu_{\min}/(1+2B^{2})$.
Noting that
$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};\cdot,\cdot)=\widehat{\mathcal{B}}({\boldsymbol{r}};\cdot,\cdot)$
on ${\bf V}$, we can show with the help of the Poincaré inequality that
$\widehat{\mathcal{B}}({\boldsymbol{r}};\cdot,\cdot)$ is also positive-
definite on ${\bf V}$ for small $c_{0}$, and consequently,
$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};\cdot,\cdot)$ is positive-definite
on ${\bf V}_{h}+{\bf V}$ for small $c_{0}$. The associated energy norm is
$\|\cdot\|_{\widehat{\mathcal{B}}_{h}}=\Big{(}\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};\cdot,\cdot)\Big{)}^{1/2}.$

###  5.1 Convergence analysis

For the error analysis, following [10], we use the projection operator
$\Pi_{h}:{\bf V}\to{\bf V}_{h}$, defined by

$\int_{E_{i}}[\Pi_{h}{\bf v}-{\bf v}]\,ds=0\iff\Pi_{h}{\bf
v}(\xi_{i})=\frac{1}{|E_{i}|}\int_{E_{i}}{\bf v}\,ds,\quad{\rm for}~{}~{}1\leq
i\leq q_{h}.$

Thus, for any $K\in\mathcal{T}_{h}$, $\nabla_{h}\cdot(\Pi_{h}{\bf
v})=|K|^{-1}\int_{K}\nabla\cdot{\bf v}\,d{\boldsymbol{x}}$ on $K$ , see [10]
for the proof. So,

$\nabla\cdot{\bf v}=0\implies\nabla_{h}\cdot(\Pi_{h}{\bf v})=0.$

By a standard Bramble–Hilbert argument, it is known that

$\|{\bf v}-\Pi_{h}{\bf v}\|+h\|\nabla_{h}({\bf v}-\Pi_{h}{\bf v})\|\leq
Ch^{2}\|{\bf v}\|_{\bf H}\,.$ (5.7)

In the next theorem, we use some ideas from [7, 10] to show the convergence of
the nonconforming finite element solution of problem (1.1) over the domain
$\Omega$.

###### Theorem 5.1

Assume that $\Omega$ is convex, ${\bf f}\in{\bf L}^{2}(\Omega)$, and
$\max\\{|\mu_{x_{1}}|,|\mu_{x_{2}}|\\}\leq c_{0}$ on $\Omega\times U$ for a
sufficiently small $c_{0}$. For every ${\boldsymbol{r}}\in U\times U,$ let
${\bf u}$ and ${\bf u}_{h}$ be the solutions of problems (1.1) and (5.4),
respectively. Under Assumption (A3), we have $\|{\bf
u}(\cdot,{\boldsymbol{r}})-{\bf
u}_{h}(\cdot,{\boldsymbol{r}})\|_{\widehat{\mathcal{B}}_{h}}\leq Ch\|{\bf
f}\|$, where $C$ depends on $\Omega$, $\mu$ and $\widehat{\lambda}$, and
linearly on $\|\nabla\mu(\cdot,{\boldsymbol{y}})\|_{L^{\infty}(\Omega)}$ and
$\|\nabla\widehat{\lambda}(\cdot,{\boldsymbol{z}})\|_{L^{\infty}(\Omega)}$,
but is independent of the Lamé parameter factor $\Lambda$.

Proof The bilinear operator
$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};\cdot,\cdot)$ is symmetric and
positive-definite on ${\bf V}_{h}+{\bf V}$, and so, by the second Strang lemma
[27], we have

$\|{\bf u}-{\bf u}_{h}\|_{\widehat{\mathcal{B}}_{h}}\leq\inf_{{\bf
v}_{h}\in{\bf V}_{h}}\|{\bf u}-{\bf
v}_{h}\|_{\widehat{\mathcal{B}}_{h}}+\sup_{{\bf v}_{h}\in{\bf
V}_{h}\backslash\\{0\\}}\frac{|\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf
u}-{\bf u}_{h},{\bf v}_{h})|}{\|{\bf v}_{h}\|_{\widehat{\mathcal{B}}_{h}}}\,.$
(5.8)

Using (5.6), (3.11), the identity $\nabla_{h}\cdot\Pi_{h}{\bf
u}^{*}=\nabla_{h}\cdot\Pi_{h}{\bf u}$ (because $\nabla_{h}\cdot{\bf
u}^{*}=\nabla_{h}\cdot{\bf u}$), the projection bound in (5.7), and the second
regularity estimate in Theorem 3.2, we reach

$\displaystyle\inf_{{\bf v}_{h}\in{\bf V}_{h}}\|{\bf u}-$ $\displaystyle{\bf
v}_{h}\|_{\widehat{\mathcal{B}}_{h}}^{2}\leq C\|{\bf u}-\Pi_{h}{\bf
u}\|_{\widehat{\mathcal{B}}_{h}}^{2}$ (5.9) $\displaystyle\leq C\|{\bf
u}-\Pi_{h}{\bf u}\|_{h}+C\Lambda\|\nabla_{h}\cdot({\bf u}^{*}-\Pi_{h}{\bf
u}^{*})\|+C\|{\bf u}-\Pi_{h}{\bf u}\|_{h}\|{\bf u}-\Pi_{h}{\bf u}\|$
$\displaystyle\leq Ch^{2}\Big{(}\|{\bf u}\|_{\bf H}+\Lambda\|{\bf
u}^{*}\|_{\bf H}\Big{)}^{2}\leq Ch^{2}\Big{(}\|{\bf u}\|_{\bf
H}+\Lambda\|\nabla\cdot{\bf u}\|_{\bf V}\Big{)}^{2}\leq Ch^{2}\|{\bf
f}\|^{2}\,.$

To estimate the second term in (5.8), we make the splitting

$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf u}-{\bf u}_{h},{\bf
v}_{h})=\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf u},{\bf
v}_{h})-\langle{\bf f},{\bf v}_{h}\rangle=J_{1}+J_{2}-J_{3},$ (5.10)

where, after applying the divergence theorem over each element,

$\displaystyle J_{1}$ $\displaystyle=\sum_{K\in\mathcal{T}_{h}}\int_{\partial
K}\mu\,(\nabla{\bf u}\,{\bf n})\cdot{\bf v}_{h}\,d{\boldsymbol{x}},\quad
J_{2}=\sum_{K\in\mathcal{T}_{h}}\int_{\partial K}(\mu+\lambda)(\nabla\cdot{\bf
u})\,{\bf n}\cdot{\bf v}_{h}\,d{\boldsymbol{x}},$ $\displaystyle{\rm
and}~{}~{}J_{3}=\sum_{K\in\mathcal{T}_{h}}\int_{\partial
K}v_{1_{h}}u_{2}\,{\bf K}_{\mu}\cdot{\bf n}\,d{\boldsymbol{x}},$

with ${\bf u}=(u_{1},u_{2})^{T}$ and ${\bf v}_{h}=(v_{1_{h}},v_{2_{h}})^{T}$,
and with ${\bf n}$ denoting the outward unit normal vector.

The task now is to estimate $J_{1}$, $J_{2}$ and $J_{3}.$ For $J_{1},$ define
the jump matrix function $[[{\bf v}_{h}]]$ across the edge $e\in{\mathcal{E}}$
between the two elements $K_{1}$ and $K_{2}$, with outward unit normal vectors
${\bf n}_{1}$ and ${\bf n}_{2}$, respectively, as $[[{\bf v}_{h}]]={\bf
v}_{h}|_{K_{1}}{\bf n}_{1}^{T}+{\bf v}_{h}|_{K_{2}}{\bf n}_{2}^{T}$. If
$e\in{\mathcal{F}}$ is a boundary edge that belongs to $K_{1}$, then $[[{\bf
v}_{h}]]={\bf v}_{h}|_{K_{1}}{\bf n}_{1}^{T}$. Thus, $J_{1}$ can be expanded
as

$J_{1}=\sum_{e\in{\mathcal{E}}\cup{\mathcal{F}}}\int_{e}\mu\,\nabla{\bf
u}:[[{\bf v}_{h}]]\,d{\boldsymbol{x}}\,.$

Due to the continuity of the function ${\bf v}_{h}$ at the midpoint of any
edge $e\in{\mathcal{E}}$ and since ${\bf v}_{h}={\bf 0}$ at the midpoint of
any edge $e\in{\mathcal{F}}$, $[[{\bf v}_{h}]]={\bf 0}$ on the midpoint of
every $e\in{\mathcal{E}}\cup{\mathcal{F}}$. Consequently, by using the
Bramble-Hilbert Lemma and homogeneity arguments, and with the help of the
Trace Theorem (see [8, Section 1.6] or [26]), we get

$J_{1}\leq Ch\Big{(}\sum_{e\in{\mathcal{E}}\cup{\mathcal{F}}}\|\mu\nabla{\bf
u}\|_{H^{1}(\bar{K}_{e})}\|^{2}\Big{)}^{1/2}\Big{(}\sum_{e\in{\mathcal{E}}\cup{\mathcal{F}}}|e|[[{\bf
v}_{h}]]\|_{{\bf L}^{2}(e)}^{2}\Big{)}^{1/2}.$

Here, $\bar{K}_{e}\in K_{e}$ where $K_{e}$ is the set of triangles in
${\mathcal{T}}_{h}$ having $e$ as an edge. Hence, using the inequality
$|e|\,\|[[{\bf v}_{h}]]\|^{2}_{{\bf L}^{2}(e)}\leq Ch^{2}\sum_{K\in
K_{e}}\|\nabla{\bf v}_{h}\|^{2}_{{\bf L}^{2}(K)}$ (because $[[{\bf
v}_{h}]]={\bf 0}$ at the midpoint of $e$) yields

$|J_{1}|\leq Ch\|\mu\nabla{\bf u}\|_{\bf V}\|{\bf
v}_{h}\|_{\mathcal{B}_{h}}\leq Ch\|{\bf u}\|_{\bf H}\|{\bf
v}_{h}\|_{\mathcal{B}_{h}}\,.$

To estimate $J_{2}$, we follow the above procedure, taking into account that
the jump scalar function $[[{\bf v}_{h}]]={\bf v}_{h}|_{K_{1}}\cdot{\bf
n}_{1}+{\bf v}_{h}|_{K_{2}}\cdot{\bf n}_{2}$ across the edge
$e\in{\mathcal{E}}$, while $[[{\bf v}_{h}]]={\bf v}_{h}|_{K_{1}}\cdot{\bf
n}_{1}$ on the boundary edge $e$ that belongs to $K_{1}$. We get $|J_{2}|\leq
Ch\|(\mu+\lambda)\nabla\cdot{\bf u}\|_{\bf V}\|{\bf
v}_{h}\|_{\mathcal{B}_{h}}.$ Again, a similar procedure can also be followed
to estimate $J_{3}$ and obtain

$|J_{3}|\leq Ch\|\nabla\mu\|_{W^{1,\infty}(\Omega)}\|u_{2}\|_{\bf
V}\|v_{1_{h}}\|_{\mathcal{B}_{h}}\leq Ch\|{\bf u}\|_{\bf V}\|{\bf
v}_{h}\|_{\mathcal{B}_{h}}\,,$

where $[[{\bf v}_{h}]]=v_{1_{h}}|_{K_{1}}{\bf n}_{1}+v_{1_{h}}|_{K_{2}}{\bf
n}_{2}$ on $e\in{\mathcal{E}}\cap\overline{K}_{1}\cap\overline{K}_{2}$ and
$[[{\bf v}_{h}]]=v_{1_{h}}|_{K_{1}}{\bf n}_{1}$ on
$e\in{\mathcal{F}}\cap\overline{K}_{1}$. Now, inserting the above estimates in
(5.10) and using Theorem 3.2 and (5.6) yield

$|\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf u}-{\bf u}_{h},{\bf
v}_{h})|\leq Ch\Big{(}\|{\bf u}\|_{\bf H}+\Lambda\|\nabla\cdot{\bf u}\|_{{\bf
V}}\Big{)}\|{\bf v}_{h}\|_{\mathcal{B}_{h}}\leq Ch\|{\bf f}\|\,\|{\bf
v}_{h}\|_{\widehat{\mathcal{B}}_{h}}\,.$ (5.11)

Using this and (5.9) together with (5.8) will complete the proof of this
theorem. $\Box$
For measuring the QMC finite element error in (2.5), we estimate
$\mathcal{L}({\bf u}-{\bf u}_{h})$ in the next theorem.

###### Theorem 5.2

In addition to the assumptions of Theorem 5.1, we assume that
$\mathcal{L}:{\bf L}^{2}(\Omega)\to\mathbb{R}$ is a bounded linear functional
($|\mathcal{L}({\bf w})|\leq\|\mathcal{L}\|\,\|{\bf w}\|$), then

$|\mathcal{L}({\bf u}(\cdot,{\boldsymbol{r}}))-\mathcal{L}({\bf
u}_{h}(\cdot,{\boldsymbol{r}}))|\leq Ch^{2}\|{\bf
f}\|\|\mathcal{L}\|,~{}~{}~{}{\rm
for~{}every}~{}~{}{\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}})\in
U\times U,$

where $C$ depends on $\Omega$, $\mu$, $\widehat{\lambda}$, and linearly on
$\|\nabla\mu(\cdot,{\boldsymbol{y}})\|_{L^{\infty}(\Omega)}$ and
$\|\nabla\widehat{\lambda}(\cdot,{\boldsymbol{z}})\|_{L^{\infty}(\Omega)}$,
but is independent of the Lamé parameter factor $\Lambda$. Furthermore, by
using the equality

$\|{\bf v}\|=\sup_{{\bf w}\in{\bf L}^{2}(\Omega),{\bf w}\neq{\bf
0}}\frac{|\mathcal{L}({\bf v})|}{\|{\bf w}\|},\quad{\rm for~{}any}~{}~{}{\bf
v}\in{\bf L}^{2}(\Omega),~{}~{}{\rm with}~{}~{}\mathcal{L}({\bf
v})=\langle{\bf w},{\bf v}\rangle,$

the above estimate leads to the following optimal convergence bound:

$\|{\bf u}(\cdot,{\boldsymbol{r}})-{\bf u}_{h}(\cdot,{\boldsymbol{r}})\|\leq
C\,h^{2}\|{\bf f}\|,~{}~{}~{}{\rm for~{}every}~{}~{}{\boldsymbol{r}}\in
U\times U\,.$

Proof We start our proof by replacing $\ell$ with $\mathcal{L}$ in (3.1), and
then, by Theorem 3.1, the new parametric variational problem: find ${\bf
u}_{\mathcal{L}}\in{\bf V}$ such that

$\mathcal{B}({\boldsymbol{r}};{\bf u}_{\mathcal{L}},{\bf v})=\mathcal{L}({\bf
v})\quad{\rm for~{}all}~{}~{}{\bf v}\in{\bf V},~{}~{}{\rm
for~{}every}~{}~{}{\boldsymbol{r}}\in U\times U,$ (5.12)

has a unique solution. Furthermore, using Theorem 3.2 (with ${\bf
u}_{\mathcal{L}}$ in place of ${\bf u}$), we deduce that

$\|{\bf u}_{\mathcal{L}}(\cdot,{\boldsymbol{r}})\|_{{\bf
V}}+\Lambda\|\nabla\cdot{\bf u}_{\mathcal{L}}(\cdot,{\boldsymbol{r}})\|\leq
C\|\mathcal{L}\|_{{\bf V}^{*}}\,.$ (5.13)

Then by proceeding as in the proof of (3.5) in Theorem 3.2 (with ${\bf
u}_{\mathcal{L}}$ in place of ${\bf u}$), and using that the linear functional
$\mathcal{L}$ is bounded in ${\bf L}^{2}(\Omega)$, we conclude that

$\|{\bf u}_{\mathcal{L}}(\cdot,{\boldsymbol{r}})\|_{\bf
H}+\|(\lambda+\mu)\nabla\cdot{\bf
u}_{\mathcal{L}}(\cdot,{\boldsymbol{r}})\|_{{\bf V}}\leq C\|\mathcal{L}\|.$
(5.14)

The nonconforming finite element solution of ${\bf u}_{\mathcal{L}}$ is
defined by: Find ${\bf\Theta}_{h}\in{\bf V}_{h}$ such that

$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf\Theta}_{h},{\bf
v}_{h})=\mathcal{L}({\bf v}_{h}),\quad\text{for all}\quad{\bf v}_{h}\in{\bf
V}_{h}\,.$ (5.15)

By following the proof of Theorem 5.1 and using (5.14), we obtain the
following error estimate

$\|{\bf
u}_{\mathcal{L}}(\cdot,{\boldsymbol{r}})-{\bf\Theta}_{h}(\cdot,{\boldsymbol{r}})\|_{\widehat{\mathcal{B}}_{h}}\leq
Ch\|{\bf u}_{\mathcal{L}}(\cdot,{\boldsymbol{r}})\|_{\bf H}\leq
Ch\|\mathcal{L}\|.$ (5.16)

From the linearity property of $\mathcal{L}$, (5.12), (5.15), and the
symmetric property of
$\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};\cdot,\cdot)$,

$\mathcal{L}({\bf u}(\cdot,{\boldsymbol{r}})-{\bf
u}_{h}(\cdot,{\boldsymbol{r}}))=\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf
u}_{\mathcal{L}},{\bf
u})-\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf\Theta}_{h},{\bf
u}_{h})=G_{1}+G_{2}+G_{3},$ (5.17)

where

$G_{1}=\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf
u}_{\mathcal{L}}-{\bf\Theta}_{h},{\bf u}-\Pi_{h}{\bf
u})+\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf u}-{\bf
u}_{h},{\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}})+\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf u}-{\bf
u}_{h},{\bf u}_{\mathcal{L}}-\Pi_{h}{\bf u}_{\mathcal{L}}),$
$G_{2}=\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf u}-{\bf
u}_{h},\Pi_{h}{\bf u}_{\mathcal{L}})\quad{\rm and}\quad
G_{3}=\widehat{\mathcal{B}}_{h}({\boldsymbol{r}};{\bf
u}_{\mathcal{L}}-{\bf\Theta}_{h},\Pi_{h}{\bf u})\,.$

By the Cauchy-Schwarz inequality, Theorem 5.1, (5.16), the bound in (5.9) in
addition to a similar bound with ${\bf u}_{\mathcal{L}}$ in place of ${\bf u}$
where the regularity estimate in (5.14) is needed, we reach

$|G_{1}|\leq\|{\bf
u}_{\mathcal{L}}-{\bf\Theta}_{h}\|_{\widehat{\mathcal{B}}_{h}}\|{\bf
u}-\Pi_{h}{\bf u}\|_{\widehat{\mathcal{B}}_{h}}+\|{\bf u}-{\bf
u}_{h}\|_{\widehat{\mathcal{B}}_{h}}(\|{\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}}\|_{\widehat{\mathcal{B}}_{h}}+\|{\bf
u}_{\mathcal{L}}-\Pi_{h}{\bf
u}_{\mathcal{L}}\|_{\widehat{\mathcal{B}}_{h}})\leq Ch^{2}\|{\bf
f}\|\,\|\mathcal{L}\|.$

To estimate $G_{2}$, we follow the approach used for bounding (5.10), that is,
the derivation between (5.10) and (5.11). To obtain a sharp bound, it is
necessary to use the identity $[[\Pi_{h}{\bf u}_{\mathcal{L}}]]=[[\Pi_{h}{\bf
u}_{\mathcal{L}}-{\bf u}_{\mathcal{L}}]]$ across the skeleton of the finite
element mesh, which is due to the continuity of ${\bf u}_{\mathcal{L}}$ on
$\Omega$ and also to the fact that ${\bf
u}_{\mathcal{L}}(\cdot,{\boldsymbol{r}})={\bf 0}$ on $\Gamma$. In conclusion,
we get

$|G_{2}|\leq Ch\Big{(}\|{\bf u}\|_{\bf H}+\Lambda\|\nabla\cdot{\bf u}\|_{{\bf
V}}\Big{)}\|\Pi_{h}{\bf u}_{\mathcal{L}}-{\bf u}_{\mathcal{L}}\|_{h}\leq
Ch^{2}\|{\bf f}\|\,\|\mathcal{L}\|\,.$

By arguing as above, we obtain a similar estimate of $G_{3}$ with ${\bf
u}_{\mathcal{L}}$ and ${\bf u}$ in place of ${\bf u}$ and ${\bf
u}_{\mathcal{L}}$, respectively, and so, $|G_{3}|\leq Ch^{2}\|{\bf
f}\|\,\|\mathcal{L}\|.$ Finally, inserting the achieved estimates of $G_{1}$,
$G_{2}$ and $G_{3}$ in (5.17) will complete the proof. $\Box$

##  6 Numerical experiments

In this section, we illustrate the theoretical convergence results using four
examples. The source code is hosted on GitHub [25].

Figure 1: The numbers in blue give the values of $j=\rho(k,l)$ if
$n=4$.$k$$l$$1$$2$$3$$4$$1$$2$$3$$4$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$ Table 1:
Errors and convergence rates for the problem in Example 1 when $\Lambda=1$.
Conforming elements
---
$h$ | DoF |  $\mathbf{L}^{2}$ error |  $\mathbf{H}^{1}$ error
0.262 | 552 | 1.22e-01 |  | 1.09e+00 |
0.131 | 2330 | 3.11e-02 | 1.976 | 5.43e-01 | 1.008
0.065 | 9570 | 7.82e-03 | 1.992 | 2.71e-01 | 1.002
0.033 | 38786 | 1.96e-03 | 1.998 | 1.36e-01 | 1.001
0.016 | 156162 | 4.90e-04 | 1.999 | 6.78e-02 | 1.000
Nonconforming elements
$h$ | DoF |  $\mathbf{L}^{2}$ error |  $\mathbf{H}^{1}$ error
0.262 | 1778 | 4.20e-02 |  | 1.07e+00 |
0.131 | 7240 | 1.06e-02 | 1.992 | 5.34e-01 | 0.997
0.065 | 29216 | 2.64e-03 | 1.997 | 2.67e-01 | 0.999
0.033 | 117376 | 6.62e-04 | 1.999 | 1.34e-01 | 1.000
0.016 | 470528 | 1.65e-04 | 2.000 | 6.68e-02 | 1.000
Table 2: Errors and convergence rates for the problem in Example 1 when
$\Lambda=1000$. Conforming elements
---
$h$ | DoF |  $\mathbf{L}^{2}$ error |  $\mathbf{H}^{1}$ error
0.262 | 552 | 1.48e+00 |  | 3.65e+00 |
0.131 | 2330 | 5.80e-01 | 1.354 | 1.62e+00 | 1.172
0.065 | 9570 | 1.81e-01 | 1.681 | 6.34e-01 | 1.355
0.033 | 38786 | 5.04e-02 | 1.842 | 2.47e-01 | 1.361
0.016 | 156162 | 1.34e-02 | 1.915 | 9.86e-02 | 1.323
Nonconforming elements
$h$ | DoF |  $\mathbf{L}^{2}$ error |  $\mathbf{H}^{1}$ error
0.262 | 1778 | 4.01e-02 |  | 1.01e+00 |
0.131 | 7240 | 1.01e-02 | 1.986 | 5.09e-01 | 0.997
0.065 | 29216 | 2.54e-03 | 1.995 | 2.54e-01 | 0.999
0.033 | 117376 | 6.35e-04 | 1.998 | 1.27e-01 | 1.000
0.016 | 470528 | 1.59e-04 | 1.999 | 6.36e-02 | 1.000

#### Example 1.

We consider a purely deterministic problem with Lamé parameters $\mu=\mu_{0}$
and $\lambda=\lambda_{0}$, in effect taking the random variables
${\boldsymbol{y}}={\boldsymbol{z}}=0$ in the expansions (1.2). Based on
Theorems 5.1 and 5.2, we expect to observe optimal order convergence rates in
both ${\bf H}^{1}(\Omega)$ and ${\bf L}^{2}(\Omega)$-norms, uniformly in
$\Lambda$, provided that all the imposed assumptions are satisfied. We chose
$\lambda_{0}(\boldsymbol{x})=\Lambda(1+\tfrac{1}{2}\sin(2x_{1}))$ and
$\mu_{0}(\boldsymbol{x})=1+x_{1}+x_{2}$ for ${\boldsymbol{x}}=(x_{1},x_{2})$
in the square $\Omega=(0,\pi)^{2}$, and chose the body force so that the exact
solution is

${\bf u}(\boldsymbol{x})=\begin{bmatrix}u_{1}(x_{1},x_{2})\\\
u_{2}(x_{1},x_{2})\end{bmatrix}=\begin{bmatrix}\bigl{(}\cos(2x_{1})-1\bigr{)}\sin(2x_{2})\\\
\bigl{(}1-\cos(2x_{2})\bigr{)}\sin(2x_{1})\end{bmatrix}+\frac{\sin(x_{1})\sin(x_{2})}{\Lambda}\begin{bmatrix}1\\\
1\end{bmatrix}.$

Notice ${\bf u}=\mathbf{0}$ on $\partial\Omega$. Table 1 compares the errors
and convergence rates of our nonconforming scheme with those using the
standard conforming, piecewise-linear Galerkin finite element method, when
$\Lambda=1$. We used a family of unstructured triangulations constructed by
uniform refinement of an initial coarse mesh. As expected, we observe
convergence of order $h^{2}$ in the $\mathbf{L}^{2}$ norm and order $h$ in the
(broken) $\mathbf{H}^{1}$ norm. For a given mesh size $h$, the nonconforming
scheme requires more than three times as many degrees of freedom as the
conforming method, but with an $\mathbf{L}^{2}$ error about a third as large.
The errors in the (broken) $\mathbf{H}^{1}$ norm are essentially the same for
both methods. However, when $\Lambda=1000$ (a nearly incompressible case), as
expected, we see in Table 2 that the accuracy of the conforming method is much
poorer, whereas the errors for the nonconforming scheme are almost unchanged.

Figure 2: Instance of $\lambda$ from Example 2 with $\alpha=2$, $\Lambda=1$
and $s_{1}=253$.

For the remaining examples, we chose the _unit_ square $\Omega=(0,1)^{2}$ as
our spatial domain. To construct random coefficients, let $\alpha>1$ and
define

$\psi_{kl}({\boldsymbol{x}})=\phi_{kl}({\boldsymbol{x}})=\frac{\sin(k\pi
x_{1})\,\sin(l\pi x_{2})}{M_{\alpha}(k+l)^{2\alpha}}\quad\text{where}\quad
M_{\alpha}=\lim_{n\to\infty}\sum_{l=1}^{n}\sum_{k=1}^{n-l}\frac{1}{(k+l)^{2\alpha}}=\sum_{j=1}^{\infty}\frac{j-1}{j^{2\alpha}}.$

Notice $M_{\alpha}=\zeta(2\alpha-1)-\zeta(2\alpha)$ in terms of the Riemann
zeta function. By defining $\rho(k,l)$ as shown in Figure 1, we can switch to
single-index labelling, putting

$z_{j}=z_{kl}\quad\text{and}\quad\psi_{j}=\phi_{j}=\psi_{kl}\quad\text{for
$j=\rho(k,l)$.}$

We find that

$1+\smash[t]{\sum_{i=1}^{m-2}i\leq\rho(k,l)\leq\sum_{i=1}^{m-1}i=\frac{m(m-1)}{2}}\quad\text{if
$k+l=m$,}$

and in particular if $j=\rho(k,l)$ then $2j\leq(k+l)^{2}$ so

$\|\psi_{j}\|_{L^{\infty}(\Omega)}\leq\smash[b]{\frac{1}{M_{\alpha}(k+l)^{2\alpha}}\leq\frac{(2j)^{-\alpha}}{M_{\alpha}}}$

and

$\|\nabla\psi_{j}\|_{L^{\infty}(\Omega)}\leq\smash[t]{\frac{\pi(k^{2}+l^{2})^{1/2}}{M_{\alpha}(k+l)^{2\alpha}}\leq\frac{\pi(k+l)^{1-2\alpha}}{M_{\alpha}}\leq\frac{\pi(2j)^{-(\alpha-\frac{1}{2})}}{M_{\alpha}}.}$

Consider the sum

$S_{n,\alpha}({\boldsymbol{x}})=\smash[t]{\sum_{l=1}^{n}\sum_{k=1}^{n-l}z_{kl}\psi_{kl}({\boldsymbol{x}})=\sum_{k=1}^{M}\sum_{l=1}^{M}a_{kl}\sin(k\pi
x_{1})\sin(l\pi x_{2})},$

where $a_{kl}=M_{\alpha}^{-1}(k+l)^{-2\alpha}$ if $1\leq k\leq n-l$ and $1\leq
l\leq n$, and zero otherwise. Since $|z_{kl}|\leq 1/2$, $-1/2\leq
S_{n,\alpha}({\boldsymbol{x}})\leq 1/2$ for all $n$ and ${\boldsymbol{x}}$, so
our assumptions (A1)–(A3) hold for $\alpha>3/2$, and (A4) holds for
$1/\alpha<p<1$ (and $q=p$). In all the computations described below,
$\alpha=2$.

Table 3: Errors and convergence rates for $\mathcal{I}_{N}\mathcal{L}_{h}$ in Example 2 for $\Lambda=1$ (left) and $\Lambda=1,000$ (right). Thus, $\mu=\mu_{0}$ is deterministic and $\lambda$ is random. $\Lambda=1$ | $\Lambda=1000$
---|---
$N$ | $\mathcal{I}_{N}\mathcal{L}_{h}$ | Error | Rate | $\mathcal{I}_{N}\mathcal{L}_{h}$ | Error | Rate
16 | -0.4024722078 | 1.74e-05 |  | -0.0031572258 | 9.85e-07 |
32 | -0.4024592212 | 4.38e-06 | 1.986 | -0.0031565029 | 2.62e-07 | 1.912
64 | -0.4024559125 | 1.07e-06 | 2.028 | -0.0031563107 | 6.95e-08 | 1.913
128 | -0.4024551004 | 2.63e-07 | 2.032 | -0.0031562566 | 1.54e-08 | 2.173

By applying a 2-dimensional fast sine transform to the array $a_{kl}$, we
obtained the values of $S_{n,\alpha}$ on an $M\times M$ grid over the unit
square $\Omega$. Bivariate cubic spline interpolation was then used to compute
$S_{n,\alpha}({\boldsymbol{x}})$ at each quadrature point ${\boldsymbol{x}}$
required by the finite element solver. We chose $M=256$ for computing
$\lambda$ and $\mu$, and $M=512$ for computing $\nabla\mu$ (the latter
requiring mixed fast sine and cosine transforms). Code profiling showed that
for each ${\boldsymbol{r}}=({\boldsymbol{y}},{\boldsymbol{z}})$, the CPU time
to compute ${\bf u}_{h}(\cdot,{\boldsymbol{r}})$ was split roughly as follows:
40% for constructing the interpolants, 40% for assembling the linear system,
and 20% for solving the linear system. We also observed that naive evaluation
of $S_{n,\alpha}$ was about 20 times slower than using fast transforms and
interpolation.

Table 4: Errors and convergence rates for the expected value of $\mathcal{L}_{h}$ in Example 3 for $\Lambda=1$ (left) and $\Lambda=1000$ (right). Thus, $\lambda=\lambda_{0}$ is deterministic and $\mu$ is random. $\Lambda=1$ | $\Lambda=1000$
---|---
$N$ | $\mathcal{I}_{N}\mathcal{L}_{h}$ | Error | Rate | $\mathcal{I}_{N}\mathcal{L}_{h}$ | Error | Rate
16 | -0.6984920721 | 1.73e-04 |  | -0.0036268972 | 1.07e-08 |
32 | -0.6983612152 | 4.21e-05 | 2.038 | -0.0036268891 | 2.67e-09 | 2.005
64 | -0.6983301483 | 1.11e-05 | 1.930 | -0.0036268871 | 6.60e-10 | 2.015
128 | -0.6983216359 | 2.54e-06 | 2.122 | -0.0036268866 | 1.59e-10 | 2.057

For the linear functional, we chose the average displacement in the
$x_{2}$-direction, that is, $\mathcal{L}({\bf
u})=\int_{\Omega}u_{2}\,d{\boldsymbol{x}}$, and for the source term, ${\bf
f}({\boldsymbol{x}})=[1-x_{2}^{2},2x_{1}-20]$. Starting from an unstructured
triangular mesh on $\Omega$, we performed three uniform refinements to obtain
a family of four meshes, the finest of which had $h=0.0190$ and $30,848$ free
nodes. For each ${\boldsymbol{r}}$ and each mesh, we computed the finite
element solution ${\bf u}_{h}(\cdot,{\boldsymbol{r}})$ and evaluated
$\mathcal{L}({\bf u}_{h}(\cdot,{\boldsymbol{r}}))=\mathcal{L}({\bf
u}(\cdot,{\boldsymbol{r}}))+O(h^{2})$. Richardson extrapolation then yielded a
value $\mathcal{L}_{h}({\boldsymbol{r}})=\mathcal{L}({\bf
u}(\cdot,{\boldsymbol{r}}))+O(h^{8})$. The finite element equations were
solved via the conjugate gradient method, preconditioned by the stiffness
matrix $\mathbf{P}$ obtained from using constant coefficients
$\lambda_{0}({\boldsymbol{x}})=\Lambda$ and $\mu_{0}({\boldsymbol{x}})=1$. In
this way, after computing the sparse Cholesky factorization of $\mathbf{P}$
once for a given mesh, only around ten iterations were needed to compute ${\bf
u}_{h}(\cdot,{\boldsymbol{r}})$ for each choice of ${\boldsymbol{r}}$.

#### Example 2.

Consider a deterministic coefficient
$\mu({\boldsymbol{x}})=\mu_{0}({\boldsymbol{x}})=1+x_{1}+x_{2}$, together with
a random coefficient $\lambda$ defined by taking
$\widehat{\lambda}_{0}({\boldsymbol{x}})=1$ in (1.2). In our computations, we
worked with the truncated expansion
$\lambda_{s_{2}}({\boldsymbol{x}},{\boldsymbol{z}})=\Lambda[1+S_{n,\alpha}({\boldsymbol{x}})]$
for $n=22$, so that $s_{2}=n(n+1)/2=253$, and with the values $\Lambda=1$ and
$\Lambda=1000$. Figure 2 shows the contour plot of $\lambda_{s_{2}}$ for a
particular choice of ${\boldsymbol{z}}$. Keeping $s_{2}$ and $h$ fixed, Table
3 shows the behaviour of the error in the expected value
$\mathcal{I}_{N}\mathcal{L}_{h}$ as we increase the number $N$ ($N=N_{2}$) of
the high-order QMC points (generated by a Python package in [17]). The error
was calculated by comparison to the reference value obtained with $N=512$.
Consistent with Theorem 2.3 (or Theorem 2.2) for $q>1/\alpha=1/2$, we observe
convergence of order $N^{-1/q}\approx N^{-2}$ for both choices of $\Lambda$.

#### Example 3.

Next, consider a deterministic
$\lambda({\boldsymbol{x}})=\lambda_{0}({\boldsymbol{x}})=\Lambda[1+\tfrac{1}{2}\sin(2\pi
x_{1})]$ and a random $\mu$, putting $\mu_{0}({\boldsymbol{x}})=1$ in (1.2).
The finite element and QMC approximations were performed as in Example 2, with
$s_{1}=253$. In Table 4, we once again observe convergence of order $N^{-2}$
($N=N_{1}$), both for $\Lambda=1$ and $\Lambda=1000$. The error was calculated
by comparison to the reference value obtained with $N=512$. Again, consistent
with Theorem 2.3 (or Theorem 2.2) for $p>1/\alpha=1/2$, we observe a
convergence of order $N^{-1/p}\approx N^{-2}$ for both choices of $\Lambda$.

Table 5: Errors and convergence rates for the expected value of $\mathcal{L}_{h}$ in Example 4 for $\Lambda=1$ (left) and $\Lambda=1000$ (right). Thus, both $\lambda$ and $\mu$ are random. $\Lambda=1$ | $\Lambda=1000$
---|---
$N$ | $\mathcal{I}_{N}\mathcal{L}_{h}$ | Error | Rate | $\mathcal{I}_{N}\mathcal{L}_{h}$ | Error | Rate
16 | -0.6942608771 | 3.39e-04 |  | -0.0031639526 | 1.01e-06 |
32 | -0.6941636764 | 2.42e-04 | 0.488 | -0.0031632154 | 2.71e-07 | 1.896
64 | -0.6939439764 | 2.19e-05 | 3.465 | -0.0031630179 | 7.34e-08 | 1.884
128 | -0.6939282911 | 6.20e-06 | 1.820 | -0.0031629609 | 1.64e-08 | 2.161

#### Example 4.

Finally, consider both $\lambda$ and $\mu$ random. This time, we chose $n=15$
so that $s_{1}=s_{2}=n(n+1)/2=120$, with QMC points in dimension
$s_{1}+s_{2}=240$. The error was calculated by comparison to the reference
value obtained with $N=1024$. The results, shown in Table 5, again exhibit
order $N^{-2}$ ($\approx N^{-\min(1/p,1/q)}=N^{-1/p}=N^{-1/q})$ convergence,
which confirms the predicted theoretical convergence results in Theorem 2.3
(or Theorem 2.2).

## References

  * [1] M. Ainsworth and C. Parker, Unlocking the secrets of locking: Finite element analysis in planar linear elasticity, Computer Methods in Applied Mechanics and Engineering, 395, 115034 (2022).
  * [2] C. Amrouche and V. Girault, On the existence and regularity of the solution of Stokes problem in arbitrary dimension, Proc. Japan Acad., 67, 171–175 (1991).
  * [3] D. N. Arnold, L. R. Scott and M. Vogelius, Regular inversion of the divergence operator with Dirichlet boundary conditions on a polygon, Ann. Scuola Norm. Sup. Pisa Cl. Sei., 15, 169–192 (1988).
  * [4] S. C. Brenner, Poincaré-Friedrichs inequalities for piecewise H1 functions, SIAM J. Numer. Anal., 41, 306– 324 (2003).
  * [5] M. Botti, D. A. Di Pietro and A. Guglielmana, A low-order nonconforming method for linear elasticity on general meshes, Computer Methods in Applied Mechanics and Engineering, 354, 96–118 (2019).
  * [6] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Elasticity Theory, Cambridge University Press, New York, 2007.
  * [7] S. C. Brenner and Li-Yeng Sung, Linear finite element methods for planar linear elasticity, Math. Comp., 59, 321–338 (1992).
  * [8] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Third Edition, Springer, 2008.
  * [9] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods. Springer-Verlag, New York, 1991.
  * [10] M. Crouzeix and P. A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations, I, RAIRO R-3, 33–75 (1973).
  * [11] J. Dick, Walsh spaces containing smooth functions and Quasi-Monte Carlo rules of arbitrary high order, SIAM J. Numer. Anal., 46, 1519–1553 (2008).
  * [12] J. Dick, F.Y. Kuo, Q.T. Le Gia, D. Nuyens and C. Schwab, Higher order QMC Petrov-Galerkin discretization for affine parametric operator equations with random field inputs, SIAM J. Numer. Anal., 52, 2676–2702 (2014).
  * [13] J. Dick, Q.T. Le Gia, K. Mustapha, and T. Tran, Quasi-Monte Carlo sparse grid Galerkin finite element methods for linear elasticity equations with uncertainties, (2023). https://arxiv.org/pdf/2310.06187.pdf
  * [14] M. Eigel, C. J. Gittelson, Ch. Schwab and E. Zander, Adaptive stochastic Galerkin FEM, Comput. Methods Appl. Mech. Eng., 270, 247–269 (2014).
  * [15] R. S. Falk, Nonconforming finite element methods for the equations of linear elasticity, Math. Comp., 57, 529–550 (1991).
  * [16] G.P. Galdi, An Introduction to the Mathematical Theory of the Navier-Stokes Equations, Second Edition, Springer Monographs in Mathematics, 2011.
  * [17] R. N. Gantner and Ch. Schwab, Computational Higher-Order Quasi-Monte Carlo Integration. Tech. Report 2014-25, Seminar for Applied Mathematics, ETH Zürich, 2014.
  * [18] T. Goda, Good interlaced polynomial lattice rules for numerical integration in weighted Walsh spaces, J. Comput. Appl. Math., 285, 279–294 (2015).
  * [19] T. Goda and J. Dick, Construction of interlaced scrambled polynomial lattice rules of arbitrary high order, Foundations of Computational Mathematics, 15, 1245–1278 (2015).
  * [20] V. H. Hoang, T. C. Nguyen and B. Xia, Polynomial approximations of a class of stochastic multiscale elasticity problems, Z. Angew. Math. Phys., 67:78, (2016).
  * [21] A. Khan, C. E. Powell and D. J. Silvester, Robust preconditioning for stochastic Galerkin formulations of parameter-dependent nearly incompressible elasticity equations, SIAM J. Sci. Comput. 41, A402–A421 (2019).
  * [22] A. Khan, A. Bespalov, C. E. Powell and D. J. Silvester, Robust a posteriori error estimation for parameter-dependent linear elasticity equations, Math. Comp., 90, 613–636 (2021).
  * [23] R. B. Kellogg and J. E. Osborn, A regularity result for the Stokes problem on a convex polygon, J. Funct. Anal., 21, 397–431 (1976).
  * [24] W. McLean, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge University Press, 2000.
  * [25] https://github.com/billmclean/ElasticityQMC.jl
  * [26] R. Scott, A sharp form of the Sobolev trace theorems. J. Funct. Anal., 25, 70–80 (1977).
  * [27] G. Strang and G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ, 1973.
  * [28] M. Vogelius, An analysis of the p-version of the ﬁnite element method for nearly incompressible materials. Uniformly valid, optimal error estimates, Numer. Math., 41, 39–53 (1983).
  * [29] B. Xia and V. H. Hoang, Best N-term GPC approximations for a class of stochastic linear elasticity equations, Mathematical Models and Methods in Applied Sciences, 24, 513–552 (2014).