Uncertainty Quantification with Lévy-type Random Fields Von der Fakultät Mathematik und Physik der Universität Stuttgart zur Erlangung der Würde eines Doktors der Naturwissenschaften (Dr. rer. nat.) genehmigte Abhandlung Vorgelegt von Andreas Stein aus Alsfeld Hauptberichterin: Prof. Dr. Andrea Barth Mitberichter: Prof. Dr. Christian Rohde Prof. Dr. Christoph Schwab Tag der mündlichen Prüfung: 18.06.2020 Institut für Angewandte Analysis und Numerische Simulation der Universität Stuttgart 2020 Acknowledgements I would like to express my sincere gratitude to my advisor Prof. Dr. Andrea Barth for her guidance, support and patience during my time at the University of Stuttgart. It is a pleasure to work with her and I could not have imagined having a better advisor and mentor for my doctoral studies. Also, I would like to thank Prof. Dr. Christian Rohde and Prof. Dr. Christoph Schwab for their willingness to be the co-examiners of this thesis. Another thank goes to my present and former colleagues from the research group of computational uncertainty quantification at the University of Stuttgart for all the in- sightful discussions throughout the years. Moreover, I would like to thank Dr. Alexan- der Grimm for his valuable suggestions that lead to a significant improvement of this thesis. During my time as doctoral student, I have been partly funded by the German Re- search Foundation (DFG) as part of the Cluster of Excellence in Simulation Technology (EXC 310/2) at the University of Stuttgart, and it is gratefully acknowledged. Furthermore, I have to thank Elke Gangl for her administrative support in the submission process of this thesis. Finally, I would like to thank all my colleagues at the IANS in Stuttgart for creating the pleasant atmosphere at the institute and making my time in Stuttgart enjoyable. iii Contents Acknowledgements iii Abstract viii Zusammenfassung x 1 Introduction 1 1.1 The contribution of this thesis . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 The articles in this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 "A study of elliptic partial differential equations with jump dif- fusion coefficients" . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.2 "Numerical analysis for time-dependent advection-diffusion prob- lems with random discontinuous coefficients" . . . . . . . . . . . 6 1.2.3 "A multilevel Monte Carlo algorithm for parabolic advection- diffusion problems with discontinuous coefficients" . . . . . . . . 7 1.2.4 "Approximation and simulation of infinite-dimensional Lévy pro- cesses" . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.5 "A stochastic transport equation with Lévy noise: Fully discrete numerical approximation" . . . . . . . . . . . . . . . . . . . . . 8 1.2.6 Overall structure . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 Discretization techniques in UQ 10 2.1 Numerical methods for PDEs . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Monte Carlo methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3 Lévy processes and stochastic integrals 19 3.1 Covariance operators and spectral expansions . . . . . . . . . . . . . . 20 3.2 Lévy random fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3 Stochastic integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 iv CONTENTS Declaration to the cumulative part 31 4 Elliptic PDEs with jump diffusion coefficient 33 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 Elliptic boundary value problems . . . . . . . . . . . . . . . . . . . . . 36 4.3 Discontinuous random elliptic problems . . . . . . . . . . . . . . . . . . 40 4.3.1 Jump-diffusion coefficients and their approximations . . . . . . . 40 4.3.2 Existence of solutions . . . . . . . . . . . . . . . . . . . . . . . . 44 4.3.3 Convergence of the approximated diffusion coefficient . . . . . . 46 4.3.4 Convergence of the approximated solution . . . . . . . . . . . . 52 4.4 Adaptive path-wise discretization . . . . . . . . . . . . . . . . . . . . . 53 4.4.1 Adaptive triangulations . . . . . . . . . . . . . . . . . . . . . . . 55 4.5 Monte Carlo methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.5.1 Monte Carlo and multilevel Monte Carlo estimators . . . . . . . 58 4.5.2 Coupled multilevel Monte Carlo . . . . . . . . . . . . . . . . . . 64 4.6 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.6.1 Numerical examples in 1D . . . . . . . . . . . . . . . . . . . . . 66 4.6.2 Numerical results in 2D . . . . . . . . . . . . . . . . . . . . . . 71 5 Random discontinuous parabolic problems 78 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.2 Initial boundary value problems . . . . . . . . . . . . . . . . . . . . . . 79 5.3 Random discontinuous problems . . . . . . . . . . . . . . . . . . . . . . 87 5.3.1 Jump-diffusion coefficients and their approximations . . . . . . . 87 5.4 Path-wise discretization schemes . . . . . . . . . . . . . . . . . . . . . . 94 5.4.1 Sample-adapted spatial discretization . . . . . . . . . . . . . . . 95 5.4.2 Temporal discretization . . . . . . . . . . . . . . . . . . . . . . . 103 5.5 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.5.1 Numerical examples in 1D . . . . . . . . . . . . . . . . . . . . . 108 5.5.2 Numerical examples in 2D . . . . . . . . . . . . . . . . . . . . . 112 6 MLMC for discontinuous parabolic PDEs 116 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 6.2 Parabolic problems with random discontinuous coefficients . . . . . . . 118 6.3 Numerical approximation of the solution . . . . . . . . . . . . . . . . . 123 6.3.1 Approximated diffusion coefficients . . . . . . . . . . . . . . . . 123 6.3.2 Semi-discretization by sample-adapted finite elements . . . . . . 125 6.3.3 Fully discrete pathwise approximation . . . . . . . . . . . . . . 129 v CONTENTS 6.4 Estimation of moments by multilevel Monte Carlo methods . . . . . . . 131 6.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 7 Infinite-dimensional Lévy processes 140 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 7.3 Simulation by Fourier inversion . . . . . . . . . . . . . . . . . . . . . . 150 7.3.1 A piecewise constant approximation . . . . . . . . . . . . . . . . 151 7.3.2 Inversion of the characteristic function . . . . . . . . . . . . . . 154 7.4 Generalized hyperbolic Lévy processes . . . . . . . . . . . . . . . . . . 169 7.5 Numerics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 7.5.1 Notes on implementation . . . . . . . . . . . . . . . . . . . . . . 184 7.5.2 Approximation of a GH field . . . . . . . . . . . . . . . . . . . . 185 7.5.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 188 8 A stochastic transport problem with Lévy noise 193 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 8.2 Stochastic partial differential equations with Lévy noise . . . . . . . . . 196 8.3 The Stochastic transport equation . . . . . . . . . . . . . . . . . . . . . 201 8.4 Temporal discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 8.5 Discontinuous Galerkin spatial discretization . . . . . . . . . . . . . . . 214 8.6 Noise approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221 8.7 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 225 8.8 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230 9 Conclusions and outlook 235 Bibliography 237 vi Abstract A countless number of models in the natural sciences, engineering and economics are based on partial differential equations (PDEs). Due to insufficient data or measurement errors, certain characteristics of the underlying PDE are subject to uncertainty, and are usually modeled by continuous and/or Gaussian random fields. Although analytically tractable, the applications of continuous or Gaussian random fields are limited: spatial and temporal discontinuities cannot be captured and Gaussian distributions notoriously underestimate the probability of rare events. To this end, the focus of this thesis is on uncertainty quantification with Lévy-type random fields, a certain class of discontinuous stochastic objects that provide a significant extension to the existing methodology. In a nutshell, this dissertation explains how to incorporate Lévy-type random fields into PDEs and how the corresponding solutions become accessible by the means of discretization and simulation. The first main contribution is the introduction of a novel type of random field, consisting of a Gaussian part and a spatially discontinuous jump field. This Lévy-type field serves as a coefficient in advection-diffusion equations and allows to model, for instance, sudden changes in the permeability of a porous medium far more realistically than state-of-the-art continuous models. In contrast to the few examples in the liter- ature, the discontinuous random coefficient in this thesis provides a unique flexibility as it is able to generate virtually any stochastic geometry. Apart from random PDEs with discontinuous coefficients, hyperbolic transport equations with Lévy noise as source term are considered. The noise processes take values in a (infinite-dimensional) Hilbert space and involve temporal discontinuities. Therefore, heavy-tailed random perturbations are introduced to the transport problem, and the resulting stochastic equation may be utilized, e.g., as a model for the dynamics in commodity forward markets. Due to the lack of tractable discretization schemes for the underlying stochastic PDE, this models have been, up to now, only of theoretical interest. This thesis paves the way to finally apply the forward model with Lévy noise in practice, as it provides the first fully discrete scheme for the corresponding stochastic viii ABSTRACT transport problem. In both cases, PDEs with random jump coefficients or with Lévy noise as source term, regularity is inherently low due to the discontinuities and state-of-the-art nu- merical algorithms are prohibitive. To remedy this issue, several advanced schemes for discontinuous random problems are introduced that outperform standard methods in terms of convergence rates and computational effort. A comprehensive numerical analysis is provided and the superior performance of the new approaches is validated by numerous numerical experiments. Moreover, an emphasis is put on the approximation of infinite-dimensional, discon- tinuous random fields, a crucial part in the discretization of stochastic PDEs. A part of this thesis covers the sampling of Hilbert space-valued Lévy processes, and introduces a sampling technique combining truncated Karhunen-Loève expansions with discrete Fourier inversion. This algorithm stands out as the arguably most flexible of a very limited number of methods for the approximation and simulation of Lévy fields. ix Zusammenfassung Unzählige Modelle in Naturwissenschaft, Technik und Ökonomie basieren of partiellen Differentialgleichungen (PDEs). Aufgrund von Messfehlern oder unvollständigen Daten sind gewisse Parameter dieser PDEs mit Unsicherheiten behaftet und werden da- her in der Regel durch stetige und/oder Gaußsche Zufallsfelder modelliert. Obwohl diese gute analytische Eigenschaften haben, sind die Anwendungen von stetigen oder Gaußschen Zufallsfeldern limitiert: Sprünge in Zeit oder Ort können nicht abgebildet werden und Gaußsche Verteilungen unterschätzen oft die Wahrscheinlichkeiten von extremen Ereignissen. Aus diesen Gründen thematisiert die vorliegende Arbeit Uncer- tainty Quantification unter Verwendung einer Klasse unstetiger Zufallsfelder, die eine signifikante Erweiterung der existierenden Theorie darstellen. Konkret wird erklärt wie man unstetigen Zufallsfelder in PDEs integriert und wie die entsprechenden Lösungen durch Diskretisierungsmethoden und Simulation approximiert werden können. Der erste wichtige Beitrag dieser Arbeit ist die Einführung eines neuartigen Zufalls- feldes, das aus einem Gaußschen Anteil, sowie einem unstetigen Sprungfeld besteht. Dieses unstetige Zufallsfeld findet dann Verwendung als Koeffizient in Advektions- Diffusionsgleichungen, beispielsweise um abrupte Änderungen der Permeabilität in porösen Medien zu simulieren. In der Literatur gibt es bisher nur wenige, sehr spezielle Beispiele von unstetigen stochastischen Koeffizienten. Im Gegensatz dazu besticht der Koeffizient in dieser Arbeit durch seine Flexibilität, da er es ermöglicht praktisch jede Art zufälliger Geometrie abzubilden. Neben zufälligen PDEs mit unstetigen Koeffizienten werden ebenfalls hyperbolis- che Transportgleichungen mit einem Lévy-Prozess als Quellterm auf der rechten Seite betrachtet. Dieser Prozess ist unstetig in der Zeit und nimmtWerte in einem unendlich- dimensionalen Hilbertraum an. Dadurch erhält die Transportgleichung Störungen mit sogenannten heavy-tails, also mit hoher Wahrscheinlichkeit für Extremereignisse, und kann zum Beispiel als Modell für Terminkontrakte in Rohstoffmärkten benutzt wer- den. Da geeignete Diskretisierungsverfahren für deren Implementierung bisher fehlten, waren diese Modelle bis jetzt eher von theoretischem Interesse. In dieser Arbeit wird x ZUSAMMENFASSUNG zum ersten Mal eine volle Diskretisierung für die zugehörigen stochastischen Transport- probleme vorgestellt, und dadurch ermöglicht die Lévy-Modelle für Terminkontrakte auch erstmals in der Praxis anzuwenden. Ein Problem, das sowohl bei PDEs mit zufälligen unstetigen Koeffizienten, als auch bei stochastischen PDEs mit Lévy-Rauschen auftritt, ist die niedrige Regularität der Lösung. Dies macht übliche numerische Verfahren zur Lösung von PDEs ineffizient und aus diesem Grund werden mehrere spezielle Methoden für Probleme mit zufälligen Unstetigkeiten entwickelt. Alle Verfahren werden vollständig numerisch analysiert und zahlreiche Experimenten belegen die besseren Performance im Vergleich zu Standard- Methoden. Ein weiterer Schwerpunkt dieser Arbeit ist die Approximation und Simulation von unendlich-dimensional Lévy-Prozessen, ein essentieller Baustein in der Diskretisierung von stochastischen PDEs. Dazu wird eine Technik vorgestellt die auf einer Kombination von abgeschnittenen Karhunen-Loève-Entwicklungen und diskreter Fourier Inversion besteht. Dieser Algorithmus ragt als die wahrscheinlich flexibelste von sehr wenigen existieren Methoden zur Simulation von Lévy-Feldern heraus. xi 1 Introduction A vast range of phenomena in the natural sciences, engineering, finance and economics are modeled by partial differential equations (PDEs). Examples include the description of porous media and subsurface flows ([68, 107])1, the propagation of aeroacoustic waves ([200]), turbulent flow models ([180]), phase separation in the mixture of fluids ([150]), the heat distribution on CPUs ([55]), optimal control problems in traffic flow networks ([95]), aquifer systems to store thermal energy ([141]), epidemic models for the spreading of diseases ([153]) and parasite populations ([103]), environmental pollution models ([206]), the metabolism of glucose in the human liver ([183]), the valuation of financial derivatives ([42, 111]) and assessment of credit risks ([41]) or models of warfare and pursuit in game theory ([117]). The quantities of interest in these applications are in general functionals of the solution to the PDE. For instance, in a subsurface flow model the solution of the underlying system describes a pressure field, and one may be interested in evaluating the pressure function at specific points in the domain or to know the average pressure in a certain area. Due to insufficient data, measurement errors or incomplete information, however, certain parameters of the PDE may not be known a-priori. In a subsurface flow model, the permeability of a medium is in general only measured at a discrete set of points and unknown in between. Another common example is the pricing of stock options, where the model parameters have to be estimated statistically based on a time series of quoted market prices. To model this effects more realistically, uncertainty quantification has become an increasingly important and popular field of research in the last decades. A common approach in uncertainty quantification is to replace certain characteristics of a PDE model by stochastic objects, for instance, by considering a random differential operator, such as in [1, 14, 15, 25, 29, 53, 59, 63, 85, 90, 96, 148, 164, 165, 190, 194], or to introduce some noise as a function-valued stochastic process, see for instance [23, 24, 28, 35, 67, 69, 86, 121, 132, 151, 152, 173, 198]. The listed references are 1All references listed in the introduction, in the articles and in the conclusion of this thesis can be found in the unified bibliography after Chapter 9 at the end of this dissertation. 1 1.1. THE CONTRIBUTION OF THIS THESIS merely a snapshot of the research activity throughout the recent years and are by far not exhaustive. Throughout the literature, it is standard practice to model the uncertainty by con- tinuous random fields (see e.g. [53, 59, 63, 85, 96, 164, 194]). Due to continuity, these objects have natural representations in terms of Fourier series that may be exploited in analysis and simulation. The major drawback of continuous random fields is, however, that applications are drastically limited as temporal and/or spatial discontinuities can- not be captured. For example, they are not suitable to model flows through porous media or composite materials where sudden changes of permeability at interfaces oc- cur. Another shortcoming is that time-continuous stochastic processes essentially follow Gaussian distributions, which notoriously underestimate the probabilities of extreme risks. For instance, when applied as a subsurface flow model, Gaussian random fields may not reflect the presence of peaks and asymmetric distributions in the permeability appropriately, see [202]. As another example, Gaussian stock return models fail to match real market data in many cases and are unable to predict the occurrence of rare events, such as financial crisis, properly (see e.g. [89]). 1.1 The contribution of this thesis To overcome the aforementioned problems, this thesis focuses on uncertainty quan- tification involving a certain type of discontinuous random objects, which is referred to as Lévy-type random fields. In a nutshell, this dissertation explains how to incor- porate random discontinuous characteristics into PDEs and how the solutions to this problems become accessible by the means of discretization and simulation. As a conse- quence, various problems can now be modeled much more realistically, since the severe restrictions imposed by continuous or Gaussian random objects have vanished. The corresponding solutions may be simulated with controlled bias and reasonable effort using the advanced numerical schemes developed in this thesis. As one of the main contributions, a novel Lévy-type random field with spatial dis- continuities is introduced. The decisive feature of this jump field is its unique flexibility, since it allows to model virtually any stochastic geometry. This is achieved by con- structing the coefficient as the sum of a Gaussian random field and a discontinuous jump term (for this reason the term Lévy-type is borrowed from stochastic analysis). Moreover, temporal Lévy processes taking values in an infinite-dimensional Hilbert space are considered. These objects are function-valued stochastic processes in the classical sense with discontinuous paths in time. Compared to a standard Gaussian noise process, a Lévy process allows to introduce more realistic heavy-tailed perturba- 2 CHAPTER 1. INTRODUCTION tions into a dynamical system. Each kind of Lévy-type random field is utilized as a stochastic parameter in a PDE, either as random coefficient or source term. A fundamental part in the discretization of stochastic PDEs is to obtain tractable closed-form approximations of the underlying random field. While this can be challeng- ing for continuous Gaussian random fields, all difficulties are amplified when dealing with more general Lévy-type random objects. As a further contribution, this thesis in- troduces a new approximation method for general Hilbert space-valued Lévy noises in Chapter 7. In fact, the proposed algorithm stands out as one among very few methods for general Lévy noises and is arguably the most flexible one. This is somewhat surpris- ing at first glance, since there is a vast literature on the numerical analysis of (parabolic) stochastic PDEs with Lévy noise, see for instance [23, 26, 33, 49, 74, 133]. Therein, however, the proposed noise approximation either suffers from strong assumptions or, as in the majority of cases, is completely neglected. Moreover, certain sampling tech- niques for the space-time Lévy source term are also applicable for the spatial random jump coefficient, as outlined in Chapters 4 and 5. Possible applications for PDEs involving Lévy-type random fields are motivated by irregular, discontinuous structures like fractured or porous media, rock strata with spherical inclusions, composite materials and alloys. To model flows, charge- or pressure distributions in this structures, the spatial Lévy-type random field may be utilized as diffusion- or advection coefficient in a PDE. The literature on PDEs with random discontinuous coefficients is sparse, with exceptions being [104, 140, 205], and the analysis is in general tailored to a very specific shape and geometry. In contrast to this, the Lévy-type random field from this thesis has the ability to represent any desired stochastic structure and therefore significantly extends the existing methodology. Another important example are transport equations with Lévy noise, a popular model for commodity- and interest futures markets (see [24, 52, 173]). So far, this models have mainly been examined from a theoretical perspective and have not been applied to valuate forward contracts, since tractable fully discrete approximations have not been available. Chapter 8 of this thesis closes this gap and stands out as the first (and to this point only) article containing a complete numerical analysis for the full discretization of this type of stochastic transport problem. Additionally, all theoretical findings are confirmed by numerical experiments in scenarios with varying regularity. Similar results have not yet been provided in the literature: in general, the much simpler case of finite-dimensional and/or Gaussian noise is investigated and only semi-discrete schemes are proposed (see e.g. [34, 130, 135, 154] for results on nonlinear hyperbolic PDES with Lévy- respectively Gaussian noise). On a further note, there seems to be a general lack of numerical experiments in the literature on SPDEs with Lévy noise. 3 1.2. THE ARTICLES IN THIS THESIS For both scenarios, PDEs with random jump coefficients as well as stochastic PDEs with Lévy noise as source term, suitable concepts of solutions are investigated to ensure well-posedness. The discontinuities entail inherently low regularity and, consequently, standard numerical algorithms converge, if at all, at very poor rates. Even if conver- gence is theoretically ensured, oscillations and instabilities may occur and the corre- sponding schemes are useless in any application. To overcome the issues of common methods, this thesis introduces a variety of novel algorithms to discretize discontinu- ous random PDE problems. A comprehensive numerical analysis is provided for each method and stability is ensured to prevent oscillations in the approximated solutions. To this end, several new techniques of proof have been developed, as applicability of existing results is rather limited. For instance, it is crucial to control for stochastic geometries or simulation biases, effects that naturally do not occur in deterministic or simpler stochastic problems. The superior performance of our novel approaches compared to state-of-the-art methods is validated in numerous experiments. At this point it needs to be emphasized that the considered Lévy-type random ob- jects provide an actual extension to the current state of research, since the aforemen- tioned continuous random fields may be recovered as a special case. The most impor- tant examples for this are log-normal diffusion coefficients as in [53, 54, 96, 98, 108, 194] and Q-Wiener processes as space-time driving noise, see for instance [132, 136, 198]. Therefore, all results and proposed algorithms of this thesis are universal and may readily be applied to a large class of important problems in the field of uncertainty quantification. 1.2 The articles in this thesis The core of this thesis consists of five independent research articles: • A study of elliptic partial differential equations with jump diffusion coefficients, • Numerical analysis for time-dependent advection-diffusion problems with random discontinuous coefficients, • A multilevel Monte Carlo algorithm for parabolic advection-diffusion problems with discontinuous coefficients, • Approximation and simulation of infinite-dimensional Lévy processes, and • A stochastic transport problem with Lévy noise: Fully discrete numerical approx- imation. 4 CHAPTER 1. INTRODUCTION As the articles have been developed partly in parallel, the order is not strictly chronolog- ical, but motivated by thematic classification. The first group consists of Chapters 4–6 that consider elliptic/parabolic problems with random Lévy-type coefficients involving spatial discontinuities. They are connected very closely as Chapter 5 and 6 are build on its predecessor(s) and were therefore also finished in the given order. It is com- mon to refer to the type of stochastic equations in Chapters 4–6 as random PDEs to emphasize that the source of randomness lies within the PDE’s coefficients, and there- fore in the corresponding differential operator. This is to obtain a distinction from so-called stochastic PDEs (SPDEs) with deterministic coefficients respectively differ- ential operators, but a stochastic source term on the right hand side of the equation. Chapters 7 and 8 constitute the second group of research articles and focus on SPDEs with Lévy noise as external source of uncertainty. The overarching theme of this thesis, however, is to introduce Lévy-type objects as stochastic PDE parameters to allow for more flexibility compared to the state-of-the-art continuous and/or Gaussian random fields. Simulation and discretization of the PDE problem at hand becomes significantly more involved due to the discontinuous random structures. Thus, the scientific litera- ture lacks tractable algorithms and this thesis closes the existing gap for a variety of problems. In uncertainty quantification, approximation schemes not only need to include the discretization of spatial and temporal domains, but also the discretization of the stochastic domain. While the first part is achieved by modifying schemes for determin- istic PDEs, the second part involves the approximation of infinite-dimensional random variables and sampling techniques like Monte Carlo algorithms. As outlined through- out the central part of this thesis, sampling becomes especially challenging for Lévy- type random objects. Moreover, depending on whether random PDEs or SPDEs are considered, some aspects of discretization may differ drastically. In the first case, path-wise error bounds (similar as for deterministic PDEs) can be exploited, where path-wise means that this estimates only hold for a specific realization of the prob- lem. To obtain Lp-estimates with p ≥ 1 and almost sure convergence, it is then still crucial to control all appearing terms with respect to the entire stochastic domain. When discretizing SPDEs as in Chapter 8, it is necessary to switch from a path-wise to a mean-square-type perspective of discretization. The Lévy-noise on the right hand side of the SPDE is given by an Itô integral with each path involving a possibly in- finite number of temporal discontinuities. This stochastic integral is in turn defined as the mean-square-limit of a sequence of simple integrals and, therefore, the mean- square theory of the corresponding Itô calculus has to be applied (details are given in Chapter 3 of this thesis). Consequently, path-wise approximation techniques are not 5 1.2. THE ARTICLES IN THIS THESIS applicable and one is usually bound to L2-type convergence results. A short overview on discretization methods with regard to uncertainty quantifica- tion is given in Chapter 2, the specific method for each problem in Chapters 4– 8 is sketched in the subsequent summary and is found in detail within the corresponding article. 1.2.1 A. Barth and A. Stein: "A study of elliptic partial differential equa- tions with jump diffusion coefficients" We consider an elliptic equation with random coefficient, which may for instance ac- count for the uncertain permeability of a given medium in a subsurface flow model. As an extension of this methodology to flows in heterogeneous\fractured\porous me- dia, we incorporate jumps in the diffusion coefficient. More precisely, we consider a second order elliptic problem where the random coefficient is given by the sum of a (continuous) Gaussian random field and a (discontinuous) jump part. To estimate moments of the solution to the resulting random partial differential equation, we use path-wise finite element approximations combined with multilevel Monte Carlo sampling. In order to account for the discontinuities and improve the convergence of the path-wise approximation, the spatial domain is decomposed with respect to the jump positions in each sample, leading to path-dependent finite element grids. Hence, it is not possible to create a sequence of grids which is suitable for each sample path a-priori. We address this issue by an adaptive multilevel Monte Carlo algorithm, where the discretization on each level is sample-dependent and fulfills given refinement conditions. 1.2.2 A. Barth and A. Stein: "Numerical analysis for time-dependent advec- tion-diffusion problems with random discontinuous coefficients" In this article we extend the elliptic diffusion problem from Chapter 4 to a time- dependent parabolic advection-diffusion equation. Specifically, a scenario with coupled advection and diffusion coefficients that are modeled as sums of continuous random fields and discontinuous jump components is considered. We employ the same adaptive, path-wise finite element discretization for the numerical approximation of the solution as in Chapter 4. In this previous work, we have merely assumed convergence in mean- square with respect to the H1(D)-norm for the discretization of the spatial domain D based on our numerical experiments. Here, however, we provide a rigorous proof on the path-wise and mean-squared convergence rate based on a set of suitable assumptions. As it turns out, this becomes surprisingly involved, since the order of convergence 6 CHAPTER 1. INTRODUCTION depends on the shape of the discontinuities. Thus, the random geometry has to be controlled for each sample to obtain bounds on the mean-squared error. To stabilize the numerical approximation and accelerate convergence, the discrete space-time grid is chosen with respect to the varying discontinuities in each sample of the coefficients, leading to a stochastic formulation of the Galerkin projection and the finite element basis. We provide several numerical experiments to verify our theoretical results and show that the regime of assumptions cannot be relaxed significantly. 1.2.3 A. Barth and A. Stein: "A multilevel Monte Carlo algorithm for para- bolic advection-diffusion problems with discontinuous coefficients" In many applications for elliptic or parabolic random PDEs, such as in the previous two articles, the aim is to estimate moments of certain quantities of interest (QoIs). Those QoIs are in general functionals of the corresponding path-wise solution, rather than the solution itself. As there are only path-wise numerical solutions at hand to which a given functional may be applied, we obtain a biased version of the true QoI. For instance, in Chapters 4 and 5 this bias is due to the (adaptive) finite element approximation. It is, however, often possible to express this bias in terms of the L2(D)-error of the spatial discretization. In this article we consider the same problem setting and assumptions as in Chapter 5 and extend the existing a-priori error bounds of the path-wise adaptive discretization from the H1(D)-norm to the L2(D)-norm. To this end, we solve a time-dependent parabolic dual problem, and derive an improved decay rate based on our results for the H1(D)-error in Chapter 5. As expected, the L2(D)-error decays twice as fast on a given spatial grid, allowing us to bound the path-wise bias for functionals of the approximated solution. Based on this result, we introduce a multilevel Monte Carlo algorithm to estimate the moments of a given QoI. A numerical example shows that the adaptive discretization is superior to a standard finite element approach when applied to a multilevel Monte Carlo estimator. 1.2.4 A. Barth and A. Stein: "Approximation and simulation of infinite- dimensional Lévy processes" The recurring motivation in the first three articles of this thesis was to replace a spa- tially continuous stochastic object by a more general, discontinuous Lévy-type random field. We now follow the same idea for time-dependent stochastic objects, and shift our attention to infinite-dimensional Lévy processes, also called (time-dependent) Lévy fields. In fact, the random coefficients in Chapters 4–6 can actually be interpreted 7 1.2. THE ARTICLES IN THIS THESIS as a spatial and stationary analogue to the classical notion of time-dependent Lévy processes from stochastic analysis. As already indicated in the introduction of this thesis, the approximation of Lévy fields is challenging: For square integrable fields beyond the Gaussian case, it is no longer given that the one-dimensional distributions in the spectral representation with respect to the covariance operator are independent. When simulated via a Karhunen- Loève expansion a set of dependent, but uncorrelated, one-dimensional Lévy processes has to be generated. The dependence structure among the one-dimensional processes ensures that the resulting field exhibits the correct point-wise marginal distributions. To approximate the respective (one-dimensional) Lévy-measures, a numerical method, called discrete Fourier inversion, is developed. For this method, Lp-convergence rates can be obtained and, under certain regularity assumptions, mean-square and Lp- convergence of the approximated field is proved (we emphasize that the Fourier inver- sion technique has also been applied to sample spatial Lévy-type random coefficients, as shown in the numerical examples from Chapters 4 and 5). Further, the class of gener- alized hyperbolic Lévy fields is introduced, where the point-wise marginal distributions are dependent, but uncorrelated, subordinated Wiener processes. For this specific class one may derive point-wise marginal distributions in closed form. Numerical examples, including hyperbolic and normal-inverse Gaussian Lévy fields, demonstrate the effi- ciency of the approach. 1.2.5 A. Barth and A. Stein: "A stochastic transport equation with Lévy noise: Fully discrete numerical approximation" Using the results from Chapter 7 as foundation, we introduce fully discrete schemes for SPDEs with Lévy noise. A particular interesting example are linear hyperbolic SPDEs, which serve as a model for the dynamics of interest rate and energy forward markets. To be more precise, the forward rate is modeled as the solution to a transport equation with a space-time stochastic process as driving noise. Again, our motivation is to capture discontinuities in time and allow for heavy-tailed distributions, thus we consider Hilbert space-valued Lévy fields as driving noise terms. The numerical discretization of the corresponding SPDE involves several difficulties: Low spatial and temporal regularity of the solution to the problem entails slow convergence rates and instabilities for space/time-discretization schemes. Even if this problem is resolved, a fully discrete approximation has yet to take the driving noise process into account. The Lévy field admits values in a possibly infinite- dimensional Hilbert space U , hence projections into a finite-dimensional subspace of 8 CHAPTER 1. INTRODUCTION U for each discrete point in time are necessary. Finally, as mentioned in the previous Chapter 7, unbiased sampling from the resulting Lévy field is not necessarily possi- ble. We introduce a fully discrete approximation scheme to address the above issues: A discontinuous Galerkin approach for the spatial approximation is coupled with a suitable time stepping scheme to avoid numerical oscillations. Moreover, we use the Fourier inversion technique combined with truncated Karhunen-Loève expansions from Chapter 7 to obtain a suitable approximation of the Lévy field. We provide a rigorous error analysis for the space-time discretization which yields together with the results from Chapter 7 a bound on the mean-squared overall discretization error. As before, we confirm our theoretical results by several numerical examples. 1.2.6 Overall structure In addition to the five articles, the second chapter of this thesis contains a prelimi- nary discussion on numerical schemes for PDEs and Monte Carlo methods with a view towards uncertainty quantification. Thereafter, Chapter 3 collects some probabilistic results on infinite-dimensional random variables, Lévy processes and stochastic inte- gration. This part provides a theoretical background on random fields and stochastic processes, especially for Chapters 7 and 8 that consider Lévy driving noise. After the two introductory parts on discretization techniques and probability theory, Chapters 4– 8 follow as the core of this thesis and contain the articles previously outlined in this section. Finally, this thesis is concluded by some remarks and perspectives for future work in Chapter 9. 9 2 Discretization techniques in uncertainty quantification This chapter gives an overview of numerical schemes for PDEs and Monte Carlo meth- ods, the discretization techniques applied in the main part of this thesis. A brief description of the key ideas and a literature review with regard to their applications in uncertainty quantification is provided. Advantages as well as possible shortcomings are highlighted and the algorithms of choice are motivated by comparing them to sev- eral alternative approaches. A more detailed and formal description for each scheme is then given in Chapters 4–8 for the specific problem at hand. 2.1 Numerical methods for PDEs This section discusses some standard numerical approximation schemes for PDEs which are applied in Chapters 4–6 and 8 to PDEs with discontinuous random characteristics. The limitations of each approach are pointed out in this context and it is indicated how this problems have eventually been overcome in the forthcoming chapters. One of the most popular methods for the numerical approximation of PDEs is the finite element (FE) method, which can be traced back to the works [8, 65, 87, 115, 184, 195]. There is an extensive amount of literature on the FE method with applications to deterministic PDEs, for instance [47, 58, 102, 177] to just name a few. In the recent years, FE-based methods have also successfully been applied in all areas of uncertainty quantification, examples are found in [15, 23, 27, 29, 54, 85]. The basic idea of the FE method may be illustrated by the stationary, elliptic boundary value problem −∇ · (a(x)∇u(x)) = f(x), x ∈ D, u(x) = 0, x ∈ ∂D, (2.1) on a bounded domain D ⊂ Rd. For any real q > 0, let Hq(D) denote the standard Sobolev space over D with the notational convention H0(D) := L2(D). The suitable 10 CHAPTER 2. DISCRETIZATION TECHNIQUES IN UQ solution space is given by H1 0 (D) := {v ∈ H1(D)| γv = 0}, with γ : H1(D) → L2(∂D) being the trace operator. Hence, the weak formulation of Problem (2.1) is to find u ∈ V := H1 0 (D) such that ∫ D a(x)∇u(x) · ∇v(x)dx = ∫ D f(x)v(x)dx, v ∈ V. (2.2) Existence and uniqueness of a weak solution is guaranteed by the Lax-Milgram lemma if a satisfies a uniform ellipticity condition and if f is an element of the dual space of V , i.e. f ∈ V ∗ = H−1(D). For the FE approximation of u, the infinite-dimensional space V is replaced by a suitable finite dimensional subspace Vh ⊂ V associated to some refinement parameter h > 0. In general, h > 0 is the maximum diameter of a tesselation Kh of D and Vh is the space of all continuous functions consisting of piecewise polynomials with respect to Kh. The (maximum) degree of the polynomials is given by p ∈ N and each basis function of Vh only has a small local support on a few elements of Kh. The discrete version of Problem (2.2) is then to find uh ∈ Vh such that ∫ D a(x)∇uh(x) · ∇vh(x)dx = ∫ D f(x)vh(x)dx, vh ∈ Vh, and the discretization error is bounded provided that u is sufficiently regular: Theorem 2.1.1. [177, Chapter 6.2] Let u be the solution to Problem (2.2) such that u ∈ Hq(D) for some q > 1. Furthermore, let Vh be the FE space containing all piecewise polynomials up to degree p ∈ N and let uh be the FE approximation of u. Then, there is a constant C > 0, independent of u and h, such that ||u− uh||Hm(D) ≤ C‖u‖Hq(D)h min(p+1,q)−m, m ∈ {0, 1}. In Chapters 4–6, we investigate a version of Problem (2.1), where the diffusion co- efficient a is a L2(D)-valued random field. More precisely, a is the sum of a continuous Gaussian part and a stochastic jump field, hence a is discontinuous on D. The corre- sponding weak solution u then is a random function with very low path-wise regularity, in general only u ∈ Hq(D) with q < 3 2 holds almost surely. Therefore, convergence with respect to h is slow and thus reasonably good approximations become computationally expensive. This problem is addressed by introducing a modified FE method, which exploits the higher piecewise regularity of u in between the interfaces of a. We show that this allows us to recover better convergence rates despite the low global regularity 11 2.1. NUMERICAL METHODS FOR PDES of u seems to permit this at first glance. There are of course several other discretization schemes for PDEs, such as fi- nite difference (FD)-, finite volume (FV)-, spectral- and collocation methods (see [14, 51, 81, 143, 144, 177]). These are, however, less suited for PDEs with discon- tinuous coefficients. The FD method is based on uniformly spaced grids, hence an adaptive approach to match the given discontinuities of a as described above is out of reach. Spectral- or collocation methods have in general superior convergence properties compared to FE, but require high regularity and periodicity of the problem. As these assumptions are too restrictive for the setting in Chapters 4–6, it can not be expected that spectral- or collocation approaches outperform the FE method. A discussion of FV methods for hyperbolic PDEs follows shortly. The class of hp-finite element methods, see [189], also needs to be mentioned. In this modification, the FE mesh and/or the polynomial degree is refined/adjusted for certain areas in the domain based on a-posteriori error estimates. This becomes particularly useful if the problem at hand is only irregular at certain areas of the domain, but rather smooth otherwise. Thus, hp-FE methods are actually a promising advancement to the methodology in Chapters 4–6. They are, however, analytically and computationally challenging, and their application to PDEs with discontinuous random coefficients is beyond the scope of this thesis (see also the discussion in Chapter 9). The FE method turns out to be a suitable spatial discretization to solve diffusion- dominated elliptic and parabolic problems. In case of advection-dominated and hy- perbolic PDEs as considered in Chapter 8, however, FE methods become inherently unstable and oscillating. To circumvent these issues, discontinuous Galerkin (DG) methods have been developed as a generalization of the FE method. Introduced in the pioneering works [142, 179] to solve a neutron transport problem, DG methods are, by now, the state-of-the-art approach for the numerical solution of (non-)linear hyperbolic PDEs. Comprehensive overviews on DG methods and their applications can be found for instance in [61, 109, 185]. Although successfully applied to deterministic PDEs, so far, there are not many applications of the DG method to problems in uncertainty quantification , notable exceptions include [90, 110, 120]. The DG approach is also based on the weak formulation of the PDE and basically the same as for the FE method, but the finite-dimensional test function space Vh now consists of all functions built from piecewise polynomials on Kh. This allows especially for discontinuous approximations with more degrees of freedom than their FE counterparts. Since left and right limits of the DG basis differ at the interfaces of Kh, it is necessary to determine a numerical flux at the discontinuities. This is the key feature of DG methods, as the numerical flux can be adjusted to the specific problem to 12 CHAPTER 2. DISCRETIZATION TECHNIQUES IN UQ increase stability of the approximation. Consequently, a DG approach is introduced to discretize the linear hyperbolic SPDE in Chapter 8 and it is shown that this is indeed advantageous to standard FE methods. Another popular discretization scheme for hyperbolic PDEs is the finite volume (FV) method, see e.g. [81, 143]. Therein, the corresponding equation is integrated over small volumes and the divergence theorem is applied to obtain integrals over the surface of each volume. This yields a piecewise constant approximation with discontinuities across the surface of each volume, and consequently a numerical flux similar as for the DG method has to be chosen. In contrast to the FD approach, FV schemes may be applied to irregular domains and non-uniform grids. Moreover, FV methods are conservative if the underlying PDE is a hyperbolic conservation law. Applied to the SPDE in Chapter 8, it turns out the FV approach may be recovered as a subclass of DG methods with piecewise constant basis functions and midpoint quadrature. Hence, the DG approach is more general and thus expected to perform at least as good as FV schemes. To obtain tractable approximations to time-dependent problems, it is still necessary to consider a discretization of the temporal domain. One possibility is to view time as an additional spatial variable and apply a FE approach to the (d + 1)-dimensional space-time domain. This results in a coupling of all spatio-temporal nodal points and the problem has to be solved simultaneously at all discrete points in time. Thus, the resulting scheme has no time-stepping or iterative character, and is not applicable to a large class of time-dependent problems. To this end, it is more practical to apply a FD approach for time integration, either based on forward-, backward or symmetric differences. The corresponding time stepping schemes are the forward respectively backward Euler method and the Crank- Nicolson scheme. In general, the latter approach has the best convergence properties, but requires more regularity of the solution, both with respect to time and space. The backward Euler scheme is unconditionally stable, while pure forward methods often lead to explosions in the numerical solution. In fact, as linear random/stochastic PDEs are considered, the FD approaches can also be interpreted as a DG discretization of the temporal dimension. For a suitably chosen numerical flux, this is equivalent to using piecewise constant functions in the Euler methods and a piecewise linear DG basis in the Crank-Nicolson scheme. The solution to the parabolic problem in Chapters 5 and 6 is smooth in time, but non-differentiable in space. Hence, Crank-Nicolson schemes entail numerical oscilla- tions and the only useful time integrator is the backward Euler method. In Chapter 8, 13 2.2. MONTE CARLO METHODS we consider stochastic integrals acting as source term on the right hand side of a lin- ear transport problem. Due to their construction as the L2-limit of simple integrals (see Section 3.3 in the next chapter), their only reasonable approximation is given by forward differences of the form ∫ ti+1 ti Ψ(s)dL(s) ≈ Ψ(ti)(L(ti+1)− L(ti)), (2.3) where L is a Lévy process and Ψ is an admissible integrand for L. Evaluating Ψ on the right hand side in Eq. (2.3) at any other t ∈ (ti, ti+1] would result in the approximation of a different type of stochastic integral for which the martingale property and the Itô isometry from Theorem 3.3.2 do not hold anymore. On the other hand, applying a forward integration in time to the differential operator is out of the question due to stability reasons. This results in a mixed backward-forward time stepping scheme, where a backward Euler or Crank-Nicolson approach is applied on the left hand side of the PDE and a forward time stepping on the right hand side. As an alternative time stepping method, exponential integrator schemes have suc- cessfully been applied to SPDEs, see [129, 198]. They require however, that the under- lying differential operator is essentially a Laplacian to exploit its spectral basis represen- tation. Clearly, this is not given in the hyperbolic setting in Chapter 8 and it turns out that exponential integrators do not perform superior compared to the other discussed schemes. Neither do higher order temporal approximations such as Runge-Kutta-DG methods([62]) promise any advantages, as they require high spatio-temporal regularity. In Chapter 8, the solution to the SPDE is only mean-square Hölder-continuous in time, at best Lipschitz-continuous in space and the temporal convergence rate is dominated by the approximation of the stochastic integral on the right hand side anyway. For further research, however, Runge-Kutta-DG methods might be useful when moving from linear to nonlinear hyperbolic problems. This is due to the fact that they rely on fully explicit, stable time stepping schemes and thus avoid solving a nonlinear system of equations in every time step. 2.2 Monte Carlo methods The quantities of interest in uncertainty quantification are ultimately statistics of the stochastic model at hand, including moments or quantiles of the solution itself as well as of functionals applied to the solution. This problem can always be reformulated to estimate the mean value of a suitably chosen random variable. Arguably, the most popular and intuitive approach to achieve this is the Monte Carlo method, dating back 14 CHAPTER 2. DISCRETIZATION TECHNIQUES IN UQ to Ulam, von Neumann and Fermi. The historical development of the Monte Carlo method is outlined in [159], a general introduction can be found in [83, 94, 160]. In short, Monte Carlo approaches estimate the expected value of a random variable X by drawing a large number of independent samples from the distribution of X and then taking the arithmetic mean over all samples. To illustrate this technique, let (Ω,F ,P) be a complete probability space and let X : Ω→ R be random variable such that E(|X|) < +∞. For some M ∈ N, independent samples X(1), . . . , X(M) are drawn from the distribution of X and the Monte Carlo estimator of E(X) is defined via EM(X) := 1 M M∑ i=1 X(i). Clearly, E(EM(X)) = E(X), and the strong law of large numbers ensures that lim M→∞ EM(X) = E(X), P-almost surely. Moreover, if X is square-integrable, i.e. E(X2) < +∞, the independence and identical distribution of the samples X(1), . . . , X(M) yield E((EM(X)− E(X))2) = 1 M2 M∑ i,j=1 E(X(i)X(j))− 2E(X) 1 M M∑ i=1 E(X(i)) + E(X)2 = 1 M E(X2) + M(M − 1)− 2M2 +M2 M2 E(X)2 = Var(X) M . In this case, the root-mean-squared error (RMSE) of the Monte Carlo estimator is RMSE := E((EM(X)− E(X))2)1/2 = Var(X)1/2 √ M . (2.4) Monte Carlo methods are easy to implement and straightforward to parallelize, very robust and mean-square-convergence is ensured as soon as X is square-integrable from Eq. (2.4). Their major disadvantage is the inherently slow convergence of order O(M−1/2) for the RMSE. For this reason, estimating E(X) by EM(X) may become prohibitive if the simulation of X is computationally expensive. For example, consider the case X = Φ(u) where u is the solution of a (elliptic) PDE as in Eq. (2.1) with a random diffusion coefficient, i.e. a is a suitable measurable mapping a : Ω → L∞(D) such that a(ω, ·) : D → (0,∞) holds for almost all ω ∈ Ω. Moreover, let Φ : V → R be a given functional on the solution space V = H1 0 (D). 15 2.2. MONTE CARLO METHODS Samples from the exact solution u are in general out of reach, it is, however, possible to obtain approximate samples uh by the FE method, where h > 0 denotes the mesh refinement parameter. Hence, samples of Xh := Φ(uh) are available and E(X) may be estimated via EM(Xh). Assuming that E(‖u‖2 H2(D)) < +∞ and that Φ ∈ V ∗, it follows by Theorem 2.1.1 E((EM(Xh)− E(X))2)1/2 ≤ |E(Xh −X)|+ E((EM(Xh))− E(Xh))2)1/2 ≤ ‖Φ‖V ∗E(‖u− uh‖H1(D)) + Var(Xh)1/2 √ M ≤ C‖Φ‖V ∗ ( E(‖u‖2 H2(D))1/2h+ E(‖u‖2 H2(D))1/2 √ M ) ≤ C ( h+ 1√ M ) . Hence, to obtain a small RMSE, it is necessary to generate many expensive samples (M ≈ h−2) on a grid with small refinement parameter h > 0. To this end, several techniques have been developed to increase efficiency, a particularly effective example is the multilevel Monte Carlo method. Invented by Heinrich in [106] and popularized by Giles in [92] for the pricing of financial derivatives, multilevel Monte Carlo has been applied to a vast range of uncertainty quantification problems in the last decade, see e.g. [28, 29, 54, 66, 120, 194], and Chapters 4 and 6 of this thesis. To illustrate the main idea behind multilevel Monte Carlo estimation, consider again the aforementioned example where X = Φ(u) is the functional of a solution to a random PDE. Rather than a fixed parameter h, we now consider a sequence of decreas- ing refinements h0 > · · · > hL > 0 with L ∈ N (e.g. h` := 2−`−1) and the associated hierarchy of approximated random variables X0 := Φ(uh0), . . . , XL := Φ(uhL). With the convention X−1 := 0, we may express the highest-accuracy approximation XL by the telescopic sum XL = L∑ `=0 X` −X`−1. Using this representation, each correctionX`−X`−1 is estimated by the standard Monte Carlo method using a level-dependent number of M` ∈ N independent samples. This yields the multilevel Monte Carlo estimator EL(XL) := L∑ `=0 EM` (X` −X`−1) = L∑ `=0 1 M` M∑̀ i=1 X (i,`) ` −X(i,`) `−1 (2.5) of E(XL). In Eq. (2.5), it is crucial to actually sample from the distribution of the 16 CHAPTER 2. DISCRETIZATION TECHNIQUES IN UQ correction term X` − X`−1, meaning the samples X(i,`) ` and X (i,`) `−1 must be based on the same set of random variables. This dependency is emphasized by the second superscript `, the sampled corrections across the levels are again independent. Since E(EL(XL)) = E(XL), the RMSE of the multilevel Monte Carlo estimator is bounded by E((EL(XL)− E(X))2)1/2 ≤ C ( hL + ( L∑ `=0 Var(X` −X`−1) M` )1/2) . Given that the variances Var(X` − X`−1) decay sufficiently fast in ` = 0, . . . , L, one may generate many inexpensive samples on the lower levels and, due to the lower variance for large `, it is sufficient to sample only a few expensive corrections on the higher levels. This yields a fast decreasing sequence M0 > · · · > ML and in the example above a RMSE of O(2−L) is achieved with L levels, refinements h` ≈ 2−` and level-dependent numbers of samples decaying from M0 ≈ 22L to ML ≈ L2, as for instance shown in Chapter 4. To obtain a similar RMSE with a singlelevel Monte Carlo estimator, M ≈ 22L samples of the high resolution approximation XL with hL ≈ 2−L are required. Thus, multilevel Monte Carlo approaches can reduce computational time by several orders of magnitudes compared to standard Monte Carlo estimation, from months to days, hours and even to minutes. The precise computational gains from multilevel Monte Carlo are stated under suitable assumptions in the, by now famous, complexity theorem in [92]. Due to their low requirements and good accessibility, multilevel Monte Carlo methods are the algorithms of choice to estimate quantities of interest or moments from PDEs with discontinuous random features. Apart from sampling-based methods, such as Monte Carlo; other approaches to dis- cretize the stochastic domain, e.g. stochastic Galerkin (SG) and stochastic collocation (SC) methods have been developed for uncertainty quantification. Here, the random PDEs are regarded as a class of high-dimensional parametrized PDEs to exploit certain structural properties of the problem. A very prominent example is the elliptic diffusion problem as in Eq. (2.1) with a (strictly) positive random field a : Ω → L∞(D) as coefficient. It can be shown that a admits the Karhunen-Loève expansion a = E(a) + ∑ i∈N ϕiξi (2.6) with respect to a basis (ϕi, i ∈ N) ⊂ L∞(D) and a sequence (ξi, i ∈ N) of centered random variables (see also Theorem 3.1.6 in Chapter 3). The corresponding parameter space is then formed by the range of (ξi, i ∈ N), for instance by [−1, 1]N if the ξi are uniformly distributed on [−1, 1]. Both, SG and SC methods, have successfully been applied to PDEs with random 17 2.2. MONTE CARLO METHODS coefficients, see for instance [14, 25, 77, 90, 165] and the references therein. In par- ticular, the articles [63, 190] need to be mentioned, where the authors show that SG methods may be superior to sampling-based algorithms, provided the underlying prob- lem is sufficiently smooth with respect to its stochastic domain. Translated to the example above, this means essentially that ‖ϕi‖L∞(D) decays fast enough with respect to the index i. Then, the series in Eq. (2.6) may be cut off after only a few terms to obtain a reasonably good representation of a. In the jump-diffusion setting from Chapters 4–6, however, a involves random jumps and hence the decay of ‖ϕi‖L∞(D) is very slow. Consequently, a large number of terms in the truncated Karhunen-Loève expansion of a are necessary and SG/SC methods are prohibitive. Moreover, it is still an open question to find suitable basis functions ϕi and random variables ξi in order to represent discontinuous coefficients as in Eq. (2.6). Another popular evolution of the Monte Carlo algorithm is the class of quasi-Monte Carlo (QMC) methods. In this approach, the pseudo-random samples are replaced by deterministic low discrepancy sequences to achieve a higher convergence rate of the RMSE than O(M−1/2) with respect to the number of samples. The estimation of moments in uncertainty quantification applications is then generally regarded as a quadrature problem in a high-dimensional parameter space and it is possible to recover improved rates of order O(M−1+ε) for arbitrary small ε > 0. A detailed introduction to QMC for integration in high dimensions can be found in [72]. QMC approaches have been applied to random PDEs with diffusion coefficients as in Eq. (2.6) (see [88, 137]) under similar regularity assumptions on the coefficients as SG/SC methods. More recently, however, QMC approaches were also used in the case that a is log- normal, see [96, 98, 108], and a reasonable truncated series-approximation of a involves a large number of terms. This development indicates that QMC techniques may also be promising for random PDEs with discontinuous coefficients, especially since they may be combined with multilevel Monte Carlo as in [108] (see also the discussion in Chapter 9). In addition, QMC methods are based on the evaluation of inverse cumulative density functions and may thus be utilized in the Fourier inversion algorithm from Chapter 7. Naturally, a further topic for research is then the application of QMC methods to SPDEs with Lévy noise as in Chapter 8. 18 3 Lévy processes and stochastic integrals in Hilbert spaces The aim of this chapter is to familiarize the reader with random fields, stochastic processes and -integration in general Hilbert spaces, concepts that are necessary to un- derstand the variety of random objects in the core of this thesis. The first part contains details about covariance operators and spectral expansions, which are recurrently used in all articles of this thesis. Thereafter, some results on Hilbert space-valued Lévy pro- cesses and stochastic integration are collected, providing the foundation for the driving noise processes in Chapters 7 and 8. Thus, the style of this chapter is naturally of more mathematical-formal and contrasts somewhat the preceding one, where the key ideas of each discretization method are roughly sketched and followed by a discussion. The basic concept of Lévy processes in one dimension dates back to Lévy and Khintchine in [127, 145, 146, 147]. Since then, there has been extensive research on Lévy processes, well-known standard works include for instance [6, 40, 187]. The stochastic integrals and the corresponding calculus in Chapter 8 is in the sense of Itô, who introduced a generalization of the Riemann-Stieltjes integral with a Brownian motion as integrator in [118, 119]. For a basic review on Itô calculus with respect to Brownian motions see [125, 168], more general Lévy processes as integrators are discussed in [6]. The stochastic calculus for the corresponding Hilbert-space valued processes is investigated in [67, 158, 173]. Lévy processes and stochastic integration have become two of the most important concepts in probability theory and have been utilized to model various problems in financial mathematics and the natural sciences, examples are given in [17, 67, 94, 173, 188] to name just a few. Throughout this chapter, let T = [0, T ] for T > 0 be a finite time interval and let (Ω,F , (Ft, t ∈ T),P) be a filtered probability space satisfying the usual conditions, i.e. F0 is P-complete and (Ft, t ∈ T) is right-continuous. For any separable metric space (E, ‖ · ‖E), the Borel σ-algebra with respect to E is defined as the smallest σ-algebra containing all open balls in E and denoted by B(E). Furthermore, let (U, (·, ·)U) be a separable Hilbert space, let L(U) be the set of all linear operators on U and let 19 3.1. COVARIANCE OPERATORS AND SPECTRAL EXPANSIONS L+ N(U) ⊂ L(U) be the set of all nonnegative, symmetric, nuclear operators on U . The space of all Bochner-integrable random variables in U with p-th moment is for any p ∈ [1,∞) given by Lp(Ω;U) := {X : Ω→ U is strongly measurable and ∫ Ω ‖X(ω)‖pU < +∞}. 3.1 Covariance operators and spectral expansions As the L2-theory of random fields is investigated in this chapter, it is natural to consider covariance operators and the associated spectral expansions. To this end, the first two moments for square-integrable, U -valued random variables need to be defined. Definition 3.1.1. For any X ∈ L1(Ω;U), the mean of X is defined by the Lebesgue- Bochner integral E(X) := ∫ Ω X(ω)dP(ω) ∈ U. If X ∈ L2(Ω;U), the covariance operator of X is the unique Q ∈ L+ N(U) such that E((X − E(X), φ)U(X − E(X), ψ)U) = (Qφ, ψ)U , φ, ψ ∈ U. The following results ensures that the covariance operator Q ∈ L+ N(U) is well- defined and unique for each square-integrable X: Lemma 3.1.2. For any X ∈ L2(Ω;U), there is a unique covariance operator Q ∈ L+ N(U). Proof. For the existence of Q see [138, Chapter 2]. Now assume that there exist two covariance operators Q, Q̃ ∈ L+ N(U) for X ∈ L2(Ω;U). By Definition 3.1.1 and the linearity of Q, Q̃ it follows ‖Q− Q̃‖2 L(U) = sup φ∈U,φ6=0 ‖(Q− Q̃)φ‖2 U ‖φ‖2 U = sup φ∈U,φ6=0 (Qφ, (Q− Q̃)φ)2 U − (Q̃φ, (Q− Q̃)φ)2 U ‖φ‖2 U = 0, and hence Q = Q̃. By the Hilbert-Schmidt theorem, there is a decreasing sequence of non-negative eigenvalues η1 ≥ η2 ≥ · · · ≥ 0 of Q with zero as only accumulation point and the 20 CHAPTER 3. LÉVY PROCESSES AND STOCHASTIC INTEGRALS corresponding eigenfunctions (ei, i ∈ N) form an orthonormal basis of U . Moreover, Q is nuclear and therefore of trace class, i.e. Tr(Q) := ∑ i∈N (Qei, ei)U = ∑ i∈N ηi < +∞. The square-root of Q is defined via Q1/2φ := ∑ i∈N √ ηi(φ, ei)Uei, φ ∈ U. Since Q1/2 is not necessarily injective, the pseudo-inverse of Q is given by Q−1/2ϕ := φ, if Q1/2φ = ϕ and ‖φ‖U = inf ϕ∈U :Q1/2ϕ=φ {‖ϕ‖U}. Definition 3.1.3. Let X ∈ L2(Ω;U) with covariance operator Q ∈ L+ N(U). Then, the set U := Q1/2(U) equipped with the scalar-product (ϕ1, ϕ2)U := (Q−1/2ϕ1, Q −1/2ϕ2)U , ϕ1, ϕ2 ∈ U , is called the reproducing kernel Hilbert space (RKHS) of X. It turns out to be more convenient to study the RKHS U rather than Q, as restrict- ing U to suitable subspaces does not change the RKHS: Theorem 3.1.4. [173, Theorem 7.4] Let X ∈ L2(Ω;U) with covariance operator Q ∈ L+ N(U) and let the Hilbert space (U, (·, ·)U) be continuously embedded into an- other Hilbert space (Ũ , (·, ·) Ũ ). If Q̃ ∈ L+ N(Ũ) denotes the covariance operator of X considered as Ũ-valued random variable X ∈ L2(Ω; Ũ), then Q1/2(U) = Q̃1/2(Ũ). Theorem 3.1.4 becomes particularly useful if X takes values in a Sobolev space and the corresponding embedding theorems are at hand. Example 3.1.5. In the important case that U = L2(D) for some compact D ⊂ Rd, it can be shown that Q and Q1/2 are integral operators with symmetric, nonnegative definite kernel functions kQ, kQ1/2 : D2 → R, i.e. Qφ = ∫ D kQ(x, ·)φ(x)dx, Q1/2φ = ∫ D kQ1/2(x, ·)φ(x)dx, φ ∈ U, see e.g. [173, Appendix A]. By Mercer’s theorem, the reproducing kernel kQ may be 21 3.1. COVARIANCE OPERATORS AND SPECTRAL EXPANSIONS represented in terms of the eigenbasis of Q by k(x, y) = ∑ i∈N ηiei(x)ei(y), x, y ∈ D, and hence kQ(x, ·) = kQ(·, x) ∈ U . Since any ϕ ∈ U is of the form ϕ = Q1/2φ = ∑ i∈N √ ηi(φ, ei)Uei for some φ ∈ U , Mercer’s theorem yields that kQ satisfies for any x ∈ D (ϕ, kQ(x, ·))U = ∑ i∈N (Q1/2φ, ηiei)Uei(x) = ∑ i∈N (Q1/2φ,Q1/2ei)U √ ηiei(x) = ∑ i∈N √ ηi(φ, ei)Uei(x) = ϕ(x), also called the reproducing kernel property. Moreover, the RKHS U is given by U = {φ ∈ L2(D)| ∑ i∈N, ηi 6=0 (φ, ei)2 U ηi < +∞}, (φ, ψ)U = ∑ i∈N, ηi 6=0 (φ, ei)U(ψ, ei)U ηi . To conclude this section, the following useful spectral representation ofX is recorded. Theorem 3.1.6 (Karhunen-Loève expansion). Let X ∈ L2(Ω;U) with covariance op- erator Q and denote by ((ηi, ei), i ∈ N) the (ordered) eigenpairs of Q. Then, X admits the expansion X = E(X) + ∑ i∈N ξiei, where the ξi are real-valued, centered random variables with Cov(ξi, ξj) = ηiδij for i, j ∈ N. The above series converges in L2(Ω;U) with truncation error bounded by ‖ ∞∑ i>N ξiei‖2 L2(Ω;U) := E ( ‖ ∞∑ i>N ξiei‖2 U ) = ∑ i>N ηi, N ∈ N. (3.1) Proof. As (ei, i ∈ N) is an orthonormal basis in U , X − E(X) is expanded via X − E(X) = ∑ i∈N ((X − E(X)), ei)Uei. 22 CHAPTER 3. LÉVY PROCESSES AND STOCHASTIC INTEGRALS Now define ξi := ((X−E(X)), ei)U . It is immediate that the ξi are centered and satisfy Cov(ξi, ξj) = E((X − E(X)), ei)U(X − E(X)), ej)U) = (Qei, ej)U = ηiδij. To show the convergence in L2(Ω;U), let XN := ∑N i=1 ξiei for any N ∈ N and note that E(‖XN‖2 U) = E ( ‖ N∑ i=1 ξiei‖2 U ) = N∑ i,j=1 E((ξiei, ξjej)U) = N∑ i=1 ηi. Since Tr(Q) = ∑ i∈N ηi < +∞, XN is a L2(Ω;U)-Cauchy-sequence and thus converges to the limit X − E(X) by the completeness of L2(Ω;U). Consequently, E(‖X −XN‖2 U) = E ( ‖ ∞∑ i>N ξiei‖2 U ) = ∑ i>N ηi. Example 3.1.7. If X in Theorem 3.1.6 is Gaussian with mean µ ∈ U and covariance operator Q ∈ L+ N(U), it holds that X = µ+ ∑ i∈N √ ηieiZi, (3.2) where (Zi, i ∈ N) is a sequence of independent, one-dimensional standard normally distributed random variables (see for instance [3, Chapter 3]). The representation in Eq. (3.2) becomes particularly useful for the simulation ofX: the sum may be truncated after a finite number of terms and one then only needs to sample from a one-dimensional standard normal distribution. For an introduction of Gaussian measures on Banach- and Hilbert spaces, see [138] or [173, Chapter 3.5]. Remark 3.1.8. Note that the independence for the random variables in the Karhunen- Loève expansion only holds for the Gaussian case in Example 3.1.7. For general X ∈ L2(Ω;U) as in Theorem 3.1.6, the sequence (ξ, i ∈ N) consists of uncorrelated, but not independent random variables (see also Proposition 3.2.9). As Q is of trace-class, the truncation error in Eq. (3.1) depends on the decay of the eigenvalues and can be made arbitrary small for sufficiently large N . A similar representation as in Eq. (3.2) is also given for Lévy processes in the next subsection, which is the fundamental observation for the approximation algorithm introduced in Chapter 7. 23 3.2. LÉVY RANDOM FIELDS 3.2 Lévy random fields Hilbert space-valued Lévy fields form the basis of the driving noise processes in Chap- ters 7 and 8. In this section, the class of Lévy processes is introduced and a fundamental decomposition theorem in terms of Gaussian and compound Poisson processes is given. Definition 3.2.1. A U -valued stochastic process L = (L(t), t ∈ T) is called Lévy process if • L has stationary and independent increments, • L(0) = 0 P-a.s. and • L is stochastically continuous, i.e. for all ε > 0 and t ∈ T it holds lim s→t, s∈T P(||L(t)− L(s)||U > ε) = 0. To obtain a clear distinction from finite-dimensional Lévy processes, L is sometimes called Lévy field if dim(U) = +∞. Lévy processes have in general discontinuous trajectories, regardless of the dimensionality of U . It can be shown, however, that every Lévy process has a modification with càdlàg paths ([173, Theorem 4.3]). The only exceptional case with continuous trajectories are Hilbert-space valued Brownian motions or Wiener processes: Definition 3.2.2. A zero mean Lévy process W on U with continuous trajectories is called a Wiener process. As subsequently outlined, every Lévy process is the sum of a continuous Wiener process and certain discontinuous jump processes. To this end, the following basic results for Hilbert space-valued Wiener processes are recorded. Theorem 3.2.3. [67, Proposition 4.1] A Wiener process W on U is a Gaussian, square-integrable process and W (t) ∼ N (0, tQ), where Q is the covariance operator of W (1). Moreover, if ((ηi, ei), i ∈ N) are the orthonormal eigenpairs of Q, then W admits the series expansion W (t) = ∑ i∈N (W (t), ei)Uei, t ∈ T, where ((W (t), ei)U , t ∈ T) are (scaled) independent, real-valued Wiener processes with covariance Cov((W (t), ei)U , (W (s), ei)U) = min(t, s)ηi. 24 CHAPTER 3. LÉVY PROCESSES AND STOCHASTIC INTEGRALS Since W is centered Gaussian, the characteristic function of W reads E(ei(W (t),φ)U ) = exp ( − t 2(Qφ, φ)U ) , φ ∈ U. (3.3) There is also a U -valued version of a compound Poisson process, which will be the other building block for the Lévy process. Definition 3.2.4. Let ν̃ be a finite measure on (U,B(U)) with ν̃({0}) = 0. A compound Poisson process with jump intensity measure ν̃ on U is a càdlàg Lévy process L on U satisfying P(L(t) ∈ Û) = e−ν̃(U)t∑ k∈N tk k! ν̃ ∗k(Û), Û ∈ B(U), t ∈ T. Above, ν̃∗k denotes the k-fold convolution of ν̃, given by ν̃∗k(Û) := (ν̃ ∗ · · · ∗ ν̃)(Û), (ν̃ ∗ ν̃)(Û) := ∫ U ν̃(Û − φ)ν̃(dφ). Theorem 3.2.5. [173, Chapter 4.3] Let L be a compound Poisson process on U with jump intensity measure ν̃. 1. The characteristic function of L is given by E(ei(L(t),φ)U ) = exp ( − t ∫ U 1− ei(ψ,φ)U ν̃(dψ) ) , φ ∈ U, t ∈ T. 2. There is a sequence of i.i.d. U-valued random variables (Ξ, i ∈ N) with law Ξi ∼ ν̃ ν̃(U) and a Poisson process N = (N(t), t ∈ T) with intensity ν̃(U), independent of (Ξ, i ∈ N), such that L(t) = N(t)∑ i=1 Ξi. 3. The compound Poisson process L is integrable if and only if ∫ U ‖φ‖U ν̃(dφ) < +∞. If L is integrable, E(L(t)) = t ∫ U φν̃(dφ), t ∈ T, and (L(t)− E(L(t)), t ∈ T) is called compensated compound Poisson process. With the notions of Wiener processes and compound Poisson processes on U , the most important representation result for Lévy processes, known as the Lévy-Khintchine decomposition, can now be stated. 25 3.2. LÉVY RANDOM FIELDS Theorem 3.2.6. [173, Theorem 4.23 (Lévy-Khintchine decomposition)] Every Lévy process L has the representation L(t) = µ0t+W (t) + ∑ k∈N LAk(t)− t ∫ Ak φν(dφ) + LA0(t), t ∈ T, (3.4) where • µ0 ∈ U is a deterministic mean function, • W is a Wiener process with covariance operator Q ∈ L+ N(U), • Ak := {φ ∈ U | rk ≤ ‖φ‖U ≤ rk−1} and A0 := {φ ∈ U | ‖φ‖U ≥ r0} for an arbitrary sequence (rk, k ∈ N) ⊂ R>0, decreasing monotonously to zero, • ν is a measure on (U,B(U)) such that ∫ U min(‖y‖2 U , 1)dν(y) < +∞ (also called Lévy measure) and (LAk , k ∈ N) are compound Poisson processes with finite intensity measures ν̃k := ν|Ak , and • all summands in Eq. (3.4) are independent stochastic processes on U . Remark 3.2.7. The Lévy-Khintchine decomposition states that every Lévy process is the sum of a Wiener process W with drift µ0, a compound Poisson process and a superposition of compensated compound Poisson processes. It needs to be emphasized that changing the cutoff treshold r0 in the Lévy-Khintchine decomposition affects the mean µ0. Hence, the above representation depends on r0 and is a-priori not unique. As seen in Theorem 3.2.8 below, uniqueness can be achieved by fixing r0 = 1, the cor- responding representation is also referred to as Lévy-Itô decomposition. If ν(U) < +∞, then L may be represented by a finite number of (compensated) compound Poisson processes. In this case, L has almost surely finitely many jumps in every bounded interval in T and L is said to be of finite activity. Due to the Lévy-Khintchine decomposition, the following is immediate: Theorem 3.2.8 (Lévy-Khintchine formula). Let L be a Lévy process on U . Then, the characteristic function of L(t) is given by E[exp(i(φ, L(t))U)] = exp(tΨL(φ)), φ ∈ U, t ∈ T, with characteristic exponent ΨL(φ) = i(µ, φ)U − 1 2(Qφ, φ)U + ∫ U ei(φ,ψ)U − 1− i(φ, ψ)U1||ψ||U<1(ψ)ν(dψ). 26 CHAPTER 3. LÉVY PROCESSES AND STOCHASTIC INTEGRALS Proof. Due to the independence of all terms in Eq. (3.4), the characteristic exponent is of additive form. The claim follows with the characteristic functions of the Wiener process and the compound Poisson processes from Eq. (3.3) and Theorem 3.2.5, re- spectively. Using r0 = 1 in Theorem 3.2.6 yields the indicator function 1||ψ||U<1 in ΨL. For a more detailed proof, see e.g. [173, Chapter 4] and the references therein. Compared with the Lévy-Khintchine decomposition, the Lévy-Khintchine formula assumes a cutoff at r0 = 1 for the non-compensated Poisson process. Hence, the above representation is unique and L is uniquely determined by (µ,Q, ν), also known as the characteristic triplet. Similar to the Wiener process, L may be expanded in terms of an orthonormal basis of U and one-dimensional Lévy processes. Proposition 3.2.9. [173, Theorem 4.39] If (φi, i ∈ N) is an orthonormal basis of U , then L has the series expansion L(t) = ∑ i∈N (L(t), φi)Uφi, t ∈ T, and the processes (L(t), φi)U are real-valued Lévy processes with càdlàg paths. The above series converges in probability uniformly in T. If L is square-integrable, the processes (L(t), φi)U are uncorrelated and indepen- dence only holds in the case that L = W is a Wiener processes on U with the marginals (L(t), φi)U becoming real-valued Brownian motions. To obtain an expansion with re- spect to the eigenbasis of the covariance operator of L, it is essential that L is actually square-integrable. As seen in the next result, the integrability of L depends solely on the moments of ν. Theorem 3.2.10. Let L be a U-valued Lévy process with characteristic triplet (µ,Q, ν). For any t ∈ T and n ∈ N it holds that E(‖L(t)‖nU) < +∞ ⇔ ∫ φ∈U,‖φ‖≥1 ‖φ‖nUν(dφ) < +∞. Moreover, if L is square-integrable, there exits µL ∈ U and QL ∈ L+ N(U) such that for any s, t ∈ T and φ, ψ ∈ U : E((L(t), φ)U) = t(µL, φ)U E((L(t)− tµ, φ)U(L(s)− sµ, ψ)U) = min(t, s)(QLφ, ψ)U E(‖L(t)− tµL‖2 U) = tTr(QL). 27 3.3. STOCHASTIC INTEGRATION Proof. The first part of the claim follows from the Lévy-Khintchine decomposition in the same fashion as for finite-dimensional Lévy processes, see e.g. [6, Theorem 2.5.2]. The second part is given in [173, Theorem 4.44]. Remark 3.2.11. In general, µL and QL above are not identical to µ resp. Q from the Lévy triplet, but merely µL = µ+µν and QL = Q+Qν , with µν , Qν stemming from the jump intensity measure ν. However, as QL has an orthonormal eigenbasis (ei, i ∈ N) with corresponding eigenvalues ηi ≥ 0, L has the series expansion L(t) = ∑ i∈N (L(t), ei)Uei, t ∈ T, where the one-dimensional marginal Lévy processes (L(t), ei)U are uncorrelated and E((L(t), ei)U) = t(µL, ei)U , and Var((L(t), ei)U) = tηi. 3.3 Stochastic integration This section provides a construction of stochastic integrals in Itô’s sense with square- integrable Lévy processes on U as integrator. The fundamental property of the corre- sponding Itô integrals, an infinite-dimensional version of the Itô isometry, is stated as the central result at the end of this section. For completeness, it needs to be mentioned that the class of integrators may be generalized to the space of all square-integrable, U -valued martingales, see [173, Chapter 8]. However, as only square-integrable Lévy processes have been considered in Chapters 7 and 8, this section is restricted to the case that L is a centered and square-integrable Lévy process on U for the sake of brevity. Let QL be the covariance operator of L, and observe that by Theorem 3.2.10 the processes L = (L(t), t ∈ T) and (‖L(t)‖2 U − t tr(QL), t ∈ T) are U - and real-valued martingales, respectively. Furthermore, let U := Q 1/2 L (U) be the RKHS associated to L and let (H, (·, ·)H) be another separable Hilbert space, not necessarily identical to U . The set of all linear operators from U to H is given by L(U,H) and LHS(U , H) denotes the space of all Hilbert-Schmidt operators from U to H. The construction of the stochastic integral is based on the following class of piece- wise constant, L(U,H)-valued integrands. Definition 3.3.1. A L(U,H)-valued process Ψ : T × Ω → L(U,H) is called a simple process if there are t0, . . . , tm ∈ T, operators Ψj ∈ L(U,H) and events Aj ∈ Ftj for 28 CHAPTER 3. LÉVY PROCESSES AND STOCHASTIC INTEGRALS j = 0, . . . ,m− 1 such that Ψ(t) = m−1∑ j=0 1Aj1[tj ,tj+1)(t)Ψj, t ∈ T. The set of all L(U,H)-valued simple processes is denoted by S(U,H). For any simple process Ψ ∈ S(U,H), define the stochastic integral IΨ(t) = ∫ t 0 Ψ(s)dL(s) := m−1∑ j=1 1AjΨj(L(min(t, tj+1))− L(min(t, tj))). It is straightforward to see that the process (IΨ(t), t ∈ T) is a H-valued, square- integrable martingale with zero mean (as L is centered). Since L has independent increments, there holds the Itô isometry for simple integrals given by ‖Ψ‖2 L,t : = E ( ‖ ∫ t 0 Ψ(s)dL(s)‖2 H ) = ∫ t 0 E(‖Ψ(s)‖2 LHS(U ,H))ds for any t ∈ T. Note that ‖ ·‖L,T defines a seminorm on S(U,H) and let L2 L,T (H) be the completion of (S(U,H), ‖ · ‖L,T ). With this at hand, the class of admissible integrands is naturally extended to L2 L,T via IΨ(t) = ∫ t 0 Ψ(s)dL(s) := lim n→∞ ∫ t 0 Ψn(s)dL(s), Ψ ∈ L2 L,T (H), where (Ψn, n ∈ N) ⊂ S(U,H) is an approximating sequence of Ψ and the limit is taken in the L2(Ω;H)-sense. To conclude this chapter, L2 L,T (H) is characterized as the space of all predictable, LHS(U , H)-valued processes and the Itô isometry is generalized to the space of all admissible integrands L2 L,T (H): Theorem 3.3.2. [173, Chapter 8.6] Let L be a square-integrable, zero-mean Lévy pro- cess on U with covariance operator QL and RKHS given by U = Q 1/2 L (U). Then, the space of admissible integrands for L is L2 L,T (H) = {Ψ : T× Ω→ LHS(U , H) |Ψ is square-integrable and predictable}. Moreover, for any Ψ ∈ L2 L,T (H), the process (IΨ(t), t ∈ T ) is a centered, square- integrable martingale, and for t ∈ T it holds that E ( ‖ ∫ t 0 Ψ(s)dL(s)‖2 H ) = ∫ t 0 E(‖Ψ(s)Q1/2 L ‖2 LHS(U,H))ds = ∫ t 0 E(‖Ψ(s)‖2 LHS(U ,H))ds. 29 Declaration to the cumulative part The central part consists of the five research articles in Chapters 4–8 and the publica- tion/review status of each article is: • Andrea Barth and Andreas Stein, "A study of elliptic partial differential equa- tions with jump diffusion coefficients": published in ASA/SIAM Journal on Un- certainty Quantification, 2018, SIAM, volume 6, issue 4, pp. 1707–1742. • Andrea Barth and Andreas Stein, "Numerical analysis for time-dependent advec- tion-diffusion problems with random discontinuous coefficients": submitted to Stochastics and Partial Differential Equations: Analysis and Computations July 2020, currently in the first stage of review, preprint available on ArXiv. • Andrea Barth and Andreas Stein, "A multilevel Monte Carlo algorithm for parabolic advection-diffusion problems with discontinuous coefficients": published in Proceedings of the 13th International Conference in Monte Carlo & Quasi- Monte Carlo Methods in Scientific Computing, 2020, Springer. • Andrea Barth and Andreas Stein, "Approximation and simulation of infinite- dimensional Lévy processes": published in Stochastics and Partial Differential Equations: Analysis and Computations, 2018, Springer, volume 6, issue 2, pp. 286–334. • Andrea Barth and Andreas Stein, "A stochastic transport problem with Lévy noise: Fully discrete numerical approximation": submitted to BIT Numerical Mathematics in October 2019, currently in the second stage of review, preprint available on ArXiv. I hereby declare that the list of co-authors is exhaustive and that I have not repro- duced, without acknowledgement, the work of another. I have majorly contributed to all research articles above, including particularly the phase of problem selection, litera- ture research, the derivation of theoretical results, conduct and evaluation of numerical experiments as well as the writing process. 31 DECLARATION TO THE CUMULATIVE PART Compared to the published/submitted version of the articles, I have modified Chap- ters 4–8 in the following way: • Some of the figures have been generated newly to obtain a better resolution and a consistent style throughout this thesis. They are, however, based on the exact same data that has been used in the corresponding research article. • As far as possible, notation has been homogenized throughout Chapters 4–8. • A short comment on the results has been added at the end of Chapters 4 and 6. • All articles have been reformatted in the style of this thesis. • A few typos and grammatical errors have been corrected. • The bibliography for all articles has been unified and is given at the end of this thesis. Except for this minor alterations, all articles in Chapters 4–8 are reproductions of their corresponding published/submitted research papers. Andreas Stein Stuttgart, August 2020 32 4 A study of elliptic partial differential equations with jump diffusion coefficients Andreas Stein and Andrea Barth First published in "ASA/SIAM Journal on Uncertainty Quantification", 2018, SIAM, volume 6, issue 4, pp. 1707–1742. URL: https: // epubs. siam. org/ doi/ abs/ 10. 1137/ 17M1148888 . Copyright c©by SIAM and ASA. Unauthorized reproduction of this article is prohibited. Abstract: As a simplified model for subsurface flows elliptic equations may be utilized. Insufficient measurements or uncertainty in those are commonly modeled by a random coefficient, which then accounts for the uncertain permeability of a given medium. As an extension of this methodology to flows in heterogeneous, fractured or porous media, we incorporate jumps in the diffusion coefficient. These discontinuities then represent transitions in the media. More precisely, we consider a second order elliptic problem where the random coefficient is given by the sum of a (continuous) Gaussian random field and a (discontinuous) jump part. To estimate moments of the solution to the resulting random partial differential equation, we use a path-wise numerical approximation combined with multilevel Monte Carlo sampling. In order to account for the discontinuities and improve the convergence of the path-wise ap- proximation, the spatial domain is decomposed with respect to the jump positions in each sample, leading to path-dependent grids. Hence, it is not possible to create a sequence of grids which is suitable for each sample path a-priori. We address this issue by an adaptive multilevel algorithm, where the discretization on each level is sample-dependent and fulfills given refinement conditions. 33 https://epubs.siam.org/doi/abs/10.1137/17M1148888 4.1. INTRODUCTION 4.1 Introduction Uncertainty quantification plays an increasingly important role in a wide range of prob- lems in the Engineering Sciences and Physics. Examples of sources of uncertainty are imprecise or insufficient measurements and noisy data. In the underlying dynami- cal system this is modeled via a stochastic operator, stochastic boundary conditions and/or stochastic data. As an example, to model subsurface flow more realistically the coefficient of an (essentially) elliptic equation is assumed to be stochastic. A common approach in the literature is to use (spatially) correlated random fields that are built from uniform distributions or colored log-normal fields. The resulting marginal dis- tributions of the field are (shifted) normally, resp. log-normally distributed. Neither choice is universal enough to accommodate all possible types of permeability, espe- cially not if fractures are incorporated (see [202]), the medium is very heterogeneous or porous. The last decade has been an active research period on elliptic equations with random data. A non-exhaustive list of publications in this field includes [1, 14, 15, 29, 59, 63, 85, 148, 165, 190, 194]. One can find various ways to approximate the distribution or moments of the solution to the elliptic equation. Next to classical Monte Carlo methods, their multilevel variants and other variance reduction techniques have been applied. The concept of multilevel Monte Carlo simulation has been developed in [106] to calculate parametric integrals and has been rediscovered in [92] to estimate the expected value of functionals of stochastic differential equations. Ever since, multilevel Monte Carlo techniques have been successfully applied to various problems, for instance in the context of elliptic random PDEs in [1, 29, 59, 148, 194] to just name a few. These sampling algorithms are fundamentally different from approaches using Polynomial Chaos. The latter suffer from the fact that one is rather restricted when it comes to possibilities to model the stochastic coefficient. While in the case of fields built from uniform distributions or colored log-normal fields these algorithms can outperform sampling strategies, approaches – like stochastic Galerkin methods – are less promising in our discontinuous setting due to the rather involved structure of the coefficient. In fact, it is even an open problem to define them in the case that a Lévy-field is used. Our main objective in this paper is to show existence and uniqueness of the solution to the elliptic equation when the coefficient is modeled as a jump-diffusion. By that we mean a field which consists of a deterministic, a Gaussian and non-continuous part. As we show in the numerical examples, this jump-diffusion coefficient can be used to model a wide array of scenarios. This generalizes the work in [148] and uses partly [113]. To approximate the expectation of the solution we develop and test, further, variants of the 34 CHAPTER 4. ELLIPTIC PDES WITH JUMP DIFFUSION COEFFICIENT multilevel Monte Carlo method which are tailored to our problem: namely adaptive 1 and coupled multilevel Monte Carlo methods. Adaptivity is needed in the jump- diffusion setting, since the coefficient is not continuous. Our analysis shows that the non-adaptivity in a multilevel setting with non-continuous coefficients entails a larger error than when an adaptive algorithm is used. This result is not surprising, since the essence of the multilevel algorithm is that many samples are calculated on coarse grids, where the distributional error from misjudging jump-locations is high. Adaptivity, however, comes to the price that in certain scenarios solving the underlying system of equation becomes computationally expensive. In these settings the advantageous time-to-error performance of adaptive methods may be worse. As a simplified version of Multifidelity Monte Carlo sampling (see [171]), we introduce a coupled multilevel Monte Carlo estimator that reuses samples across levels and is preferred when sampling from a certain distribution is computationally expensive. The coupled algorithms outperform, in general, algorithms with a standard sampling strategy of multilevel Monte Carlo, as it actually reduces the mean square error. In Section 4.2 we introduce the model problem, define path-wise weak solutions of random partial differential equations (PDEs) and show almost sure existence and uniqueness under relatively weak assumptions on the model parameters. The main contribution of this section is the existence and uniqueness result in Theorem 4.2.5, which is then readily transferred to the special case of a jump-diffusion problem in the subsequent sections of this article. In Section 4.3 we define the jump-diffusion coefficient and construct suitable approximations. Both stochastic parts of the jump- diffusion coefficient are approximated: The Gaussian one by a standard truncation of its basis representation, and the jump part (if direct sampling from the jump distribution is not possible) by a technique based on Fourier Inversion. We show Lp-type convergence for all existing moments of the approximation. From this result convergence of the approximated solution follows immediately. The approximated solution has still to be discretized to actually estimate moments of it (Section 4.4). The Galerkin-type discretization is directly furthered into an adaptive scheme. Section 4.5 then introduces the sampling methods that are used, namely Monte Carlo, multilevel Monte Carlo and a coupled variant of it. An extensive discussion of numerical examples in one and two dimensions, in Section 4.6, concludes the paper. 1Throughout this chapter, the term "adaptive" refers to finite element grids that are aligned a- priori to the discontinuities of the diffusion coefficient. This is contrast to the "adaptive finite element method" which is based on a-posteriori error estimates and in general includes several stages of mesh refinement. To avoid confusion, the terminology has been changed to "sample-adapted grids" instead of "adaptive grids" in the next two chapters. For this chapter, however, I decided to keep the terminology according to the original article and inserted this footnote to clarify. 35 4.2. ELLIPTIC BOUNDARY VALUE PROBLEMS 4.2 Elliptic boundary value problems and existence of solutions We consider the following random elliptic equation in a general setting before we specify in Section 4.3 our precise choice of coefficient function. Let (Ω,F ,P) be a complete probability space and D ⊂ Rd, for some d ∈ N, be a bounded and connected Lipschitz domain. In this paper we consider the linear, random elliptic problem −∇ · (a(ω, x)∇u(ω, x)) = f(ω, x) in Ω×D, (4.1) where a : Ω × D → R is a stochastic jump-diffusion coefficient and f : Ω × D → R is a random source function. The Lipschitz boundary ∂D consists of open (d − 1)- dimensional manifolds which are grouped into two disjoint subsets Γ1 and Γ2 such that Γ1 6= ∅ and ∂D = Γ1 ∪ Γ2. We impose mixed Dirichlet-Neumann boundary conditions u(ω, x) = 0 on Ω× Γ1, a(ω, x) #»n · ∇u(ω, x) = g(ω, x) on Ω× Γ2, (4.2) on Eq. (4.1), where #»n is the outward unit normal vector to Γ2 and g : Ω × Γ2 → R, assuming that the exterior normal derivative #»n · ∇u on Γ2 is well-defined for any u ∈ C1(D). To obtain a path-wise variational formulation of this problem, we use the standard Sobolev space H1(D) equipped with the norm ‖v‖H1(D) = (∫ D |v|2 + ‖∇v‖2 2dx )1/2 for v ∈ H1(D), where ‖ · ‖2 denotes the Euclidean norm on Rd. On the Lipschitz domain D, the existence of a bounded, linear operator γ : H1(D)→ H1/2(∂D) with γ : H1(D) ∩ C∞(D)→ H1/2(∂D), v 7→ v|∂D and ‖γv‖H1/2(∂D) ≤ CD‖v‖H1(D) (4.3) for v ∈ H1(D) and some constant CD > 0, dependent only on D, is ensured by the trace theorem, see for example [73]. At this point, one might argue that the trace operator γ needs to be defined path-wise for any ω ∈ Ω, since the Neumann part of the boundary conditions in Eq. (4.2) may contain a random function g : Ω× Γ2 → R. 36 CHAPTER 4. ELLIPTIC PDES WITH JUMP DIFFUSION COEFFICIENT This is true if one works with γ on Γ2 ⊂ ∂D, as the trace γv then has to match the boundary condition given by g(ω, ·) on Γ2 for P-almost all ω ∈ Ω and v ∈ H1(D). In our case, for simplicity, we may treat γ independently of ω ∈ Ω, since we only consider the trace operator on the homogeneous boundary part Γ1 to define V as follows: The subspace of H1(D) with zero trace on Γ1 is then V := {v ∈ H1(D)| γv|Γ1 = 0}, with norm ‖v‖V := (∫ D |v|2 + ‖∇v(x)‖2 2dx )1/2 . Remark 4.2.1. The condition Γ1 6= ∅ implies that V is a closed linear subspace of H1(D). We may as well work with non-homogeneous boundary conditions on the Dirichlet part, i.e. u(ω, x) = g1(ω, x) for g1 : Ω × Γ1 → R. The corresponding trace operator γ is still well defined if, for P-a.e. ω ∈ Ω, g1(ω, ·) can be extended to a function g̃1(ω, ·) ∈ H1(D). Then, we consider for P-a.e. ω ∈ Ω the problem −∇ · (a(ω, x)∇((u− g̃1)(ω, x))) = f +∇ · (a(ω, x)∇g̃1(ω, x)) on D, (u− g̃1)(ω, x) = 0 on Γ1 and a(ω, x) #»n · ∇((u− g̃1)(ω, x)) = g(ω, x)− a(ω, x) #»n · ∇g̃1(ω, x) on Γ2. But this is in fact a version of Problem (4.1) equipped with Eqs. (4.2) where the source term and Neumann-data have been changed (see also [80, p. 317]). As the coefficient and the boundary conditions are given by random functions, the solution u is also a random function. Besides path-wise properties, u may also have certain integrability properties with respect to the underlying probability measure. To this end, we introduce the space of Bochner integrable random variables resp. random functions (see [67] for an overview). Definition 4.2.2. Let (X , ‖ · ‖X ) be a Banach space and define the norm ‖ · ‖Lp(Ω;X ) for a X -valued random variable ϕ : Ω→ X as ‖ϕ‖Lp(Ω,X ) := E(‖ϕ‖pX )1/p for 1 ≤ p < +∞ esssupω∈Ω‖ϕ‖X for p = +∞ . The corresponding space of Bochner-integrable random variables is then given by Lp(Ω;X ) := {ϕ : Ω→ X is strongly measurable and ‖ϕ‖Lp(Ω;X ) < +∞}. 37 4.2. ELLIPTIC BOUNDARY VALUE PROBLEMS The following set of assumptions on a, f and g allows us to show existence and uniqueness of the solution to Eq. (4.1). Consequently, we denote by V ′ the topological dual of a vector space V . Assumption 4.2.3. Let H := L2(D). For P-almost all ω ∈ Ω it holds that: • a−(ω) := infx∈D a(ω, x) > 0 and a+(ω) := supx∈D a(ω, x) < +∞. • 1/a− ∈ Lp(Ω;R), f ∈ Lq(Ω;H) and g ∈ Lq(Ω;L2(Γ2)) for some p, q ∈ [1,∞] such that r := (1/p+ 1/q)−1 ≥ 1. Remark 4.2.4. In Assumption 4.2.3, we did not establish a uniform elliptic bound on a in Ω (see for example [148]), neither did we assume a certain spatial regularity as in [54, 194]. The relatively weak assumptions are natural in our context, since in Section 4.3 we model a as a jump-diffusion coefficient and uniform bounds or assumptions on Hölder- continuity are too restrictive. For the investigation of problem (4.1) with piecewise Hölder-continuous coefficients we refer to [194]. We may identify H with its dual and work on the Gelfand triplet V ⊂ H ' H ′ ⊂ V ′. Hence, Assumption 4.2.3 guarantees that f(ω, ·) ∈ V ′, and, similarly, g(ω, ·) ∈ H−1/2(Γ2) for P-almost all ω ∈ Ω. For fixed ω ∈ Ω, multiplying the random PDE (4.1) with a test function v ∈ V and integrating by parts yields the integral equation ∫ D a(ω, x)∇u(ω, x) · ∇v(x)dx = ∫ D f(ω, x)v(x)dx+ ∫ Γ2 g(ω, x)[Tv](x)dx. (4.4) Consider the bilinear form Ba(ω) Ba(ω) : V × V → R, (u, v) 7→ ∫ D a(ω, x)∇u(x) · ∇v(x)dx and Fω : V → R, v 7→ ∫ D f(ω, x)v(x)dx+ ∫ Γ2 g(ω, x)[Tv](x)dx, where the integrals in Fω are understood as the duality pairings ∫ D f(ω, x)v(x)dx = V ′〈f(ω, ·), v〉V and∫ Γ2 g(ω, x)[Tv](x)dx = H−1/2(Γ2)〈g(ω, ·), T v〉H1/2(Γ2). Equation (4.4) then leads to the path-wise variational formulation of Problem (4.1): For P-almost all ω ∈ Ω, given f(ω, ·) ∈ V ′ and g(ω, ·) ∈ H−1/2(Γ2), find u(ω, ·) ∈ V such that Ba(ω)(u(ω, ·), v) = Fω(v) (4.5) 38 CHAPTER 4. ELLIPTIC PDES WITH JUMP DIFFUSION COEFFICIENT for all v ∈ V . A function u(ω, ·) ∈ V that fulfills the path-wise variational formulation is then called path-wise weak solution to Problem (4.1). Theorem 4.2.5. If Assumption 4.2.3 holds, then there exists a unique path-wise weak solution u(ω, ·) ∈ V to Problem (4.5) for P-almost all ω ∈ Ω. Furthermore, u ∈ Lr(Ω;V ) and ‖u‖Lr(Ω;V ) ≤ C(a−,D, p)(‖f‖Lq(Ω;H) + ‖g‖Lq(Ω;L2(Γ2))), where C(a−,D, p) > 0 is a constant depending only on the indicated parameters. Proof. Choose ω ∈ Ω such that Assumption 4.2.3 is fulfilled. For all u, v ∈ V , we obtain by the Cauchy-Schwarz inequality |Ba(ω)(u, v)| ≤ (∫ D (a(ω, x))2‖∇u(x)‖2 2dx ∫ D ‖∇v(x)‖2 2dx )1/2 ≤ a+(ω)‖u‖V ‖v‖V . On the other hand, Ba(ω)(u, u) ≥ a−(ω) ∫ D ‖∇u(x)‖2 2dx = a−(ω) 2 (‖∇u‖2 L2(D) + ‖∇u‖2 L2(D)) ≥ a−(ω) 2 (‖∇u‖2 L2(D) + C−2 |D|‖u‖ 2 L2(D)) ≥ a−(ω) 2 min(1, C−2 |D|)‖u‖ 2 V , where C2 |D| > 0 stems from the constant in Poincaré’s inequality, C|D|, and only depends on |D|. Hence the bilinear form Ba(ω) : V ×V → R is continuous and coercive. We use that H ' H ′ ⊂ V ′ and the trace theorem (Equation (4.3)) to bound Fω by Fω(v) ≤ ‖f(ω, ·)‖V ′‖v‖V + ‖g(ω, ·)‖H−1/2(Γ2)‖Tv‖H1/2(Γ2) ≤ (‖f(ω, ·)‖H + CD‖g(ω, ·)‖L2(Γ2))‖v‖V This shows that Fω is a bounded linear functional on V (and therefore continuous) for almost all ω ∈ Ω. The existence of a unique path-wise weak solution u(ω, ·) is then guaranteed by the Lax-Milgram lemma P-almost surely. If u(ω, ·) is a solution of 39 4.3. DISCONTINUOUS RANDOM ELLIPTIC PROBLEMS Eq. (4.5) for given f(ω, ·) ∈ H and g(ω, ·) ∈ L2(Γ2), then a−(ω) 2 min(1, C−2 |D|)‖u(ω, ·)‖2 V ≤ Ba(ω)(u(ω, ·), u(ω, ·)) = Fω(u(ω, ·)) ≤ (‖f(ω, ·)‖H + CD‖g(ω, ·)‖L2(Γ2))‖u(ω, ·)‖V . Using Hölder’s and Minkowski’s inequality together with r = (1/p+ 1/q)−1 ≥ 1 yields ‖u‖Lr(Ω;V ) ≤ 2 max(1, CD) min(1, C−2 |D|) E ( a−p− )1/p E ( (‖f‖H + ‖g‖L2(Γ2))q )1/q ≤ 2 max(1, CD) min(1, C−2 |D|) ‖1/a−‖Lp(Ω;R)︸ ︷︷ ︸ :=C(a−,D,p) (‖f‖Lq(Ω;H) + ‖g‖Lq(Ω;L2(Γ2)) < +∞. In the next section, we introduce the diffusion coefficient a, which allows us to incorporate discontinuities at random points or areas in D. We show the existence and uniqueness of a weak solution to the disc