Numerical Comparison of Leja and Clenshaw-Curtis
Dimension-Adaptive Collocation for Stochastic Parametric
Electromagnetic Field Problems

arXiv:1712.07223v1 [] 19 Dec 2017

Dimitrios Loukrezis1,2 , Ulrich Römer3 , and Herbert De Gersem1,2
1

Institut für Theorie Elektromagnetischer Felder, Technische Universität Darmstadt,
Schlossgartenstraße 8, 64289 Darmstadt, Germany

2

Graduate School for Computational Engineering, Technische Universität Darmstadt,
Dolivostraße 15, 64293 Darmstadt, Germany
3

Institut für Dynamik und Schwingungen, Technische Universität Braunschweig,
Schleinitzstraße 20, 38106 Braunschweig, Germany

Abstract
We consider the problem of approximating the output of a parametric electromagnetic field
model in the presence of a large number of uncertain input parameters. Given a sufficiently
smooth output with respect to the input parameters, such problems are often tackled with
interpolation-based approaches, such as the stochastic collocation method on tensor-product or
isotropic sparse grids. Due to the so-called curse of dimensionality, those approaches result in
increased or even forbidding computational costs. In order to reduce the growth in complexity with the number of dimensions, we employ a dimension-adaptive, hierarchical interpolation
scheme, based on nested univariate interpolation nodes. Clenshaw-Curtis and Leja nodes satisfy the nestedness property and have been found to provide accurate interpolations when the
parameters follow uniform distributions. The dimension-adaptive algorithm constructs the approximation based on the observation that not all parameters or interactions among them are
equally important regarding their impact on the model’s output. Our goal is to exploit this
anisotropy in order to construct accurate polynomial surrogate models at a reduced computational cost compared to isotropic sparse grids. We apply the stochastic collocation method
to two electromagnetic field models with medium- to high-dimensional input uncertainty. The
performances of isotropic and adaptively constructed, anisotropic sparse grids based on both
Clenshaw-Curtis and Leja interpolation nodes are examined. All considered approaches are
compared with one another regarding the surrogate models’ approximation accuracies using a
cross-validation error metric.
keywords– dimension adaptivity, Clenshaw-Curtis, computational electromagnetics, electromagnetic field simulations, hierarchical interpolation, Leja, sparse grids, stochastic collocation,
uncertainty quantification.

1

Introduction

More often than ever before, the design phase of electric and electronic devices, e.g. electric machines
or high-frequency (HF) components, incorporates parameter studies in order to predict the device’s
1

behavior under uncertainty. This uncertainty, e.g. with respect to the device’s geometry, materials or
other, often stems from tolerances during the manufacturing process. As part of those uncertainty
quantification (UQ) studies, one typically investigates a specific output of the device, called the
quantity of interest (QoI), and tries to estimate statistical moments or sensitivities, with the goal
of reducing the risk of malfunction, misfire or other type of failure. Another common goal is the
construction of a sufficiently accurate surrogate model, which could serve as a computationally
inexpensive substitute to the original model. This latter goal will be the focus of the present paper.
Given a QoI which depends smoothly on a set of input parameters, UQ problems can be efficiently
resolved with the use of polynomial approximations. Most commonly, stochastic Galerkin [2, 12, 22],
stochastic collocation [1, 3, 6, 35] or point collocation [5, 24, 25] methods are employed. The
stochastic Galerkin method is often labeled as “intrusive”, due to the fact that dedicated solvers
have to be developed in order to tackle the stochastic problem at hand. The additional programming
effort is usually regarded as a major disadvantage, especially in the case of complex computational
models whose software and underlying solvers are difficult to be accessed, modified or otherwise
manipulated. Therefore, and despite the fact that stochastic Galerkin methods have nice properties
for error analysis and estimation, collocation methods are generally preferred in this context, as they
allow for a non-intrusive, black-box use of the original computational models. It must be noted that
the separation of methods into intrusive and non-intrusive is an ongoing topic of discussion, see e.g.
[13]. In the context of the present paper we shall retain the usual distinction. Comparisons between
stochastic and point collocation methods, see e.g. [9], indicate that the former tends to provide
superior accuracies and convergence rates. However, since these approaches differ significantly, a
fair comparison between the two is still an open research topic, as also indicated in [24].
A common bottleneck of all aforementioned methods is the so-called “curse of dimensionality” [4],
i.e. convergence rates deteriorate and computational costs increase with the number of considered
input parameters, often exponentially. As a possible remedy, state-of-the-art methods employ sparse,
adaptively constructed polynomial approximations, see e.g. [7, 26, 29, 32] for adaptive stochastic
collocation methods and [5, 23] for adaptive point collocation methods. In the case of adaptive
stochastic collocation methods, sparse approximations are based on nested interpolation nodes.
When the parameters follow uniform distributions, a typical choice is Clenshaw-Curtis nodes [29, 32].
More recently, Leja interpolation nodes have been proposed as a possible alternative [7, 26]. While
generally not free of the curse of dimensionality, adaptive methods exploit possible anisotropies
among the input parameters and their interactions regarding their impact upon the QoI. Assuming
that such anisotropies exist, adaptive approximations may enable studies with a comparably large
number of input parameters. Besides adaptive approximations, alternative methods aiming to tackle
the curse of dimensionality are also available, e.g. tensor decompositions (see [15] and the references
therein) or active subspaces (see [8] and the references therein).
In the present paper, we consider medium- to high-dimensional UQ electromagnetic field (EMF)
problems. Our goal is to construct surrogate models which could reliably replace the original ones
in computationally demanding tasks, e.g. failure probability estimations or robust optimization.
Therefore, the surrogate model’s accuracy is a critical factor. For that purpose, we investigate the
application of the sparse grids stochastic collocation method in the context of EMF simulations.
Specifically, we compare isotropic and adaptive, anisotropic sparse collocation schemes employing
Clenshaw-Curtis and Leja interpolation nodes with respect to accuracy and computational work.
Isotropic sparse grids are constructed with Smolyak’s formula [33]. The adaptive construction of
sparse grids is based on a dimension-adaptive algorithm, similar to the ones presented in [11] for
adaptive quadrature and in [20] for adaptive interpolation. In the search of an acceptable compromise
between computational work and approximation accuracy, adaptivity is a topic which is receiving
increasing attention in uncertainty quantification. In the same context, the choice of collocation
2

nodes employed in adaptive schemes remains an active field of research. However, comparisons of
different strategies have rarely been addressed. Available literature is scarce and so far considers
problems in a purely mathematical context, e.g. approximations of multivariate functions or parametric elliptic partial differential equations [7, 26, 27]. To our knowledge, such comparisons have not
been performed in the context of concrete engineering applications, such as EMF models and simulations. Similarly, Leja nodes remain practically unknown as an option for UQ or other parametric
problems in complex engineering applications.
The rest of the paper is organized as follows. In Section 2 we offer a general description of
the approximation problem at hand. In Section 3 we present the approximation method of choice,
i.e. the interpolation-based stochastic collocation method on sparse grids. Univariate and tensorproduct collocation, which are used as “building blocks” for sparse grids collocation, are presented
in Sections 3.1 and 3.2, respectively. The construction of isotropic and adaptive, anisotropic sparse
grids is presented in Sections 3.3.1 and 3.3.2, respectively. Section 4 presents the EMF applications
and the respective numerical results. All considered collocation methods are first tested on an
academic benchmark example in Section 4.1. A real-world application is presented in Section 4.2.

2

Problem Setting

Let (Θ, Σ, P ) be a probability space, θ ∈ Θ a random event and Y = (Y1 , Y2 , . . . , YN ) a vector of
N independent random variables (RVs) defined on (Θ, Σ, P ). We denote with y = Y (θ) a random
realization of the input RVs and with %Y their joint probability density function (PDF), such that
Y : Θ → Γ ⊂ RN and %Y : Γ → R+ . In the context of the present work, Γ is an N -dimensional
hyper-rectangle Γ = Γ1 × Γ2 × · · · × ΓN . A uniform univariate PDF %n (yn ) corresponds to each Γn ,
n = 1, 2, . . . , N . Due to the statistical independence of the RVs, the joint PDF is given by
%Y (y) =

N
Y

%n (yn ) .

(1)

n=1

We are interested in approximating an output of a computational model which depends on the
RVs Y. This output is typically referred to as the quantity of interest (QoI) and is here denoted
with q. Preferably, the dependence of the QoI on the input RVs is smooth, ideally analytic. We seek
an approximation in the form of a polynomial surrogate model, such that
X
q (y) ≈ q̃ (y) =
si Ψi (y) ,
(2)
i

where Ψi are multivariate polynomials and si are real or complex coefficients. Multivariate polynomials Ψi can be conveniently chosen as tensor-product polynomials
Ψi (y) =

N
Y

ψin (yn ) ,

(3)

n=1

where ψin refers to a univariate polynomial in dimension n. Each combination of local indices
i1 , . . . , iN in (3) is uniquely associated to a global index i.
We aim at reducing the error of the surrogate model uniformly in Γ, which is relevant in the
context of failure probabilities. We assume that the model which provides the considered QoI is
complex, in the sense that its software and underlying solvers are difficult to be accessed and modified. Therefore, we rely on non-intrusive methods, which use the original computational model in
3

a black-box fashion, i.e. without any modifications. Moreover, we assume that the computational
cost of one call to the original computational model, e.g. in order to compute the QoI for a given
set of parameter values, is significantly higher than all other numerical operations during the surrogate model’s construction. Therefore, the total cost of the approximation can be predominantly
attributed to the required number of original model evaluations, i.e. calls to the original model.
Consequently, we are interested in minimizing the number of model evaluations while still achieving
desired accuracies.

3

Sparse Grids Stochastic Collocation

Using the stochastic collocation method, the polynomial approximation (2) is constructed by interpolating specific values of the QoI. Those values are provided by evaluating the QoI on a set
of realizations
of the input RVs, called collocation points. We denote the set of collocation points
n
o
(i)
with Ξ = y
and its cardinality with #Ξ. The choice of collocation points depends on the PDF
%Y which characterizes the input RVs. Lagrange polynomials are typically employed for the interpolation. Since the QoI has to be evaluated on all collocation points and because each evaluation
requires a call to the original deterministic model, the computational cost of the method depends
predominantly on #Ξ.

3.1

Univariate collocation

Univariate interpolation rules are used as building blocks for stochastic collocation in multiple dimensions, respectively, multiple random parameters. Therefore, let us first consider the case of only
one such parameter, Yn .
We introduce the non-negative integer `n ∈ N0 , called the interpolation level. The corresponding
univariate grid of interpolation nodes, respectively, collocation points, is denoted with Ξ`n . The
number of univariate collocation points #Ξ`n is associated with the level `n through a monotonically
increasing “level-to-nodes” function mn : N0 → N, such that #Ξ`n = mn (`n ). The choice of
collocation points depends on the PDF %n . In the case of uniform distributions, common choices are
Gauss-Legendre, Clenshaw-Curtis or Leja interpolation nodes. Gauss-Legendre nodes correspond to
zeros of Legendre polynomials, Clenshaw-Curtis nodes to extrema of Chebyshev polynomials and
Leja nodes are sequences of points such that
y (K) = arg max

y∈[a,b]

K−1
Y
k=0

y − y (k) ,

(4)

where a and b are the limits of the corresponding uniform distribution.
n
omn (`n )
(i )
Given a set of univariate collocation points Ξ`n = yn n
, the nodal Lagrange basis
in =1

polynomials of degree mn (`n ) − 1 are defined as

mn (`n )

li`nn (yn ) =

Y

k=1,k6=in

l10 (yn ) = 1.

4

(k)

yn − yn
(i )

(k)

yn n − yn

,

(5)

Assuming now that the QoI depends only on RV Yn and denoting the univariate Lagrange operator
with I`n , the interpolation formula reads
mn (`n )

I`n q (yn ) =

3.2

X

in =1



q yn(in ) li`nn (yn ) .

(6)

Tensor-product collocation

In its simplest form, multivariate collocation consists of tensor-product combinations of univariate
interpolation grids and operators. We introduce the multi-index ` = (`1 , `2 , . . . , `N ) ∈ NN
0 , which
contains the interpolation level for each RV. Generally, indices `1 , `2 , . . . , `N can have different values
from one another. The special case where `n = `, ∀n = 1, 2, . . . , N , is called “isotropic” tensorproduct collocation.
The multivariate collocation points are given as the tensor grid
Ξ` = Ξ`1 × Ξ`2 × · · · × Ξ`N
n
om1 (`1 ) n
om2 (`2 )
n
omN (`N )
(i )
(i )
(i )
= y1 1
× y2 2
× · · · × yNN
,
i1 =1

i2 =1

(7)

iN =1

with cardinality #Ξ` = #Ξ1 #Ξ2 · · · #ΞN . Every multivariate collocation
point y(i) ∈ 
Ξ` is uniquely

(i1 ) (i2 )
(i )
(i)
defined by the multi-index i = (i1 , i2 , . . . , iN ), such that y = y1 , y2 , . . . , yNN . The corre-

sponding multivariate nodal Lagrange polynomials, L`i are defined as
L`i (y) =

N
Y

li`nn (yn ) .

(8)

n=1

Denoting the multivariate Lagrange operator corresponding to multi-index ` with I` , the interpolation formula reads
I` [q] (y) = (I`1 ⊗ I`2 ⊗ · · · ⊗ I`N ) q (y)
m1 (`1 ) m2 (`2 )

=

X

i1 =1

=

X

X

i2 =1

i : y(i) ∈Ξ`

mN (`N )

···

X

iN =1

N

Y
(i ) (i )
(i )
q y1 1 , y2 2 , . . . , yNN
li`nn (yn )



q y(i) L`i (y) .

(9)

n=1

While simple in its conception and construction, tensor-product stochastic collocation becomes
intractable for a large number N of RVs. The curse of dimensionality is particularly evident in the
N
case of isotropic tensor grids, where #Ξ` = m (`) , i.e. the required computational work increases
exponentially with respect to N . Therefore, the use of tensor-product stochastic collocation is
restricted to low-dimensional settings.

3.3

Sparse Grids Collocation

In high-dimensional settings, collocation on sparse grids is typically employed as a way to mitigate
the computational cost of the full tensor-product collocation. Sparse grids are based on specific
combinations of tensor-product grids, resulting in a reduced number of collocation points and, respectively, model evaluations. Depending on the grid’s construction process, we distinguish between
isotropic and adaptive, anisotropic sparse grids and corresponding collocation schemes.
5

3.3.1

Isotropic sparse grids

We first refer to the construction of isotropic sparse grids [1, 3, 6, 28]. We define the approximation
level ` ∈ N0 and enforce the restriction
|`| =

N
X

n=1

`n ≤ `.

(10)

The sparse grid is then constructed as a linear combination of all tensor-product grids Ξ` , as defined in (7), whose multi-index ` satisfies the admissibility condition (10). Further, introducing the
difference operators
∆`n := I`n − I`n −1 ,

∆` := ∆`1 ⊗ ∆`2 ⊗ · · · ⊗ ∆`N ,

where I−1 is the null operator, the Smolyak interpolation formula [33] reads


X
X
N −1
S` q (y) =
∆` q (y) =
(−1)`−|`|
I` q (y) .
` − |`|
|`|≤`

(11)
(12)

(13)

`−N +1≤|`|≤`

We should note that formula (13) is in general not interpolatory, except for the case of nested
univariate collocation points, for which holds that Ξ`n −1 ⊂ Ξ`n [3]. The corresponding sparse grid
of collocation points is given by
[
ΞS` =
Ξ` .
(14)
`−N +1≤|`|≤`

While a closed-form formula for the computation
of the exact number
of collocation points #ΞS`


N −1
is not available, the growth complexity is O m (`) (log m (`))
[6], which is much milder than


N
the O m (`)
for the full tensor-product collocation. Therefore, even though isotropic sparse grids
are obviously not free of the curse of dimensionality, they are usable up to a moderate number of
RVs.
3.3.2

Adaptive anisotropic sparse grids

It is often the case that certain parameters or parameter combinations and interactions have a more
significant impact on the QoI than others. This parameter anisotropy can be exploited by using
anisotropic sparse grids [29]. In this way, the sparse grid is constructed only with those multiindices and collocation points which add significant contributions to the approximation. Therefore,
we expect to save computational effort, while obtaining similar or better accuracies. In this section,
we present a greedy, dimension-adaptive strategy to achieve this goal. Our approach is based on
the methods and algorithms presented in [11] for quadrature purposes and in [20] for interpolation
purposes. While Leja nodes have not been considered in those works, the extension of the presented
methods and algorithms to Leja nodes is straightforward. Similar dimension-adaptive approaches
are suggested in [7, 26] for Leja nodes and in [17, 32] for Clenshaw-Curtis nodes.
We first re-visit the univariate collocation described in Section 3.1 and restrict the choice of
interpolation nodes to nested sequences Ξ`n −1 ⊂ Ξ`n . The nestedness constraint is satisfied by
Clenshaw-Curtis nodes with the level-to-nodes function mn (`n ) = 2`n + 1. As can be observed in
(4), Leja nodes are nested by construction and the corresponding level-to-nodes function is mn (`n ) =
6

`n + 1. We should note that, while in this work we only consider uniformly distributed parameters,
the ideas and algorithms are applicable to any distribution, assuming that a corresponding family
of nested nodes is available. E.g., in the case of Gaussian distributions, Gauss-Kronrod-Patterson
nodes satisfy the nestedness constraint [32]. Moreover, the use of Leja nodes with non-uniform
probability measures, particularly the ones relative to the Wiener-Askey polynomial chaos [36], is
discussed in [26].
Employing nestedness, higher order approximations can be hierarchically constructed by a stepby-step refinement of the interpolation level and re-use of the lower level approximations at each
step [7, 20, 32]. Using the difference operator defined in (11) the hierarchical interpolation operator
is defined as
`n
X
I `n =
∆k .
(15)
k=0

We define the hierarchical Lagrange basis polynomials h`inn , such that
(
(i )
yn n ∈ Ξ`n \ Ξ`n −1 ,
l`n (yn ) ,
`n
hin (yn ) = i`nn −1
(i )
hin (yn ) , yn n 6∈ Ξ`n \ Ξ`n −1 ,

(16)

where li`nn is the nodal Lagrange basis polynomial defined in (5). Then, the hierarchical univariate
Lagrange interpolation operator reads
 



X
(17)
I`n q (y) = I`n −1 q (y) +
q yn(in ) − I`n −1 q yn(in ) h`inn (y) .
{z
}
|
(in )
in : yn ∈Ξ` \Ξ` −1
n

n

n
hierarchical surplus s`in

(i )

The hierarchical surplus s`inn determines the contribution of the new interpolation node yn n ∈
Ξ`n \ Ξ`n −1 compared to the already available approximation I`n −1 .
In the multivariate case, we assume that the multi-indices which are employed in the interpolation
form the multi-index set Λ. The corresponding interpolation operator is denoted with IΛ and the
corresponding sparse grid with ΞΛ . We enforce upon the multi-index set Λ a monotonicity constraint,
such that
∀` ∈ Λ ⇒ ` − en ∈ Λ, ∀n = 1, 2, . . . , N, with `n > 0,
(18)
where en = (δmn )1≤m≤N is the n-th unit vector and δmn denotes the Kronecker delta. Monotone
sets are also known as downward-closed sets. Then, the multivariate interpolation is defined as
X
IΛ q (y) =
∆` ,
(19)
`∈Λ

where the multivariate difference operators ∆` are defined as in (12). It must be noted that the
isotropic sparse grids described in Section 3.3.1 satisfy the monotonicity constraint by construction.
Assuming a multi-index ` 6∈ Λ such that the set Λ ∪ ` is monotone, the hierarchical interpolation
scheme reads
 



X
(20)
IΛ∪` q (y) = IΛ q (y) +
q y(i) − IΛ q y(i) Hi` (y) ,
{z
}
i : y(i) ∈ΞΛ∪` \ΞΛ |
hierarchical surplus s`i

where Hi` are multivariate hierarchical Lagrange polynomials, defined as
Hi` (y) =

N
Y

n=1

7

h`inn (yn ) .

(21)

Data: QoI q, number of RVs N , monotone multi-index set Λ, tolerance , budget B
Result: sparse grid ΞΛ∪A(Λ) , approximation IΛ∪A(Λ)
repeat
`new = arg max`∈A(Λ) η`
update Λ and IΛ as in (20)
update A (Λ) as in (23)
compute s`i , ∀` ∈ A (Λ) as in (20)
compute η` , ∀` ∈ A (Λ) as in (24)
until termination[A (Λ) , , B];
Algorithm 1: Dimension-adaptive collocation algorithm
Regarding the adaptive refinement of the polynomial approximation and the construction of the
corresponding sparse grid, i.e. the choice of the multi-index ` to be added to set Λ in (20), several
approaches are available in literature, see e.g. [7, 11, 14, 17, 26, 27, 29, 32]. A short description of
the dimension-adaptive algorithm employed in the context of this work is provided in Algorithm 1.
A general overview of the approach follows.
We assume that a monotone multi-index set Λ is available. All possible refinements of the
multi-indices comprising set Λ form the refinement set
R (Λ) = {` + en , ∀` ∈ Λ, ∀n = 1, 2, . . . , N } .

(22)

The admissible multi-indices to be added to set Λ form the admissible set
A (Λ) = {` ∈ R (Λ) : ` 6∈ Λ and Λ ∪ ` is monotone} .

(23)

We note that, when an approximation and its corresponding monotone multi-index set are not
already available, Algorithm 1 is initialized with Λ = {(0, 0, . . . , 0)}. Each one of the admissible
indices in A (Λ) defines a set of collocation points to be possibly added to ΞΛ . Using (20), we
compute the hierarchical surpluses which correspond to the new collocation points. The dimensionadaptivity is based on error indicators η` , ` ∈ A (Λ), computed after the corresponding hierarchical
surpluses. In the context of the present work, error indicators are given by
X
1
s`i ,
(24)
η` =
# (ΞΛ∪` \ ΞΛ )
(i)
i:y

∈ΞΛ∪` \ΞΛ
`∈A(Λ)

as proposed in [11]. Other suggestions on computing local error indicators can be found in [17, 26, 32].
The admissible multi-index to be added to Λ is the one with the maximum contribution compared to
the already available approximation, i.e. the one corresponding to the maximum η` . The algorithm
terminates either when a pre-defined simulation budget B, i.e. number of model evaluations, is
reached or when the total contribution of set A (Λ) is below a specified tolerance  [11, 17, 26, 20, 32].
The termination conditions can be formally formulated as
#ΞΛ∪A(Λ) ≥ B,
X
η` ≤ .

(25a)
(25b)

`∈A(Λ)

Since the hierarchical surpluses and collocation points for the admissible multi-indices are already
available, the sparse grid and the corresponding approximation after the algorithm’s termination is
constructed with all multi-indices in the set Λ ∪ A (Λ).
8

4

Numerical Experiments

We test four variants of the sparse grids stochastic collocation method, namely isotropic and adaptive, anisotropic sparse grids collocation based on Clenshaw-Curtis and Leja interpolation nodes.
In the presented results, “CC” and “LJ” refer to Clenshaw-Curtis and Leja nodes, respectively.
Accordingly, subscript “iso” refers to isotropic sparse grids and subscript “adapt” to anisotropic,
adaptively constructed grids. The Sparse Grids MATLAB Kit [34] has been employed for the implementation of isotropic sparse grids. Results of the Clenshaw-Curtis isotropic sparse grids have
been validated with the Sparse Grids Interpolation Toolbox [19, 20], which has also been used for
implementing the Clenshaw-Curtis adaptive collocation. We note that the authors of [19, 20] use the
term “Clenshaw-Curtis” for nodes given by the Newton-Cotes rule and the term “Chebyshev” for the
extrema of Chebyshev polynomials which are typically known as Clenshaw-Curtis nodes. Therefore,
the “Chebyshev” option should be used in the toolbox. An in-house, Python-based implementation
has been developed for the adaptive Leja collocation. As part of our implementation, univariate
Leja nodes are computed with the Python package Chaospy [10]. All collocation variants are applied to two EMF applications, namely a dielectric slab waveguide and a Stern-Gerlach magnet.
Both computational models employ in-house developed MATLAB/Octave and Python softwares.
The accuracy of the surrogate models is measured with a cross-validation error metric. Using a
validation set with M random realizations of the input parameters, the cross-validation error cv is
computed as




cv =
max
qe y(m) − q y(m) ,
(26)
m=1,2,...,M

where qe and q denote the surrogate and the original model, respectively.

4.1

Dielectric Slab Waveguide

We consider an academic example from the field of high-frequency computational electromagnetics.
In particular, we consider a three-dimensional, rectangular, dielectric slab waveguide, as the one
illustrated in Fig. 1. Such structures provide us with simple models, typically used to study
wave-confinement mechanisms. The waveguide operates in the frequency range [fmin , fmax ]. We
consider uncertainties with respect to the waveguide’s geometry and material parameters. We shall
construct approximations with respect to the waveguide’s reflection coefficient, also known as the
S11 parameter, which quantifies the reflection of the incoming wave at port 1 of the waveguide. For
this simple model, an analytical solution for S11 exists and is used in this first example in order to
avoid the consideration of discretization errors.
4.1.1

Model

As can be seen in Fig. 1, the waveguide is filled with vacuum (blue area) and has a dielectric material
in its middle (yellow area). Its geometry is defined by its width w along the x-axis, its height h along
the y-axis and the vacuum offset d and dielectric filling length l along the z-axis. The red planes
denote the waveguide’s input and output ports, respectively port 1 and 2. With the exception of the
ports, the walls of the waveguide are considered to be perfect electrical conductors (PEC). Using

9

Figure 1: 3D waveguide model with dielectric filling. The yellow area denotes the dielectric filling
and the blue area the vacuum. The red planes denote the waveguide ports.
the Cartesian coordinate system, the waveguide’s computational domain is formally defined as
Ω = [0, w] × [0, h] × [0, 2d + l] ,

(27a)

ΓPEC = {(x, y, z) ∈ ∂Ω : z 6= 0 ∧ z 6= 2d + l},
Γin = {(x, y, z) ∈ ∂Ω : z = 0},

Γout = {(x, y, z) ∈ ∂Ω : z = 2d + l},

E × n = 0,

n × curlE + γn × (n × E) = 0,

(27c)
(27d)

such that ∂Ω = ΓPEC ∪ Γin ∪ Γout . Maxwell’s source problem for the electric field E reads

curl µ−1 curlE − ω 2 εE = 0,
in Ω,
n × curlE + γn × (n × E) = Uinc ,

(27b)

(28a)

on ΓPEC ,

(28b)

on Γin ,

(28c)

on Γout ,

(28d)

where γ = jk inc and k inc refer to the wavenumber of the incoming wave represented by Uinc . See
[18] for details. A Debye relaxation model is used for the description of the dielectric material’s
behavior, such that the material’s electric permittivity ε and magnetic permeability µ depend on
the waveguide’s angular frequency ω as follows:
!
K
X
εs,k − ε∞
ε = ε0 ε∞ +
,
(29a)
1 + jωτε,k
k=1
!
K
X
µs,k − µ∞
.
(29b)
µ = µ0 µ∞ +
1 + jωτµ,k
k=1

In (29), K denotes the relaxation model’s order, index 0 refers to the material properties in vacuum,
index ∞ to high frequencies and index s to the static case. For both material properties, τ denotes
the characteristic relaxation time of the medium. We assume that the incoming field coincides with
10

Table 1: Nominal values of the waveguide’s geometrical, material and source parameters.
Parameter
operating frequency
width
height
filling length
vacuum offset
vacuum permittivity
high-frequency permittivity
static permittivity 1
static permittivity 2
vacuum permeability
high-frequency permeability
static permeability 1
static permeability 2
permittivity relaxation time 1
permittivity relaxation time 2
permeability relaxation time 1
permeability relaxation time 2

Symbol
f
w
h
l
d
ε0
ε∞
εs,1
εs,2
µ0
µ∞
µs,1
µs,2
τε,1
τε,2
τµ,1
τµ,2

Value, Model 1
18-22
30
3
7
5
8.85 · 10−12
1
2
4π · 10−7
1
2.4
1/ (40π)
1.1/ (40π)
-

Value, Model 2
18-22
30
3
7
5
8.85 · 10−12
1
2
2.2
4π · 10−7
1
2
3
1/ (40π)
1.1/ (40π)
1/ (40π)
2/ (40π)

Units
[GHz]
[mm]
[mm]
[mm]
[mm]
[F/m]
[F/m]
[F/m]
[F/m]
[H/m]
[H/m]
[H/m]
[H/m]
[ns]
[ns]
[ns]
[ns]

the fundamental TE10 mode e10 and that higher order modes are quickly attenuated in the structure.
In this case, the reflection coefficient at port 1 is an important output quantity, given as
Z
q = S11 = C inc
E · e10 dx,
(30)
Γin

where C inc is a normalization constant.
In the most general case, we may consider uncertainties with respect to geometry, material or
source parameters. Then, the parametric counterpart of (28) reads


−1
2
curl µ (y) curlE − ω (y) ε(y)E = 0,
in Ω(y),
(31a)
E × n = 0,

n × curlE + γn × (n × E) = U
n × curlE + γn × (n × E) = 0,

inc

(y) ,

on ΓPEC (y),

(31b)

on Γin (y),

(31c)

on Γout (y),

(31d)

ρY -almost everywhere in Γ. A common approach in uncertainty quantification with random geometries is to pull-back the parametric equations to a fixed reference domain. This ensures the
tensor-product structure of the solution space. Here, we do not approximate the solution itself, but
only a scalar quantity of interest and hence, this transformation is not required. The parametric
quantity of interest is given as
Z
q(y) = C inc
E(y) · e10 (y) dx.
(32)
Γin (y)

We consider uncertainties with respect to the waveguide’s geometry and material properties.
The corresponding RVs follow uniform distributions centered around the respective nominal values,
presented in Table 1. For the geometrical parameters we allow a 10% variation around the nominal
11

10−2

collocation points

105
104
cv error

10−5

103

10−8

102

iso/LJ
iso/CC

101

adapt/CC

10−11

iso/LJ

0

2

4
6
8
10
approximation level

12

12

adapt/LJ

101

(a) level-to-nodes relation.

max. interpol. level

iso/CC

102
103
104
model evaluations

105

(b) Cross-validation errors.

12

12

adapt/CC
adapt/LJ

10
8
5

4

4

4

0 0

0
f

w

8

3

3

εs,1

µs,1

0 0

h

l

d

(c) Maximum interpolation level per parameter.

Figure 2: Numerical results for the dielectric slab waveguide with 1st order Debye relaxation and 7
parameters.
value, while an 1% variation is assumed for material parameters. Moreover, the waveguide’s operating frequency may take any value in the frequency range [fmin , fmax ]. Frequency is therefore added
to the parameter vector with the assumption that it is uniformly distributed in the frequency range.
Two waveguide models are considered, corresponding to medium- and high-dimensional uncertainty, respectively. The first model employs 1st order Debye relaxation and a total of 7 parameters,
namely f , w, h, l, d, εs,1 , µs,1 . The second model employs 2nd order Debye relaxation and a total
of 15 parameters, namely f , w, h, l, d, εs,1 , εs,2 , ε∞ , τε,1 , τε,2 , µs,1 , µs,2 , µ∞ , τµ,1 , τµ,2 .
4.1.2

Numerical results

Results for the 1st waveguide model with 7 parameters are presented in Fig. 2. Results for the
2nd waveguide model with 15 parameters are shown in Fig. 3. For the isotropic cases, Fig. 2(a)
and 3(a) show the relation between the approximation level and the required computational work,
equivalently number of collocation points. As already discussed in Section 3, each collocation point
12

104

cv error

collocation points

10−3

105

103
10

10−4
10−5

iso/CC
iso/LJ

2
iso/CC

adapt/CC

10−6

iso/LJ

101

2
4
6
approximation level

adapt/LJ

101

(a) level-to-nodes relation.

6

max. interpol. level

6

6

102
103
104
model evaluations

105

(b) Cross-validation errors.

6

adapt/CC
adapt/LJ

4

4

4
3
2

2

2

2

2

3

00
f

w

h

3

3

3

3
2

2
1

0

3

1

1

1

1

1

1

1

00
l

d

εs,1 εs,2 ε∞ τε,1 τε,2 µs,1 µs,2 µ∞ τµ,1 τµ,2

(c) Maximum interpolation level per parameter.

Figure 3: Numerical results for the dielectric slab waveguide with 2nd order Debye relaxation and
15 parameters.
corresponds to a single model evaluation. The curse of dimensionality is particularly evident in
Fig. 3(a). Cross-validation error results based on a random sample with M = 104 entries are shown
in Fig. 2(b) and 3(b). For the most accurate adaptive approximations available, we present the
maximum interpolation levels per parameter in Fig. 2(c) and 3(c) as an indication to the anisotropy
among the input parameters.
In the case of the 1st waveguide model, Leja nodes generally outperform Clenshaw-Curtis for
both isotropic and dimension-adaptive collocation. However, the advantage of Leja nodes is not
significant. Overall, we may say that both choices of collocation nodes perform comparably well,
with Leja nodes having a slight edge. As shown in Fig. 2(c), parameters h and d are completely
ignored by the dimension-adaptive algorithm. The most refined parameters are f , w and l, followed
be the two material parameters. Due to its level-to-nodes function mn (`n ) = `n + 1, the most
accurate Leja-based surrogate model employs polynomials of lower degrees than its Clenshaw-Curtis
counterpart, where mn (`n ) = 2`n + 1.
In the case of the high-dimensional 2nd waveguide model, Clenshaw-Curtis nodes outperform
13

8
7

y [mm]

6
4
3
1
(1)

Ωp

Ωbeam Ω(2)
p

(3)

0-8

Ωp

(a) 3D magnet model.

(3)

Ωp

5

2

dΩp

(2)

Ωp = Ωair

(1)

Ωp

-6

Ωbeam
-4

-2
0
x [mm]

2

4

(b) Zoom in magnet’s pole region.

Figure 4: Right: 3D model of one-half of a Rabi-type Stern-Gerlach magnet. Left: Zoom in the
magnet’s pole region. Modified pictures from [30].
Leja in the case of isotropic sparse grids. The opposite is true for the dimension-adaptive collocation, where Leja nodes perform consistently better than Clenshaw-Curtis. For the most part, the
advantage of Leja nodes is significant, particularly for errors below 10−5 . As in the case of model
1, parameters h and d are disregarded by the dimension-adaptive algorithm, which again refines
parameters f , w and l the most, followed by the material parameters. In this case, the resulting
polynomial degrees for both choices of interpolation nodes are comparable.

4.2

Stern-Gerlach Magnet

We consider a real-world application, in particular a Rabi-type Stern-Gerlach magnet (see Fig.
4(a)), similar to the one described in [21] and further studied in [30, 31]. This magnet is currently
in use at KU Leuven. Stern-Gerlach magnets are typically employed for the magnetic separation
of atom beams or clusters. A key design requirement is a homogeneous magnetic field with a
strong gradient. Due to design and manufacturing imperfections, the pole region might suffer from
geometrical uncertainties, which in their turn affect the field homogeneity and gradient. The aim of
this study is to apply the aforementioned UQ methods in order to quantify the impact of geometrical
uncertainties onto the QoI, which is here the average magnetic field gradient in the magnet’s beam
area.
4.2.1

Model

All computations are performed using a linearized two-dimensional model of the magnet’s crosssection, as in [30]. The magnet’s pole region, denoted with Ωp , is the only domain which is spatially
(i)
resolved. Domain Ωp is decomposed into distinct patches Ωp , i = 1, . . . , 3, with numbering from left
(1)
(2)
(3)
(2)
to right, such that Ωp = Ωp ∪ Ωp ∪ Ωp . Region Ωp refers to the air gap inside the magnet’s pole
(1)
(3)
region, while regions Ωp and Ωp to the regions on the left and right of the air gap, respectively,
as in Fig. 4(b). The contributions of the remaining yoke part and the coils are taken into account
by a field-circuit coupling and a magnetic equivalent circuit [30]. More precisely, in a first step, the
14

magnetic vector potential and the magnetic flux through the iron yoke are computed for the entire
geometry. The values are denoted as A0z and Φ0 , respectively. Then, the coupling is realized by
imposing
Az = AΦ = Φ/Φ0 A0z , on ∂Ωp ,
(33)
where Φ is recomputed for a different geometry using magnetic circuit theory. Let I and Nc represent
the coil current and the number of turns in the winding, respectively. Then this relation can be
abstractly written as F (Az , Φ) = Nc I, where F refers to the magnetomotive force. For details, the
reader is referred to [30]. In summary, the field-circuit coupled problem reads
div (ν gradAz ) = 0,
Az − AΦ = 0,

in Ωp ,

(34a)

on ∂Ωp ,

(34b)

F (Az , Φ) = Nc I,

(34c)

~ is obtained as B
~ =
where ν refers to the magnetic reluctivity. The magnetic flux density B
(∂y Az , −∂x Az , 0). The magnet’s beam area is denoted with Ωbeam and lies inside the air gap of
~

(2)

B|
the pole region, denoted with Ωair , where Ωair = Ωp (see Fig. 4(b)). Denoting with τ (x, y) = ∂|∂x
the magnetic field gradient in the x-direction, the average field gradient in the beam area is given
by
Z
1
τ (x, y) dΩ.
(35)
q = τavg =
|Ωbeam | Ωbeam

Isogeometric analysis (IGA) is employed for the spatial discretization [16]. In IGA, both the
solution variable Az and the geometry are described in terms of non-uniform rational B-splines
(NURBS). A generic NURBS curve reads
~
R(ξ)
=

N
X

P~i Nip (ξ),

i=1

ξ ∈ [0, 1],

(36)

where P~i and Nip refer to a control point and a NURBS basis function of degree p, respectively.
NURBS basis functions are defined as
wi B p (ξ)
Nip (ξ) = PN i p ,
j=1 wj Bj (ξ)

(37)

with weights wi and B-spline basis functions Bip , respectively.
(i)
The original NURBS curves defining the three patches Ωp , i = 1, 2, 3, are depicted in Fig. 5 in
blue, while a geometry variation is given in red. We introduce random geometry deformations by
regarding the control points and weights of these curves as uncertain. More precisely, we introduce
10 RVs corresponding to the x and y coordinates of 4 control points and 2 weights, as presented in
Table 2. Then, we obtain a random reluctivity as
ν(y) = νiron 1Ω(1) (y) + νair 1Ω(2) (y) + νiron 1Ω(3) (y) ,
p

p

p

(38)

with 1Ω(i) denoting the characteristic function of patch i and νiron and νair denoting the reluctivity
p
of iron and air, respectively. Accordingly, the parametric field-circuit coupled problem reads
div (ν(y) gradAz ) = 0,
Az − AΦ = 0,

in Ωp ,

(39a)

on ∂Ωp ,

(39b)

F (Az , Φ) = Nc I,
15

(39c)

8
7
6

Original geometry
Original control points
Geometry variation

(2)

Ωp = Ωair

(3)

Ωp

y [mm]

5

4
2

4
3
2
1
0
-8

1

3

(1)
Ωp

-6

Ωbeam
-4

-2
x [mm]

0

2

4

Figure 5: Original and varied geometry of the magnet’s pole region. Modified picture from [30].
with ρY -almost everywhere in Γ. Assuming that problem (39) is well-posed, the quantity of interest
is itself a random variable.
Table 2: Control point uncertainties and uniform distribution limits.
Parameter
x1
y1
x2
y2
w2
x3
y3
x4
y4
w4

4.2.2

Bottom Limit
−3.00
1.66
−5.0
2.50
0.35
1.50
1.89
−2.0
4.0
0.37

Upper Limit
−2.00
2.0
−2.0
4.50
2.85
2.50
2.5
2.0
5.50
2.87

Units
[mm]
[mm]
[mm]
[mm]
[-]
[mm]
[mm]
[mm]
[mm]
[-]

Numerical results

All results are presented in Fig. 6. Fig. 6(a) shows the relation between approximation level
and number of collocation points of the isotropic sparse grids for the case of 10 parameters. The
Stern-Gerlach magnet simulation is based on an expensive numerical model, therefore approximation
levels ` = 5 for the Clenshaw-Curtis case and or ` = 7 for the Leja case already result in hampering
16

105
104
cv error

collocation points

101

3

10

102

100
10−1

iso/LJ
iso/CC

10

iso/LJ

101

2
4
6
approximation level

8

max. interpol. level

adapt/LJ

101

102
103
model evaluations

104

(b) Cross-validation errors.

9

adapt/CC

9
8

adapt/LJ

6

6
4

adapt/CC

−2

(a) level-to-nodes relation.

8

iso/CC

5

5
4

3

5

3

4

3

3

2

2
x1

y1

5

4
3

4
3

2
w1

x2

y2

w2

x3

y3

x4

y4

(c) Maximum interpolation level per parameter.

Figure 6: Numerical results for the Stern-Gerlach magnet with 10 parameters.
computational costs. Because of the high computational cost of the magnet model, cross-validation
errors are now computed using a sample with M = 103 entries. Error results are presented in Fig.
6(b). Maximum interpolation levels per parameter for the most accurate approximations are shown
in Fig. 6(c).
We observe that isotropic Clenshaw-Curtis collocation outperforms significantly its Leja-based
counterpart. Both dimension-adaptive schemes outperform isotropic collocation. In this case, computational gains are not as pronounced as in the waveguide example. This could be attributed to the
considered QoI and possible difficulties in its approximation. Another possible explanation could
be that strong anisotropy among the parameters is not present in this problem. This explanation is
supported by the results presented in Fig. 6(c) for the Clenshaw-Curtis case, where all parameters
are refined with levels either 3 or 4. Stronger parameter anisotropy can be observed in the Leja
case, however this can be attributed to the more granular level-to-nodes function of Leja nodes.
It is possible that the approximation levels become comparable as the approximation accuracy increases. In general, the dimension-adaptive algorithm seems to perform similarly for both choices of
interpolation nodes. For the most part, Leja nodes have a small advantage over Clenshaw-Curtis,
17

however both families of interpolation nodes perform comparably. The algorithm seems to recognize
deformations in the y-direction as more important than deformations in the x-direction, hence coordinates y2 , y3 and y4 are refined more than all other parameters. For all parameters, polynomials
employed in the Clenshaw-Curtis approximation are of significantly higher orders than in the Leja
case.

5

Summary and Conclusions

We have investigated the performance of the sparse grids stochastic collocation method, when used
for the construction of accurate polynomial surrogates of computational EMF models. Two models suffering from medium- to high-dimensional input uncertainty have been used for that purpose.
We have employed both isotropic and adaptively constructed, anisotropic sparse grids based on
Clenshaw-Curtis and Leja interpolation nodes, under the assumption that all input parameters are
statistically independent and uniformly distributed. In all examples, dimension-adaptive collocation outperforms the isotropic variant. In cases characterized by strong parameter anisotropy, the
computational gains are even impressive. In all examples, Leja-based dimension-adaptive collocation either has an advantage over its Clenshaw-Curtis-based counterpart with respect to surrogate
model accuracy as measured with the cross-validation error metric, or performs comparably well.
These findings are in agreement with the available literature results, see e.g. [26]. Follow-up studies based on the present work should consider the use of sufficiently accurate surrogate models for
computationally demanding tasks, e.g. surrogate-based robust optimization or failure probability
estimations. Moreover, while the use of Leja nodes in UQ problems with non-uniformly distributed
parameters is suggested in [26], the authors provide numerical examples only for the uniform and
Gaussian case. Therefore, the use of Leja points with different probability distributions should be
further explored.
Acknowledgments
All authors would like to thank our colleague Andreas Pels for providing the Stern-Gerlach magnet
model. The first and third author would like to acknowledge the support of the Graduate School of
Computational Engineering, Technische Universität Darmstadt.

References
[1] I. Babuska, F. Nobile, and R. Tempone, A Stochastic Collocation Method for Elliptic
Partial Differential Equations with Random Input Data, SIAM Review, (2010), pp. 317–355.
[2] I. Babuska, R. Tempone, and G. E. Zouraris, Galerkin Finite Element Approximations
of Stochastic Elliptic Partial Differential Equations, SIAM J. Numerical Analysis, 42 (2004),
pp. 800–825.
[3] V. Barthelmann, E. Novak, and K. Ritter, High Dimensional Polynomial Interpolation
on Sparse Grids, Adv. Comput. Math., 12 (2000), pp. 273–288.
[4] D. Bellman, Dynamic Programming, Princeton University Press, 1957.
[5] G. Blatman and B. Sudret, Adaptive Sparse Polynomial Chaos Expansion Based on Least
Angle Regression, J. Comput. Physics, 230 (2011), pp. 2345–2367.

18

[6] H.-J. Bungartz and M. Griebel, Sparse Grids, Acta Numerica, 13 (2004), pp. 147–269.
[7] A. Chkifa, A. Cohen, and C. Schwab, High-Dimensional Adaptive Sparse Polynomial Interpolation and Applications to Parametric PDEs, Foundations of Computational Mathematics,
14 (2014), pp. 601–633.
[8] P. G. Constantine, Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter
Studies, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2015.
[9] M. S. Eldred, Comparison of Non-Intrusive Polynomial Chaos and Stochastic Collocation
Methods for Uncertaitny Quantification, In: AIAA, (2009).
[10] J. Feinberg and H. P. Langtangen, Chaospy: An Open Source Tool for Designing Methods
of Uncertainty Quantification, J. Comput. Science, 11 (2015), pp. 46–57.
[11] T. Gerstner and M. Griebel, Dimension-Adaptive Tensor-Product Quadrature, Computing,
71 (2003), pp. 65–87.
[12] R. G. Ghanem and P. D. Spanos, Stochastic Finite Elements: A Spectral Approach, SpringerVerlag New York, Inc., New York, NY, USA, 1991.
[13] L. Giraldi, A. Litvinenko, D. Liu, H. G. Matthies, and A. Nouy, To Be or Not to Be
Intrusive? The Solution of Parametric and Stochastic Equations - the ”Plain Vanilla” Galerkin
Case, SIAM J. Scientific Computing, 36 (2014).
[14] M. Griebel, Adaptive Sparse Grid Multilevel Methods for Elliptic PDEs Based on Finite
Differences, Computing, 61 (1998), pp. 151–180.
[15] W. Hackbusch, Numerical Tensor Calculus, Acta Numerica, 23 (2014), p. 651–742.
[16] T. J. Hughes, J. A. Cottrell, and Y. Bazilevs, Isogeometric Analysis: CAD, Finite
Elements, NURBS, Exact Geometry and Mesh Refinement, Computer methods in applied mechanics and engineering, 194 (2005), pp. 4135–4195.
[17] J. D. Jakeman and T. Wildey, Enhancing Adaptive Sparse Grid Approximations and Improving Refinement Strategies Using Adjoint-Based A Posteriori Error Estimates, J. Comput.
Physics, 280 (2015), pp. 54–71.
[18] J.-M. Jin, The Finite Element Method in Electromagnetics, John Wiley & Sons, 2015.
[19] A. Klimke, Sparse Grid Interpolation Toolbox – User’s Guide, Tech. Rep. IANS 2007/017,
University of Stuttgart, 2007.
[20] A. Klimke and B. I. Wohlmuth, Algorithm 847: Spinterp: Piecewise Multilinear Hierarchical Sparse Grid Interpolation in MATLAB, ACM Trans. Math. Softw., 31 (2005), pp. 561–579.
[21] B. Masschaele, T. Roggen, H. De Gersem, E. Janssens, and T. T. Nguyen, Design
of a Strong Gradient Magnet for the Deflection of Nanoclusters, IEEE Transactions on Applied
Superconductivity, 22 (2012).
[22] H. G. Matthies and A. Keese, Galerkin Methods for Linear and Nonlinear Elliptic Stochastic
Partial Differential Equations, Informatik-Berichte der Technischen Universität Braunschweig,
2003-08.

19

[23] G. Migliorati, Adaptive Polynomial Approximation by Means of Random Discrete Least
Squares, in ENUMATH, A. Abdulle, S. Deparis, D. Kressner, F. Nobile, and M. Picasso, eds.,
vol. 103 of Lecture Notes in Computational Science and Engineering, Springer, 2013, pp. 547–
554.
[24] G. Migliorati, F. Nobile, E. von Schwerin, and R. Tempone, Approximation of Quantities of Interest in Stochastic PDEs by the Random Discrete L2 Projection on Polynomial
Spaces, SIAM J. Scientific Computing, 35 (2013).
[25]

, Analysis of Discrete L2 Projection on Polynomial Spaces with Random Evaluations,
Foundations of Computational Mathematics, 14 (2014), pp. 419–456.

[26] A. Narayan and J. D. Jakeman, Adaptive Leja Sparse Grid Constructions for Stochastic
Collocation and High-Dimensional Approximation, SIAM J. Scientific Computing, 36 (2014).
[27] F. Nobile, L. Tamellini, and R. Tempone, Comparison of Clenshaw–Curtis and Leja
Quasi-Optimal Sparse Grids for the Approximation of Random PDEs, Springer International
Publishing, Cham, 2015, pp. 475–482.
[28] F. Nobile, R. Tempone, and C. G. Webster, A Sparse Grid Stochastic Collocation Method
for Partial Differential Equations with Random Input Data, SIAM J. Numerical Analysis, 46
(2008), pp. 2309–2345.
[29]

, An Anisotropic Sparse Grid Stochastic Collocation Method for Partial Differential Equations with Random Input Data, SIAM J. Numerical Analysis, 46 (2008), pp. 2411–2442.

[30] A. Pels, Z. Bontinck, J. Corno, H. De Gersem, and S. Schöps, Optimization of a
Stern-Gerlach Magnet by Magnetic Field-Circuit Coupling and Isogeometric Analysis, IEEE
Transactions on Magnetics, 51 (2015).
[31] U. Römer, S. Schöps, and H. De Gersem, A Defect Corrected Finite Element Approach
for the Accurate Evaluation of Magnetic Fields on Unstructured Grids, J. Comput. Physics, 335
(2017), pp. 688–699.
[32] B. Schieche, Unsteady Adaptive Stochastic Collocation on Sparse Grids, PhD Thesis, TU
Darmstadt, 2012.
[33] S. A. Smolyak, Quadrature and Interpolation Formulas for Tensor Products of Certain Classes
of Functions, Dokl. Acad. Nauk SSSR, 4 (1963), pp. 240–243.
[34] L. Tamellini and F. Nobile, Sparse Grids MATLAB Kit. URL: http://csqi.epfl.ch. Accessed:
15.08.2017.
[35] D. Xiu and J. S. Hesthaven, High-Order Collocation Methods for Differential Equations with
Random Inputs, SIAM J. Scientific Computing, 27 (2005), pp. 1118–1139.
[36] D. Xiu and G. E. Karniadakis, The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations, SIAM J. Scientific Computing, 24 (2002), pp. 619–644.

20

