# Random set solutions to elliptic and hyperbolic partial differential
equations

Jelena Karakašević Unit of Engineering Mathematics, University of Innsbruck,
Technikerstraße 13, 6020 Innsbruck, Austria,
<EMAIL_ADDRESS>Michael Oberguggenberger Unit of
Engineering Mathematics, University of Innsbruck, Technikerstraße 13, 6020
Innsbruck, Austria<EMAIL_ADDRESS>
###### Abstract

The past decades have seen increasing interest in modelling uncertainty by
heterogeneous methods, combining probability and interval analysis, especially
for assessing parameter uncertainty in engineering models. A unifying
mathematical framework admitting the combination of a wide range of such
methods is the theory of random sets, describing input and output of a
structural model by set-valued random variables. The purpose of this paper is
to highlight the mathematics behind this approach. The modelling and
computational implications are discussed and demonstrated with the help of
prototypical partial differential equations—a scalar elliptic equation from
elastostatics and hyperbolic systems arising in elastodynamics.

Keywords. Partial differential equations; random fields; interval parameters;
random sets; elastostatics; elastodynamics

##  1 Introduction

Many models in structural engineering are cast in the form of (systems of)
ordinary or partial differential equations. Static problems, e.g. from
elasticity theory, give rise to elliptic partial differential equations.
Dynamic problems (e.g. vibrations of elastic bodies or acoustic wave
propagation) typically are formulated in terms of hyperbolic partial
differential equations or systems. The quantities of interest are solutions to
(elliptic) equilibrium equations or to (hyperbolic) equations of motion. These
equations contain coefficients that—as a rule—depend on further parameters,
such as material constants or geometric dimensions. In addition, the
parameters may vary in space and time. An important task in engineering design
is the analysis of the uncertainty of the computed output. This uncertainty
traces back to uncertainty in the input parameters, which may be due to many
sources: random fluctuations, lack of information, conflicting assessments,
random measurement errors, systematic measurement errors, fluctuations due to
spatial inhomogeneity, errors made by assigning parameter status to state
variables, model insufficiency, just to name a few. Uncertainty in the model
itself is important, but will be left aside in this paper.

Available information on parameter variability may consist in frequency
distributions obtained from large samples, values from small samples or single
measurements, interval bounds, experts’ point estimates, or educated guesses
from experience.

Accordingly, an increasing desire has risen in the engineering community to
introduce methods of quantifying the uncertainty beyond probability, such as
interval analysis, set-valued models, fuzzy sets, evidence theory, random
sets, sets of probability measures, imprecise probability, lower and upper
previsions, info-gap-analysis, etc. This paper focusses on combinations of
probability and interval analysis. This includes random variables, random
fields, intervals, interval fields, random variables with interval parameters
and random fields with interval parameters (such as mean, variance,
correlation length).

As already noted in [1, 50, 51], all these types and their combinations can be
accommodated in the framework of _random sets_. Random sets are set-valued
random variables, whereby the values may be sets of real numbers, subsets of
finite dimensional vector spaces, or subsets of function spaces. For example,
random variables are single-valued random sets and intervals are random sets
on a one-point probability space. The purpose of this paper is to highlight
some important mathematical modelling questions that arise when using random
sets. The considerations will be exemplified with the help of a prototypical
elliptic equation and a class of hyperbolic systems.

When uncertainties in the coefficient functions of ordinary or partial
differential equations are modelled as random sets, random fields with
interval parameters, or (finite) interval fields, the solutions of these
equations are no longer single-valued, but rather random sets with values in
function spaces. Evaluating the solution at a fixed point in space or time
will result in a real-valued random set. However, here an important issue
arises: in order to properly qualify as a random set, the set-valued solution
has to satisfy a certain measurability condition (to be detailed in Section
2.3). It turns out that the essential ingredient for this condition to hold is
that the solution should depend continuously on the coefficients of the
equation, and in turn also on the hyperparameters of the random variables or
random fields defining the coefficients. This continuous dependence is already
needed in pure interval modelling to prove that the values of the solutions
are intervals. It is also needed in probabilistic modelling to infer that the
solutions are random variables or random fields. In these two (classical)
cases, the arguments are rather straightforward. However, in models combining
interval and probabilistic inputs, additional subtleties are involved. One of
the main goals of the present article is to formulate a rigorous mathematical
framework that allows one to infer that the solutions are random sets. The
elliptic and the hyperbolic case require different arguments, the essence of
which will be worked out in the corresponding sections.

Here is a short review of relevant literature on uncertainty quantification in
engineering models, combining interval methods with probabilistic methods.
Adopting the language of [32], the structural model can be cast in the form
$Z=g(X,Y)$ where $X$ and $Y$ refer to different types of inputs in a hybrid
probability-interval uncertainty analysis; $Z$ is the system response or
output, $g$ is the performance function or input-output map. The map $g$ can
be given, for example, by an ordinary or partial differential equation, a
finite element model, or a black box algorithm. Here $X$ refers to the
probabilistic part, $Y$ to the interval part—both may be multi-dimensional. If
both parts are present, the output is necessarily both random and set-valued.
There are three subcases. In the first case, $X$ and $Y$ are considered as
separate, non-interacting entities. The set-valued output arises through their
combination in the input-output map $g$. In the second case, $X=X(Y)$, that
is, $X$ is a random variable or random field with interval hyperparameters.
This is often referred to as an imprecise random variable or an imprecise
random field. The third possibility is that $Y=Y(X)$, the interval bounds are
random variables. This is usually referred to as a random interval, a special
case of a random set.

Typical representative articles concerning the first case are [42, 63], where
$n$-DOF structures with interval material parameters under stochastic
excitation have been studied. Structures described by stochastic ordinary
differential equations with random interval parameters have been addressed in
[60], while wave equations with interval propagation speed under space-time
white noise excitation have been addressed in [53]. Combinations of intervals
and random vectors as input parameters in beams and truss structures have been
considered in [22]. The paper [38] merges intervals and random variables by
means of evidence theory (finite random sets) in an application to squeal
analysis of a brake system.

The second case extends from imprecise random variables as material parameters
in frame and grid structures [64] to finite element models with imprecise
random field inputs [11, 16], employing a Karhunen-Loève expansion with
interval coefficients. The applications in Sections 4 and 5 of this paper are
in the same spirit. A combination of inputs consisting of imprecise random
variables and imprecise random fields can be found in [50].

As mentioned, the third case concerns random interval input and is addressed
in the finite dimensional version in [22, 51], in combination with stochastic
differential equations in the quoted paper [60] and in [55, 56].

There are various papers which contrast and compare interval modelling and
probabilistic modelling side by side in beam and plate structures [43, 62, 64,
65]. The question of identification and model calibration under hybrid
uncertainty has been addressed in the literature as well. A procedure for
identifying interval models from measurements has been proposed in [15].
Constructing random sets from statistical data has been considered in [51]. In
the recent paper [36] a Bayesian-type updating procedure for random variables
with interval constraints has been presented, while a whole range of
calibration methods has been combined in [27] for quantification of mixed
uncertainties, reliability analysis and design of structures.

Survey papers on modelling uncertainty in engineering by hybrid probability-
interval approaches are [6, 48, 32]. A recent review paper on numerical
methods for uncertainty propagation in probability-box models is [14]. Further
details and references on numerical algorithms can be found in Subsection 6.2.
For the theory of random sets, their interpretations and applications the
reader is referred to the monographs [9, 41, 46], notably [7] for engineering
applications.

The plan of the paper is as follows. Section 2 will be devoted to the required
mathematical set-up. These facts are known, but it appears useful to elaborate
them in some detail for the later applications. Section 3 singles out the
principal mathematical argument that allows one to prove that the output
quantity is a random set. The combination of arguments giving a full proof of
this result seems to be new. Section 4 details the considerations for elliptic
equations as an application of the findings of Section 3, including a simple
explicit numerical example. Section 5 deals with hyperbolic systems and
elaborates the continuity argument in the hyperbolic case, also having a short
look at hyperbolic equations and systems in higher space dimensions. Section 6
contains a comparison with the parametric approach and numerical aspects, as
well as some remarks on modelling correlations. Detailed proofs can be found
in the dissertation of the first author [45].

The present paper is an extended and updated version of the proceedings
contribution [34].

##  2 Mathematical preliminaries

This section specifies some mathematical notions which are required for a
rigorous treatment of random sets. Recall that a subset $S$ of a finite
dimensional vector space such as $n$-dimensional Euclidean space
$\mathbb{R}^{d}$ is _open_ , if for each point $x$ in $S$, each ball of the
form $\\{y\in\mathbb{R}^{d}:\|y-x\|<\eta\\}$ with center in $x$ and
sufficiently small radius $\eta$ is also contained in $S$. The set $S$ is
_closed_ , if with each converging sequence whose terms belong to $S$, the
limit also belongs to $S$. The subset $S$ is _bounded_ , if the norms $\|x\|$
of the elements of $S$ are all smaller than a fixed constant. Further, finite
dimensional Euclidean spaces are _separable_ , that is, they contain a
countable dense subset (a sequence of points such that every ball in
$\mathbb{R}^{d}$ contains at least one element of the sequence), and they are
_complete_ (that is, if $x_{k}$ is a sequence such that the sum of its norms
$\sum_{k=1}^{\infty}\|x_{k}\|$ converges, then the sum
$\sum_{k=1}^{\infty}x_{k}$ also converges).

All these notions can be literally generalized to infinite dimensional normed
spaces. A typical example is the space of continuous functions on a bounded,
closed interval $I$, denoted by ${\mathcal{C}}(I)$, with norm
$\|f\|=\max_{x\in I}|f(x)|$. This space is separable and complete.

###  2.1 Random variables

A real-valued random variable $a$ is characterized by its distribution
function $F_{a}(x)$ which defines the probabilities $P$ of basic events of the
form $\\{a\leq x\\}$ via $P(a\leq x)=F_{a}(x)$. This point of view suffices
for most purposes. However, in the sequel _families of random variables_ will
be needed, and it will be necessary to specify the domains of the random
variables under consideration. The usual mathematical setting is a probability
space $(\Omega,\Sigma,P)$ where $\Omega$ is a set (the collection of all
elementary events), $\Sigma$ is the family of measurable subsets of $\Omega$,
and $P$ is a probability measure (that assigns to every set belonging to
$\Sigma$ a number between 0 and 1, satisfying the usual axioms). In this point
of view, a random variable is a mapping from $\Omega$ into the real line such
that the events $\\{\omega\in\Omega:a(\omega)\leq x\\}$ are measurable (i.e.,
belong to $\Sigma$). Thus a probability
$P\big{(}\\{\omega\in\Omega:a(\omega)\leq x\\}\big{)}$ can be assigned to such
events, denoted again by $P(a\leq x)$, and these probabilities in turn define
the distribution function of the random variable $a$. From the probabilities
$P(a\leq x)$, the probabilities of open and closed intervals, their
complements, their countable unions and intersections etc. can be computed.
The family of subsets of the real line which can be obtained by indefinite
repetition of these set theoretic operations defines the so-called _Borel
$\sigma$-algebra_. The Borel $\sigma$-algebra can also be characterized as the
smallest $\sigma$-algebra containing all open subsets of $\mathbb{R}$. A
random variable has the property that all sets of the form
$\\{\omega\in\Omega:a(\omega)\in B\\}$, where $B$ is any Borel set, are
measurable. Accordingly, random variables are often referred to as measurable
functions.

If the probability space is finite, say $\Omega=\\{1,\ldots,m\\}$, one may
take the power set of $\Omega$ as $\sigma$-algebra, and every real-valued
function on $\Omega$ is measurable. A typical infinite probability space is
standard Gaussian space with $\Omega=\mathbb{R}$, $\Sigma$ the Borel
measurable subsets of $\mathbb{R}$, and the probability of an event $B$ given
by the Gaussian integral $P(B)=(2\pi)^{-1/2}\int_{B}{\mathrm{e}}\hskip
1.0pt^{-\omega^{2}/2}\hskip 1.0pt{\mathrm{d}}\omega;$ the probability measure
$P$ is given through the Gaussian density ${\mathrm{e}}\hskip
1.0pt^{-\omega^{2}/2}/\sqrt{2\pi}$. On standard Gaussian space, every
continuous real-valued function on $\Omega$ is a random variable.

These notions can be easily extended to finite dimensional random variables,
i.e., with values in $\mathbb{R}^{d}$ (then the joint probability
distributions have to be specified). Measurable functions with values in an
infinite dimensional normed space can be defined in the same way.

###  2.2 Random fields

A random field is a process that assigns a random variable $q(x)$ to every
point $x$ in a region $D$ in space. To define the probabilistic properties of
the field, the joint distributions of the values at any finite number of
points $q(x_{1}),\ldots q(x_{n})$ should be specified. If the random field is
Gaussian, it is completely specified by the mean value
$\mu_{q}={\mathrm{E}}(q(x))$ and the second moments, i.e., the covariance
$C(x,y)={\mathrm{COV}}(q(x),q(y))$ for any two points $x,y$. If it is (weakly)
homogeneous and isotropic, the covariance depends only on the distance
$\rho=\|x-y\|$ of the points and is of the form $C(x,y)=\sigma^{2}c(\rho)$
with the field variance $\sigma^{2}$ and the autocorrelation function
$c(\rho)$. A typical autocorrelation function is of the form

$c(\rho)=\exp\Big{(}-\frac{|\rho|}{\ell}\Big{)},$ (1)

where $\ell$ is the so-called correlation length. The standard method of
simulation of a random field is based on the _Karhunen-Loève expansion_. If
$q(x)$ is any mean zero random field with finite second moments given by the
covariance $C(x,y)$, its Karhunen-Loève expansion is obtained by solving the
eigenvalue problem

$\int_{D}C(x,y)\varphi_{k}(y)\hskip 1.0pt{\mathrm{d}}y=c_{k}\varphi_{k}(x)$
(2)

which has a sequence of positive eigenvalues $c_{k}$ and orthonormal
eigenfunctions $\varphi_{k}(x)$ (orthonormality in mean square). Then

$q(x)=\sum_{k=1}^{\infty}\sqrt{c_{k}}\xi_{k}\varphi_{k}(x)$ (3)

where the $\xi_{k}$ are uncorrelated random variables with unit variance. If
the process is Gaussian, the $\xi_{k}$ are independent and distributed
according to ${\mathcal{N}}(0,1)$.

For the numerical simulation, the spatial region is discretized by a grid and
the $\varphi_{k}$ are taken, e.g., piecewise constant on the grid elements.
The eigenvalue problem becomes a matrix eigenvalue problem, and the series
(3), with approximate eigenvalues and eigenfunctions, truncated after a finite
number $M$ of terms, can be used for Monte Carlo simulations of the field
trajectories (i.e., realizations of the field).

In case $q(x)$ is a mean zero homogeneous Gaussian random field on
$D=\mathbb{R}$ with autocorrelation function (1) and variance $\sigma^{2}$, an
alternative to the Karhunen-Loève expansion is the representation as an
Ornstein-Uhlenbeck process, namely as solution to the Langevin stochastic
differential equation

$\hskip 1.0pt{\mathrm{d}}q=-\frac{1}{\ell}q\,\hskip
1.0pt{\mathrm{d}}x+\sqrt{\frac{2}{\ell}}\,\sigma\,\hskip
1.0pt{\mathrm{d}}W,\quad q|_{x=0}\sim{\mathcal{N}}(0,\sigma^{2}),$ (4)

see e.g. [5]. Here $W=W(x,\omega)$ denotes Wiener process on the real line. In
this case, realizations of the random field can be conveniently generated by
simulating solutions to the Langevin equation. An additional advantage of the
representation through (4) is that the dependence of $q(x)$ on the field
parameters $\sigma,\ell$ is more explicit than through the Karhunen-Loève
expansion (where $\sigma,\ell$ enter in the coefficients $c_{k}$ and the
eigenfunctions $\varphi_{k}$). The representation (4) allows one to directly
prove that the assignment $(\sigma,\ell)\to q(x,\omega)$ is continuous at
fixed $x$ and $\omega$, an important ingredient in the random set case
discussed in Section 3.

###  2.3 Random sets

In general, a random set is a set-valued random variable satisfying certain
measurability conditions, to be detailed below. The simplest case arises when
the underlying probability space is finite. In this case, one speaks of finite
random sets or Dempster-Shafer structures. Such a structure is given by
finitely many closed subsets $A_{i},i=1,\dots,m$ of a target space
${\mathcal{A}}$, usually Euclidean space $\mathbb{R}^{d}$, called the focal
elements, each of which comes with a probability weight $p_{i}=p(A_{i})$,
$\sum p_{i}=1$.

Viewed as a random set, a Dempster-Shafer structure is given by an $m$-point
probability space $\Omega=\\{1,2,\ldots,m\\}$ with probability masses
$\\{p_{1},p_{2},\ldots,p_{m}\\}$; the assignment $i\to A_{i}$ is the defining
set-valued random variable $A$. However, starting with Subsection 2.4 and
particularly in Section 3, random sets defined on infinite probability spaces
will be needed. Also, the target space can be infinite-dimensional, for
example, a function space containing the system outputs.

In the sequel, $(\Omega,\Sigma,P)$ denotes a general probability space. The
target space will be either Euclidean space $\mathbb{R}^{d}$ or more generally
a normed, complete and separable space ${\mathcal{A}}$. In analogy to the
Borel $\sigma$-algebra on $\mathbb{R}$ introduced in Subsection 2.1, the Borel
$\sigma$-algebra on ${\mathcal{A}}$ is defined as the smallest
$\sigma$-algebra containing all open subsets of ${\mathcal{A}}$. Given a set-
valued map $\omega\to A(\omega)$ where each $A(\omega)$ is a subset of the
target space ${\mathcal{A}}$, upper and lower inverses of subsets $B$ of the
target space ${\mathcal{A}}$ are defined by

$A^{-}(B)=\\{\omega\in\Omega:A(\omega)\cap B\neq\emptyset\\},\qquad
A_{-}(B)=\\{\omega\in\Omega:A(\omega)\subset B\\}.$ (5)

The requirements that $\omega\to A(\omega)$ is a _random set_ are

  * •

each $A(\omega)$ is closed;

  * •

the upper inverses $A^{-}(B)$ are measurable subsets of $\Omega$ for every
closed subset $B$ of the target space ${\mathcal{A}}$.

In the later applications, one would like to admit also open sets $B$ or more
generally Borel sets $B$ in the second condition. In addition, it will be
useful to relate it to the existence of certain measurable selections. Recall
that a _measurable selection_ of a random set $A$ is a random variable $a$
such that $a(\omega)\in A(\omega)$ for all $\omega\in\Omega$. A _Castaing
representation_ is a sequence of measurable selections $a_{k}(\omega)\in
A(\omega)$ such that for all $\omega$, $\\{a_{k}(\omega):k=1,2,3,\ldots\\}$ is
dense in $A(\omega)$.

In this respect, an important technical tool is the _fundamental measurability
theorem_. In order to obtain all subsequent equivalences, one needs to assume
that $(\Omega,\Sigma,P)$ is a _complete probability space_ , meaning that
every subset of a set of probability zero belongs to $\Sigma$. This can always
be achieved by enlarging $\Sigma$ so as to include these subsets.

The fundamental measurability theorem states that, assuming the first
condition (each $A(\omega)$ closed) that the following conditions are
equivalent:

  * •

the upper inverses $A^{-}(B)$ are measurable subsets of $\Omega$ for every
closed subset $B$ of the target space ${\mathcal{A}}$;

  * •

the upper inverses $A^{-}(B)$ are measurable subsets of $\Omega$ for every
open subset $B$ of the target space ${\mathcal{A}}$;

  * •

the upper inverses $A^{-}(B)$ are measurable subsets of $\Omega$ for every
Borel subset $B$ of the target space ${\mathcal{A}}$;

  * •

$\omega\to A(\omega)$ admits a Castaing representation.

For the rather elaborate proof, the reader is referred to [9, 41]. Here are
some first consequences of the fundamental measurability theorem. First, the
lower inverses of events $B$ can be written as
$A_{-}(B)=\\{\omega\in\Omega:A(\omega)\cap B^{\sf c}\neq\emptyset\\}^{\sf c}$,
where $B^{\sf c}$ denotes the complement of $B$, and hence are measurable due
to the fundamental measurability theorem and the measurability of $A^{-}(B)$.
Second, thanks to the required measurability properties, it is legitimate to
introduce upper and lower probabilities

$\displaystyle\overline{P}(B)$ $\displaystyle=$ $\displaystyle
P\big{(}A^{-}(B))\ =\ P(\\{\omega\in\Omega:A(\omega)\cap B\neq\emptyset\\})$
(6) $\displaystyle\underline{P}(B)$ $\displaystyle=$ $\displaystyle
P\big{(}A_{-}(B))\ =\ P(\\{\omega\in\Omega:A(\omega)\subset B\\})$ (7)

for any Borel set $B\subset{\mathcal{A}}$.

The concept of random sets encompasses many other well-established methods of
uncertainty modelling. Clearly, every interval defines a random set (as a
random variable on a one-point probability space). Similarly, every random
variable $a$ on a probability space $(\Omega,\Sigma,P)$ defines a (single-
valued) random set $\omega\to A(\omega)=\\{a(\omega)\\}$. The upper and lower
probability of an event $B$ coincide with its probability $P(B)$. Further,
every normalized fuzzy set can be viewed as a random set. Here the probability
space is the interval $\Omega=[0,1]$, equipped with the uniform probability
distribution, and the focal elements $A(\omega)$ are just the $\omega$-level
sets. It is not difficult to prove that the possibility measure $\pi(B)$ of a
subset $B$ of the real line coincides with its upper probability [25]. Indeed,
and rather obviously,

$\pi(B)=\sup_{\omega\in[0,1]}\\{A(\omega)\cap
B\neq\emptyset\\}=P(\\{\omega\in\Omega:A(\omega)\cap
B\neq\emptyset\\})=\overline{P}(B).$

In the special case of finite random sets, the measurability condition holds
automatically; the only remaining condition is that the focal elements $A_{i}$
are closed.

###  2.4 Random sets generated by families of random variables

To set the stage, consider random variables $a_{\lambda}$, defined on the same
probability space $(\Omega,\Sigma,P)$, and depending on a set of parameters
$\lambda\in\Lambda$. Suppose the random variables have values in a target
space ${\mathcal{A}}$. Then one can define a set-valued map $\omega\to
A(\omega)=\\{a_{\lambda}(\omega):\lambda\in\Lambda\\}.$ The question is
whether this assignment defines a random set. Before answering this question,
consider the example of an imprecise Gaussian family

$a_{\lambda}(\omega)=\mu+\sigma\Phi^{-1}(\omega)$ (8)

with interval parameters
$\lambda=(\mu,\sigma)\in\Lambda=[\underline{\mu},\overline{\mu}]\times[\underline{\sigma},\overline{\sigma}].$
Here the underlying probability space is the unit interval $\Omega=(0,1)$ with
the uniform distribution, and $\Phi$ is the standard normal distribution
function. It follows that each random variable $a_{\lambda}$ is distributed
according to ${\mathcal{N}}(\mu,\sigma^{2})$. The corresponding subset of the
real line is

$A(\omega)=[\underline{a}(\omega),\overline{a}(\omega)]$

with

$\underline{a}(\omega)=\min\big{(}a_{\lambda}(\omega):\lambda\in\Lambda\big{)},\quad\overline{a}(\omega)=\max\big{(}a_{\lambda}(\omega):\lambda\in\Lambda\big{)}.$

Due to the fact that the parameter set $\Lambda$ is a (two-dimensional) closed
and bounded interval and the fact that the assignment
$(\mu,\sigma)\to\mu+\sigma\Phi^{-1}(\omega)$ is continuous at fixed $\omega$,
it follows that each $A(\omega)$ is a closed and bounded interval, indeed.
Some members of the family and the resulting random intervals are depicted in
Figure 1 (left).

Figure 1: Some members of the family $a_{\lambda}$ (left), random set with
focal element (right). The parameters $\lambda=(\mu,\sigma)$ were taken as
$-1\leq\mu\leq 1$, $1\leq\sigma\leq 2$.

Anticipating that $\omega\to A(\omega)$ is a random set, one may define the
upper and lower distribution functions

$\underline{F}(b)=P(\omega:\overline{a}(\omega)\leq
b),\quad\overline{F}(b)=P(\omega:\underline{a}(\omega)\leq b).$ (9)

Anticipating also that the functions $\underline{a}$, $\overline{a}$ bounding
the random intervals are random variables, one obtains simply that
$\underline{F}$ is the distribution function of $\overline{a}$ and
$\overline{F}$ is the distribution function of $\underline{a}$.

Here is the argument why $\omega\to A(\omega)$ is a random set. First, one
constructs a Castaing representation by taking a countable dense subset of
parameter values $\lambda_{k}=(\mu_{k},\sigma_{k})$ in
$\Lambda=[\underline{\mu},\overline{\mu}]\times[\underline{\sigma},\overline{\sigma}]$,
for example, all points with rational coordinates. Again from the continuity
of the assignment $\lambda\to a_{\lambda}(\omega)$ at fixed $\omega$, i.e.,
$(\mu,\sigma)\to\mu+\sigma\Phi^{-1}(\omega)$, see (8), it follows that the set
$\\{a_{\lambda_{k}}(\omega):k=1,2,3,\ldots\\}$ is dense in $A(\omega)$ for
every fixed $\omega$. Let $B$ be an open set. Then

$\begin{array}[]{rcl}A^{-}(B)&=&\\{\omega\in\Omega:A(\omega)\cap
B\neq\emptyset\\}\\\ &=&\\{\omega\in\Omega:\mbox{ there\ is\ }k\mbox{\ such\
that\ }a_{\lambda_{k}}(\omega)\in B\\}\\\
&=&\bigcup_{k}\\{\omega\in\Omega:a_{\lambda_{k}}(\omega)\in B\\}\end{array}$
(10)

is measurable as a countable union of measurable sets; note that the
individual sets are measurable, because the maps $\omega\to
a_{\lambda}(\omega)$ are measurable at each fixed $\lambda$. The fundamental
measurability theorem implies that $A^{-}(B)$ is measurable for every Borel
set $B$, in particular, for every closed set $B$. This is the required
property for $\omega\to A(\omega)$ to be a random set.

Using again that $\\{a_{\lambda_{k}}(\omega):k=1,2,3,\ldots\\}$ is dense in
$A(\omega)$, it follows that
$\underline{a}(\omega)=\inf\\{a_{\lambda_{k}}(\omega):k=1,2,3,\ldots\\}$ and
$\overline{a}(\omega)=\sup\\{a_{\lambda_{k}}(\omega):k=1,2,3,\ldots\\}$ are
measurable as infima and suprema of countably many random variables.

_Comparison with the concept of $p$-boxes._ Given two cumulative distribution
functions $\underline{F}\leq\overline{F}$ on the real line, a _probability
box_ or _$p$ -box_, as introduced in [17], is a collection of cumulative
distribution functions $F$ such that $\underline{F}(b)\leq
F(b)\leq\overline{F}(b)$ for all $b\in\mathbb{R}$. The random set $\omega\to
A(\omega)$ above can be viewed as a $p$-box through (9), but it collects only
the cumulative distribution functions corresponding to the random variables
$a_{\lambda}$, $\lambda\in\Lambda$. Such $p$-boxes are often referred to as
_distributional_ or _parametric $p$-boxes_. The argument above shows that
distributional $p$-boxes are random sets, provided the underlying family of
random variables depends continuously on the parameter $\lambda$. On the other
hand, $p$-boxes that are understood as containing all cumulative distribution
functions between $\underline{F}(b)$ and $\overline{F}(b)$ (the original
concept of [17]) may be referred to as _distribution-free_ or _non-parametric
$p$-boxes_. One may define a corresponding set-valued map by

$\omega\to[\overline{F}^{(-1)}(\omega),\underline{F}^{(-1)}(\omega)]$

employing the quasi-inverses of $\underline{F}^{(-1)}$ and
$\overline{F}^{(-1)}$. It has been shown in [1] that this assignment defines a
random set as well; see also the discussion in [4]. The distributional and the
distribution-free $p$-boxes corresponding to $\underline{F}(b)$ and
$\overline{F}(b)$ in (9) look the same, but the distribution-free version
contains many probability distribution functions other than the one coming
from the Gaussian family (indicated in Figure 1, right).

##  3 Random sets generated by parametrized random fields

This central section lays down the basic structure that will allow one to
prove that set-valued maps obtained as solutions to partial differential
equations whose coefficients are given as intervals, random fields, or
parametrized random fields, are indeed random sets.

###  3.1 The general measurability result

For the present purpose, the most general situation to be considered is the
following.

  * •

$(\Omega,\Sigma,P)$ is a complete probability space;

  * •

the target space ${\mathcal{A}}$ is a normed, complete and separable space;

  * •

$\Lambda$ is a bounded and closed subset of some Euclidean space;

  * •

$a_{\lambda}(\omega)$, $\lambda\in\Lambda$ is a family of random variables
with values in ${\mathcal{A}}$, that is, the maps $\omega\to
a_{\lambda}(\omega)$ are measurable for each $\lambda$;

  * •

the random variables depend continuously on $\lambda$, that is, the maps
$\Lambda\to{\mathcal{A}}:\lambda\to a_{\lambda}(\omega)$ are continuous for
each $\omega$.

Then the following assertions hold: The set-valued map

$\omega\to A(\omega)=\\{a_{\lambda}(\omega):\lambda\in\Lambda\\}$

defines a random set. All focal elements $A(\omega)$ are bounded and closed
subsets of ${\mathcal{A}}$. In addition, if $\Lambda$ is a (multi-dimensional)
interval and ${\mathcal{A}}=\mathbb{R}$, the sets $A(\omega)$ are intervals.

_Proof of assertion._ First, it is clear from the continuity assumption that
all $A(\omega)$ are bounded and closed subsets of ${\mathcal{A}}$. The
assertion that they are intervals under the additional specializing hypothesis
follows from the continuity as well. To prove the random set property, one has
to show that the upper inverses $A^{-}(B)$ are measurable for every open
subset $B$ of ${\mathcal{A}}$. This follows exactly by the argument given in
equation (10) and applying the fundamental measurability theorem.

_Example._ Consider a function $u=u(\lambda,a)$ that depends on an interval
parameter $\lambda$, say $\lambda\in[\underline{\lambda},\overline{\lambda}]$,
and a fixed random variable $a$. Such a situation arises, for example, when
$u$ is the response of a system with mixed interval and random uncertainty.
Both $\lambda$ and $a$ can be multi-dimensional, so this covers the case of
discretized interval and random fields. Assume that the dependence of $u$ on
$\lambda$ and $a$ is continuous. The functions
$a_{\lambda}(\omega)=u(\lambda,a(\omega))$ form a family of random variables
exactly of the type discussed here. Collecting the responses in
$A(\omega)=\\{u(\lambda,a(\omega)):\lambda\in\Lambda\\}$ therefore leads to a
random set, that is actually a random interval
$\omega\to[\underline{a}(\omega),\overline{a}(\omega)]$ if the response is
one-dimensional.

###  3.2 Continuous dependence of random fields on their parameters

A more involved instance of a random set is obtained by taking a random field
on the real line ($x\in\mathbb{R}$) whose parameters are intervals. For
example, consider mean zero random fields $q(x,\omega,\ell)$ with fixed
variance $\sigma^{2}$, but correlation length $\ell$ varying in an interval
$[\underline{\ell},\overline{\ell}]$. A common probability space for these
random fields with different correlation lengths is needed. A convenient way
of generating such random fields, in the case of the autocorrelation function
(1), is by means of the Langevin equation (4). The common probability space is
Wiener space $(\Omega_{W},\Sigma_{W},P_{W})$ (see Appendix).

Due to the fact that the random fields $q(x,\omega,\ell)$ are solutions to
(4), the following properties hold; see e.g. [55, 56]:

  * •

For fixed $x\in\mathbb{R}$ and $\ell\in[\underline{\ell},\overline{\ell}]$,
the map $\omega\to q(x,\omega,\ell)$ is measurable.

  * •

For fixed $x\in\mathbb{R}$ and $\omega\in\Omega_{W}$, the map $\ell\to
q(x,\omega,\ell)$ is continuous.

Thus at every chosen point $x$, the assignment $\omega\to
Q(x,\omega)=\\{q(x,\omega,\ell):\ell\in[\underline{\ell},\overline{\ell}]\\}$
defines a (scalar) random set, as can be seen by the same arguments as in the
example above. The trajectories of the random field are interval-valued
curves, given by $x\to[\min(q(x,\omega,\ell),\max(q(x,\omega,\ell)]$, where
the minimum and maximum are obtained by letting $\ell$ run through
$[\underline{\ell},\overline{\ell}]$.

Alternatively, the Karhunen-Loève expansion (3) can be used to place the
random fields $q(x,\omega,\ell)$ in a common probability space $\Omega$, not
depending on $\ell$. Formally, this can be done by considering an infinite
product of standard Gaussian spaces with random elements
$\omega=(\omega_{1},\omega_{2},\omega_{3},\ldots)$. The required sequence of
Gaussian variables $\xi_{k}$ is then given by $\xi_{k}(\omega)=\omega_{k}$
(see Appendix), and the Karhunen-Loève expansion reads

$q(x,\omega,\ell)=\sum_{k=1}^{\infty}\sqrt{c_{k,\ell}}\xi_{k}(\omega)\varphi_{k,\ell}(x).$
(11)

In general, there is no way to assert that the eigenvalues $c_{k,\ell}$ and
eigenfunctions $\varphi_{k,\ell}(x)$ depend continuously on the correlation
length $\ell$, even after discretization and reduction to a matrix eigenvalue
problem. Thus one has to analyze the situation in each case individually. The
case of the autocovariance function (1) on a finite interval can be treated
explicitly, thanks to the results in [24]. Without restriction of generality,
one may assume that the finite interval is $D=[-1,1]$. Then the eigenvalue
problem (2) becomes

$\int_{-1}^{1}{{\mathrm{e}}\hskip 1.0pt^{-|x-y|/\ell}\varphi(y)}\,\hskip
1.0pt{\mathrm{d}}y=c\varphi(x).$ (12)

Following [24], problem (12) can be solved in the following way. First, the
two equations

$\frac{1}{\ell}-\alpha\tan{\alpha}=0,\quad{\rm
respectively}\quad\alpha+\frac{1}{\ell}\tan{\alpha}=0$ (13)

have two sequences of positive solutions denoted by
$\alpha_{k}=\alpha_{k}(\ell)$ and $\alpha^{\ast}_{k}=\alpha^{\ast}_{k}(\ell)$,
for even and odd $k$, respectively. The resulting eigenvalues and
eigenfunctions are

$c_{k,\ell}=\frac{2\ell}{1+\ell^{2}\alpha_{k}^{2}(\ell)},\quad
c_{k,\ell}^{\ast}=\frac{2\ell}{1+\ell^{2}\alpha^{\ast 2}_{k}(\ell)}$ (14)
$\varphi_{k,\ell}(x)=\frac{\cos{\alpha_{k}(\ell)x}}{\sqrt{1+\frac{\sin{2\alpha_{k}(\ell)}}{2\alpha_{k}(\ell)}}},\quad\varphi^{\ast}_{k,\ell}(x)=\frac{\sin{\alpha^{\ast}_{k}(\ell)x}}{\sqrt{1-\frac{\sin{2\alpha^{\ast}_{k}(\ell)}}{2\alpha^{\ast}_{k}(\ell)}}}$
(15)

for even and odd $k$, respectively. Therefore, the one-dimensional random
field $q(x,\omega,\ell)$ with covariance function given by the equation (1)
has the Karhunen-Loève expansion

$q(x,\omega,\ell)=\lim_{M\to\infty}\sum\limits_{k=1}^{M}\Big{(}\sqrt{c_{k,\ell}}\xi_{k}(\omega)\varphi_{k,\ell}(x)+\sqrt{c_{k,\ell}^{\ast}}\xi^{\ast}_{k}(\omega)\varphi^{\ast}_{k,\ell}(x)\Big{)}.$
(16)

For the truncated Karhunen-Loève expansion, one may take a $2M$-dimensional
standard Gaussian space as probability space $\Omega$, on which the random
variables $\xi_{k}$ and $\xi^{\ast}_{k}$ are defined. It is clear that the
solutions $\alpha_{k}$ and $\alpha_{k}^{\ast}$ of the deterministic equations
(13) depend continuously on $\ell$. The same is true for the eigenvalues as
well as eigenfunctions shown in (14), (15). Thus truncating the Karhunen-Loève
expansion (16) at a finite $M$, the resulting random field depends
continuously on the correlation length $\ell$ at fixed $\omega$.

###  3.3 Random fields as random variables valued in function spaces

Consider the random field $q(x,\omega,\ell)$ as in Subsection 3.2. The _paths
or trajectories (in space)_ of the random field are the maps $x\to
q(x,\omega,\ell)$, with fixed $\omega$ and $\ell$. It is well-known from the
theory of stochastic differential equations that the paths are continuous [5].
This also follows from the Kolmogorov-Chentsov theorem [33], because the
autocovariance function (1) satisfies the hypotheses of this theorem. However,
for the purposes of the present paper, more is needed.

First, jointly continuous dependence on $x$ and $\ell$ is required. This is
easy to achieve, because one can consider the _paths in space and parameter
set_ as the maps $(x,\ell)\to q(x,\omega,\ell)$, with fixed $\omega$ and apply
the previous arguments. Indeed, when using a truncated Karhunen-Loève
expansion, the argument was just given after (16). When the Langevin equation
is used, one may resort to the results of [55, 56] which assert the continuous
dependence, jointly in $x$ and $\ell$.

Second, the random fields have to be interpreted as random variables in the
space of continuous functions. This is because they enter as variable
coefficients in the partial differential equations to be solved, and the
solution depends on the whole trajectories, not only the values of the random
fields at individual points $x$. Let $x$ belong to an interval $I$ and $\ell$
to an interval $L$, both closed and bounded. Continuity of the paths in $x$
and $\ell$ can be stated by saying that the map $q(\omega):(x,\ell)\to
q(x,\omega,\ell)$ belongs to the space of continuous function
${\mathcal{C}}(I\times L)$ for fixed $\omega$. Thus the random field induces a
map $\Omega\to{\mathcal{C}}(I\times L)$, $\omega\to q(\omega)$. Standard
arguments from the theory of stochastic processes show that this map is
measurable, that is, a random variable with values in ${\mathcal{C}}(I\times
L)$. For details of the proof, see [45, Section 1.2]. The argument also shows
that

  * •

for fixed $\ell\in L$, the map $\Omega\to{\mathcal{C}}(I)$, $\omega\to
q(\omega,\ell)$ is measurable;

  * •

for fixed $\omega\in\Omega$, the map $L\to{\mathcal{C}}(I)$, $\ell\to
q(\omega,\ell)$ is continuous.

Note that the same random field $q$ may be considered in three different ways
as a measurable function of $\omega$: as an element of the function space
${\mathcal{C}}(I\times L)$, denoted by $q(\omega)$, as an element of the
function space ${\mathcal{C}}(I)$ at fixed correlation length $\ell$, denoted
by $q(\omega,\ell)$, and as a real number at each fixed $x$ and $\ell$,
denoted by $q(x,\omega,\ell)$. Whenever its values are collected in a random
set, the corresponding capital letter $Q$ will be used.

###  3.4 The mean values of a parametrized random field as an interval field

Let $\omega\to A(\omega)$ be a random set whose focal elements
$A(\omega)=[\underline{a}(\omega),\overline{a}(\omega)]$ are closed and
bounded intervals. Let $\\{a_{k}(\omega)$, k = 1,2,3,…} be a Castaing
representation. The _Aumann expectation_ ${\mathrm{E}}(A)$ is defined as the
closure of the set of all expectation values of the measurable selections
$a_{k}$ for $k=1,2,3,\ldots$, that is, ${\mathrm{E}}(A)={\rm
cl}\\{{\mathrm{E}}(a_{k}):k=1,2,3,\ldots\\}$, see [41]. It is clear that in
the considered case, the Aumann expectation is simply the interval
${\mathrm{E}}(A)=[{\mathrm{E}}(\underline{a}),{\mathrm{E}}(\overline{a})].$
Next, consider a random set generated by a parametrized random field
$\omega\to
Q(x,\omega)=\\{q(x,\omega,\ell):\ell\in[\underline{\ell},\overline{\ell}]\\}$
as in Subsection 3.2. Thanks to the continuous dependence on the parameter
$\ell$ discussed in Subsection 3.2 each $Q(x,\omega)$ is a closed and bounded
interval at fixed $x$. One may form the Aumann expectation
${\mathrm{E}}(Q(x))$, which is a deterministic interval at fixed $x$. Thus the
assignment $x\to{\mathrm{E}}(Q(x))$ defines an interval field.

##  4 Application to elliptic equations

A prototypical stochastic elliptic partial differential equation, arising e.g.
in elastostatics, is the boundary value problem

$\begin{array}[]{rclrcl}-\,{\rm div\,}\big{(}a(x)\,{\rm
grad\,}u(x)\big{)}&=&f(x),&x&{\rm in}&D\\\ u(x)&=&0,&x&{\rm in}&\partial D\\\
\end{array}$ (17)

as discussed, e.g., in [28]. Here $D$ is a bounded open domain in
$\mathbb{R}^{d}$ and $\partial D$ is its boundary. For example, in two space
dimensions $u$ might be the transversal displacement of a non-uniform membrane
under transversal load $f$. The spatially varying elastic properties are
subsumed in the coefficient $a(x)$. Alternatively, in three space dimensions,
$u(x)$ could be the pressure of a fluid in a porous medium with permeability
$a(x)$ and source rate $f(x)$ [39].

This section serves to demonstrate how the random field methods can be applied
in practice. The required solution theory is collected, and a simple numerical
example is elaborated.

In the foundations of the finite element method, problem (17) is commonly
solved by variational methods. Suitable function spaces are defined as
follows: $L^{2}(D)$ is the space of square integrable functions; $H^{1}(D)$ is
the space of square integrable functions all whose first order partial
derivatives are square integrable as well; $H^{1}_{0}(D)$ is the subspace of
functions that vanish on the boundary $\partial D$. Both $L^{2}(D)$ and
$H^{1}_{0}(D)$ are normed spaces; the squares of the norms are given by

$\|f\|_{L^{2}(D)}^{2}=\int_{D}|f(x)|^{2}\hskip
1.0pt{\mathrm{d}}x,\quad\|u\|_{H^{1}(D)}^{2}=\|u\|_{L^{2}(D)}^{2}+\|{\rm
grad\,u}\|_{L^{2}(D)}^{2}.$

The spaces are complete and separable. The variational formulation of equation
(17) then is:

> Find $u\in H^{1}_{0}(D)$ such that
>
> $\int_{D}a(x)\,{\rm grad\,}u(x)\cdot{\rm grad\,}v(x)\hskip
> 1.0pt{\mathrm{d}}x+\int_{D}f(x)v(x)\hskip 1.0pt{\mathrm{d}}x=0$ (18)
>
> for all $v\in H^{1}_{0}(D)$.

_The deterministic case._ Assume that $a$ is a continuous function on $D$
which is bounded from above and below, that is, $a\in{\mathcal{C}}(D)$ and
$\alpha\leq a(x)\leq\beta$ for some constants $\alpha,\beta>0$ and all $x$ in
$D$. Further, assume that $f$ is square integrable on $D$. Under these
conditions, it is well-known that the variational problem (18) has a unique
solution $u$ in $H^{1}_{0}(D)$, [44, 70]. For the purpose of this example, the
source term $f$ will be fixed. The solution $u$ depends on the coefficient
function $a$; this dependence will be made explicit in the notation $u=u_{a}$.
Denote the set of functions which are continuous on $D$ up to the boundary
with values between $\alpha$ and $\beta$ by ${\mathcal{C}}(D;\alpha,\beta)$.
It is equipped with the norm $\|a\|=\max_{x\in D}|a(x)|$. The crucial
assertion is that the map

${\mathcal{C}}(D;\alpha,\beta)\to H^{1}_{0}(D):a\to u_{a}$

is continuous. This non-trivial result can be proved by invoking the
variational formulation (18), as it was done in [45]. A different proof, for
more general partial differential equations, can be found in [44].

_The random field case._ To start with, let $a(x,\omega)$, $\omega\in\Omega$,
be a random field with continuous paths lying between $\alpha$ and $\beta$. As
noted in Subsection 3.3, it can be viewed as a random variable with values in
${\mathcal{C}}(D;\alpha,\beta)$. Due to the continuity of the above map, the
assignment $\omega\to u_{a(\omega)}$ is measurable, thus the solution $u_{a}$
is a random variable with values in $H^{1}_{0}(D)$. Intuitively, this means
that $u_{a}$ is a random field whose paths belong to $H^{1}_{0}(D)$. In order
to interpret it as a classical random field, one should be able to assign
values $u_{a(\omega)}(x)$ at each point $x$ in $D$. In one space dimension,
this is easy because $H^{1}_{0}(D)$ is continuously imbedded in
${\mathcal{C}}(D)$, and so the pointwise evaluations

${\mathcal{C}}(D;\alpha,\beta)\to\mathbb{R}:a\to u_{a}(x)$

are continuous as well. Thus the maps $\omega\to u_{a(\omega)}(x)$ are random
variables for each $x$, and the paths $x\to u_{a(\omega)}(x)$ are continuous
functions. In space dimensions two and three, more theory is needed. First,
the input random field $a(x,\omega)$ should have Lipschitz continuous paths,
that is, satisfying an inequality of the form $|a(x,\omega)-a(y,\omega)|\leq
C|x-y|$ for some constant $C$ and all $x,y$. This is not true of the random
field generated by the Langevin equation, but it is true of the random field
generated by the truncated Karhunen-Loève expansion. Then one can invoke
regularity theory saying that the solution to problem (17) is in
$H^{2}(D^{\prime})$ in any open subregion $D^{\prime}$ of $D$ and depends
continuously on $a$ in this space [44]. Again, $H^{2}(D^{\prime})$ is
continuously imbedded in ${\mathcal{C}}(D^{\prime})$, and so one obtains a
classical random field $u_{a(\omega)}(x)$, defined for all $x$ in the interior
of $D$.

_The random set case._ Let $a(x)$ be a parametrized random field. More
specifically, it will be constructed as follows. In two space dimensions, take
a closed and bounded interval $I$ such that $D\subset I\times I$ and a random
field of the form $q(x,\omega,\ell)$ as in Subsection 3.2 and set

$a(x,\omega,\ell)=\mu(x)+\chi\big{(}q(x_{1},\omega,\ell)q(x_{2},\omega,\ell)\big{)},\quad{\rm
for}\quad x=(x_{1},x_{2}){\ \rm in\ }D.$ (19)

Here $\mu(x)>0$ is the mean field (assumed fixed and deterministic) and $\chi$
is a cut-off function; this is needed because the individual realizations of a
Gaussian random field can be arbitrarily large with non-zero probability. In
the end, the cut-off has to be chosen so as to guarantee that $\alpha\leq
a(x,\omega,\ell)\leq\beta$.

One could also use different correlation lengths in the two variables and/or
multiply the fields with different standard deviations
$\sigma_{1},\sigma_{2}$. The construction in one or three space dimensions is
analogous.

When $\ell$ varies in an interval $L$, a set-valued solution to (17) can be
defined as $U(\omega)=\\{u_{a(\omega,\ell)}:\ell\in L\\}$. Combining the
continuity of the map $a\to u_{a}$ with the statement at the end of Subsection
3.3, one obtains that

  * •

for fixed $\ell\in L$, the map $\Omega\to H^{1}_{0}(D)$, $\omega\to
u_{a(\omega,\ell)}$ is measurable;

  * •

for fixed $\omega\in\Omega$, the map $L\to H^{1}_{0}(D)$, $\ell\to
u_{a(\omega,\ell)}$ is continuous.

It follows from Subsection 3.1 that $U(\omega)$ is a random set in the target
space $H^{1}_{0}(D)$. Again, in one space dimension one may form the set-
valued functions $U(\omega,x)=\\{u_{a(\omega,\ell)}(x):\ell\in L\\}$ at each
$x$ in $D$. By the same argument, the $U(\omega,x)$ are random sets in
$\mathbb{R}$; they form the point evaluations of the random set solution
$U(\omega)$ at the points $x$ in $D$. In space dimensions two and three, a
modified argument as in the random field case has to be applied.

_Numerical example._ The purpose of this paragraph is to illustrate the method
in a simple example, the displacement of an L-shaped membrane with non-
uniform, random elastic properties. This has the character of a text book
example, and so all units are taken dimensionless. The L-shaped domain $D$ can
be seen in Figure 2 (left). The equation for the displacement is (17). A
constant load $f(x)\equiv 1$ is assumed, and the random field $a(x)$ is taken
of the form (19) with constant mean $\mu(x)\equiv 1$, autocovariance function
of the form (1), and field variance identically equal to one as well. In this
elliptic problem, it suffices to guarantee that the realizations of the random
field $a(x,\omega,\ell)$ stay above zero. Thus a lower cut-off function $\chi$
was chosen at the value 0.1. To produce the random set solution, the
correlation length was taken from the interval $0.5\leq\ell\leq 1.5$.

Numerically, the L-shaped region was subdivided in a grid of size $18\times
18$, and problem (17) was solved by the finite element method using a Matlab-
program. To generate the random field, the Karhunen-Loève expansion truncated
at $M=130$ was used. The value of $a(x,\omega,\ell)$ entered in each finite
element was chosen as the average of the values at its four nodes. A Monte
Carlo sample size of 500 was used. At each realization of the underlying
$2M$-dimensional Gaussian variable, see (16), the finite element program was
executed with eleven values of $\ell$ between 0.5 and 1.5, stepsize 0.1.

As outlined above, this procedure results in an approximation to the random
set solution $U(\omega)$. The solution is a random set in $H_{0}^{1}(D)$ and
in $H^{2}(D^{\prime})$ in any open subregion, actually everywhere except near
the sharp inside corner [29]. Consequently, each horizontal or vertical slice
through the solution is a random set in the space of continuous functions. At
fixed $x=(x_{1},x_{2})$, each realization of the solution is an interval. The
mean (Aumann expectation) is an interval field. Figure 2 (left) shows a
realization of the displacement random field at fixed correlation length
$\ell=0.5$. Realizations of the full random set solution can be visualized by
plotting slices at fixed $x_{1}$ or $x_{2}$. Such a cut at $x_{2}=0.4444$ is
shown in Figure 2 (center). Above each fixed $x_{1}$, the realization is an
interval. Finally, the mean value field can also be depicted by slices. Figure
2 (right) shows a slice of this interval field at $x_{2}=0.4444$. The upper
and lower bounding curves of the mean displacement are shown together with the
mean solutions at correlation lengths 0.5, 0.6, …, 1.5.

Figure 2: Some aspects of the random set solution. Trajectory of displacement
random field at fixed correlation length $\ell=0.5$ (left); slice at
$x_{2}=0.4444$ through a realization of the full random set solution (center);
slice at $x_{2}=0.4444$ through the mean value interval field, showing
individual mean fields for varying correlation lengths inside.

It is clear that the strategy outlined in the simple example (17) can be
extended to more general elliptic equations and systems, using the variational
formulation (see [44, 70] for second order scalar elliptic equations and [12]
for more general elliptic systems of elasticity theory).

##  5 Application to hyperbolic systems

Many equations from one-dimensional elastodynamics (and other fields, such as
acoustics or hydrodynamics) can be written as a hyperbolic system

$\begin{array}[]{rcl}\left(\partial_{t}+a_{i}(x,t)\partial_{x}\right)u_{i}(x,t)&=&\sum_{j=1}^{n}f_{ij}(x,t)u_{j}(x,t)+g_{i}(x,t),\\\
u_{i}(x,0)&=&u_{0i}(x),\quad{\rm for}\ i=1,\ldots,n,\end{array}$ (20)

where $x\in\mathbb{R}$ and $t\in\mathbb{R}$ denote the spatial and time
variable, respectively. The solution vector is $(u_{1},\ldots,u_{n})$, the
driving term is $(g_{1},\ldots,g_{n})$, the coefficient matrix is
$(f_{ij})_{i,j=1,\ldots,n}$. The coefficients $a_{i}$ are assumed to be real-
valued and signify the propagation speeds of the solution components. The
focus of this section will be on the case where the coefficients as well as
the initial data are random fields/random sets on a probability space
$(\Omega,\Sigma,P)$.

A typical example from elastodynamics concerns longitudinal waves in thin,
non-homogeneous rods, described by the wave equation

$\rho(x)\partial_{t}^{2}u(x,t)-\partial_{x}(E(x)\partial_{x}u(x,t))=q(x,t)$
(21)

where $u$ is the longitudinal displacement, $\rho$ the density, $E$ the
modulus of elasticity, and $q$ the body force [26, Chapter 2]. Equation (21)
is readily cast in the form of (20)

$\begin{array}[]{rcl}\left(\partial_{t}+a(x)\partial_{x}\right)u_{1}(x,t)&=&f(x,t)(u_{2}(x,t)-u_{1}(x,t))+g(x,t),\\\
\left(\partial_{t}-a(x)\partial_{x}\right)u_{2}(x,t)&=&f(x,t)(u_{2}(x,t)-u_{1}(x,t))+g(x,t)\end{array}$
(22)

by setting $u_{1}(x,t)=(\partial_{t}-a(x)\partial_{x})u(x,t)$,
$u_{2}(x,t)=(\partial_{t}+a(x)\partial_{x})u(x,t)$ with $a=\sqrt{E/\rho}$,
$f=\partial_{x}E/\rho-\partial_{x}a/2$ and $g=q/\rho$.

A second example is provided by various scalar transport or advection-reaction
equations in the form

$\begin{array}[]{rcl}\partial_{t}u(x,t)+a(x,t)\partial_{x}u(x,t)&=&f(x,t)u(x,t)+g(x,t),\\\
u(x,0)&=&u_{0}(x).\end{array}$ (23)

The construction of classical, deterministic solutions to the hyperbolic
system (20) is a vector-valued version of the construction for the scalar
equation (23). Thus it will suffice to briefly indicate how (23) can be
treated. The systems case can be found, e.g., in [47, Section 13]. The
solution can be constructed by the method of characteristics, which consists
in solving the system of differential equations

$\begin{array}[]{rcl}\frac{\partial}{\partial\tau}\gamma(x,t,\tau)&=&a\left(\gamma(x,t,\tau),\tau\right),\\\
\gamma(x,t,t)&=&x.\end{array}$ (24)

This results in the characteristic curves $\tau\to\gamma(x,t,\tau)$, passing
through the point $x$ at time $\tau=t$. Equation (23) is equivalent to the
integral equation

$u(x,t)=u_{0}(\gamma(x,t,0))+\int_{0}^{t}\big{(}f(\gamma(x,t,\tau),\tau)u(\gamma(x,t,\tau),\tau)+g(\gamma(x,t,\tau),\tau)\big{)}\,d\tau,$
(25)

which in turn can be solved by successive iteration.

The issue to be addressed here is when $a$ and $f$, $g$ (and possibly $u_{0}$)
are random fields or parametrized random fields on the probability space
$(\Omega,\Sigma,P)$. In order to obtain the solution as a random field or as a
random set, additional properties of the coefficients are required.

_The random field case._ The starting observation is that equation (24) turns
into the random differential equation

$\frac{\partial}{\partial\tau}\gamma(x,t,\tau,\omega)=a\left(\gamma(x,t,\tau,\omega),\tau,\omega\right),\quad\gamma(x,t,t,\omega)=x.$

Following the theory of random differential equations [8, 30], a continuous
random field solution $\gamma(x,t,\tau,\omega)$ can be obtained provided the
transport coefficient $a(x,t,\omega)$ is pathwise Lipschitz continuous and
globally bounded. The next issue is the composition of the random fields $f$,
$g$, $u_{0}$ with $\gamma$, for example, $(x,t,\tau,\omega)\to
f(\gamma(x,t,\tau,\omega),\tau,\omega)$. This composition is jointly
measurable (in all four variables) and pathwise continuous, if both $\gamma$
and $f$ are pathwise continuous random fields. Finally, the integral equation
(25), now with random field input, has to be solved by successive
approximation (on suitable bounded domains of determinacy), taking account of
the measurability at each iteration step. The precise requirements are:

  * •

For fixed $\omega\in\Omega$, $a(x,t,\omega)$ is locally Lipschitz in $(x,t)$;

  * •

$|a(x,t,\omega)|$ has a uniform global bound independently of $(x,t,\omega)$;

  * •

$|f(x,t,\omega)|$ is locally bounded in $(x,t)$ independently of $\omega$;

  * •

$f$, $g$ and $u_{0}$ are pathwise continuous random fields.

Under these assumptions, the integral equation (25) has a pathwise continuous
random field solution $u(x,t,\omega)$. Note that the assumptions on $a$ and
$f$ can be enforced by a cut-off procedure similar to (19). If, in addition,
the random fields $a$, $f$, $g$, $u_{0}$ are pathwise differentiable,
$u(x,t,\omega)$ is a pathwise classical solution to the differential equation
(23). When solving, e.g., the wave equation (21), one has to take care that
the fields $\rho$ and $E$ are sufficiently regular so that the coefficients in
system (22) have the required properties.

Another approach to hyperbolic systems with random field coefficients, based
on weak solutions, can be found in [40]. It applies in any space dimension and
leads to a different set of requirements on the coefficients and data.

_The random set case._ Assume that some or all of the random fields $a$, $f$,
$g$, $u_{0}$ depend on interval parameters. As in the previous sections, this
leads to a set-valued solution $\omega\to U(\omega)$ to (23) in a suitable
space to be determined. In order to prove that the set-valued solution is a
random set, one may follow exactly the same strategy as in Section 4. Only the
continuous dependence of the solution on the coefficients and data requires
different arguments, as outlined next. The continuity of the map
$(a,f,g,u_{0})\to u=u_{(a,f,g,u_{0})}$ can be split up into the continuity of
the maps

$a\to\gamma\quad{\rm and}\quad(\gamma,f,g,u_{0})\to u.$

To be specific, start with a bounded and closed initial interval
$K_{0}=[-\kappa,\kappa]$ and a time horizon $T>0$. Denote by $c$ the uniform
global bound of $|a(x,t,\omega)|$, as listed under the requirements. The
region

$K_{T}=\\{(x,t)\in\mathbb{R}^{2}:-T\leq t\leq T,|x|\leq\kappa-c|t|\\}$

is a domain of determinacy for $K_{0}$. That is, when $(x,t)\in K_{T}$, all
characteristic curves $\tau\to\gamma(x,t,\tau,\omega)$ joining $(x,t)$ with
the $x$-axis lie also in $K_{T}$. Denote by ${\rm Lip}(K_{T})$ the space of
Lipschitz continuous functions on $K_{T}$ and by $\Gamma(K_{T}\times[-T,T])$
the subset of ${\mathcal{C}}(K_{T}\times[-T,T])$ which comprises all
characteristic curves $\gamma$ arising from elements $a\in{\rm Lip}(K_{T})$.
The first step is to show that the map

${\rm Lip}(K_{T})\to\Gamma(K_{T}\times[-T,T]):a\to\gamma$

is continuous, when the spaces are equipped with the supremum norm. The second
step is to prove that the map

$\Gamma(K_{T}\times[-T,T])\to{\mathcal{C}}(K_{T}):(\gamma,f,g,u_{0})\to u$

is continuous in the supremum norm. The latter is a rather standard argument
in the theory of differential equations using the Gronwall lemma. Details are
elaborated in [45]. Since the initial choice of $K_{0}$ was arbitrary, the
continuity of the map $(a,f,g,u_{0})\to u=u_{(a,f,g,u_{0})}$ with values in
${\mathcal{C}}(\mathbb{R}^{2})$ follows. Thus if the random fields
$(a,f,g,u_{0})$ depend continuously on a (possibly multidimensional) parameter
$\lambda\in\Lambda$, the same arguments as in Sections 3 and 4 show that

$U(\omega)=\\{u_{(a(\omega,\lambda),f(\omega,\lambda),g(\omega,\lambda),u_{0}(\omega,\lambda))}:\lambda\in\Lambda\\}$

is a random set in ${\mathcal{C}}(\mathbb{R}^{2})$.

_The multidimensional case:_ The generalization to the transport equation in
higher space dimensions,

$\partial_{t}u(x,t)+a(x,t)\cdot{\rm grad\,}u(x,t)=f(x,t)u(x,t)+g(x,t)$

is straightforward, using the method of characteristics. Here
$x\in\mathbb{R}^{d}$, $a=(a_{1},\ldots,a_{d})$, $u$, $f$ and $g$ are still
scalar functions. In the multidimensional analogues of (20), (21), it is more
challenging to establish continuous dependence of the solution on the data.
One approach could be to formulate the problem in terms of semigroup theory
and then invoke Kato’s perturbation results [35]. For the system of linear
acoustics this has been done, e.g., in [23, Section 6]. Another approach takes
recourse to the Fourier integral operator representation of the solution [52].

##  6 Further results, numerical aspects

###  6.1 The parametric point of view

Given a parametrized family $a_{\lambda},\lambda\in\Lambda$, of random
variables or random fields, and a map $a_{\lambda}\to u_{a_{\lambda}}$
describing the response of a system, the approach favored in this paper has
been to first construct the random set

$U(\omega)=\\{u_{a_{\lambda}(\omega)}:\lambda\in\Lambda\\}$ (26)

and then to define lower and upper probabilities of events $B$ by

$\underline{P}(B)=P(U_{-}(B)),\quad\overline{P}(B)=P(U^{-}(B)).$

An alternative approach, the _parametric viewpoint_ , would consist in
skipping the random set step and defining lower and upper probabilities
induced by the family $a_{\lambda},\lambda\in\Lambda$, directly:

$P_{\rm low}(B)=\inf_{\lambda\in\Lambda}P(u_{a_{\lambda}}(B)),\quad P_{\rm
upp}(B)=\sup_{\lambda\in\Lambda}P(u_{a_{\lambda}}(B)).$

This approach gives tighter bounds:

$\underline{P}(B)\leq P_{\rm low}(B)\leq P_{\rm upp}(B)\leq\overline{P}(B),$

but deprives one of using the concepts of the theory of random sets. The
numerical computation of the two types of upper and lower probabilities are
different, but of comparable effort. For a comparison of the two approaches,
from the conceptual and numerical points of view, see [3, 21].

###  6.2 Numerical aspects

As is well-known, the computational effort when implementing heterogeneous
uncertainty models is high: it involves both Monte Carlo sampling and the
computation of interval bounds. Suppose the quantity of interest, given
through a model response function $a_{\lambda}\to u_{a_{\lambda}}$ as above,
is one-dimensional. This could be the model response at a single point $x$ or
the value of a performance criterion. Further, suppose the lower and upper
distribution functions $\underline{F}(b)$, $\overline{F}(b)$ of this quantity
are sought.

The procedure in the random set framework is as follows: If the hypotheses of
the previous sections are met, the quantity of interest (26) is a random
interval $U(\omega)=[\underline{u}(\omega),\overline{u}(\omega)]$, and the
lower and upper distribution functions are simply given by
$\underline{F}(b)=P(\overline{u}\leq b)$, $\overline{F}(b)=P(\underline{u}\leq
b)$. The following algorithm produces a numerical approximation to these lower
and upper distribution functions. From there, further probabilistic indicators
such as bounds on exceedence probabilities or moment bounds can be computed.
It is assumed that the parameter set $\Lambda$ indexing the family of random
variables $a_{\lambda}$ is a $d$-dimensional interval (or a sufficiently
regular subset of such an interval). It is also assumed that the random
variables $a_{\lambda}$ are defined on a common probability space
$(\Omega,\Sigma,P)$. The following steps are required:

  * •

Discretize $\Lambda$ in a grid $\lambda_{1},\ldots,\lambda_{M}$.

  * •

Generate a sample $\omega_{1},\ldots,\omega_{N}$ distributed according to the
underlying probability $P$ on the common probability space $\Omega$.

  * •

For each $\omega_{k}$, estimate

$\underline{u}(\omega_{k})=\min_{i=1,\ldots,M}u_{a_{\lambda_{i}}}(\omega_{k}),\quad\overline{u}(\omega_{k})=\max_{i=1,\ldots,M}u_{a_{\lambda_{i}}}(\omega_{k}).$

Each $[\underline{u}(\omega_{k}),\overline{u}(\omega_{k})]$ is a realization
of the random interval $[\underline{u},\overline{u}]$.

  * •

Then $\underline{u}(\omega_{k})$ and $\overline{u}(\omega_{k})$,
$k=1,\ldots,N,$ are samples of $\underline{u}$ and $\overline{u}$,
respectively. Their empirical cumulative distribution functions provide
estimates for $\underline{F}(b)$ and $\overline{F}(b)$.

In the parametric point of view, the assumptions on the index set $\Lambda$
are the same. It is not necessary to require that the random variables
$a_{\lambda}$ are defined on the same probability space. Denote by
$P_{\lambda}$ the probability distribution of $a_{\lambda}$. The algorithm is
as follows:

  * •

Discretize $\Lambda$ in a grid $\lambda_{1},\ldots,\lambda_{M}$.

  * •

For each $\lambda_{i}$, generate a sample
${a_{\lambda_{i}}}(1),\ldots,{a_{\lambda_{i}}}(N)$ of the random variable
${a_{\lambda_{i}}}$ (that is, generate a sample of size $N$ from the
distribution $P_{\lambda_{i}}$).

  * •

Estimate the probability $P_{\lambda_{i}}((-\infty,b])$ using the sample, for
example by means of the Monte Carlo estimator
$\frac{1}{N}\sum_{j=1}^{N}\chi_{(-\infty,b]}(a_{\lambda_{i}}(j))$ where
$\chi_{(-\infty,b]}$ denotes the characteristic function of the interval
$(-\infty,b]$.

  * •

Estimate

$F_{\rm low}(b)\approx\min_{i=1,\ldots,M}P_{\lambda_{i}}((-\infty,b]),\quad
F_{\rm upp}(b)\approx\max_{i=1,\ldots,M}P_{\lambda_{i}}((-\infty,b]).$

In these crude forms, both algorithms require $M\cdot N$ evaluations of the
(generally expensive) model function $u$. Accordingly, many approaches have
been devised to reduce computational cost. In the random set case, a
stochastic response surface, like a polynomial chaos expansion [24], reduces
the number $M$ to a level required for the accuracy of the response surface.
Evaluating a Monte Carlo sample of size $N$, given the response surface, is
cheap [50]. For further sampling concepts for random sets, see [1].

Here is a short list of methods currently being developed for the reduction of
computational cost in imprecise probability models, in particular, for
imprecise random fields: Approximations by _inner and outer bounds_ [2, 3, 4,
68, 77]; propagation methods for _$p$ -boxes_ [13, 14, 27, 54, 75, 76];
methods employing _polynomial chaos expansion_ [37, 49, 61]; _probability
plots_ , an enhancement of the first order reliability method (FORM) proposed
by [31]. A very efficient method appears to be _advanced line sampling_ [10,
66] that intertwines the two required loops and can reduce $M$ and $N$
simultaneously. In the parametric point of view, a promising method has become
_importance sampling_ , or rather _reweighting_ a single importance Monte
Carlo sample, which means $M=1$. Still $N$ model evaluations are required [18,
19, 71, 72, 73, 78]. This can of course be aided by a deterministic response
surface.

###  6.3 Modelling correlations

The question of how to model correlations in hybrid probability-interval
uncertainty analysis has occurred implicitly at various places in this
article. The purpose of this subsection is to spell out more explicitly what
kind of approaches are possible. Essentially, there are three ways to
formulate correlations.

_Correlations introduced by the underlying probabilistic model._ In the
approach to random sets based on families of random variables, correlations
can be defined through their joint distributions. These correlations are
carried along when replacing distributional hyperparameters by intervals. This
is generally the case in random field models with interval correlation lengths
as e.g in [50] as well as in Section 4. More generally, random fields with
families of correlation structures have been considered in [16], especially
analyzing the corresponding Karhunen-Loève expansions with interval valued
eigenvalues and eigenfunctions. In models with interval parameters and random
field excitation, correlations are introduced through the random field, but
become interval valued in output and/or state variables. Such a situation has
been investigated in connection with spectral finite elements (finite elements
combined with polynomial chaos expansions) in [11].

_Copulas and multivariate capacities._ A very interesting and practical
approach is due to [1, 2], who introduced the notion of _random sets of
indexable type_. These random sets are defined on $\Omega=(0,1]^{d}$ and have
values in a target product space
${\mathcal{A}}={\mathcal{A}}_{1}\times\cdots\times{\mathcal{A}}_{d}$. The
focal elements are of product form $A_{1}(\omega_{1})\times\cdots\times
A_{d}(\omega_{d})$ where $\omega=(\omega_{1},\ldots,\omega_{d})\in\Omega$ and
each $A_{j}(\omega_{j})$ belongs to a suitable family of subsets of
${\mathcal{A}}_{j}$ (closed subsets for the purpose of the present paper). A
correlation structure can be imposed on a random set of indexable type by
prescribing the joint probability measure on $\Omega$ with the help of a
$d$-dimensional copula. The approach has been further developed in connection
with simulation methods in [4].

Alternatively, copulas can also be introduced so as to act on the target
space, rather than on the underlying probability space $\Omega$. A bit of
terminology is in order. The upper probability of a Borel subset
$B\subset{\mathcal{A}}$ defined by (6) can be viewed as a so-called _capacity
functional_ [41], say on the family ${\mathcal{K}}({\mathcal{A}})$ of compact
subsets of ${\mathcal{A}}$, while the lower probability defines a _containment
functional_. In the multidimensional case, there are two possible ways of
extending capacity functionals, following [59]. The first one is to define a
capacity functional on the family of sets belonging to
${\mathcal{K}}({\mathcal{A}}_{1})\times\cdots\times{\mathcal{K}}({\mathcal{A}}_{d})$
directly, in which case it is termed a _multivariate capacity functional_. The
second choice is to define a capacity functional on the larger family of sets
${\mathcal{K}}({\mathcal{A}}_{1}\times\cdots\times{\mathcal{A}}_{d})$,
referred to as a _joint capacity functional_. It has been shown in [58] that
multivariate containment functionals can be reconstructed from their
univariate marginals by means of a family of copulas (but not by a single
copula, in general). On the other hand, joint capacity functionals can always
be reconstructed from their marginal capacities by means of a single copula,
provided the marginal capacities are maxitive [59] (this is the case if the
marginals arise from normalized fuzzy sets as indicated in Subsection 2.3).
The relation of this approach to the one of [2] and others has been discussed
in [57].

_Structured focal elements._ The third possibility to model dependence is
through the structure of the focal elements. The geometric approach would
replace, say in the case of two variables, the rectangles
$A_{1}(\omega_{1})\times A_{2}(\omega_{2})$ defining the joint focal sets by
rhombic shapes [7, 69]. Alternatively, any random set delivers a family of
probability distributions, see e.g. [48]. One may keep the focal elements as
rectangles but restrict the admissible joint probability distributions [20].
In the case of interval fields, one may use explicitly defined basis functions
multiplied by interval factors [74] or _spatial dependency functions_ [62],
which are the counterpart to Karhunen-Loève expansions in the interval
setting.

##  7 Summary and conclusions

The paper addressed mathematical and modelling issues when heterogeneous
uncertainty is entered in engineering models. The model inputs were assumed to
consist of random variables or random fields, combined with interval
parameters. This led to model outputs which were both random and set-valued.
This set-up is not new; it actually has been widely adopted in the engineering
literature [11, 22, 38, 42, 50, 51, 60, 63, 64, 66], see also the literature
review in the introduction. However, few papers have addressed the question
whether the set-valued model output satisfies the measurability conditions
which make it a random set in its proper sense [3, 53, 55, 56]. The main
contribution of this paper was to show how the theoretical measurability
questions can be addressed, and that the continuous dependence of the solution
and the involved random and interval quantities on their parameters is the
essential ingredient. This is important in as much as computing upper and
lower probabilities in the sense of (6), (7) is only justified when the model
output is a random set.

It was assumed here that the underlying probability spaces were infinite and
that the quantity of interest was the model output itself. If one is only
interested in upper and lower bounds on certain probabilities or moments, one
may avoid collecting the output in a random set. Rather, one may work directly
with the families of random variables defining the input, propagating them
through the model and computing the desired probability bounds. In this
parametric approach, which has been adopted and discussed e.g. in [3, 16, 18,
72], the measurability issue is much simpler: it reduces to proving that the
propagated variables are measurable (that is, still random variables), which
is guaranteed if the input-output map is continuous.

Another simple situation in which the measurability question does not arise is
the case of finite random sets, because on a finite probability space (with
the power set as the $\sigma$-algebra) every map is measurable.

The focus was on model outputs that were generated through solutions to
elliptic and hyperbolic partial differential equations. The most difficult
part was the proof that the model output depended continuously on the input
parameters and hyperparameters. There is no general method of proof; different
arguments are needed for each type of equations. It is expected that new
arguments will be needed when considering other types of equations, such as
parabolic equations. As the paper was oriented toward the theoretical
background of a widely adopted framework, only a simple illustrative example
was given. Many larger scale applications can be found in the papers quoted
above. For the same reason, computational aspects have only been briefly
addressed. For more extensive analyses of the computational effort, the reader
is referred to [4, 14, 18, 21].

An important future research question is the modelling of correlations. It is
expected that the framework of [59] can encompass most other approaches. In
particular, it will be interesting to investigate how it applies in the case
of interval fields.

## Acknowledgements

This work was supported by the FWF doctoral program Computational
Interdisciplinary Modelling W1227 of the Austrian Science Fund as well as a
doctoral scholarship of the University of Innsbruck. Thanks are due to Thomas
Fetz for providing the Matlab-code solving the deterministic FE-problem (17).
The authors are grateful to the two referees for valuable comments. Thanks to
the first referee, the important issue of correlation modelling has been
included. The detailed and critical comments of the second referee led to many
clarifications and to substantial improvements of the presentation.

## Appendix: Explicit construction of infinite probability spaces

In Subsection 3.2 it was necessary to introduce common probability spaces for
random fields with different variances and correlation lengths, first for
random fields defined through the Langevin equation, and second for the
infinitely many Gaussian random variables arising in the Karhunen-Loève
expansion. The appendix briefly explains the mathematical details.

_Wiener space $(\Omega_{W},\Sigma_{W},P_{W})$._ Here
$\Omega_{W}={\mathcal{C}}([0,\infty))$ is the space of continuous functions
$\omega=\omega(x)$ on the half-line. It is a metric space; the metric is
generated by the semi-norms $\|\omega\|_{a,b}=\max_{a\leq x\leq
b}|\omega(x)|,$ $0\leq a<b<\infty$. Then $\Sigma_{W}$ is the corresponding
Borel $\sigma$-algebra. The probability measure $P_{W}$ is the measure on
$\Omega_{W}$ which is generated by the transition probabilities of Wiener
process [67]. In this setting, Wiener process
$W:[0,\infty)\times\Omega_{W}\to\mathbb{R}$ is defined through the point
evaluations $W(x,\omega)=\omega(x)$. Note that the realizations of the process
are identified with the random elements, as is usual practice when
constructing probability spaces on which stochastic processes live. The one-
point and two-point probabilities are, for points $0<x<y$,

$P_{W}(a\leq\omega(x)\leq b)=\int_{a}^{b}p(x,z)\hskip
1.0pt{\mathrm{d}}z\quad\mbox{with}\quad p(x,z)=\frac{1}{\sqrt{2\pi
x}}\exp\big{(}-\frac{z^{2}}{2x}\big{)}$ $P_{W}(a\leq\omega(x)\leq
b,c\leq\omega(y)\leq
d)=\int_{a}^{b}p(x,z_{1})\int_{c}^{d}p(y-x,z_{2}-z_{1})\hskip
1.0pt{\mathrm{d}}z_{2}\hskip 1.0pt{\mathrm{d}}z_{1}$

and similarly for the joint probabilities at finitely many points.

_Infinite products of Gaussian spaces._ For $k=1,2,3,\ldots$, let
$\Omega_{k}=\mathbb{R}$ with the Borel $\sigma$-algebra $\Sigma_{k}$ and the
standard Gaussian measure $P_{k}$,

$P_{k}(B_{k})=\frac{1}{\sqrt{2\pi}}\int_{B_{k}}\exp(-\frac{1}{2}\omega_{k}^{2})\hskip
1.0pt{\mathrm{d}}\omega_{k}$

for any Borel subset $B_{k}\subset\mathbb{R}$. The spaces
$(\Omega_{k},\Sigma_{k},P_{k})$ are all identical copies of standard Gaussian
space, but the index is needed to form their product. The standard Gaussian
variable on $\Omega_{k}=\mathbb{R}$ is defined through the identity map ${\rm
id}(\omega_{k})=\omega_{k}$. Again, the realizations of the random variable
are identified with the random elements.

Let $\Omega$ be the infinite product of the spaces $\Omega_{k}$. The random
elements are sequences $\omega=(\omega_{1},\omega_{2},\omega_{3},\ldots)$. The
probability measure on $\Omega$ is the product measure

$P(B_{1}\times B_{2}\times\ldots\times
B_{n})=\prod_{k=1}^{n}P_{k}(B_{k})=(2\pi)^{-n/2}\prod_{k=1}^{n}\int_{B_{k}}\exp(-\frac{1}{2}\omega_{k}^{2})\hskip
1.0pt{\mathrm{d}}\omega_{k}$

for any finite collection of Borel sets $B_{1},B_{2},\ldots,B_{n}$. The
infinite dimensional random variable $\xi:\Omega\to\mathbb{R}$,
$\xi(\omega)=(\xi_{1}(\omega),\xi_{2}(\omega),\ldots)$ is given through the
projections to $\Omega_{k}$, that is,
$\xi_{k}(\omega)=\xi_{k}(\omega_{1},\omega_{2},\ldots)=\omega_{k}$.

## References

  * [1] D.A. Alvarez,  On the calculation of the bounds of probability of events using infinite random sets,  International Journal of Approximate Reasoning, 43 (2006) 241–267. https://doi.org/10.1016/j.ijar.2006.04.005.
  * [2] D.A. Alvarez,  A Monte Carlo-based method for the estimation of lower and upper probabilities of events using infinite random sets of indexable type,  Fuzzy Sets and Systems, 160 (2009) 384–401. https://doi.org/10.1016/j.fss.2008.08.006.
  * [3] D.A. Alvarez, J.E. Hurtado, J. Ramírez,  Tighter bounds on the probability of failure than those provided by random set theory,  Computers & Structures, 189 (2017) 101–113. http://dx.doi.org/10.1016/j.compstruc.2017.04.006.
  * [4] D.A. Alvarez, F. Uribe, J.E. Hurtado,  Estimation of the lower and upper bounds on the probability of failure using subset simulation and random set theory,  Mechanical Systems and Signal Processing, 100 (2018) 782–801. https://doi.org/10.1016/j.ymssp.2017.07.040.
  * [5] L. Arnold,  Stochastic Differential Equations: Theory and Applications,  Wiley-Interscience, New York-London-Sydney, 1974.
  * [6] M. Beer, S. Ferson, V. Kreinovich,  Imprecise probabilities in engineering analysis,  Mechanical Systems and Signal Processing, 37 (2013) 4–29. http://dx.doi.org/10.1016/j.ymssp.2013.01.024.
  * [7] A. Bernardini, F. Tonon,  Bounding Uncertainty in Civil Engineering – Theoretical Background,  Springer, Berlin, 2010.
  * [8] H. Bunke,  Gewöhnliche Differentialgleichungen mit zufälligen Parametern,  Akademie-Verlag, Berlin, 1972.
  * [9] C. Castaing, M. Valadier,  Convex analysis and measurable multifunctions,  Springer-Verlag, Berlin-New York, 1977.
  * [10] M. de Angelis, E. Patelli, M. Beer,  Advanced Line Sampling for efficient robust reliability analysis,  Structural Safety, 52 (2015) 170–182. http://dx.doi.org/10.1016/j.strusafe.2014.10.002.
  * [11] M.D. Do, W. Gao, C. Song, M. Beer,  Interval spectral stochastic finite element analysis of structures with aggregation of random field and bounded parameters,  International Journal for Numerical Methods in Engineering, 108 (2016) 1198–1229. https://doi.org/10.1002/nme.5251.
  * [12] G. Duvaut, J.-L. Lions (1976),  Inequalities in Mechanics and Physics,  Springer-Verlag, Berlin-New York, 1976.
  * [13] M. Faes, M. Broggi, K.-K. Phoon, M. Beer,  Distribution-free P-box processes: definition and simulation,  in: A. Sofi, G. Muscolino, R.L. Muhanna (Eds.), REC 2021: Proceedings of the 9th International Workshop on Reliable Engineering Computing, 2021, pp. 380–394. http://ww2new.unime.it/REC2021/proceedings/REC2021_Proceedings.pdf.
  * [14] M. Faes, M. Daub, S. Marelli, E. Patelli, M. Beer,  Engineering analysis with probability boxes: A review on computational methods,  Structural Safety, 93 (2021) 102092. https://doi.org/10.1016/j.strusafe.2021.102092.
  * [15] M. Faes, D. Moens,  Identification and quantification of spatial interval uncertainty in numerical models,  Computers and Structures, 192 (2017) 16–33. http://dx.doi.org/10.1016/j.compstruc.2017.07.006.
  * [16] M. Faes, D. Moens,  Imprecise random field analysis with parametrized kernel functions,  Mechanical Systems and Signal Processing, 134 (2019) 106334. https://doi.org/10.1016/j.ymssp.2019.106334.
  * [17] S. Ferson, L. Ginzburg, V. Kreinovich, D. Myers, K. Sentz,  Constructing probability boxes and Dempster-Shafer structures,  Technical Report SAND2002-4015, Sandia National Laboratories, Albuqueque and Livermore, 2003. https://doi.org/10.2172/809606.
  * [18] T. Fetz,  Efficient Computation of upper probabilities of failure,  in: C. Bucher, B.R. Ellingwood, D.M. Frangopol (Eds.), Proceedings of the 12th International Conference on Structural Safety and Reliability (ICOSSAR 2017), TU-MV Media Verlag GmbH, Wien, 2017, pp. 493–502. https://owncloud.tuwien.ac.at/index.php/s/DnQ2sWRcnhT8Zfq#pdfviewer.
  * [19] T. Fetz,  Improving the convergence of iterative importance sampling for computing upper and lower expectations,  Proceedings of Machine Learning Research, 103 (2019) 185–193. https://proceedings.mlr.press/v103/fetz19a.html.
  * [20] T. Fetz, M. Oberguggenberger,  Propagation of uncertainty through multivariate functions in the framework of sets of probability measures,  Reliability Engineering & System Safety, 85 (2004) 73–87. https://doi.org/10.1016/j.ress.2004.03.004.
  * [21] T. Fetz, M. Oberguggenberger,  Imprecise random variables, random sets, and Monte Carlo simulation,  International Journal of Approximate Reasoning, 78 (2016) 252–264. https://doi.org/10.1016/j.ijar.2016.06.012.
  * [22] W. Gao, C. Song, F. Tin-Loi,  Probabilistic interval analysis for structures with uncertainty,  Structural Safety, 32 (2010) 191–199. https://doi.org/10.1016/j.strusafe.2010.01.002.
  * [23] C. Garetto, M. Oberguggenberger,  Symmetrisers and generalised solutions for strictly hyperbolic systems with singular coefficients,  Mathematische Nachrichten, 288 (2015) 185–205. https://doi.org/10.1002/mana.201400192.
  * [24] R. Ghanem, P. Spanos,  Stochastic Finite Elements: A Spectral Approach,  Springer, New York, 1991.
  * [25] I.R. Goodman, H.T. Nguyen,  Fuzziness and randomness,  in: C. Bertoluzza, M. Gil, D.A. Ralescu (Eds.), Statistical Modeling, Analysis and Management of Fuzzy Data, Physica-Verlag, Heidelberg, 2002, pp. 3–21.
  * [26] K.F. Graff,  Wave Motion in Elastic Solids,  Oxford University Press, Oxford, 1975.
  * [27] A. Gray, A. Wimbush, M. de Angelis, P.O. Hristov, D. Calleja, E. Miralles-Dolz, R. Rocchetta,  From inference to design: A comprehensive framework for uncertainty quantification in engineering with limited information,  Mechanical Systems and Signal Processing, 165 (2022) 108210. https://doi.org/10.1016/j.ymssp.2021.108210.
  * [28] M. Grigoriu,  Probabilistic models for stochastic elliptic partial differential equations,  Journal of Computational Physics, 229 (2010) 8406–8429. https://doi:10.1016/j.jcp.2010.07.023.
  * [29] P. Grisvard,  Elliptic Problems in Nonsmooth Domains.  Pitman Advanced Publishing Program, Boston, 1985.
  * [30] X. Han, P.E. Kloeden,  Random Ordinary Differential Equations and Their Numerical Solution,  Springer, Singapore, 2017.
  * [31] J.E. Hurtado, D.A. Alvarez, J.A. Paredes,  Interval reliability analysis under the specification of statistical information on the input variables,  Structural Safety, 65 (2017) 35–48. https://doi.org/10.1016/j.strusafe.2016.12.005.
  * [32] C. Jiang, J. Zheng, X. Han,  Probability-interval hybrid uncertainty analysis for structures with both aleatory and epistemic uncertainties: a review,  Structural and Multidisciplinary Optimization, 57 (2018) 2485–2502. https://doi.org/10.1007/s00158-017-1864-4.
  * [33] O. Kallenberg,  Foundations of Modern Probability,  Springer, New York, 1997.
  * [34] J. Karakašević, M. Oberguggenberger,  Random set solutions to partial differential equations,  in: A. Sofi, G. Muscolino, R.L. Muhanna (Eds.) REC 2021: Proceedings of the 9th International Workshop on Reliable Engineering Computing, 2021, pp. 279–294. http://ww2new.unime.it/REC2021/proceedings/REC2021_Proceedings.pdf.
  * [35] T. Kato,  Perturbation Theory for Linear Operators,  Springer-Verlag, Berlin-New York, 1976.
  * [36] M. Kitahara, S. Bi, M. Broggi, M. Beer,  Nonparametric Bayesian stochastic model updating with hybrid uncertainties,  Mechanical Systems and Signal Processing, 163 (2022) 108195. https://doi.org/10.1016/j.ymssp.2021.108195.
  * [37] H.B. Liu, C. Jiang, Z. Xiao,  Efficient uncertainty propagation for parameterized p-box using sparse-decomposition-based polynomial chaos expansion,  Mechanical Systems and Signal Processing, 138 (2020) 106589. https://doi.org/10.1016/j.ymssp.2019.106589.
  * [38] H. Lü, W.-B. Shangguan, D. Yu,  A new hybrid uncertainty analysis method and its application to squeal analysis with random and interval variables,  Probabilistic Engineering Mechanics, 51 (2018) 1–10. https://doi.org/10.1016/j.probengmech.2017.11.001.
  * [39] H.G. Matthies,  Stochastic finite elements: computational approaches to stochastic partial differential equations,  ZAMM. Zeitschrift für Angewandte Mathematik und Mechanik, 88 (2008) 849–873. https://doi.org/10.1002/zamm.200800095.
  * [40] S. Mishra, Ch. Schwab, J. Šukys,  Multi-level Monte Carlo finite volume methods for uncertainty quantification of acoustic wave propagation in random heterogeneous layered medium,  Journal of Computational Physics, 312 (2016) 192–217. https://doi.org/10.1016/j.jcp.2016.02.014.
  * [41] I. Molchanov,  Theory of Random Sets,  Springer-Verlag, Berlin, 2005.
  * [42] G. Muscolino, R. Santoro, A. Sofi,  Reliability analysis of structures with interval uncertainties under stationary stochastic excitations,  Computer Methods in Applied Mechanics and Engineering, 300 (2016) 47–69. http://dx.doi.org/10.1016/j.cma.2015.10.023.
  * [43] G. Muscolino, A. Sofi, M. Zingales,  One-dimensional heterogeneous solids with uncertain elastic modulus in presence of long-range interactions: Interval versus stochastic analysis,  Computers & Structures, 122 (2013) 217–229. http://dx.doi.org/10.1016/j.compstruc.2013.03.005.
  * [44] J. Nečas,  Les Méthodes Directes en Théorie des Équations Elliptiques,  Masson et Cie, Paris; Academia, Prague, 1967.
  * [45] J. Nedeljković,  Partial differential equations with random set and random field coefficients,  PhD thesis, Universität Innsbruck, Austria, 2020. https://diglib.uibk.ac.at/urn:nbn:at:at-ubi:1-68791.
  * [46] H.T. Nguyen,  An Introduction to Random Sets,  Chapman & Hall, Boca Raton, 2006.
  * [47] M. Oberguggenberger,  Multiplication of Distributions and Applications to Partial Differential Equations,  Longman, Harlow UK, 1992.
  * [48] M. Oberguggenberger,  Combined methods in nondeterministic mechanics,  in: I. Elishakoff, C. Soize (Eds.), Nondeterministic Mechanics, Springer, Wien, 2013, pp. 263–356.
  * [49] M. Oberguggenberger,  Stochastic response surfaces, interval analysis, and the reliability of structures,  in: M. Beer, S.-K. Au, J.W. Hall (Eds.), Vulnerability, Uncertainty, and Risk: Quantification, Mitigation and Management, American Society of Civil Engineers, Reston VA, 2014, pp. 864–875. http://dx.doi.org/10.1061/9780784413609.
  * [50] M. Oberguggenberger,  Analysis and computation with hybrid random set stochastic models,  Structural Safety, 52 (2015) 233–243. http://dx.doi.org/10.1016/j.strusafe.2014.05.008.
  * [51] M. Oberguggenberger, W. Fellin,  Reliability bounds through random sets: Non-parametric methods and geotechnical applications,  Computers & Structures, 86 (2008) 1093–1101. https://doi.org/10.1016/j.compstruc.2007.05.040.
  * [52] M. Oberguggenberger, M. Schwarz,  On the measurability of stochastic Fourier integral operators,  in: P. Boggiatto, M. Cappiello, E. Cordero, S. Coriasco, G. Garello, A. Oliaro, J. Seiler (Eds.), Advances in Microlocal and Time-Frequency Analysis, Birkhäuser/Springer, Cham, 2020, pp. 383–401. .
  * [53] M. Oberguggenberger, L. Wurzer,  Random set solutions to stochastic wave equations,  Proceedings of Machine Learning Research, 103 (2019) 314–323. https://proceedings.mlr.press/v103/oberguggenberger19a.html.
  * [54] J. Sadeghi, M. de Angelis, E. Patelli,  Robust propagation of probability boxes by interval predictor models,  Structural Safety, 82 (2020) 101889. https://doi.org/10.1016/j.strusafe.2019.101889.
  * [55] B. Schmelzer,  On solutions to stochastic differential equations with parameters modeled by random sets,  International Journal of Approximate Reasoning, 51 (2010) 1159–1171. https://doi.org/10.1016/j.ijar.2010.08.006.
  * [56] B. Schmelzer,  Set-valued assessments of solutions to stochastic differential equations with random set parameters,  Journal of Mathematical Analysis and Applications, 400 (2013) 425–438. https://doi.org/10.1016/j.jmaa.2012.11.019.
  * [57] B. Schmelzer,  Sklar’s theorem for minitive belief functions,  International Journal of Approximate Reasoning, 63 (2015) 48–61. https://doi.org/10.1016/j.ijar.2015.05.010.
  * [58] B. Schmelzer,  Joint distributions of random sets and their relation to copulas,  International Journal of Approximate Reasoning, 65 (2015) 59–69. https://doi.org/10.1016/j.ijar.2015.01.007.
  * [59] B. Schmelzer,  Multivariate capacity functionals vs. capacity functionals on product spaces,  Fuzzy Sets and Systems, 364 (2019) 1–35. https://doi.org/10.1016/j.fss.2018.07.005.
  * [60] B. Schmelzer, M. Oberguggenberger, C. Adam,  Efficiency of tuned mass dampers with uncertain parameters on the performance of structures under stochastic excitation,  Journal of Risk and Reliability, 224 (2010) 297–308. https://doi.org/10.1243/1748006XJRR310.
  * [61] R. Schöbi, B. Sudret,  Uncertainty propagation of p-boxes using sparse polynomial chaos expansions,  Journal of Computational Physics, 339 (2017) 307–327. https://doi.org/10.1016/j.jcp.2017.03.021.
  * [62] A. Sofi,  Structural response variability under spatially dependent uncertainty: Stochastic versus interval model,  Probabilistic Engineering Mechanics, 42 (2015) 78–86. http://dx.doi.org/10.1016/j.probengmech.2015.09.001
  * [63] A. Sofi, F. Giunta, G. Muscolino,  Reliability analysis of randomly excited FE modelled structures with interval mass and stiffness via sensitivity analysis,  Mechanical Systems and Signal Processing, 163 (2022) 107990. https://doi.org/10.1016/j.ymssp.2021.107990.
  * [64] A. Sofi, G. Muscolino, F. Giunta,  Propagation of uncertain structural properties described by imprecise Probability Density Functions via response surface method,  Probabilistic Engineering Mechanics, 60 (2020) 103020. https://doi.org/10.1016/j.probengmech.2020.103020.
  * [65] A. Sofi, E. Romeo,  A unified response surface framework for the interval and stochastic finite element analysis of structures with uncertain parameters,  Probabilistic Engineering Mechanics, 54 (2018) 25–36. http://dx.doi.org/10.1016/j.probengmech.2017.06.004.
  * [66] J. Song, M. Valdebenito, P. Wei, M. Beer, Z. Lu,  Non-intrusive imprecise stochastic simulation by line sampling,  Structural Safety, 84 (2020) 101936. https://doi.org/10.1016/j.strusafe.2020.101936.
  * [67] D.W. Stroock,  Lectures on Stochastic Analysis: Diffusion Theory,  Cambridge University Press, Cambridge, 1987.
  * [68] F. Tonon,  Using random set theory to propagate epistemic uncertainty through a mechanical system,  Reliability Engineering & System Safety, 85 (2004) 169–181. https://doi.org/10.1016/j.ress.2004.03.010.
  * [69] F. Tonon, C.L. Pettit,  Toward a definition and understanding of correlation for variables constrained by random relations,  International Journal of General Systems, 39 (2010) 577–604. https://doi.org/10.1080/03081070903541273.
  * [70] F. Trèves,  Basic Linear Partial Differential Equations,  Academic Press, New York, 1975.
  * [71] M.C.M. Troffaes,  A note on imprecise Monte Carlo over credal sets via importance sampling,  Proceedings of Machine Learning Research, 62 (2017) 325–332. https://proceedings.mlr.press/v62/troffaes17a.html.
  * [72] M.C.M. Troffaes,  Imprecise Monte Carlo simulation and iterative importance sampling for the estimation of lower previsions,  Mechanical Systems and Signal Processing, 101 (2018) 31–48. https://doi.org/10.1016/j.ijar.2018.06.009.
  * [73] M.C.M. Troffaes, T. Fetz, M. Oberguggenberger,  Iterative importance sampling for estimating expectation bounds under partial probability specifations,  in: M. De Angelis (Ed.), REC 2018: Proceedings of the 8th International Workshop on Reliable Engineering Computing, Liverpool University Press, Liverpool, 2002, pp. 147–154. http://www.rec2018.uk/papers/proceedings/proceedings.pdf.
  * [74] W. Verhaeghe, W. Desmet, D. Vandepitte, D. Moens,  Interval fields to represent uncertainty on the output side of a static FE analysis,  Computer Methods in Applied Mechanics and Engineering, 260 (2013) 50–62. https://doi.org/10.1016/j.cma.2013.03.021.
  * [75] P. Wei, J. Song, S. Bi, M. Broggi, M. Beer, Z. Lu, Z. Yue,  Non-intrusive stochastic analysis with parameterized imprecise probability models: I. Performance estimation,  Mechanical Systems and Signal Processing, 124 (2019) 349–368. https://doi.org/10.1016/j.ymssp.2019.01.058.
  * [76] P. Wei, J. Song, S. Bi, M. Broggi, M. Beer, Z. Lu, Z. Yue,  Non-intrusive stochastic analysis with parameterized imprecise probability models: II. Reliability and rare events analysis,  Mechanical Systems and Signal Processing, 126 (2019) 227–247. https://doi.org/10.1016/j.ymssp.2019.02.015.
  * [77] H. Zhang,  Interval importance sampling method for finite element-based structural reliability assessment under parameter uncertainties,  Structural Safety, 38 (2012) 1–10. https://doi.org/10.1016/j.strusafe.2012.01.003.
  * [78] J. Zhang, M.D. Shields,  On the quantification and efficient propagation of imprecise probabilities resulting from small datasets,  Mechanical Systems and Signal Processing, 98 (2018) 465–483. http://dx.doi.org/10.1016/j.ymssp.2017.04.042.