# Quasi-Monte Carlo sparse grid Galerkin finite element methods for linear
elasticity equations with uncertainties††thanks: This work was supported by
the Australian Research Council grant DP220101811.

J. Dick, Q. T. Le Gia, K. Mustapha and T. Tran School of Mathematics and
Statistics, University of New South Wales, Sydney, Australia

###### Abstract

We explore a linear inhomogeneous elasticity equation with random Lamé
parameters. The latter are parameterized by a countably infinite number of
terms in separated expansions. The main aim of this work is to estimate
expected values (considered as an infinite dimensional integral on the
parametric space corresponding to the random coefficients) of linear
functionals acting on the solution of the elasticity equation. To achieve
this, the expansions of the random parameters are truncated, a high-order
quasi-Monte Carlo (QMC) is combined with a sparse grid approach to approximate
the high dimensional integral, and a Galerkin finite element method (FEM) is
introduced to approximate the solution of the elasticity equation over the
physical domain. The error estimates from (1) truncating the infinite
expansion, (2) the Galerkin FEM, and (3) the QMC sparse grid quadrature rule
are all studied. For this purpose, we show certain required regularity
properties of the continuous solution with respect to both the parametric and
physical variables. To achieve our theoretical regularity and convergence
results, some reasonable assumptions on the expansions of the random
coefficients are imposed. Finally, some numerical results are delivered.

##  1 Introduction

In this work we investigate and analyze the application of the high-order
quasi-Monte Carlo (QMC) sparse grid method combined with the conforming
Galerkin finite element methods (FEMs) to solve a linear elastic model with
uncertainties (see [20] for an interesting overview of how to incorporate
uncertainty into material parameters in linear elasticity problems). More
specifically, we consider the case where the properties of the elastic
inhomogeneous material are varying spatially in an uncertain way by using
random Lamé parameters which are parametrized by a countably infinite number
of parameters. This leads to randomness in both the Young modulus ($E$) and
the Poisson ratio ($\nu$). We intend to measure theoretically the efficiency
of our numerical algorithm through the expectation of the random solution over
the random field.

The equation governing small elastic deformations of a body $\Omega$ in
$\mathbb{R}^{d}$ ($d\in\\{2,3\\}$) with polyhedral boundary can be written as

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

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

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

with ${\bf u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})$ being the displacement
vector field of dimension $d$, and the symmetric strain-rate tensor
${\boldsymbol{\varepsilon}}({\bf u}):=(\nabla{\bf u}+(\nabla{\bf
u})^{T})/2\,.$ Here, ${\bf f}\in[L^{2}(\Omega)]^{d}$ 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 Lamé elasticity parameter $\lambda$
is related to the compressibility of the material, and the shear modulus $\mu$
is related to how the material behaves under normal stress and shear stress,
for the material $\Omega$ containing uncertainties. To parametrize these
uncertainties, we assume that $\mu=\mu({\boldsymbol{x}},{\boldsymbol{y}})$ and
$\lambda=\lambda({\boldsymbol{x}},{\boldsymbol{z}})$ can be 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_{0}({\boldsymbol{x}})+\sum_{j=1}^{\infty}z_{j}\phi_{j}({\boldsymbol{x}}),\quad{\boldsymbol{x}}\in\Omega,$
(1.2)

where $\\{\psi_{j}\\}$ and $\\{\phi_{j}\\}$ are orthogonal basis functions for
the $L^{2}(\Omega)$ space. The parameter vectors
${\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}}$, consist of a countable
number of parameters $y_{j}$ and $z_{k}$, respectively, which are assumed to
be i.i.d. uniformly distributed. Using independent random fields for $\lambda$
and $\mu$ in our model assumes that the compressibility and behaviour under
stress of the material are, in the range of the parameters of the random
$\lambda$ and $\mu$, independent.

The model problem (1.1) is similar to the one in [15, 24], in which a priori
analysis for so-called best $N$-term approximations of standard two-field
mixed formulations was investigated. For the well-posedness of (1.1), we
assume that there are some positive numbers $\mu_{\min},$ $\mu_{\max}$ and
$\lambda_{\max}$ so that

$0<\mu_{\min}\leq\mu({\boldsymbol{x}},{\boldsymbol{y}})\leq\mu_{\max}~{}~{}{\rm
and}~{}~{}0\leq\lambda({\boldsymbol{x}},{\boldsymbol{z}})\leq\lambda_{\max},\quad{\boldsymbol{x}}\in\Omega,\;\;{\boldsymbol{y}},{\boldsymbol{z}}\in
U.$ (A1)

Due to (A1), the values of the Poisson ratio of the elastic material
$\nu=\frac{\lambda}{2(\lambda+\mu)}$ are ranging between $0$ and $1/2$ which
is the case in most materials. Indeed, 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$ is constant. This case was studied in [17,
18] where the authors introduced a three-field PDE model with a parameter-
dependent $E$ which is amenable to discretization by stochastic Galerkin mixed
FEMs. The focus in [17] was on efficient linear algebra, while an a posteriori
error estimation was detailed in [18]. In relation, the authors in [9]
presented a framework for residual-based a posteriori error estimation and
adaptive stochastic Galerkin approximation.

If $\nu$ approaches $1/2$, then we are dealing with an elastic material that
becomes nearly incompressible. In this case, and with constant Lamé
parameters, the convergence rate of the conforming piecewise quadratic or
cubic Galerkin FEMs is short by one from being optimal [22, 23]. However, the
piecewise linear Galerkin FEM runs into trouble with the phenomenon of locking
where the convergence rates may deteriorate as $\lambda$ becomes too large.
This is owing to their inability to represent non-trivial divergence-free
displacement fields. Locking can be avoided by using a nonconforming Galerkin
FEM [23, 2, 10] or by using mixed FEMs. These and several other methods were
studied very extensively in the existing literature for the case of constant
Lamé parameters, we refer the reader to the following books [1, 3, 4].
Investigating the nearly incompressible case with stochastic Lamé parameters
is beyond the scope of this paper; it is a topic of future research.

Outlines of the paper. The next section is devoted to the statement of the
main results of this work. In Section 3, we derive the variational formulation
of the parametric model problem (1.1), and prove the existence and uniqueness
of the weak solution. We also investigate certain regularity properties of the
continuous solution ${\bf u}$ with respect to the random parameters
${\boldsymbol{y}}$ and ${\boldsymbol{z}}$, and the physical parameter
${\boldsymbol{x}}.$ These results are needed to guarantee the convergence of
the errors from both the QMC integration and the Galerkin finite element
discretization. The error from truncating the infinite series expansion in
(1.2) is investigated in Section 4. In Section 5, for every
${\boldsymbol{y}},{\boldsymbol{z}}\in U,$ we approximate the parametric
solution ${\bf u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})$ of (1.1) over the
physical domain $\Omega$ using the conforming Galerkin FEM, and discuss the
stability and error estimates. In Section 6, we investigate the high-order QMC
error from estimating the expected value of a given function over a high
dimensional field. More precisely, we use a high-order QMC rule for the random
coefficients arising from the expansion for $\lambda$ and another QMC rule for
the random coefficients arising from the expansion for $\mu$. In principle it
would be possible to use one high-order QMC rule to simulate both $\lambda$
and $\mu$ simultaneously, but this leads to a complicated design of the high-
order QMC rule. Using separate high-order rules allows us to use existing
constructions from [8]. We study two ways of combining the QMC rules, one is a
tensor product structure (Theorem 6.1) and the other one is a sparse grid
combination (Theorem 6.2). As one would expect, the QMC sparse grid
combination leads to a better rate of convergence. We end the paper with some
numerical simulations in Section 7. In a sample of four different examples, we
illustrate numerically the achieved theoretical finite element QMC convergence
results.

##  2 Main results

We start this section by introducing the following vector function spaces and
the associated norms, which we will be using throughout the remainder of the
paper. Let ${\bf V}:=[H^{1}_{0}(\Omega)]^{d}$ and the associated norm be
$\|{\bf w}\|_{{\bf
V}}:=\Big{(}\sum_{i=1}^{d}\|w_{i}\|_{H^{1}(\Omega)}^{2}\Big{)}^{1/2}$ with the
$w_{i}$’s being the components of the vector function ${\bf w}.$ For
$J=0,1,2,\cdots,$ the norm on the vector Sobolev space ${{\bf
H}^{J}}:=[H^{J}(\Omega)]^{d}$, denoted by $\|\cdot\|_{{{\bf H}^{J}}},$ is
defined in a similar fashion with $\|w_{i}\|_{H^{J}(\Omega)}$ in place of
$\|w_{i}\|_{H^{1}(\Omega)}.$ We dropout ${\bf H}^{J}$ from the norm notation
on the space ${\bf H}^{0}={\bf L}^{2}(\Omega):=[L^{2}(\Omega)]^{d}$ for $d\geq
1.$ Here, $H^{1}_{0}(\Omega)$ and $H^{J}(\Omega)$ are the usual Sobolev
spaces. Finally, ${\bf V}^{*}$ is denoted the dual space of ${\bf V}$ with
respect to the ${\bf L}^{2}(\Omega)$ inner product, with norm denoted by
$\|\cdot\|_{{\bf V}^{\ast}}$.

As mentioned earlier, and more precisely, we are interested in efficient
approximation of the expected value of the function $\mathcal{L}({\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})),$ for a certain linear functional
$\mathcal{L}:{\bf L}^{2}(\Omega)\to\mathbb{R}$, with respect to the random
variables ${\boldsymbol{y}}$ and ${\boldsymbol{z}},$ where ${\bf u}$ is the
solution of (1.1). In other words, we seek to approximate

$\Xi_{\bf u}:=\int_{U}\int_{U}\mathcal{L}\big{(}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\big{)}\,d{\boldsymbol{y}}\,d{\boldsymbol{z}}\,,$
(2.1)

where $d{\boldsymbol{y}}$ and $d{\boldsymbol{z}}$ are the uniform probability
measures on $U$. As a practical example, we may choose $\mathcal{L}$ to be a
local continuous average on some domain $\Omega_{0}\subset\Omega.$

To accomplish the above task, and for the practical implementation, the
occurred infinite sums in (1.2) must be truncated. Then, we approximate ${\bf
u}$ by ${\bf u}_{\bf s}$ which is the solution of (1.1) obtained by truncating
the infinite expansions in (1.2) where ${\bf s}=(s_{1},s_{2})$, that is,
assuming that $y_{j}=z_{k}=0$ for $j>s_{1}$ and $k>s_{2}$. Then, with
$U_{i}=[0,1]^{s_{i}}$ for $i=1,2,$ being of (finite) fixed dimension $s_{i}$,
we estimate the expected value of $\mathcal{L}({\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))$ by approximating

$\Xi_{{\bf s},{\bf u}_{\bf
s}}:=\int_{U_{2}}\int_{U_{1}}\mathcal{L}\Big{(}{{\bf u}_{{\bf
s}}}\Big{(}\cdot,{\boldsymbol{y}}-{\bf\frac{1}{2}},{\boldsymbol{z}}-{\bf\frac{1}{2}}\Big{)}\Big{)}d{\boldsymbol{y}}\,d{\boldsymbol{z}}\,.$
(2.2)

In the above finite dimensional integral, $d{\boldsymbol{y}}$ and
$d{\boldsymbol{z}}$ are the uniform probability measures on $U_{1}$ and
$U_{2}$, respectively. 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.$ We approximate such $(s_{1}+s_{2})$-dimensional integrals using a
high-order QMC quadrature. Preceding this, we intend to solve the truncated
problem over the physical domain $\Omega$ numerically via a continuous
Galerkin FEM. So, for every ${\boldsymbol{y}}\in U_{1}$ and
${\boldsymbol{z}}\in U_{2},$ we approximate the truncated solution ${\bf
u}_{\bf
s}\Big{(}\cdot,{\boldsymbol{y}}-{\bf\frac{1}{2}},{\boldsymbol{z}}-{\bf\frac{1}{2}}\Big{)}$
by the parametric spatial Galerkin finite element solution ${\bf u}_{{\bf
s}_{h}}\Big{(}\cdot,{\boldsymbol{y}}-{\bf\frac{1}{2}},{\boldsymbol{z}}-{\bf\frac{1}{2}}\Big{)}\in{\bf
V}_{h}\subset{\bf V}$ (see Section 5 for the definition of the finite element
space ${\bf V}_{h}$) with the sums in (1.2) truncated to $s_{1}$ and $s_{2}$
terms, respectively. In the third step, we estimate the expectation of the
approximation using first a tensor product of two high-order QMC methods and
second, an efficient high-order QMC sparse grid combination. In summary, we
approximate the expected value in (2.1) by the following truncated QMC
Galerkin finite element rule

$\Xi_{{\bf u}_{{\bf
s}_{h}},Q}:=\frac{1}{N_{1}\,N_{2}}\sum_{j=0}^{N_{1}-1}\sum_{k=0}^{N_{2}-1}\mathcal{L}\Big{(}{\bf
u}_{{\bf
s}_{h}}\Big{(}\cdot,{\boldsymbol{y}}_{j}-{\bf\frac{1}{2}},{\boldsymbol{z}}_{k}-{\bf\frac{1}{2}}\Big{)}\Big{)}\,,$
(2.3)

where 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}$.
Therefore, we have three sources of error: a dimension truncation error
depending on $s_{1}$ and $s_{2}$, a Galerkin discretization error depending on
the maximum finite element mesh diameter $h$ of the physical domain $\Omega$,
and a QMC quadrature error which depends on $N_{1}$ and $N_{2}$. We split the
error as:

$|\Xi_{\bf u}-\Xi_{{\bf u}_{{\bf s}_{h}},Q}|\leq|\Xi_{\bf u}-\Xi_{{\bf s},{\bf
u}_{\bf s}}|+|\Xi_{{\bf s},{\bf u}_{\bf s}}-\Xi_{{\bf s},{\bf u}_{{\bf
s}_{h}}}|+|\Xi_{{\bf s},{\bf u}_{{\bf s}_{h}}}-\Xi_{{\bf u}_{{\bf
s}_{h}},Q}|.$ (2.4)

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

$\Xi_{\bf u}-\Xi_{{\bf s},{\bf u}_{\bf
s}}=\int_{U}\int_{U}\mathcal{L}\Big{(}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})-{\bf u}_{{\bf
s}}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})\Big{)}\,d{\boldsymbol{y}}\,d{\boldsymbol{z}},$

where ${\boldsymbol{y}}=(y_{j})_{j\geq 1},{\boldsymbol{z}}=(z_{k})_{k\geq
1}\in U$, and the truncated vectors ${\boldsymbol{y}}^{s_{1}}$ and
${\boldsymbol{z}}^{s_{2}}$ are $(y_{1},y_{2},\cdots,y_{s_{1}},0,0,\cdots)$ and
$(z_{1},z_{2},\cdots,z_{s_{2}},0,0,\cdots)$, respectively. To estimate this
term, we refer to the dimension truncation error which is analyzed in Theorem
4.1. To reduce the errors from such a truncation, which is necessary from a
practical point of view, we assume that the $L^{2}(\Omega)$ orthogonal basis
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\,.$ (A2)

For the convergence from the series truncation, we assume that
$\mu_{0},\,\lambda_{0}\in L^{\infty}(\Omega)$, and

$\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\leq 1\,.$ (A3)

When $p=1$ and$/$or $q=1,$ it is essential to have

$\sum_{j=s_{1}+1}^{\infty}\|\psi_{j}\|_{L^{\infty}(\Omega)}\leq
Cs_{1}^{1-1/\varrho_{1}}~{}~{}{\rm
and/or}\sum_{j=s_{2}+1}^{\infty}\|\phi_{j}\|_{L^{\infty}(\Omega)}\leq
Cs_{2}^{1-1/\varrho_{2}},$ (A4)

for some $0<\varrho_{1},\varrho_{2}<1$. The second term on the right-hand side
of (2.4) is the finite dimensional integral of the linear functional
$\mathcal{L}$ acting on the difference between the truncated solution ${\bf
u}_{\bf s}$ and its approximation ${\bf u}_{{\bf s}_{h}}$. This can be deduced
from Theorem 5.1 by replacing the vectors ${\boldsymbol{y}}$ and
${\boldsymbol{z}}$ with ${\boldsymbol{y}}^{s_{1}}$ and
${\boldsymbol{z}}^{s_{2}}$, respectively, and using the fact that ${\bf
u}_{\bf s}$ satisfies the regularity properties in Theorem 3.2. The Lamé
parameters $\mu(\cdot,{\boldsymbol{y}})$ and $\lambda(\cdot,{\boldsymbol{z}})$
are required to be in the Sobolev space $W^{\theta,\infty}(\Omega)$ for every
${\boldsymbol{y}},\,{\boldsymbol{z}}\in U$ and for some integer
$1\leq\theta\leq r$ ($r$ is the degree of the finite element solution, so
$\theta=1$ in the case of piecewise linear Galerkin FEM). To have this, we
assume that

$\mu_{0},\,\lambda_{0}\in
W^{\theta,\infty}(\Omega),~{}~{}\sum_{j=1}^{\infty}\|\psi_{j}\|_{W^{\theta,\infty}(\Omega)}~{}~{}{\rm
and}~{}~{}\sum_{j=1}^{\infty}\|\phi_{j}\|_{W^{\theta,\infty}(\Omega)}~{}{\rm
are~{}finite}\,.$ (A5)

It is clear from Theorem 5.1 that for every
$({\boldsymbol{y}},{\boldsymbol{z}})\in U_{1}\times U_{2}$, the sequence
$|\mathcal{L}({\bf u}_{\bf s}-{\bf u}_{{\bf s}_{h}})|$ converges faster than
$\|{\bf u}_{\bf s}-{\bf u}_{{\bf s}_{h}}\|_{\bf V}$ provided that the linear
functional $\mathcal{L}$ is bounded in the $L^{2}(\Omega)$ sense (that is,
$|\mathcal{L}({\bf w})|\leq\|\mathcal{L}\|\|{\bf w}\|$ for any ${\bf w}\in{\bf
L}^{2}(\Omega)$), which is not always guaranteed.

The third term in (2.4) is the QMC quadrature error which can be estimated by
applying Theorem 6.1 and Theorem 6.2 with
${F}({\boldsymbol{y}},{\boldsymbol{z}}):=\mathcal{L}\Big{(}{\bf u}_{{\bf
s}_{h}}\Big{(}\cdot,{\boldsymbol{y}}-{\bf\frac{1}{2}},{\boldsymbol{z}}-{\bf\frac{1}{2}}\Big{)}\Big{)}.$
Noting that (the mixed derivative)
$|\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}{F}({\boldsymbol{y}},{\boldsymbol{z}})|\leq\|\mathcal{L}\|_{{\bf
V}*}\,\|\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}{\bf
u}_{{\bf
s}_{h}}\big{(}\cdot,{\boldsymbol{y}}-{\bf\frac{1}{2}},{\boldsymbol{z}}-{\bf\frac{1}{2}}\big{)}\|_{{\bf
V}}$, and then, Theorem 3.3 can be applied to verify the regularity conditions
in (6.3) and (6.7) which are necessary for the QMC error results in Theorem
6.1 and the QMC sparse grid error results in Theorem 6.2, respectively. Here,
Assumptions (A2) and (A3) (for $p=q=1$ only) are needed.

We summarize the combined error estimate in the next theorem. In addition to
Assumptions (A1)–(A5), we assume that the physical domain $\Omega$ is
${\mathcal{C}}^{\theta,1}$ or the boundary of $\Omega$ is of class
${\mathcal{C}}^{\theta+1}$ for some integer $\theta\geq 1$ (for $\theta=1$,
$\Omega$ can be convex instead) and the body force vector function ${\bf f}$
belongs to ${\bf H}^{\theta-1}(\Omega)$ (recall that ${\bf H}^{0}(\Omega)={\bf
L}^{2}(\Omega)$). These additional assumptions are needed to guarantee that
the strong solution ${\bf u}$ of (1.1) is in the space ${\bf
H}^{\theta+1}(\Omega)$, see Theorem 3.2. This property is essential for the
optimal convergence of the finite element solution of (1.1) over $\Omega.$
Throughout the paper, $C$ is a generic constant that is independent of $h$,
the number of QMC points $N_{i}$, and the dimension $s_{i}$ for $i=1,\,2,$ but
may depend on the physical domain $\Omega$ and other parameters that will be
mentioned accordingly.

###### Theorem 2.1

Let ${\bf u}$ be the solution of problem (1.1) and ${\bf u}_{{\bf s}_{h}}$ be
the Galerkin finite element solution of degree $\leq r$ (with $r\geq 1$)
defined as in (5.2) with $y_{j}=z_{k}=0$ for $j>s_{1}\geq 1$ and $k>s_{2}\geq
1.$ For $i=1,\,2$, let $N_{i}=b^{m_{i}}$ with $m_{i}$ being positive integers
and $b$ being prime. Then one can construct two interlaced polynomial lattice
rules of order $\alpha\,:=\,\lfloor 1/p\rfloor+1$ with $N_{1}$ points, and of
order $\beta:=\,\lfloor 1/q\rfloor+1$ with $N_{2}$ points where
$|m_{1}q-m_{2}p|<1$, and $p$ and $q$ are those in (A3), such that the
following QMC Galerkin finite element error bound holds: for $1\leq\theta\leq
r,$

$|\Xi_{\bf u}-\Xi_{{\bf u}_{{\bf s}_{h}},Q}|\leq C\,h^{\theta+1}\|{\bf
f}\|_{{\bf H}^{\theta-1}}\|\mathcal{L}\|\\\
+C\,\Big{(}s_{1}^{1-\max(1/p,1/\varrho_{1})}+s_{2}^{1-\max(1/q,1/\varrho_{2})}+N^{-\frac{1}{p+q}}\Big{)}\|{\bf
f}\|_{{\bf V}^{*}}\|\mathcal{L}\|_{{\bf V}^{*}}\,,$

where $N=N_{1}N_{2}$ is the total number of QMC quadrature points. The
constant $C$ depends on $b,p,q,\lambda,$ and $\mu$, but is independent of
$s_{i}$ and $m_{i}$ for $i\in\\{1,2\\},$ and $h$.

Further, there exists a combined QMC sparse grid approximation given by (6.6),
such that with the assumptions in Theorem 6.2, the above error bound remains
valid with $M^{-\min(1/p,1/q)}$ in place of $N^{-\frac{1}{p+q}}$, where $M$ is
the total number of quadrature points in the QMC sparse grid approximation.
The constant $C$ in the new bound depends on $b,p,q,\lambda,$ and $\mu$, but
is independent of $s_{1},s_{2}$, $M$, and $h$.

##  3 Weak formulation and regularity

This section is devoted to deriving the weak formulation of the parametric
elasticity equation (1.1) for each value of the parameter
${\boldsymbol{y}},\,{\boldsymbol{z}}\in U.$ Then we show 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}}$. Preceding this, we establish the existence and uniqueness
of the weak solution.

For the weak formulation of (1.1), for every
${\boldsymbol{y}},{\boldsymbol{z}}\in U$, we multiply both sides of (1.1) by a
test function ${\bf v}\in{\bf V}$, and use Green’s formula and the given
homogeneous Dirichlet boundary conditions after integrating over the physical
domain $\Omega.$ Then, the usage of the identities

$\boldsymbol{\sigma}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u}):\nabla{\bf
v}=\boldsymbol{\sigma}({\boldsymbol{y}},{\boldsymbol{z}};{\bf
u}):{\boldsymbol{\varepsilon}}({\bf v})=\lambda\nabla\cdot{\bf
u}\nabla\cdot{\bf v}+2\mu{\boldsymbol{\varepsilon}}({\bf
u}):{\boldsymbol{\varepsilon}}({\bf v})$

(the colon operator is the inner product between tensors) results in the
following parameter-dependent weak formulation: for every
${\boldsymbol{y}},\,{\boldsymbol{z}}\in U$, find ${\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\in{\bf V}$ satisfying

$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u},{\bf v})=\ell({\bf
v}),\quad\text{for all}\quad{\bf v}\in{\bf V},$ (3.1)

where the bilinear form
$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\cdot,\cdot)$ and the linear
functional $\ell(\cdot)$ are defined by

$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u},{\bf
v}):=\int_{\Omega}[2\mu\,{\boldsymbol{\varepsilon}}({\bf
u}):{\boldsymbol{\varepsilon}}({\bf v})+\lambda\nabla\cdot{\bf
u}\nabla\cdot{\bf v}]\,d{\boldsymbol{x}}~{}~{}{\rm and}~{}~{}\ell({\bf
v}):=\langle{\bf f},{\bf v}\rangle:=\int_{\Omega}{\bf f}\cdot{\bf
v}\,d{\boldsymbol{x}}.$ (3.2)

The next theorem shows the existence and uniqueness of the solution of (3.1).

###### Theorem 3.1

Assume that (A1) is satisfied. Then, for every $f\in{\bf V}^{*}$ and
${\boldsymbol{y}},\,{\boldsymbol{z}}\in U$, the parametric weak formulation
problem (3.1) has a unique solution.

Proof. Using assumption (A1) and applying the Cauchy-Schwarz inequality leads
to

$|\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf v},{\bf w})|\leq
2\mu_{\max}\|{\boldsymbol{\varepsilon}}({\bf
v})\|\,\|{\boldsymbol{\varepsilon}}({\bf w})\|+\lambda_{\max}\|\nabla\cdot{\bf
v}\|\,\|\nabla\cdot{\bf w}\|\,.$

Hence, using the inequalities $\|\nabla\cdot{\bf
w}\|\leq\sqrt{d}\,\|\nabla{\bf w}\|$ and $\|{\boldsymbol{\varepsilon}}({\bf
w})\|\leq\|\nabla{\bf w}\|$, we have

$|\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf v},{\bf
w})|\leq(d\lambda_{\max}+2\mu_{\max})\|\nabla{\bf v}\|\,\|\nabla{\bf
w}\|\leq(d\lambda_{\max}+2\mu_{\max})\|{\bf v}\|_{{\bf V}}\,\|{\bf w}\|_{{\bf
V}},$ (3.3)

for any ${\bf v},{\bf w}\in{\bf V}$. So, the bilinear form
$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\cdot,\cdot)$ is bounded on
${\bf V}\times{\bf V}$. For the coercivity property of
$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\cdot,\cdot)$, we use again
assumption (A1) in addition to Korn’s inequality to obtain

$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\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.4)

Since $\ell(\cdot)$ is a bounded linear functional on ${\bf V}$, an
application of the Lax-Milgram theorem completes the proof. $\quad\Box$

For the finite element error analysis, we discuss next some required
regularity properties of the parametric solution of (3.1). For the nearly
incompressible case (which is beyond the scope of this work), one has to be
more specific about the constant $\widetilde{C}$ in the following theorem.

###### Theorem 3.2

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

$\|{\bf u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf
V}}\leq\frac{C}{\mu_{\min}}\|{\bf f}\|_{{\bf V}^{*}}\,.$ (3.5)

If $\Omega$ is ${\mathcal{C}}^{\theta,1}$ (or the boundary of $\Omega$ is of
class ${\mathcal{C}}^{\theta+1}$) for some integer $\theta\geq 1$ (for
$\theta=1$, we may assume that $\Omega$ is convex instead), then ${\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\in{\bf V}\cap{\bf
H}^{\theta+1}(\Omega)$ (that is, it is a strong solution of problem (1.1))
provided that (A5) is satisfied and ${\bf f}\in{\bf H}^{\theta-1}(\Omega)$.
Furthermore,

$\|{\bf u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf
H}^{\theta+1}}\leq\widetilde{C}\,\|{\bf f}\|_{{\bf H}^{\theta-1}},\quad{\rm
for~{}every}~{}~{}{\boldsymbol{y}},\,{\boldsymbol{z}}\in U\,.$ (3.6)

The constant $\widetilde{C}$ depends on $\Omega$, $\mu$, $\lambda$, including
$\|\mu(\cdot,{\boldsymbol{y}})\|_{W^{\theta,\infty}(\Omega)}$ and
$\|\lambda(\cdot,{\boldsymbol{z}})\|_{W^{\theta,\infty}(\Omega)}$.

Proof. From the coercivity property in (3.4) 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{y}},{\boldsymbol{z}};{\bf u},{\bf
u})=\ell({\bf u})\leq\|{\bf f}\|_{{\bf V}^{*}}\|{\bf u}\|_{{\bf V}},$

for every ${\boldsymbol{y}},\,{\boldsymbol{z}}\in U$. Thus, the proof of the
regularity estimate in (3.5) is completed.

For every ${\boldsymbol{y}},\,{\boldsymbol{z}}\in U$, the operator
$\nabla\cdot\boldsymbol{\sigma}$ in the elasticity equation (1.1) is strongly
elliptic because the bilinear operator
$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\cdot,\cdot)$ is coercive on
${\bf V}$. Thus, due to the imposed assumptions on $\Omega$, the solution
${\bf u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})$ of (3.1) is in the space
${\bf H}^{\theta+1}$ and satisfies the regularity property in (3.6). See [5,
Theorem 6.3-6], [14, Theorems 2.4.2.5, 2.5.1.1, and 3.2.1.3] and [21, Chapter
4] for more details about strongly elliptic system and elliptic regularity.
$\quad\Box$

In the QMC (see Theorem 6.1) and QMC sparse grids (see Theorem 6.2) error
estimations, we need to bound the mixed first partial derivatives of the
parametric displacement ${\bf u}$ with respect to the random variables $y_{j}$
and $z_{k}$. This will be the topic of the next theorem. For convenience, we
introduce ${\mathcal{S}}$ to be 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.$ That is,
sequences of nonnegative integers for which only finitely many entries are
nonzero. For $\boldsymbol{\alpha}=(\alpha_{j})_{j\geq 1}$ and
$\boldsymbol{\beta}=(\beta_{j})_{j\geq 1}$ belonging to $\mathcal{S},$ the
mixed partial derivative
$\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}$
is defined by

$\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}:=\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}.$

It reduces to $\partial_{{\boldsymbol{y}}}^{\boldsymbol{\beta}}$ and
$\partial_{{\boldsymbol{z}}}^{\boldsymbol{\alpha}}$ when
$|\boldsymbol{\beta}|=0$ and $|\boldsymbol{\alpha}|=0,$ respectively.

###### Theorem 3.3

Assume that (A2) and (A3) (for $p=q=1$) are satisfied. Then, for every ${\bf
f}\in{\bf V}^{*}$, ${\boldsymbol{y}},{\boldsymbol{z}}\in U$, and
$\boldsymbol{\alpha},\boldsymbol{\beta}\in\mathcal{S}$, the solution ${\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})$ of the parametric weak problem
(3.1) satisfies

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

where

$\widetilde{\bf b}^{\boldsymbol{\alpha}}=\prod_{i\geq
1}(\widetilde{b}_{i})^{\alpha_{i}}~{}~{}{\rm and}~{}~{}\widehat{\bf
b}^{\boldsymbol{\beta}}=\prod_{i\geq
1}(\widehat{b}_{i})^{\beta_{i}},~{}~{}{\rm
with}~{}~{}\widetilde{b}_{j}=\frac{\|\psi_{j}\|_{L^{\infty}(\Omega)}}{\mu_{\min}}~{}~{}{\rm
and}~{}~{}\widehat{b}_{j}=\frac{d}{2}\,\frac{\|\phi_{j}\|_{L^{\infty}(\Omega)}}{\mu_{\min}}\,.$

Consequently,

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

where the constant $C$ depends on $\Omega$ only.

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

$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}{\bf
u},{\bf
v})=-2\sum_{\boldsymbol{\alpha}}\alpha_{j}\int_{\Omega}\psi_{j}{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}-{\bf
e}_{j}}{\bf u}):{\boldsymbol{\varepsilon}}({\bf v})\,d{\boldsymbol{x}}\\\
-\sum_{\boldsymbol{\beta}}\beta_{k}\int_{\Omega}\phi_{k}\nabla\cdot(\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta}-{\bf
e}_{k},\boldsymbol{\alpha}}{\bf u})\nabla\cdot{\bf
v}\,d{\boldsymbol{x}},\quad\text{for all}\quad{\bf v}\in{\bf V},$ (3.9)

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.

Choosing ${\bf
v}=\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}{\bf
u}$ in (3.9), then using the inequality $\|\nabla\cdot{\bf
v}\|\leq\sqrt{d}\,\|{\boldsymbol{\varepsilon}}({\bf v})\|$ and the coercivity
property of $\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\cdot,\cdot)$ in
(3.4), after some simplifications, we obtain

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

and consequently,

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

When $\boldsymbol{\beta}={\bf 0},$ the above inequality reduces to

$\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}}}^{\boldsymbol{\alpha}}{\bf
u})\|\leq\sum_{\boldsymbol{\alpha}}\alpha_{j}\widetilde{b}_{j}\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}}}^{\boldsymbol{\alpha}-{\bf
e}_{j}}{\bf u})\|\,.$ (3.11)

Recursively, we obtain

$\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{y}}}^{\boldsymbol{\alpha}}{\bf
u})\|\leq[\,|\boldsymbol{\alpha}|\,(|\boldsymbol{\alpha}|-1)\,\cdots
1]\,\prod_{i\geq
1}(\widetilde{b}_{i})^{\alpha_{i}}\|{\boldsymbol{\varepsilon}}({\bf
u})\|=|\boldsymbol{\alpha}|!\prod_{i\geq
1}(\widetilde{b}_{i})^{\alpha_{i}}\|{\boldsymbol{\varepsilon}}({\bf u})\|,$

and hence, (3.7) holds true in this case. In a similar fashion, we can show
(3.7) when $\boldsymbol{\alpha}={\bf 0}.$

When $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are both not identically
zero vectors, the above approach can be extended, however it is not easy to
follow. Owing to this, following [6], we use instead the induction hypothesis
on $n:=|\boldsymbol{\alpha}+\boldsymbol{\beta}|$. From the above contribution,
it is clear that (3.7) holds true when
$|\boldsymbol{\alpha}+\boldsymbol{\beta}|=1$. Now, assume that (3.7) is true
for $|\boldsymbol{\alpha}+\boldsymbol{\beta}|=n$, and the task is to claim it
for $|\boldsymbol{\alpha}+\boldsymbol{\beta}|=n+1.$

From (3.10) and the induction hypothesis, we have

$\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}{\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}}+\sum_{\boldsymbol{\beta}}\beta_{k}\,\widehat{b}_{k}\widetilde{\bf
b}^{\boldsymbol{\alpha}}\widehat{\bf b}^{\boldsymbol{\beta}-{\bf
e}_{k}}\Big{)}\|{\boldsymbol{\varepsilon}}({\bf u})\|\\\ =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})\|\,.$

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

Finally, since
$\|{\boldsymbol{\varepsilon}}(\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))\|\geq
C\,\|\partial_{{\boldsymbol{z}},{\boldsymbol{y}}}^{\boldsymbol{\beta},\boldsymbol{\alpha}}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf V}}$ (by Korn’s
inequality) and since $\|{\boldsymbol{\varepsilon}}({\bf u})\|\leq\|\nabla{\bf
u}\|\leq\frac{C}{\mu_{\min}}\|{\bf f}\|_{{\bf V}^{*}}$ (by (3.5)), we derive
(3.8) from (3.7). $\quad\Box$

##  4 A truncated problem

This section is dedicated to investigating the error from truncating the first
and second sums in (1.2) at $s_{1}$ and $s_{2}$ terms, respectively, from some
$s_{1},s_{2}\in\mathbb{N}.$ In other words, we set $y_{j}=0$ and $z_{k}=0$ for
$j>s_{1}$ and $k>s_{2}$, respectively. We start by defining the truncated weak
formulation problem: for every
${\boldsymbol{y}}^{s_{1}},\,{\boldsymbol{z}}^{s_{2}}\in U,$ find ${\bf u}_{\bf
s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})\in{\bf V}$, with
${\bf s}=(s_{1},s_{2})$, such that

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

Thanks to Theorem 3.1, the truncated problem (4.1) has a unique solution.
Estimating the truncation error, which is needed for measuring the QMC finite
element error in (2.4), is the topic of the next theorem. For brevity, we let

$\mu_{c}({\boldsymbol{x}},{\boldsymbol{y}})=\sum_{j=s_{1}+1}^{\infty}y_{j}\psi_{j}({\boldsymbol{x}})~{}~{}{\rm
and}~{}~{}\lambda_{c}({\boldsymbol{x}},{\boldsymbol{z}})=\sum_{j=s_{2}+1}^{\infty}z_{j}\phi_{j}({\boldsymbol{x}}),\quad{\boldsymbol{x}}\in\Omega,\
{\boldsymbol{y}},{\boldsymbol{z}}\in U.$

###### Theorem 4.1

Under Assumption (A1), for every ${\bf f}\in{\bf V}^{*}$,
${\boldsymbol{y}},{\boldsymbol{z}}\in U$, and ${\bf
s}=(s_{1},s_{2})\in\mathbb{N}^{2}$, the solution ${\bf u}_{\bf s}$ of the
truncated parametric weak problem (4.1) satisfies

$\|{\bf u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})-{\bf u}_{\bf
s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})\|_{{\bf V}}\leq
C\,\widehat{C}\|{\bf f}\|_{{\bf V}^{*}},~{}~{}{\rm
with}~{}~{}\widehat{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 (A2)–(A4) are satisfied, and if $\mathcal{L}:{\bf
V}\to\mathbb{R}$ is a bounded linear functional, (that is, $|\mathcal{L}({\bf
w})|\leq\|\mathcal{L}\|_{{\bf V}^{*}}\|{\bf w}\|_{\bf V}$ for all ${\bf
w}\in{\bf V}$), then for every ${\boldsymbol{y}},\,{\boldsymbol{z}}\in U,$ we
have

$\displaystyle|\mathcal{L}({\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))-\mathcal{L}$ $\displaystyle({\bf
u}_{\bf s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}}))|\leq
C\,\Big{(}s_{1}^{1-\max(1/p,1/\varrho_{1})}+s_{2}^{1-\max(1/q,1/\varrho_{2})}\Big{)}\,\|{\bf
f}\|_{{\bf V}^{*}}\|\mathcal{L}\|_{{\bf V}^{*}},$ (4.2)

for some $0<p\,,q\leq 1$ (see (A3)) and for some
$0<\varrho_{1}\,,\varrho_{2}<1$ (see (A4)) . Here, the (generic) constant $C$
depends on $\Omega$, $\mu_{\max}$, $\mu_{\min}$, and $\lambda_{\max}$.

Proof. From the variational formulations in (3.1) and (4.1), we notice that

$\mathcal{B}({\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}};{\bf u}_{\bf
s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})-{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}),{\bf
v})=\mathcal{B}({\boldsymbol{y}}-{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}-{\boldsymbol{z}}^{s_{2}};{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}),{\bf v}).$

Following the steps in (3.3) and using the achieved estimate in (3.5), we have

$\displaystyle|\mathcal{B}({\boldsymbol{y}}-{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}-$
$\displaystyle{\boldsymbol{z}}^{s_{2}};{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}),{\bf v})|$ $\displaystyle\leq
C\max_{{\boldsymbol{x}}\in\Omega,\ {\boldsymbol{y}},{\boldsymbol{z}}\in
U}(|\lambda_{c}({\boldsymbol{x}},{\boldsymbol{y}})|+|\mu_{c}({\boldsymbol{x}},{\boldsymbol{z}})|)\,\|{\bf
u}\|_{{\bf V}}\,\|{\bf v}\|_{{\bf V}}\leq C\,\widehat{C}\|{\bf f}\|_{{\bf
V}^{*}}\|{\bf v}\|_{{\bf V}}\,.$

On the other hand, by the coercivity property in (3.4), we have

$\mathcal{B}({\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}};{\bf u}_{\bf
s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})-{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}),{\bf u}_{\bf
s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})-{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))\\\ \geq C\mu_{\min}\|{\bf u}_{\bf
s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})-{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf V}}^{2}\,.$

Combining the above equations, then the first desired result follows after
simplifying by similar terms. To show (4.2), we simply use the imposed
assumption on $\mathcal{L}$ and the first achieved estimate to obtain

$\displaystyle\big{|}\mathcal{L}\big{(}{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})-{\bf u}_{\bf
s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})\big{)}\big{|}$
$\displaystyle\leq\|\mathcal{L}\|_{{\bf V}^{*}}\|{\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})-{\bf u}_{\bf
s}(\cdot,{\boldsymbol{y}}^{s_{1}},{\boldsymbol{z}}^{s_{2}})\|_{{\bf V}}\leq
C\,\widehat{C}\,\|{\bf f}\|_{{\bf V}^{*}}\|\mathcal{L}\|_{{\bf V}^{*}}.$

Hence, (4.2) is a direct consequence of the above estimate, the Stechkin
inequality

$\sum_{j\geq s+1}b_{j}\leq
C_{\varsigma}\,s^{1-\frac{1}{\varsigma}}\Big{(}\sum_{j\geq
1}b_{j}^{\varsigma}\Big{)}^{\frac{1}{\varsigma}},\quad{\rm
for}~{}~{}0<\varsigma<1,$

where $\\{b_{j}\\}_{j\geq 1}$ is a nonincreasing sequence of positive numbers,
and assumptions (A2)–(A4). $\quad\Box$

##  5 Finite element approximation

This section is devoted to introducing the Galerkin FEM of degree at most $r$
($r\geq 1$) for the approximation of the solution to the model problem (3.1)
over $\Omega$, and consequently, to problem (1.1). Stability and error
estimates are investigated. The achieved results in this section are needed
for measuring the QMC finite element error in (2.4).

We introduce a family of regular triangulation (made of simplexes)
$\mathcal{T}_{h}$ of the domain $\overline{\Omega}$ and set
$h=\max_{K\in\mathcal{T}_{h}}(h_{K})$, where $h_{K}$ denotes the diameter of
the element $K$. Let $V_{h}\subset H^{1}_{0}(\Omega)$ denote the usual
conforming finite element space of continuous, piecewise polynomial functions
of degree at most $r$ on $\mathcal{T}_{h}$ that vanish on $\partial\Omega$.
Let ${\bf V}_{h}=[V_{h}]^{d}$ be the finite element vector space. Then there
exists a constant $C$ (depending on $\Omega$) such that,

$\inf_{{\bf v}_{h}\in{\bf V}_{h}}\|{\bf v}-{\bf v}_{h}\|_{{\bf V}}\leq
Ch^{\theta}\|{\bf v}\|_{{\bf H}^{\theta+1}},\quad{\rm
for}~{}~{}1\leq\theta\leq r.$ (5.1)

Motivated by the weak formulation in (3.1), we define the parametric finite
element approximate solution as: find ${\bf
u}_{h}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\in{\bf V}_{h}$ such that

$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u}_{h},{\bf
v}_{h})=\ell({\bf v}_{h}),\quad\text{for all}\quad{\bf v}_{h}\in{\bf
V}_{h},\quad{\rm for~{}every}~{}~{}{\boldsymbol{y}},\,{\boldsymbol{z}}\in U.$
(5.2)

Assuming that (A1) is satisfied, then, for every ${\bf f}\in{\bf V}^{*}$ and
every ${\boldsymbol{y}},\,{\boldsymbol{z}}\in U$, the finite element scheme
defined in (5.2) has a unique parametric solution ${\bf
u}_{h}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\in{\bf V}_{h}$. This can be
shown by mimicking the proof of Theorem 3.1 because ${\bf V}_{h}\subset{\bf
V}.$ Furthermore, the finite element solution is also stable; the bound in
(3.5) remains valid with ${\bf u}_{h}$ in place of ${\bf u}$, that is,

$\|{\bf u}_{h}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf
V}}\leq\frac{C}{\mu_{\min}}\|{\bf f}\|_{{\bf V}^{*}}\,.$ (5.3)

In the next theorem, we discuss the ${\bf V}$-norm error estimate from the
finite element discretization. Then, and as in Theorem 4.1, for measuring the
QMC finite element error in (2.4), we derive an estimate that involves a
linear functional acting on the difference ${\bf u}-{\bf u}_{h}.$

###### Theorem 5.1

For every ${\boldsymbol{y}},\,{\boldsymbol{z}}\in U,$ let ${\bf u}$ and ${\bf
u}_{h}$ be the solutions of problems (1.1) and (5.2), respectively. Assuming
that ${\bf u}$ satisfies the regularity properties in (3.6) for some integer
$1\leq\theta\leq r.$ Under Assumption (A1) and (A5), we have

$\|{\bf u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})-{\bf
u}_{h}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf V}}\leq
C\,h^{\theta}\|{\bf f}\|_{{\bf H}^{\theta-1}},$ (5.4)

when ${\bf f}\in{\bf H}^{\theta-1}(\Omega)$. Moreover, if $\mathcal{L}:{\bf
L}^{2}(\Omega)\to\mathbb{R}$ is a bounded linear functional (that is,
$|\mathcal{L}({\bf w})|\leq\|\mathcal{L}\|\,\|{\bf w}\|$), then

$|\mathcal{L}({\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))-\mathcal{L}({\bf
u}_{h}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))|\leq C\,h^{\theta+1}\|{\bf
f}\|_{{\bf H}^{\theta-1}}\|\mathcal{L}\|,\quad{\rm
for~{}every}~{}~{}{\boldsymbol{y}},\,{\boldsymbol{z}}\in U.$ (5.5)

The constant $C$ depends on $\Omega$, $\mu_{\max}$, $\mu_{\min}$, and
$\lambda_{\max}$, but not on $h$.

Proof. The proof of (5.4) follows a standard argument for finite element
approximations and is included here for completeness. From the weak
formulation in (3.1) and the finite element scheme in (5.2), we have the
following orthogonality property

$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u}-{\bf u}_{h},{\bf
v}_{h})=0,\quad\text{for all}\quad{\bf v}_{h}\in{\bf V}_{h}\,.$ (5.6)

By using the above equation, the coercivity and boundedness of
$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\cdot,\cdot)$, we obtain,

$\displaystyle\|{\bf u}-{\bf u}_{h}\|_{{\bf V}}^{2}$ $\displaystyle\leq
C\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u}-{\bf u}_{h},{\bf
u}-{\bf u}_{h})=C\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u}-{\bf
u}_{h},{\bf u}-{\bf w}_{h})\leq C\|{\bf u}-{\bf u}_{h}\|_{{\bf V}}\,\|{\bf
u}-{\bf w}_{h}\|_{{\bf V}},$

for all ${\bf w}_{h}\in{\bf V}_{h}$. This implies $\|{\bf u}-{\bf
u}_{h}\|_{{\bf V}}\leq C\|{\bf u}-{\bf w}_{h}\|_{{\bf V}}$ for all ${\bf
w}_{h}\in{\bf V}_{h}$. Thus, (5.4) follows from (5.1) and the regularity
estimate in (3.6).

To show (5.5), we use the so-called Nitsche trick. We first replace $\ell$ in
(3.1) by $\mathcal{L}$, and consider a new parametric variational problem:
find ${\bf u}_{\mathcal{L}}\in{\bf V}$ such that

$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u}_{\mathcal{L}},{\bf
v})=\mathcal{L}({\bf v}),\quad\text{for all}\quad{\bf v}\in{\bf V}.$ (5.7)

By Theorem 3.1, this problem has a unique solution for every
${\boldsymbol{y}},{\boldsymbol{z}}\in U$. Hence, using Theorem 3.2 (with ${\bf
u}_{\mathcal{L}}$ in place of ${\bf u}$) and the given assumption on
$\mathcal{L}$, we conclude that $\|{\bf u}_{\mathcal{L}}\|_{{\bf H}^{2}}\leq
C\|\mathcal{L}\|.$ Therefore, by repeating the above argument, we deduce

$\|{\bf u}_{\mathcal{L}}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})-{\bf
u}_{\mathcal{L},h}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})\|_{{\bf V}}\leq
Ch\|{\bf u}_{\mathcal{L}}\|_{{\bf H}^{2}}\leq Ch\|\mathcal{L}\|,$ (5.8)

where ${\bf u}_{\mathcal{L},h}\in{\bf V}_{h}$ is the finite element
approximation of ${\bf u}_{\mathcal{L}}$. By using successively the linearity
of $\mathcal{L}$, equation (5.7), the symmetry of
$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\cdot,\cdot)$, the Galerkin
orthogonality (5.6), and the boundedness of
$\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};\cdot,\cdot)$, we obtain

$\displaystyle|\mathcal{L}({\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))-\mathcal{L}({\bf
u}_{h}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))|$
$\displaystyle=|\mathcal{L}({\bf
u}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}})-{\bf
u}_{h}(\cdot,{\boldsymbol{y}},{\boldsymbol{z}}))|$
$\displaystyle=|\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u}-{\bf
u}_{h},{\bf
u}_{\mathcal{L}})|=|\mathcal{B}({\boldsymbol{y}},{\boldsymbol{z}};{\bf u}-{\bf
u}_{h},{\bf u}_{\mathcal{L}}-{\bf u}_{\mathcal{L},h})|$ $\displaystyle\leq
C\|{\bf u}-{\bf u}_{h}\|_{{\bf V}}\,\|{\bf u}_{\mathcal{L}}-{\bf
u}_{\mathcal{L},h}\|_{{\bf V}}.$

The required estimate (5.5) now follows from (5.4) and (5.8) . $\quad\Box$

##  6 QMC method and sparse grids

Our aim is to measure the QMC finite element error which occurs in the third
term on the right hand side of (2.4). To serve this purpose, the current
section is dedicated to investigate both the high-order QMC and the high-order
QMC sparse grid errors from estimating the finite dimensional integral

${\mathcal{I}}_{\bf
s}{F}:=\int_{U_{2}}\int_{U_{1}}{F}({\boldsymbol{y}},{\boldsymbol{z}})\,d{\boldsymbol{y}}\,d{\boldsymbol{z}}\,.$
(6.1)

Recall that $U_{i}=[0,1]^{s_{i}}$ are of fixed dimensions $s_{i}$ for $i=1,2,$
and ${\bf s}=(s_{1},s_{2})$. We approximate ${\mathcal{I}}_{\bf s}{F}$ via an
equal-weight rule of the form:

$Q_{\bf s,\bf
N}[{F}]:=\frac{1}{N_{1}\,N_{2}}\sum_{k=0}^{N_{2}-1}\sum_{j=0}^{N_{1}-1}{F}({\boldsymbol{y}}_{j},{\boldsymbol{z}}_{k}),\quad{\rm
with}~{}~{}{\bf N}=(N_{1},N_{2}),$ (6.2)

where $N_{i}=b^{m_{i}}$ for a given prime $b$ and a given positive integer
$m_{i}$, with $i=1,2.$ The QMC points
$\\{{\boldsymbol{y}}_{0},\ldots,{\boldsymbol{y}}_{N_{1}-1}\\}$ belong to
$U_{1}$ and $\\{{\boldsymbol{z}}_{0},\ldots,{\boldsymbol{z}}_{N_{2}-1}\\}$
belong to $U_{2}$. We shall analyze, in particular, $Q_{\bf s,\bf N}$ being
deterministic, interlaced high-order polynomial lattice rules as introduced in
[7] and as considered for affine-parametric operator equations in [8]. To this
end, to generate a polynomial lattice rule in base $b$ with $N_{1}$ points in
$U_{1}$, we need a _generating vector_ of polynomials ${\bf
g}(x)=(g_{1}(x),\ldots,g_{s_{1}}(x))\in[P_{m_{1}}({\mathbb{Z}}_{b})]^{s_{1}}$,
where $P_{m_{1}}({\mathbb{Z}}_{b})$ is the space of polynomials of degree less
than $m_{1}$ in $x$ with coefficients taken from a finite field
${\mathbb{Z}}_{b}$.

For each integer $0\leq n\leq b^{m_{1}}-1$, we associate $n$ with the
polynomial

$n(x)=\sum_{i=1}^{m_{1}}\eta_{i-1}x^{i-1}\quad\in{\mathbb{Z}}_{b}[x],$

where $(\eta_{m_{1}-1},\ldots,\eta_{0})$ is the $b$-adic expansion of $n$,
that is $n=\sum_{i=1}^{m_{1}}\eta_{i-1}\,b^{i-1}\,.$ We also need a map
$v_{m_{1}}$ which maps elements in ${\mathbb{Z}}_{b}(x^{-1})$ to the interval
$[0,1)$, defined for any integer $w$ by

$v_{m_{1}}\left(\sum_{\ell=w}^{\infty}t_{\ell}x^{-\ell}\right)=\sum_{\ell=\max(1,w)}^{m_{1}}t_{\ell}b^{-\ell}.$

Let $P\in{\mathbb{Z}}_{b}[x]$ be an irreducible polynomial with degree
$m_{1}$. The classical polynomial lattice rule
$\mathcal{S}_{P,b,m_{1},s_{1}}({\bf g})$ associated with $P$ and the
generating vector ${\bf g}$ is comprised of the quadrature points

${\boldsymbol{y}}_{n}=\left(v_{m_{1}}\Big{(}\frac{n(x)g_{j}(x)}{P(x)}\Big{)}\right)_{1\leq
j\leq s_{1}}\in[0,1)^{s_{1}},\quad{\rm for}~{}~{}n=0,\ldots,N_{1}-1.$

In a similar fashion, we define the quadrature points
${\boldsymbol{z}}_{n}\in[0,1)^{s_{2}}$ for $n=0,\ldots,N_{2}-1$. In this case,
the _generating vector_ of polynomials is ${\bf g}=(g_{1},\ldots,g_{s_{2}})$.

Classical polynomial lattice rules give almost first order of convergence for
integrands of bounded variation. To obtain high-order of convergence, an
interlacing procedure described as follows is needed. Following [12, 13], the
_digit interlacing function_ with digit interlacing factor
$\alpha\in\mathbb{N}$, $\mathscr{D}_{\alpha}:[0,1)^{\alpha}\to[0,1)$, is
defined by

$\mathscr{D}_{\alpha}(x_{1},\ldots,x_{\alpha})=\sum_{i=1}^{\infty}\sum_{j=1}^{\alpha}\xi_{j,i}b^{-j-(i-1)\alpha}\;,$

where $x_{j}=\sum_{i\geq 1}\xi_{j,i}\,b^{-i}$ for $1\leq j\leq\alpha$. For
vectors, we set $\mathscr{D}^{s}_{\alpha}:[0,1)^{\alpha s}\to[0,1)^{s}$ with

$\mathscr{D}^{s}_{\alpha}(x_{1},\ldots,x_{\alpha
s})=(\mathscr{D}_{\alpha}(x_{1},\ldots,x_{\alpha}),\ldots,\mathscr{D}_{\alpha}(x_{(s-1)\alpha+1},\ldots,x_{s\alpha}))\;.$

Then, an interlaced polynomial lattice rule of order $\alpha$ with $b^{m}$
points in $s$ dimensions is a QMC rule using
$\mathscr{D}_{\alpha}(\mathcal{S}_{P,b,m,\alpha s}({\bf g}))$ as quadrature
points, for some given modulus $P$ and generating vector ${\bf g}$.

Next, we derive the error from approximating the integral ${\mathcal{I}}_{\bf
s}{F}$ in (6.1) by the QMC quadrature formula $Q_{\bf s,\bf N}[{F}]$ in (6.2).
The proof mainly relies on [8, Theorem 3.1].

###### 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 with $\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}}$ and
$\boldsymbol{\varphi}_{s_{1}}=(\varphi_{j})_{1\leq j\leq s_{2}},$ and let
$\alpha\,:=\,\lfloor 1/p\rfloor+1$ and $\beta\,:=\,\lfloor 1/q\rfloor+1$.
Assume that $F$ satisfies the following regularity properties: for any
${\boldsymbol{y}}\in U_{1}$, ${\boldsymbol{z}}\in U_{2}$,
$\boldsymbol{\alpha}\in\\{0,1,\ldots,\alpha\\}^{s_{1}}$, and
$\boldsymbol{\beta}\in\\{0,1,\ldots,\beta\\}^{s_{2}}$, the following
inequalities hold

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

where the constant $c$ is independent of ${\boldsymbol{y}}$,
${\boldsymbol{z}}$, $s_{1}$, $s_{2}$, and of $p$ and $q.$ Then one can
construct two interlaced polynomial lattice rules of order $\alpha$ with
$N_{1}$ points, and of order $\beta$ with $N_{2}$ points, using a fast
component-by-component algorithm, with cost
$\mathcal{O}(\alpha\,s_{1}N_{1}(\log N_{1}+\alpha\,s_{1}))$ and
$\mathcal{O}(\beta\,s_{2}N_{2}(\log N_{2}+\beta\,s_{2}))$ operations,
respectively, so that the following error bound holds

$|{\mathcal{I}}_{\bf s}{F}-Q_{\bf s,\bf N}[{F}]|\leq
C\,\Big{(}N_{1}^{-1/p}+N_{2}^{-1/q}\Big{)}\,,$

where for $i\in\\{1,2\\},$ the generic constant $C$ depends on $b,p$ and $q$,
but is independent of $s_{i}$ and $m_{i}.$

By choosing $m_{1},m_{2}\in\mathbb{N}$ such that $|m_{1}q-m_{2}p|<1$, we
obtain that the total number of QMC points is $N=N_{1}N_{2}=b^{m_{1}+m_{2}}$
and

$|{\mathcal{I}}_{\bf s}{F}-Q_{\bf s,\bf N}[{F}]|\leq C\,N^{-\frac{1}{p+q}}\,.$

Proof. By adding and subtracting
$Q_{s_{1},N_{1}}[{F}(\cdot,{\boldsymbol{z}})]:=\frac{1}{N_{1}}\sum_{j=0}^{N_{1}-1}{F}({\boldsymbol{y}}_{j},{\boldsymbol{z}}),$
the QMC error can be decomposed as

$|{\mathcal{I}}_{\bf s}{F}-Q_{\bf s,\bf
N}[{F}]|\leq\int_{U_{2}}\Big{|}\int_{U_{1}}{F}({\boldsymbol{y}},{\boldsymbol{z}})\,d{\boldsymbol{y}}-Q_{s_{1},N_{1}}[{F}(\cdot,{\boldsymbol{z}})]\Big{|}\,d{\boldsymbol{z}}\\\
+\frac{1}{N_{1}}\sum_{j=0}^{N_{1}-1}\Big{|}\int_{U_{2}}{F}({\boldsymbol{y}}_{j},{\boldsymbol{z}})\,d{\boldsymbol{z}}-Q_{s_{2},N_{2}}[{F}({\boldsymbol{y}}_{j},\cdot)]\Big{|}\,,$

where
$Q_{s_{2},N_{2}}[{F}({\boldsymbol{y}}_{j},\cdot)]:=\frac{1}{N_{2}}\sum_{k=0}^{N_{2}-1}{F}({\boldsymbol{y}}_{j},{\boldsymbol{z}}_{k})$.
By using [8, Theorem 3.1] and the regularity assumptions in (6.3), we have

$\Big{|}\int_{U_{1}}{F}({\boldsymbol{y}},{\boldsymbol{z}})\,d{\boldsymbol{y}}-Q_{s_{1},N_{1}}[{F}(\cdot,{\boldsymbol{z}})]\Big{|}\leq
C\,N_{1}^{-1/p},$

and

$\Big{|}\int_{U_{2}}{F}({\boldsymbol{y}}_{j},{\boldsymbol{z}})\,d{\boldsymbol{z}}-Q_{s_{2},N_{2}}[{F}({\boldsymbol{y}}_{j},\cdot)]\Big{|}\leq
CN_{2}^{-1/q}\,.$

Combining the above equations, we immediately deduce the first desired
results.

Now, using $N=b^{m_{1}+m_{2}}$ and the conditions $|m_{1}q-m_{2}p|<1,$ we have

$N_{1}^{-1/p}=N^{-1/(p+q)}b^{(m_{1}+m_{2})/(p+q)-m_{1}/p}\\\
=N^{-1/(p+q)}b^{(m_{2}p-m_{1}q)/(p(p+q))}\leq N^{-1/(p+q)}b^{1/(p(p+q))}\leq
CN^{-1/(p+q)}\,.$

Similarly,

$N_{2}^{-1/q}\leq N^{-1/(p+q)}b^{1/(q(p+q))}\leq CN^{-1/(p+q)}\,,$

and therefore, the proof of the second desired estimate is completed.
$\quad\Box$

In order to reduce the computational cost (and thereby improving the
convergence rate), we next discuss a combination of the QMC rules with a
sparse grid approach. Let $\\{N_{i}^{(j)}\\}_{j\geq 1}$ be increasing
sequences of positive values for $i=1,2.$ Then

$\mathcal{I}_{\bf s}{F}=\lim_{j,k\to\infty}Q_{\mathbf{s},{\bf
N}^{j,k}}[{F}],\quad{\rm with}~{}~{}~{}{\bf
N}^{j,k}=(N_{1}^{(j)},N_{2}^{(k)}).$

We can write this as a telescoping sum

$\mathcal{I}_{\bf s}{F}=\sum_{j,k=1}^{\infty}a_{jk},~{}{\rm
with}~{}a_{jk}=\Big{(}Q_{\mathbf{s},{\bf N}^{j,k}}-Q_{\mathbf{s},{\bf
N}^{j-1,k}}-Q_{\mathbf{s},{\bf N}^{j,k-1}}+Q_{\mathbf{s},{\bf
N}^{j-1,k-1}}\Big{)}[{F}]$ (6.4)

where

$Q_{\mathbf{s},{\bf N}^{j,k}}=0~{}~{}{\rm
if}~{}~{}(j,k)\in\\{(\zeta,0),(0,\omega):\zeta,\omega\in\mathbb{N}\cup\\{0\\}\\}\,.$
(6.5)

To get a computable quantity, we need to truncate the infinite sums in (6.4).
This can be done in different ways; we choose to truncate the tensor grid
along the main diagonal of indexed points for each combination of the QMC
rules ${\bf N}^{j,k}$ where $L$ is the “level” of the sparse grid rule.
Explicitly, we truncate as follows:

${\mathcal{I}}_{{\bf s},L}[{F}]=\sum_{\begin{subarray}{c}j,k=1\\\ j+k\leq
L\end{subarray}}a_{jk}=\sum_{k=1}^{L-1}\left(Q_{\mathbf{s},{\bf
N}^{L-k,k}}[{F}]-Q_{\mathbf{s},{\bf N}^{L-k,k-1}}[{F}]\right),\quad$ (6.6)

where in the second equality we used $\sum_{\begin{subarray}{c}j,k=1\\\
j+k\leq L\end{subarray}}=\sum_{k=1}^{L-1}\sum_{j=1}^{L-k}$ and the condition
in (6.5).111In principle one could use a QMC rule
$\frac{1}{N}\sum_{n=0}^{N-1}F(\boldsymbol{x}_{n})$ in dimension $s_{1}+s_{2}$
directly, i.e. ${\boldsymbol{x}}_{n}\in[0,1]^{s_{1}+s_{2}}$, without combining
QMC with sparse grids. In this approach we would have to combine the different
weights arising from simulating $\mu$ and $\lambda$ and we would have to use
the same number of points for each part. The sparse grid approach gives us
more flexibility in that sense, and since the weights for both parts are of a
similar form as in other problems [8], we can reuse existing constructions of
higher order polynomial lattice rules. We prove next that the quadrature error
incurred on this QMC sparse grid is relatively small for a sufficiently large
$L$.

In the next theorem, for some $\vartheta>0$ such that
$p\vartheta,q\vartheta\geq 1$, we choose $N_{1}^{(j)}=b^{\lceil
jp\vartheta\rceil}$ and $N_{2}^{(j)}=b^{\lceil jq\vartheta\rceil}$ for $j\geq
1$. The purpose of $\vartheta>0$ is to avoid a situation where
$N_{i}^{(j)}=N_{i}^{(j+1)}$ for some admissible $i,j$. Choosing $\vartheta$
such that $p\vartheta,q\vartheta\geq 1$ guarantees that this cannot happen. On
the other hand, since the constant $C$ increases with $\vartheta$, we consider
$\vartheta$ as a constant. In other words, in order to reduce the error in
Theorem 6.2 one increases $L$ and therefore $M$, but keeps $\vartheta$ fixed.

###### Theorem 6.2

In addition to the assumptions of Theorem 6.1, we assume that

$|\partial^{\boldsymbol{\beta},\boldsymbol{\alpha}}_{{\boldsymbol{z}},{\boldsymbol{y}}}{F}({\boldsymbol{y}},{\boldsymbol{z}})|\leq
c(|\boldsymbol{\alpha}|+|\boldsymbol{\beta}|)!\boldsymbol{\chi}_{s_{1}}^{\boldsymbol{\alpha}}\boldsymbol{\varphi}_{s_{2}}^{\boldsymbol{\beta}}\,.$
(6.7)

Then we have

$|{\mathcal{I}}_{\bf s}[{F}]-{\mathcal{I}}_{{\bf s},L}[{F}]|\leq
CL\Big{(}N_{1}^{(L)}\Big{)}^{-1/p}+CL\Big{(}N_{2}^{(L)}\Big{)}^{-1/q},$

where the constant $C$ depends on $b,p,q,\vartheta$, but is independent of
$s_{1}$, $s_{2}$ and $L$.

Let $M$ denote the total number of quadature points used in the QMC sparse
grid rule. Then

$|{\mathcal{I}}_{\bf s}[{F}]-{\mathcal{I}}_{{\bf s},L}[{F}]|\leq C(\log
M)^{1+1_{p=q}/p}M^{-\min(1/p,1/q)},$

where $1_{p=q}$ is $1$ if $p=q$ and $0$ otherwise, and where the constant $C$
depends on $b,p,q,\vartheta$, but is independent of $s_{1}$, $s_{2}$ and $M$.

Proof. From (6.4) and (6.6),

$|{\mathcal{I}}_{\bf s}[{F}]-{\mathcal{I}}_{{\bf
s},L}[{F}]|\leq\sum_{\begin{subarray}{c}j,k=1\\\ j+k\geq
L+1\end{subarray}}^{\infty}|a_{jk}|\leq\sum_{k=L}^{\infty}\sum_{j=1}^{\infty}|a_{jk}|+\sum_{k=1}^{L-1}\,\,\sum_{j=L-k+1}^{\infty}|a_{jk}|\,.$
(6.8)

To estimate the first term on the right hand side of the above equation, and
for brevity, we let ${\bf
g}_{j}({\boldsymbol{z}})=Q_{s_{1},N_{1}^{(j)}}[{F}(\cdot,{\boldsymbol{z}})]-Q_{s_{1},N_{1}^{(j-1)}}[{F}(\cdot,{\boldsymbol{z}})]$.
Then,

$|a_{jk}|=\left|Q_{s_{2},N_{2}^{(k)}}[\mathbf{g}_{j}]-Q_{s_{2},N_{2}^{(k-1)}}[\mathbf{g}_{j}]\right|\,.$

Adding and subtracting $\int_{U_{2}}{\bf
g}_{j}({\boldsymbol{z}})\,d{\boldsymbol{z}}$ and using [8, Theorem 3.1], we
obtain

$|a_{jk}|\leq\sum_{\ell=k-1}^{k}\Big{|}Q_{s_{2},N_{2}^{(\ell)}}[\mathbf{g}_{j}]-\int_{U_{2}}{\bf
g}_{j}({\boldsymbol{z}})\,d{\boldsymbol{z}}\Big{|}\leq
C\sup_{{\boldsymbol{z}}\in
U_{2}}|\partial^{\boldsymbol{\beta}}_{\boldsymbol{z}}{\bf
g}_{j}({\boldsymbol{z}})|\,\Big{(}N_{2}^{(k-1)}\Big{)}^{-1/q}\,.$

Another application of [8, Theorem 3.1] but on the region $U_{1}$, where the
regularity assumption in (6.7) is needed here, yields

$|\partial^{\boldsymbol{\beta}}_{\boldsymbol{z}}{\bf
g}_{j}({\boldsymbol{z}})|=\Big{|}Q_{s_{1},N_{1}^{(j)}}[\partial^{\boldsymbol{\beta}}_{\boldsymbol{z}}{F}(\cdot,{\boldsymbol{z}})]-Q_{s_{1},N_{1}^{(j-1)}}[\partial^{\boldsymbol{\beta}}_{\boldsymbol{z}}{F}(\cdot,{\boldsymbol{z}})]\Big{|}\\\
\leq\sum_{\ell=j-1}^{j}\Big{|}Q_{s_{1},N_{1}^{(\ell)}}[\partial^{\boldsymbol{\beta}}_{\boldsymbol{z}}{F}(\cdot,{\boldsymbol{z}})]-\mathcal{I}_{U_{1}}[\partial^{\boldsymbol{\beta}}_{\boldsymbol{z}}{F}(\cdot,{\boldsymbol{z}})]\Big{|}\leq
C\,\Big{(}N_{1}^{(j-1)}\Big{)}^{-1/p}\,.$

Inserting this in the previous equation leads to $|a_{jk}|\leq
C\,\Big{(}N_{2}^{(k-1)}\Big{)}^{-1/q}\Big{(}N_{1}^{(j-1)}\Big{)}^{-1/p},$ and
consequently,

$|{\mathcal{I}}_{\bf s}[{F}]-{\mathcal{I}}_{{\bf s},L}[{F}]|\leq
C\sum_{k=L}^{\infty}\Big{(}N_{2}^{(k-1)}\Big{)}^{-1/q}\sum_{j=1}^{\infty}\Big{(}N_{1}^{(j-1)}\Big{)}^{-1/p}\\\
+C\,\sum_{k=1}^{L-1}\Big{(}N_{2}^{(k-1)}\Big{)}^{-1/q}\sum_{j=L-k+1}^{\infty}\Big{(}N_{1}^{(j-1)}\Big{)}^{-1/p}\,.$

Using $N_{1}^{(j)}=b^{\lceil jp\vartheta\rceil}$ and $N_{2}^{(j)}=b^{\lceil
jq\vartheta\rceil}$, we notice that

$\displaystyle|{\mathcal{I}}_{\bf s}[{F}]-{\mathcal{I}}_{{\bf s},L}[{F}]|\leq$
$\displaystyle
C\,\sum_{k=L-1}^{\infty}b^{-k\vartheta}\sum_{j=0}^{\infty}b^{-j\vartheta}+C\,\sum_{k=0}^{L-2}b^{-k\vartheta}\sum_{j=L-k}^{\infty}b^{-j\vartheta}$
$\displaystyle\leq$ $\displaystyle
Cb^{-L\vartheta}\left(\frac{b^{\vartheta}}{(1-b^{-\vartheta})^{2}}+\frac{L}{1-b^{-\vartheta}}\right).$

Hence the desired QMC sparse grid estimate is obtained.

The total number of quadrature points used in the combined QMC sparse grid
approach is

$M=\sum_{k=1}^{L-1}b^{\lceil(L-k)p\vartheta\rceil}b^{\lceil
kq\vartheta\rceil}\leq b^{2+Lp\vartheta}\sum_{k=1}^{L-1}b^{k(q-p)\vartheta}.$

For $q=p$ we have $M\leq b^{2+Lp\vartheta}(L-1),$ and for $q\neq p$ we have

$\displaystyle M\leq
b^{2+Lp\vartheta}\frac{b^{(q-p)\vartheta}-b^{L(q-p)\vartheta}}{1-b^{(q-p)\vartheta}}\leq\frac{b^{2+L\max(p,q)\vartheta-|p-q|\vartheta}}{1-b^{-|p-q|\vartheta}}\leq
C\,b^{L\max(p,q)\vartheta}\,.$

Since the error is of order $Lb^{-L\vartheta}$ we have

$|{\mathcal{I}}_{\bf s}[{F}]-{\mathcal{I}}_{{\bf s},L}[{F}]|\leq
CLb^{-L\vartheta}\leq CL^{1+1_{p=q}/p}M^{-\min(1/p,1/q)},$

where $1_{p=q}$ is $1$ if $p=q$ and $0$ otherwise. Since $M\geq
b^{(L-1)p\vartheta}$ we have $L\leq C\log M$, hence the result follows.
$\quad\Box$

##  7 Numerical experiments

In this section, we illustrate numerically the theoretical finding in Theorem
2.1. In all experiments, the physical domain $\Omega$ is chosen to be the unit
square $[0,1]^{2}$, and $\mathcal{T}_{h}$ is a family of uniform triangular
meshes with diameter $\sqrt{2}h$ obtained from uniform $J$-by-$J$ square
meshes by cutting each mesh square into two congruent triangles with $h=1/J.$
In all numerical experiments we set the base of the polynomial lattice rules
$b=2$.

Figure 1: A comparison between the body force ${\bf
f}=(x_{1}-1/2,x_{2}(x_{2}-1/2))$ (left) and the displacement ${\bf u}_{h}$
(right) with $J=15.$

Example 1: In this example, we corroborate the finite element errors and
convergence rates when $r=1$ (piecewise linear Galerkin FEM) the Lamé
parameters $\mu$ and $\lambda$ are variable but deterministic. We choose
$\mu(x_{1},x_{2})=x_{1}+x_{2}+1$ and $\lambda(x_{1},x_{2})=\sin(2\pi
x_{1})+2.$ In Figure 1, we compare graphically between the body force load
${\bf f}$ and the approximate displacement ${\bf u}_{h}$ for the case of
homogeneous Dirichlet boundary conditions. To illustrate the FEM convergence
order in Theorem 5.1 (or Theorem 2.1), we choose the body force ${\bf f}$ so
that the exact solution is

${\bf u}(x_{1},x_{2})=\left[\begin{matrix}u_{1}(x_{1},x_{2})\\\
u_{2}(x_{1},x_{2})\end{matrix}\right]=\left[\begin{matrix}2(\cos(2\pi
x_{1})-1)\sin(2\pi x_{2})\\\ (1-\cos(2\pi x_{2}))\sin(2\pi
x_{1})\end{matrix}\right]\,.$

Motivated by the equality

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

we define, for some fixed (but arbitrary) ${\bf w}\in{\bf L}^{2}(\Omega)$, the
functional $\mathcal{L}$ by: $\mathcal{L}({\bf v})=\mathcal{L}_{\bf w}({\bf
v}):=\langle{\bf w},{\bf v}\rangle$. Then, by using the convergence estimate
(5.5) in Theorem 5.1, we have $|\langle{\bf u}-{\bf u}_{h},{\bf w}\rangle|\leq
C\,h^{2}\|{\bf f}\|\|{\bf w}\|$ for ${\bf w}\in{\bf L}^{2}(\Omega)\,.$
Consequently, the equality in (7.1) leads to the following optimal ${\bf
L}^{2}(\Omega)$ estimate: $E_{h}:=\|{\bf u}-{\bf u}_{h}\|\leq C\,h^{2}\|{\bf
f}\|\,.$ To demonstrate this numerically, we compute $E_{h}$ by approximating
the ${\bf L}^{2}$-norm ($\|\cdot\|$) using the centroids of the elements in
the mesh $\mathcal{T}_{h}$. The empirical convergence rate (CR) is calculated
by halving $h$, and thus, $\text{CR}=\log_{2}(E_{h}/E_{h/2}).$

If we choose ${\bf w}={\bf 1}$ (the unitary constant vector function), then
with ${\bf v}=[v_{1}~{}~{}v_{2}]^{T},$

$\mathcal{L}({\bf v})=\mathcal{L}_{\bf 1}({\bf v})=\int_{\Omega}{\bf
v}({\boldsymbol{x}})\cdot{\bf
1}\,d{\boldsymbol{x}}=\int_{\Omega}[v_{1}({\boldsymbol{x}})+v_{2}({\boldsymbol{x}})]\,d{\boldsymbol{x}}\,,$

which is the mean of ${\bf v}$ over $\Omega=[0,1]^{2}.$ Since
$\|\mathcal{L}_{\bf 1}\|\leq 1,$ by Theorem 5.1,

$|\mathcal{L}({\bf u}-{\bf u}_{h})|=|\mathcal{L}_{\bf 1}({\bf u}-{\bf
u}_{h})|=\Big{|}\sum_{i=1}^{2}\int_{\Omega}(u_{i}-u_{i_{h}})\,d{\boldsymbol{x}}\Big{|}\leq
C\,h^{2}\|{\bf f}\|\,.$

Again, we use the centroids of the elements in the mesh $\mathcal{T}_{h}$ to
approximate the above integral. The reported numerical (empirical) convergence
rates in Table 1 illustrate the expected second order of accuracy. For a
graphical illustration of the efficiency of the approximate solution over the
global domain $\Omega$, we highlight the pointwise nodal displacement errors
in Figure 2 for $J=60.$

$J$ |  $\|{\bf u}-{\bf u}_{h}\|$ CR |  $|\mathcal{L}_{\bf 1}({\bf u}-{\bf u}_{h})|$ CR
---|---|---
8 | 3.8533e-01 |  | 1.1697e-02 |
16 | 1.1163e-01 | 1.7873 | 3.7017e-03 | 1.6599
32 | 2.9204e-02 | 1.9345 | 9.8934e-04 | 1.9037
64 | 7.3903e-03 | 1.9825 | 2.5179e-04 | 1.9743
128 | 1.8533e-03 | 1.9955 | 6.3238e-05 | 1.9933
Table 1: Example 1, Errors and empirical convergence rates for different
values of $J.$

Figure 2: Pointwise nodal errors in the displacement, $|u_{1}-u_{1_{h}}|$ on
right and $|u_{2}-u_{2_{h}}|$ on left.

Example 2: This example is devoted to confirm the QMC theoretical convergence
results when $\mu$ is random and $\lambda$ is constant. More precisely,
$\lambda=1$ and

$\mu({\boldsymbol{x}},{\boldsymbol{y}})=\frac{1}{10}\Big{(}1+\sum_{j=1}^{\infty}y_{j}\psi_{j}\Big{)},\quad{\rm
with}~{}~{}\psi_{j}=\frac{1}{j^{2}}\sin(j\pi x_{1})\sin((2j-1)\pi x_{2}),$

and for $y_{j}\in[-1/2,1/2].$ Since $\|\psi_{j}\|\leq{1}/{(2j^{2})},$
$\sum_{j=1}^{\infty}\|\psi_{j}\|^{p}$ is convergent for $p>1/2$ and
$\sum_{j=s_{1}+1}^{\infty}\|\psi_{j}\|\leq Cs_{1}^{-1}\,.$ Thus, (A3) and (A4)
are satisfied when $p>1/2$ and $\varrho_{1}=1/2,$ respectively. We discretize
on the physical domain using the quadratic FEM. Therefore, according to
Theorem 2.1, we expect the truncated QMC Galerkin finite element error to be
of order $\mathcal{O}(s_{1}^{-1}+\log(s_{1})(N_{1}^{-2}+h^{3}))$. The
appearance of the logarithmic factor $\log(s_{1})$ in front of $N_{1}^{-2}$
and $h^{3}$ is due to the facts that $\sum_{j=1}^{s_{1}}\|\psi_{j}\|^{1/2}\leq
C\log(s_{1})$ and that $\sum_{j=1}^{s_{1}}\|\nabla\psi_{j}\|\leq
C\log(s_{1})$, respectively. For measuring the error, and since the exact
solution is unknown, we rely on the reference solution $\Xi_{{\bf u}_{h}^{*}}$
which is computed using $s_{1}=256$, $J=128,$ and $1024$ QMC points. Hence, by
ignoring the logarithmic factor $\log(s_{1})$, and in the absence of the
truncated series error, we anticipate $\mathcal{O}(N_{1}^{-2})$-rates of
convergence for $N_{1}\leq J^{\frac{3}{2}}$, with $N_{1}\ll 1024$. This is
illustrated numerically in Table 2 and graphically in Figure 3 for different
values of $N_{1}$, and for fixed $s_{1}=256$ and $J=128$, with
$\mathcal{L}=\mathcal{L}_{\bf 1}$ and ${\bf f}=(2x_{1}+10,x_{2}-3).$ Note that
the middle column of Table 2 displays $|\Xi_{{\bf u}^{*}_{h}}-\Xi_{{{\bf
u}_{h}},Q}|$ where $Q$ is a quadrature for $\mathcal{L}=\mathcal{L}_{1}$ with
$N_{2}=1$ and $N_{1}$ varying, see (2.3).

$N_{1}$ | $|\Xi_{{\bf u}^{*}_{h}}-\Xi_{{{\bf u}_{h}},Q}|$ | CR
---|---|---
8 | 3.1045e-01 |
16 | 6.1906e-02 | 2.3262
32 | 1.5191e-02 | 2.0269
64 | 4.3387e-03 | 1.8079
128 | 9.3546e-04 | 2.2135
256 | 2.5232e-04 | 1.8904
Table 2: Example 2, errors and convergence rates for different values of
$N_{1}$. Figure 3: Numerical errors (err${}_{N_{1}}$) vs. $N_{1}^{-2}$ for
Example 2

Example 3: In this example, we focus on the randomness in $\lambda$ while
$\mu=1$. We choose

$\lambda({\boldsymbol{x}},{\boldsymbol{z}})=1+\sum_{j=1}^{\infty}\frac{z_{j}}{j^{2}}\sin(j\pi
x_{1})\sin((2j-1)\pi x_{2}),\quad z_{j}\in[-1/2,1/2].$

By arguing as in the preceding example, based on Theorem 2.1, we fix
$s_{2}=256$, $J=128$ and $r=2,$ then the QMC Galerkin finite element error is
expected to be of order $\mathcal{O}(N_{2}^{-2})$ whenever $N_{2}\leq J$,
where the logarithmic factor $\log(s_{2})$ is ignored. We rely again on the
reference solution $\Xi_{{\bf u}_{h}^{*}}$, which is computed as in the
previous example, in measuring the errors, and consequently, the convergence
rates. As expected, an $\mathcal{O}(N_{2}^{-2})$ convergence rate is
illustrated tabularly and graphically for different values of $N_{2}$ in Table
3 and Figure 4, respectively, for fixed $s_{2}=J=256,$ with
$\mathcal{L}=\mathcal{L}_{\bf 1}$ and ${\bf f}=(2x_{1}+10,x_{2}-3).$ Note that
the middle column of Table 2 displays $|\Xi_{{\bf u}^{*}_{h}}-\Xi_{{{\bf
u}_{h}},Q}|$ where $Q$ is a quadrature for $\mathcal{L}=\mathcal{L}_{1}$ with
$N_{1}=1$ and $N_{2}$ varying, see (2.3).

$N_{2}$ | $|\Xi_{{\bf u}^{*}_{h}}-\Xi_{{{\bf u}_{h}},Q}|$ | CR
---|---|---
8 | 7.5520e-04 |
16 | 2.0011e-04 | 1.9161
32 | 4.5223e-05 | 2.1457
64 | 1.1630e-05 | 1.9592
128 | 2.7057e-06 | 2.1038
256 | 6.5202e-07 | 2.0530
Table 3: Example 3, errors and convergence rates for different values of
$N_{2}$. Figure 4: Numerical errors (err${}_{N_{2}}$) vs. $N_{2}^{-2}$ for
Example 3

Example 4: The aim behind this example is to illustrate numerically the
achieved QMC sparse grid convergence results in Theorems 2.1 (second part) and
6.2. As before, we set the body force ${\bf f}=(2x_{1}+10,x_{2}-3)$ but now
both coefficients $\lambda$ and $\mu$ are random. We choose

$\displaystyle\mu({\boldsymbol{x}},{\boldsymbol{y}})$
$\displaystyle=1+\sum_{j=1}^{\infty}\frac{y_{j}}{j^{2}}\sin(j\pi
x_{1})\sin((2j-1)\pi x_{2}),\quad y_{j}\in[-1/2,1/2],$
$\displaystyle\lambda({\boldsymbol{x}},{\boldsymbol{z}})$
$\displaystyle=1+\sum_{j=1}^{\infty}\frac{z_{j}}{j^{2}}\sin(j\pi
x_{1})\sin((2j-1)\pi x_{2}),\quad z_{j}\in[-1/2,1/2],$

and so, $p=q=1/2$ (note that strictly speaking we have $p=q=1/2+\varepsilon$
for an arbitrary $\varepsilon>0$; in order to simplify the computation we
ignore this technicality in the following). We fix the truncation degree
$s_{1}=s_{2}=256$, the spatial mesh element size $J=128,$ and the degree of
the Galerkin FEM $r=2.$ The reference solution $\Xi_{{\bf u}_{h}^{*}}$ is
generated using a full grid of $2048\times 2048$ (that is, $b=2$ and
$m_{1}=m_{2}=11$) high-order QMC points (generated by a Python package in
[11]). The combined QMC sparse grid algorithm (6.6) (with ${\bf
N}^{L-k,k}=(2^{L-k},2^{k})$, that is, $\vartheta=2$) is implemented to compute
$\Xi_{{\bf u}_{h},Q_{L}}$. The errors between approximation $\Xi_{{\bf
u}_{h},Q_{L}}$ and the reference solution for different values of $L$ are
given in the second column of Table 4. The fourth column of the table gives
the QMC sparse grid upper error bounds $4(\log(M))^{3}M^{-2}$ predicted by
Theorem 6.2 (where the constant $C$ in the error bound is ignored). Note that
the PDE solvers can be run in parallel for distinct QMC points. To speed up
the computation, finite element PDE solvers based on examples in the FEniCS
package [19] are used on the high-performance computing platform Katana [16]
provided by UNSW, Sydney. The Python code used in the numerical experiments
together with the PBS scripts is available at
https://github.com/qlegia/Elasticity-HigherOrder-QMC.

$L$ | $M=(L-1)2^{L}$ | $|\Xi_{{\bf u}_{h},Q_{L}}-\Xi_{{\bf u}_{h}^{*}}|$ | $4(\log(M))^{3}M^{-2}$
---|---|---|---
9 | 4096 | 3.5687e-06 | 1.3720e-04
10 | 9216 | 3.4606e-06 | 3.5826e-05
11 | 20480 | 7.4782e-06 | 9.3300e-06
12 | 45056 | 1.1456e-06 | 2.4244e-06
13 | 98304 | 6.8586e-08 | 6.2884e-07
14 | 212992 | 1.4453e-07 | 1.6284e-07
15 | 458752 | 5.4749e-08 | 4.2108e-08
Table 4: Example 4, numerical and theoretical error results for QMC sparse
grid algorithm.

## References

  * [1] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Elasticity Theory, Cambridge University Press, New York, 2007.
  * [2] S. C. Brenner and Li-Y. Sung, Linear finite element methods for planar linear elasticity, Math. Comp., 59, 321–338 (1992).
  * [3] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Third Edition, Springer, 2008.
  * [4] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods. Springer-Verlag, New York, 1991.
  * [5] P. G. Ciarlet, Mathematical Elasticity, Volume I: Three-Dimensional Elasticity, North-Holland, Amsterdam, 1988.
  * [6] A. Cohen, R. De Vore and Ch. Schwab, Convergence rates of best N-terrn Galerkin ap proximations for a class of elliptic sPDEs, Found. Comput. Math., 10, 615–646 (2010).
  * [7] J. Dick, Walsh spaces containing smooth functions and Quasi-Monte Carlo rules of arbitrary high order, SIAM J. Numer. Anal., 46, 1519–1553 (2008).
  * [8] 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).
  * [9] M. Eigel, C. J. Gittelson, Ch. Schwab and E. Zander, Adaptive stochastic Galerkin FEM, Comput. Methods Appl. Mech. Eng., 270, 247–269 (2014).
  * [10] R. S. Falk, Nonconforming finite element methods for the equations of linear elasticity, Math. Comp., 57, 529–550 (1991).
  * [11] 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.
  * [12] T. Goda, Good interlaced polynomial lattice rules for numerical integration in weighted Walsh spaces, J. Comput. Appl. Math., 285, 279–294 (2015).
  * [13] T. Goda and J. Dick, Construction of interlaced scrambled polynomial lattice rules of arbitrary high order, Foundations of Computational Mathematics, 15, 1245–1278 (2015).
  * [14] P. Grisvard, Elliptic Problems in Nonsmooth Domains. Monographs and Studies in Mathematics, 24. Pitman (Advanced Publishing Program), Boston, MA, 1985.
  * [15] V. H. Hoang, T. C. Nguyen and B. Xia, Polynomial approximations of a class of stochastic multiscale elasticity problems, Z. Angew. Math. Phys., 67, 67–78 (2016).
  * [16] D. Smith and L. Betbeder-Matibet, Katana, 2010,https://doi.org/doi:10.26190/669X-A286, https://doi.org/10.26190/669x-a286.
  * [17] 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).
  * [18] 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).
  * [19] H. P. Langtangen and A. Logg, Solving PDEs in Python, Springer, 2017.
  * [20] H. G. Matthies, C. Brenner, C. Bucher, and C. G. Soares, Uncertainties in probabilistic numerical analysis of structures and solid–stochastic ﬁnite elements, Struct. Safety, 19, 283–336 (1997).
  * [21] W. McLean, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge University Press, 2000.
  * [22] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials, Math. Mod. Numer. Anal., 19, 111–143 (1985).
  * [23] 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).
  * [24] B. Xia and V. H. Hoang, Best N-term GPC approximations for a class of stochastic linear elasticity equations, Math. Models and Methods in Applied Sciences, 24, 513–552 (2014).