Effective Approximations of Stochastic Partial Differential Equations based on Wiener Chaos expansions and the Malliavin Calculus by Chia Ying Lee B. S., University of Michigan, Ann Arbor, 2005 Sc. M., Brown University, 2007 A Dissertation submitted in partial fulfillment of the requirements for the Degree of Doctor of Philosophy in the Division of Applied Mathematics at Brown University Providence, Rhode Island May 2011 c Copyright 2011 by Chia Ying Lee This dissertation by Chia Ying Lee is accepted in its present form by the Division of Applied Mathematics as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Boris L. Rozovsky, Director Recommended to the Graduate Council Date George Karniadakis, Reader Date Kavita Ramanan, Reader Approved by the Graduate Council Date Peter M. Weber, Dean of the Graduate School iii Curriculum Vitæ Chia Ying Lee was born in Penang, Malaysia on October 16, 1983, and lived and received her basic education in Singapore. She received a Bachelor of Science with High Distinction in Mathematics and Music from the University of Michigan at Ann Arbor in August 2005. From August 2005-2006, she worked at the Bioinformatics Institute, part of the Agency for Science, Technology and Research, in Singapore. She began her graduate studies in the Division of Applied Mathematics at Brown University in September 2006, where she worked under the supervision of Professor Boris Rozovsky. While at Brown, she received a Masters of Science in Applied Mathematics in May 2007. She was also awarded the Stella Dafermos Award in May 2011. iv Dedicated to My Parents, Family and Percy v Acknowledgements I would like to express my deepest gratitude to my thesis advisor Professor Boris Ro- zovsky, who has provided invaluable guidance throughout my graduate career, and has helped to point me in the right directions of my research. I am thankful for the freedom he has allowed me to discover my own interests and for the wide opportunities which, under his tutelage, was opened to me. Thus, as my academic father, I am indebted to him. I would also like to thank all the professors and fellow students in the Division and at Brown who have helped make my graduate experience at Brown an immensely fruitful and enjoyable period of professional and personal growth. Special thanks goes to my committee members Professors George Karniadakis and Kavita Ramanan, who have kindly devoted their time to referee my thesis. I must also include to thank Professors Karniadakis, Chi- Wang Shu, David Gottlieb, Hao-Min Zhou and Bjorn Sandstede, among many others, who all played instrumental roles, both indirectly and direcly, in shaping the various aspects of my research and learning. Last but not least, I am eternally grateful to my parents, family and my fiance, Percy, who have all given me incredible support, and though far away have been a constant source of motivation and strength to pursue my dreams. vi Contents Curriculum Vitæ iv Acknowledgements vi List of Tables ix List of Figures x Chapter 1. Introduction 1 Chapter 2. The Wiener Chaos Expansion and the Malliavin Calculus Framework 5 1. Gaussian white noise and the Wiener Chaos expansion 5 2. Generalized Malliavin calculus and the Wick product 11 3. A basic application of the Wiener chaos expansion to solving SPDEs via the propagator system 14 Chapter 3. Error Analysis for the Stochastic Finite Element Method 25 1. Review of the error estimates for the finite element method for deterministic PDE 27 2. The stochastic finite element method formulation 32 3. Error analysis for SPDE with time independent operators 34 4. The SFEM for SPDE with time-dependent operators 51 5. Numerical Simulations 59 Chapter 4. Unbiased perturbations of the Navier-Stokes equations 66 1. Functional analysis framework 68 2. Stationary QSNS 71 3. The time-dependent QSNS (4.1) 76 4. Long time convergence to the stationary solution 83 5. Finite Approximation by Wiener Chaos Expansions 89 vii 6. The Catalan numbers method 91 Chapter 5. Randomization of Incoherent Forcing for Improvement of Energy Approximations 94 1. Introduction 94 2. Change of Wiener chaos basis 97 3. Comparative error analysis and 1st order improvement 103 4. Examples and simulations 114 Bibliography 121 viii List of Tables 1 Absolute errors (and convergence orders) from the finite element part under the weights qk ∼ k −10 . (Values are squared of the error norm.) 63 2 Truncation error (and convergence order) for each fixed p. (Values are squared of the error norm.) 63 1 Improvement nP /nE in the number of basis elements required to attain 5% error. 118 ix List of Figures 1 Order of convergence (in N ) for the JN,p part, as the power decay of the weights qk ∼ k −s varies. (Orders are computed for the square of the error norm.) 64 1 ¯ (n) ] when the system is truncated to n coefficients, Relative errors incurred R[u under the point forcing basis (dotted line) and the eigenfunction (cosine) basis (solid line). The convection-diffusion equation was used to produce this data. 104 2 (a) Relative errors on log-log axes for increasing values of N , under the cosine basis for the heat equation. (b) Relative errors for two values of diffusion coefficients ǫ = 0.1, 0.01. The graph for ǫ = 0.01 lies above the graph for ǫ = 0.1. 104 3 Log scale plots of the relative errors incurred by the truncation of the convection- diffusion system, as well as the pure diffusion and the pure convection systems. The convection and diffusion coefficients are b = 6b0 and ǫ = 0.1, respectively. 120 x Abstract of “Effective Approximations of Stochastic Partial Differential Equations based on Wiener Chaos expansions and the Malliavin Calculus” by Chia Ying Lee, Ph.D., Brown University, May 2011 This thesis studies the application of the Wiener chaos expansion in the analysis of stochastic partial differential equations (SPDEs). Specifically, linear parabolic SPDEs and the quantized stochastic Navier-Stokes equations are considered, under the framework of the Malliavin calculus. Especially for these highly singular SPDEs, the Wiener chaos expansion is a useful tool for our study of the basic questions of solvability, regularity and dynamical behaviour, and it enables us to study approximations of the solutions of SPDEs and to quantify the errors of approximation. For the quantized stochastic Navier-Stokes equations, we use the Malliavin calculus to formulate a random perturbation of the Navier-Stokes equations that is unbiased, and we will show the existence and uniqueness of steady and time-dependent solutions, as well as the convergence to steady solution, in a stochastic weighted space. We also study a stochastic finite element method for numerical simulation of the solution of linear parabolic SPDEs and derive error estimates for the numerical solution. Finally, we show how one basis of the Wiener chaos expansion can be more efficient than another for approximating the energy of the solution, so that computational efficiency can be increased when applied to some physical applications. CHAPTER 1 Introduction In this thesis, we present analyses and numerical analyses of stochastic partial differential equations using the Malliavin calculus and Wiener chaos expansions, for two classes of SPDEs, linear parabolic SPDEs and the quantized stochastic Navier-Stokes equation. The motivation for choosing the Wiener chaos expansion and Malliavin calculus as a tool of analysis comes in large part from the type of SPDE considered. Since the discovery of the Itˆo integral spurred the development of stochastic analysis and the Itˆo calculus, the study of stochastic models in a myriad of physical, biological and economic applications has caught the wave of Itˆo calculus. However, many SPDE arising in physical and mathematical models, such as those in Uncertainty Quantification, reveal several obvious limitations of the Itˆo calculus, not least the fact that it requires a notion of adaptedness. Uncertainty Quantification frequently deals with models of phenomena for which coefficients or param- eters are not known to full certainty. Rather than neglecting the uncertainty in the model parameters, one can turn to stochastic models as a way to incorporate the uncertainty into the equations. In the simplest case, the uncertainty in a parameter may take the form of being a single random variable of a known distribution—a uniform distribution being a common choice. More complex real world examples of stochastic modelling include the stochastic pressure equations or models of flow in heterogeneous porous media, where the permeability of the medium is difficult to measure at all locations, and is instead modelled as a random field. This is a case of a noncausal system without a natural notion of a filtration, and it would be unwise to confine it into the Itˆo calculus framework. We are thus prompted to appeal to a more general stochastic calculus in order to formulate models outside of the Itˆo calculus framework. This is where the Malliavin calculus comes into the picture. In fact, the Malliavin calculus is not such a far flung idea, because it is an extension of the Itˆo calculus. The Skorokhod integral, one of the important constructs of the Malliavin calculus, extends the Itˆo integral to non-adapted integrands, and coincides with the Itˆo integral for adapted integrands. For this reason, a modeler, when considering 1 to introduce stochasticity into a model, may find the Skorokhod integral a viable modelling choice. Further reasons to work in a more general framework come from the need for more general solution concepts. A desired model for a random perturbation can, and often does, lead to an immediate difficulty of non-square integrable solutions. Early works, such as Walsh’s [64], had already shown that certain equations commonly encountered do not possess solutions with finite variance. The models with uncertain coefficients, involving multiplicative noise that acts on the highest order partial differential operator, are prime examples of equations without square integrable solutions. Lacking a “usual” solution, we are forced to broaden our notion of a solution to include solutions with infinite variance in a larger space of random elements. These spaces are the so-called weighted stochastic spaces, which include the Hida spaces and Kondratiev spaces among others, and whose elements are characterized by their Wiener chaos expansion. The Wiener chaos expansion is a classical orthogonal expansion theory for random func- tions that was first introduced by Cameron and Martin [9]. It representations a square inte- grable random element in an orthogonal expansion of the stochastic variable with respect to a basis derived from the Hermite polynomials. In our case of non-square integrable solutions or random elements, the Wiener chaos expansion is especially pertinent for representing the random elements. The application of the Wiener chaos expansion here is two-fold: to analyze solutions of SPDE through approximate solutions derived from the Wiener chaos expansion; and to quantify the efficiency of approximations of the solutions. The major theme underlying all the analysis in this thesis is the transformation of the single SPDE, via the Wiener chaos expansion, into the related propagator system of PDE. The utility of this transformation is no different from that of the Fourier transformation— the solution is understood as a collection of its expansion coefficients, and the analysis of the SPDE is achieved through the analysis of the propagator system of equations. As noted above, the Wiener chaos expansion is also key to obtaining finite approximations of solutions of SPDEs. The finite approximations, here specifically Galerkin approximations, render the analysis of the SPDE more tractable to analysis. In analogy to deterministic theory, creating approximate solutions is the first step in formulating energy estimates which are then used 2 to deduce the existence of a solution. In fact, because the conversion to the propagator system separates out the stochastic variable and leaves behind a deterministic system of equations, a large part of the stochastic analysis is founded on deterministic theory. Thus, the strength of the stochastic results depend on the strength of deterministic PDE theory. Through the use of the Wiener chaos expansion, the difficulty in the stochastic analysis is greatly reduced to understanding results from deterministic theory. Seen in a different light, the stochastic theory is in fact a generalization of the deterministic theory to the sto- chastic setting. One supporting argument for this is that the Malliavin calculus, historically developed to build a stochastic theory based on the integration by parts formula, accords us with an arsenal of conceptual tools familiar from deterministic theory—an integration by parts formula, adjoint operators and the possibility of defining weak or variational so- lutions by action on test functions. To give an example of such stochastic analogues, we will subsequently encounter the use two of the main constructs of the Malliavin calculus, the Malliavin derivative and the Malliavin divergence operator, which are adjoints of each other under the Gaussian measure. The latter is, in fact, a stochastic convolution and is equivalent to the Skorokhod integral. Related to the Malliavin divergence operator is the Wick product (see e.g. [27,31,35,42]). The Wick product is generally considered a suitable replacement of the usual product, especially when the product is between two generalized random elements for which the usual product is not well defined. Interestingly, although the Wick product was introduced independently in the seemingly unrelated field of quantum field theory, it turns out that the Wick product and the Malliavin divergence operator are closely related [46]. Both are stochastic convolutions between two random elements, and moreover, the two concepts coincide in some cases, so that it is possible in these cases to formulate stochastic equations using one or the other framework. Thus, we will see the use of both the Malliavin divergence operator and the Wick product. The thesis is organized as follows. Chapter 2 introduces the mathematical framework of the Wiener chaos expansion, Malliavin calculus and Wick product. With these tools, we then present a basic technique of applying the propagator system to deduce the solvability of a stochastic parabolic equation under the Malliavin calculus approach. In Chapter 3, we discuss a numerical algorithm based on the finite element method for solving the stochastic parabolic equation, and quantify the error incurred by the numerical solutions by deriving 3 a priori error estimates. We will consider a newly proposed, unbiased, random perturbation of the deterministic Navier-Stokes equations, called the quantized stochastic Navier-Stokes equation, in Chapter 4. We will analyze the solvability of the stationary equations, as well as the long time convergence of a time dependent solution to the steady solution. As an application of the Wiener chaos expansion, we will study in Chapter 5 how certain physical models of incoherent forcing sources can exploit a simple change of Wiener chaos expansion basis to drastically reduce computational cost. As the latter three topics are considerably distinct in their nature, we will leave further introduction to each topic to the start of their respective chapters. 4 CHAPTER 2 The Wiener Chaos Expansion and the Malliavin Calculus Framework In this chapter, we present the tools and techniques which form the framework for the analysis of solutions of SPDEs, as studied in this thesis. We will first discuss the general construction of Gaussian noise, which is the main source of stochasticity in the random perturbations of stochastic equations. This will lead to the definition of the Wiener chaos expansion and the weighted stochastic spaces. The weighted stochastic spaces are introduced as spaces that our solutions live in, and the Wiener chaos expansion is the way our solutions are represented. We then define the main operators in Malliavin calculus and the Wick product, to be used as the stochastic models for the actual incorporation of the stochasticity into an equation. Finally, with the Wiener chaos expansion and the Malliavin calculus approach, we elucidate a technique for analyzing the solvability of a parabolic SPDE by considering the equivalent propagator system for the chaos modes of the solution. 1. Gaussian white noise and the Wiener Chaos expansion Let ξ = {ξk }k≥1 be a collection of i.i.d. N (0, 1) random variables on a probability space (Ω, F, P), where F is the σ-algebra generated by {ξk }. Let U be a real separable Hilbert space with complete orthonormal basis {uk }k≥1 . Definition 1.1. The Gaussian white noise on U is the formal series ∞ X (2.1) ˙ := W ξk uk . k=1 Note that the white noise is not an element of L2 (Ω; U ), because ∞ X ˙ k2 = EkW kuk k2U = ∞ U k=1 We list some special examples of Gaussian noise commonly encountered: 5 ˙ (t), we take U = L2 (0, T ). The basis {uk } 1. For the standard 1-dimensional white noise W may, for example, be taken to be the cosine basis in U , but this is not a unique choice ˙ (t) may be understood as the formal derivative of of basis. If t represents time, then W a 1-dimensional Brownian motion W (t). 2. For stationary or spatial Gaussian white noise on a domain D ⊂ Rd , we take U = L2 (D). ˙ (x) is spatially uncorrelated; that is, for x, y ∈ D, Then the definition yields that W ˙ (x)W E[W ˙ (y)] = δx (y) where δx is the Dirac delta function at x. Indeed, for smooth φ, DX ∞ ∞ X E ∞ X ˙ (x)W hE[W ˙ (y)], φi = uk (x)uk (·), φk uk (·) = φk uk (x) = φ(x) k=1 k=1 k=1 3. If U = L2 (0, T ; H) for some separable Hilbert space H, then the cylindrical Wiener process with values in H is ∞ X (2.2) W (t) := uk Wk (t), k=1 where Wk (t), k = 1, 2, . . . , are independent standard Wiener processes. If H = L2 (D), ˙ (t, x) is called space-time white noise. the formal derivative W ˙ Q (x) with covariance operator Q2 . Let Q be 4. Correlated (or weighted) Gaussian noise W an operator on U defined by Quk = σk uk , for k = 1, 2, . . . . where {σk , k ≥ 1} are non-negative real numbers. The Gaussian noise with covariance operator Q2 is defined by ∞ X (2.3) ˙ Q (x) := W σ k uk ξk . k=1 P∞ If Q2 is nuclear or trace class, i.e., 2 k=1 σk < ∞, then it is defined by the covariance function ∞ X ˙ Q (x)W q(x, y) := E[W ˙ Q (y)] = σk2 uk (x)uk (y) k=1 6 R as the operator Q2 f (x) = D q(x, y)f (y) dy, for f ∈ L2 (D). The expansion (2.3) is the ˙ Q. Karhunen–Lo`eve expansion for W ˙ 1, W 5. Let W ˙ 2 be two white noises on U1 , U2 respectively, ∞ X ˙ i := (i) (i) W ξk uk , for i = 1, 2. k=1 ˙ i to be independent, or correlated in some way. We can define an We may assume W ˙ to accommodate both noises W abstract white noise W ˙ i into a single term, by defining (1) (2) (1) (2) u2k−1 = uk , u2k = uk , and ξ2k−1 = ξk , ξ2k = ξk . Then ∞ X ˙ := W ξk uk . k=1 The formulation of the abstract noise is useful for studying equations that are driven by two or more distinct white noises, or equations whose input data (initial or boundary conditions, or forcing terms) are measurable with respect to a different Gaussian noise than the noise driving the equation. In order to develop the L2 theory of F-measurable random elements and the Wiener chaos expansions, we first introduce some housekeeping tools. Let J = {α = (α1 , α2 , . . . ) : P αk ∈ N0 } be the set of multi-indices of finite length, |α| := k≥1 αk < ∞. Denote dim(α) = P min{k : ακ = 0 for κ > k} and d(α) = ∞ k=1 1αk >0 . We denote the zero multi-index by (0) = (0, 0, . . . ), and the unit multi-index with 1 in the kth entry by ǫk . For α, β ∈ J , Y α + β = (α1 + β1 , α2 + β2 , . . . ), α! := αk ! k≥1     α+β (α + β)! |α| |α|! = , = α α!β! α α! Q αk For a sequence of nonnegative real numbers q = (q1 , q2 , . . . ), define q α = k≥1 qk . A multi-index α can be uniquely characterized by its characteristic set Kα . Let n = |α|, and denote the ordered n-tuple Kα = (k1 , . . . , kn ) where k1 ≤ k2 ≤ · · · ≤ kn , with ki defined as follows. Let κ1 < · · · < κd(α) be the indices of α for which ακi 6= 0. Then i−1 X i X kl = κi if α κj < l ≤ α κj . j=1 j=1 7 In words, the first ακ1 entries of Kα are k1 = · · · = kακ1 = κ1 , followed by the next ακ2 entries of Kα being kακ1 +1 = · · · = kακ1 +ακ1 = κ2 , etc. The definitions of the multi-indices and their characteristic sets give rise to some useful combinatorial results which will come in handy later. We state two of these results here. Lemma 1.2. (A multinomial sum in infinite dimensions) Let ρ ~ = (ρ1 , ρ2 , . . . ) with P ρk > 0, and let ρ¯ = k≥1 ρk . Then for any n ∈ N0 , X ρα ρ¯n = . α! n! |α|=n Proof. Fix n and |α| = n. We identify α with its characteristic set Kα = (k1 , . . . , kn ). Since there are n!/α! distinct permutations of {k1 , . . . , kn }, Qn Qn X ρα X j=1 ρkj (n!/α!) X j=1 ρkj 1 = · = · α! α! (n!/α!) α! (n!/α!) |α|=n k1 ≤···≤kn k1 ,...,kn where we have multiplied by 1 and rearranged the sum over non-decreasing indices into a sum over all unordered indices. Finally, from the formula for the multinomial expansion X ρα n  n 1 X Y 1 X = ρkj = ρk α! n! n! |α|=n k1 ,...,kn j=1 k  Lemma 1.3. For all α, β ∈ J , |β|! |α − β|! |α|! ≤ . β! (α − β)! α! |α|! Proof. Let Kα = (k1 , . . . , k|α| ) be the characteristic set of α. On the RHS, α! is the number of distinct permutations of Kα . On the LHS, we partition Kα into the two subsets corresponding to Kβ and K(α−β) . Then, the number of distinct permutations of Kβ times that of K(α−β) cannot exceed the number of distinct permutations of Kα .  1.1. The Wiener chaos expansion and weighted Wiener chaos spaces. The Wiener chaos expansion is an orthogonal expansion for random elements that are measurable ˙ . Due to the Gaussian assumption, the Wiener with respect to the Gaussian white noise W chaos expansion is necessarily an expansion in the Hermite polynomials. Recall the Hermite 8 polynomial Hn (x) of degree n, x2 dn − x 2 Hn (x) = (−1)n e 2 e 2. dxn For each α ∈ J , define the random variables Y Hα (ξk ) ξα = √k . k≥1 αk ! Theorem 1.4. (Cameron and Martin [9]) The collection Ξ = {ξα , α ∈ J } is an or- thonormal basis of L2 (Ω). Ξ is referred to as the Cameron-Martin basis. Given a real separable Hilbert space X with norm |·|X , let L2 (Ω; X) be the Hilbert space of square integrable F-measurable random elements with values in X. Then the Cameron- Martin theorem provides that any square integrable random element ζ ∈ L2 (Ω; X) has the Wiener Chaos expansion with respect to the Cameron-Martin basis, X ζ= ζα ξα α∈J where ζα = E[ζξα ], and Parseval’s identity holds, X kζkL2 (Ω;X) ≡ E|ζ|2X = |ζα |2X . α∈J We will frequently encounter random elements that are not square integrable, and thus we describe a construction analogous to the construction of Sobolev scales. Define the test function space X D = {ζ = ζα ξα : ζα ∈ R and only finite number of ζα are non-zero}. α Definition 1.5. A generalized random element f with values in X is a formal series X (2.4) f= fα ξα , α∈J where fα ∈ X. f is identified with the sequence {fα , α ∈ J }. The expansion (2.4) is also called the Wiener chaos expansion of f . 9 The space D′ (X) of generalized random elements in X is the dual space of D with respect to L2 (Ω), with duality pairing X hhf, ζii = ζα fα α The space D′ is a very large space. Its elements have Wiener chaos expansions that may exhibit severe blow-up. Next, we introduce the weighted Wiener Chaos spaces that quantify the asymptotic behaviour of the Wiener chaos modes. Let R be a bounded linear operator on L2 (Ω) defined by Rξα = rα ξα for every α ∈ J , where the weights {rα , α ∈ J } are positive real numbers. Note that R is bounded if and only if the weights rα are uniformly bounded from above, that is, rα < C for all α ∈ J , for some constant C. Define the norm X kf k2RL2 (Ω;X) := |fα |2X rα2 α∈J P for f = α∈J fα ξα . The space RL2 (Ω; X) of random elements in X, is defined as the closure of L2 (Ω; X) under the norm k · kRL2 (Ω;X) ; in other words, the elements of RL2 (Ω; X) are P identified with a formal series α∈J fα ξα , where kf k2RL2 (Ω;X) < ∞. Clearly, RL2 (Ω; X) is a Hilbert space with respect to k · kRL2 (Ω;X) . The operator R−1 that is inverse to R is defined by R−1 ξα = rα−1 ξα . Let X ֒→ Y ֒→ X ′ be a normal triple of Hilbert spaces with duality pairing h·, ·iX ′ ,X . We define the space R−1 L2 (Ω; X) as the dual of RL2 (Ω; X ′ ) relative to the inner product in the space L2 (Ω; Y ). The duality pairing is given by X hhf, giiRL2 (Ω;X ′ ),R−1 L2 (Ω;X) := EhRfα , R−1 gα iX ′ ,X = hfα , gα iX ′ ,X α∈J for f ∈ RL2 (Ω; X ′ ) and g ∈ R−1 L2 (Ω; X). Similarly, R−1 L2 (Ω; X ′ ) is defined as the dual of RL2 (Ω; X) relative to the inner product in L2 (Ω; Y ). We may leave out notating the dual spaces in hh·, ·ii where it is either obvious or inconsequential. There are several classes of weights in the literature [28, 31, 37, 47]. We list a few. (1) In Section 3 and Chapter 3, we consider only admissible weights of the form qα rα2 = , |α|! 10 where q = (q1 , q2 , . . . ) is a decreasing sequence of positive real numbers. This class of weights arises naturally, for example, when studying equations where the driving white noise acts on the term of the second order partial differential operator [44, 48, 49].  (2) Kondratiev spaces. Denote the sequence (2N)−q := (2k)−q k=1,2,... . The Kon- dratiev space S−1,−q (X) is a weighted space with weights α (2N)−q rα2 = α! The Kondratiev spaces have been widely used to study various classes of SPDE (see e.g., [10, 29, 31]). 2. Generalized Malliavin calculus and the Wick product In this section, we discuss a generalized form of the Malliavin calculus, and also its relation to the Wick product. Roughly speaking, the traditional development of the subject of the Malliavin calculus begins with defining the Malliavin derivative D W˙ with respect to Gaussian white noise, N X ∂F D W˙ F (W (h1 ), . . . , W (hN )) = (W (h1 ), . . . , W (hN ))hi , ∂xi i=1 RT for a smooth function F and hi ∈ U , i = 1, . . . , N , and where W (hi ) = 0 hi dW (t). The Malliavin derivative maps elements in L2 (Ω) into L2 (Ω; U ). The Skorokhod integral δ W˙ is then a map L2 (Ω; U ) into L2 (Ω), defined via the adjoint property, E[δ W˙ (f )φ] = E[(f, D W˙ φ)U ]. (See [50, 56] for details.) The generalized form of the Malliavin calculus retains the properties of the traditional Malliavin calculus, including the adjoint property, but accords more flexibility when it comes to defining the Malliavin derivative and Malliavin divergence operator with respect to other random elements besides Gaussian white noise. For any ξk , define the Malliavin derivative 11 D ξk and Malliavin divergence operator δ ξk by1 √ √ D ξk (ξα ) := αk ξα−ǫk , and δ ξk (ξα ) := αk + 1ξα+ǫk . The Malliavin operators D ξk , δ ξk can be extended to any Cameron-Martin basis element ξβ by s  s  α α+β D ξβ (ξα ) := ξα−β , and δ ξβ (ξα ) := ξα+β . β β An important relationship between the Malliavin derivative and Malliavin divergence oper- ator is the adjoint property: for α, α′ , β ∈ J , (2.5) hhδ ξβ (ξα ), ξα′ ii = hhξα , D ξβ (ξα′ )ii By bilinearity, D u (v) and δ u (f ) can be defined for random elements u on U , v on X, and f on X ⊗ U [46]. Elementary computations with the Wiener chaos expansion yields P explicit formulas for D u (v) and δ u (f ), as follows. Let u = α uα ξα with uα ∈ U , and P P v = α vα ξα with vα ∈ X, and f = α fα ξα with fα ∈ X ⊗ U . Then  s  X X α + β  D u (v) =  vα+β ⊗ uβ  ξα , α β β  s  X X α δ u (f ) =  (fβ , uα−β )U  ξα . α β β≤α The Malliavin divergence operator is closely related to the Wick product. The Wick product is defined as s  α+β ξα ⋄ ξβ = ξα+β β and is extended by linearity to f ⋄ η, where f is a generalized X-valued random element and η is a generalized real-valued random element. Clearly, f ⋄ η = δ η (f ) in this case (with U = R). The difference between the Malliavin divergence operator and the Wick product lies in the fact that the Wick product is a point-wise product between X-valued and R-valued generalized random elements, whereas the Malliavin divergence operator, being a stochastic integral, is a convolution between U -valued and U ⊗ X-valued generalized random elements. 1Recall the multi-index notation on page 7. 12 Thus, the Wick product is a symmetric operator, whereas the Malliavin divergence operator is not symmetric (see [46]). Next, we describe some situations in which we will encounter the Malliavin derivative, Malliavin divergence operator and the Wick product: ˙ , (1) We can define the Malliavin divergence operator with respect to white noise W which is a random element on U . For f ∈ RL2 (Ω; X ⊗ U ), δ W˙ (f ) is the unique element of RL2 (Ω; X) with the property that δ W˙ (f ), ϕ RL2 (Ω;X),R−1 L2 (Ω;X ′ ) = f, D W˙ (ϕ) RL2 (Ω;X⊗U ),R−1 L2 (Ω;X ′ ⊗U ) for every ϕ ∈ R−1 L2 (Ω; X ′ ) such that D W˙ (ϕ) ∈ R−1 L2 (Ω; X ′ ⊗ U ). In the case of time white noise on a finite time interval, i.e., U = L2 (0, T ), it can be shown that the Malliavin divergence operator coincides with the Itˆo integral under the assumption of adaptedness [56]. That is, Z T δ W˙ (u) = u(t)dW (t), 0 provided u(t) is a suitable random element that is adapted to the filtration gener- ated by the Brownian motion W (t). (2) For a random element g ∈ RL2 (Ω; X), it will often arise in modelling problems ˙ . Strictly to consider the multiplication, or convolution, of g with white noise W speaking, the term δ W˙ (g) is not well defined. But, by abuse of notation, we interpret X√ [δ W˙ (g)]α = αk gk,α−ǫk k≥1 where gk,α = uk ⊗ gα . ˙ (x) is a white noise on L2 (D) with orthonormal basis {uk (x)}, and When W when g(x) is a generalized random function in x, this interpretation coincides with ˙ (x) by taking the Wick product pointwise in x. the Wick product model g(x) ⋄ W To see this, choose gk,α (x) = uk (x)gα (x), then   X√ ˙ (x) g(x) ⋄ W = αl gα−ǫl (x)ul (x) = δ W˙ (g) α l 13 If {uk } is chosen such that Dγ uk ∈ L∞ for all k and all γ = (γ1 , . . . , γd ) with ˙ ∈ D′ (W l,p ) provided g ∈ D′ (W l,p ). |γ| ≤ l, then g ⋄ W (3) In Chapter 4, we will consider a nonlinearity of the form ui ⋄ ∂xi u. Direct compu- tation gives that X q  α (2.6) (ui ⋄ ∂xi u)α = γ (uγ , ∇)uα−γ . 0≤γ≤α Each chaos mode of the Wick product is determined by a convolution among be- tween the lower order chaos modes. This observation is important, as it suggests a connection to the Catalan numbers which are themselves characterized recursively as convolutions. 3. A basic application of the Wiener chaos expansion to solving SPDEs via the propagator system The Wiener chaos expansion is used as a separation of variables technique for stochastic ordinary or partial differential equations. Just like how the Fourier expansion is used to solve for the Fourier modes of a solution of a deterministic PDE, the Wiener chaos expansion is used in an analogous way to find the Wiener chaos modes of the solution of an SPDE. The separation of the independent variable of randomness leaves behind a deterministic system of equations for the Wiener chaos modes of the solution, termed the propagator system of equations. The analysis of the SPDE is thus reduced to the analysis of a deterministic PDE system for which deterministic theory can be applied. In this section, we show a basic application of the Wiener chaos expansion to study the parabolic SPDE on a bounded domain D ⊂ Rd , du dt + Au + δ W˙ (Mu + g) = f on D × (0, T ] (2.7) u|∂D = 0 u|t=0 = v where W ˙ is a Gaussian white noise which may depend on space or time, A is a second P order elliptic operator from H01 (D) onto H −1 (D), and Mu := k Mk u ⊗ uk where Mk , k = 1, 2, . . . are bounded operators from H01 (D) into H −1 (D). We will assume that the boundary ∂D and the coefficients of A, Mk are sufficiently smooth in the space and time 14 variables. The input data v, f, g are allowed to be generalized random elements, and we ˙ . recall that they can be measurable with respect to a white noise different from W 3.1. Some notation and constants. Before we proceed, we state our notation for various constants that will show up often in the rest of this chapter and the next chapter. We assume throughout that A(t) is uniformly elliptic on (0, T ] and is coercive and bounded, coerc kuk2 , hA(t)u, ui ≥ CA H1 (2.8) 0 hA(t)u, vi ≤ b kuk CA H01 kvkH01 . ellip coerc )−1 to be the constant in Denote CA = (CA ellip (2.9) kwkH01 ≤ CA kf kH −1 for the solution of the zero Dirichlet problem A(t)w = f , for any t ∈ (0, T ]. Also denote CA to be the constant in (2.10) kwkL2 (0,T ;H01 (D)) ≤ CA (kw0 kL2 (D) + kf kL2 (0,T ;H −1 (D)) ) dw for the weak solution w of the zero Dirichlet problem dt + A(t)w = f with w(0) = w0 . (r) For the Sobolev spaces H r (D), r ≥ 1, let λk be the constants in (r) (2.11) kMk (t)wkH r−2 (D) ≤ λk kwkH r (D) , ∀w ∈ H r (D), t ∈ (0, T ] (1) ellip For brevity, we write λk = λk . Observe that Ck := λk CA are the constants defined by kA−1 Mk vkH 1 ≤ Ck kvkH 1 for all v ∈ H0X 1 . 0X 0X (r) (r) Finally, let µk , r = −1, 0, 1, . . . be the constant arising in kgk kHXr ≤ µk kgkHXr . We (−1) will write µk = µk . −1 We will use shorthand to denote the spaces: for example, we will write RΩ L2T HX to denote RL2 (Ω; L2 ((0, T ); H −1 (D))). Also, H0X 1 denotes H 1 (D). 0 3.2. Existence and uniqueness of solutions. We begin by defining the notion of a weak solution of (2.7). 15 −1 Definition 3.1. A weak solution of (2.7), with f, g ∈ RΩ L2T HX and v ∈ RΩ L2X , is a 1 such that for every φ ∈ R−1 with D φ ∈ R−1 U , process u ∈ RΩ L2T H0X Ω ˙ W Ω Z t Z t (2.12) hhu(t), φii = hhv, φii − hhAu + δ W˙ (Mu + g), φiids + hhf, φiids 0 0 −1 with equality in L2T HX . The existence and uniqueness result for a weak solution of (2.7) for v, f deterministic and g ≡ 0 has been shown in [49]. The proof relies on the Equivalence Theorem 3.2 that relates the weak solution to the propagator system (2.13). P Theorem 3.2. The process u = α u α ξα ∈ RΩ L2T H0X 1 is a solution of (2.7), if and only if, for each α ∈ J , Z t X√ Z t (2.13) uα (t) = vα − Auα (s) + αk (Mk uα−ǫk + gk,α−ǫk ) ds + fα (s)ds 0 k≥1 0 −1 holds in HX for a.e. t ∈ [0, T ]. Proof. See [49].  Using the techniques from [47] (Theorem 9.4) or [45] (Proposition 4.2), we can extend the existence and uniqueness result to the case when v, f, g are random, and determine the conditions for the weighted spaces that u may belong to, in terms of the spaces that the input data belong to. qα Theorem 3.3. Let the weights R, with rα2 = |α|! , satisfy X 2 2 (2.14) qk C A λk < 1. k≥1 −1 (1) If the input data v ∈ L2X and f, g ∈ L2T HX are deterministic, then there exists a 1 , and unique weak solution u ∈ RΩ L2T H0X   kukRΩ L2 H 1 ≤ C kvkL2 + kf kL2 H −1 + kgkL2 H −1 T 0X X T X T X where C depends only on R, A, M and T . 16 ¯ Ω L2 and f, g ∈ R ¯ Ω L2 H −1 for some r¯α2 = ρα (2) Assume v ∈ R X T X |α|! . Also assume, in addition to (2.14), that qk are chosen to satisfy X qk (2.15) <1 ρk k≥1 Then there exists a unique solution u ∈ RΩ L2T H0X 1 , and   kukRΩ L2 H 1 ≤ C kvkR¯ Ω L2 + kf kR¯ Ω L2 H −1 + kgkR¯ Ω L2 H −1 T 0X X T X T X ¯ A, M and T . where C depends only on R, R, Proof. Step 1. Assume v, f, g are non-random. This case has been studied in [49] for g = 0. The proof here is essentially the same. The propagator system is Z t u(0) (t) = v + Au(0) (s) + f (s)ds 0 Z t  uǫk (t) = Auǫk (s) + Mk u(0) (s) + gk (s) ds 0 Z t X√ uα (t) = Auα (s) + αk Mk uα−ǫk (s)ds, |α| ≥ 2 0 k Let Φt = eA t be the semigroup generated by A. Then for α = (0), Z t u(0) (t) = Φt v + Φt−s f (s)ds 0 and from the deterministic parabolic estimates, ku(0) kL2 H 1 ≤ CA (kvkL2 + kf kL2 H −1 ) T 0X X T X For α = ǫk , Z t  uǫk (t) = Φt−s1 Mk u(0) (s1 ) + gk ds1 0 17 and again applying the deterministic parabolic estimates,   kuǫk kL2 H 1 ≤ CA kMk u(0) kL2 H −1 + kgk kL2 H −1 T 0X T X T X   ≤ CA λk CA ku(0) kL2 H 1 + µk kgkL2 H −1 T 0X T X  ≤ CA M ~λCA (kvkL2 + kf kL2 H −1 + kgkL2 H −1 ) X T X T X µk where M = supk (1 ∨ λk CA ). For |α| = n ≥ 2, with characteristic set Kα = (k1 , . . . , kn ), it can be shown by induction that Z Z Z s2 1 X t sn  uα (t) = √ ... Φt−sn Mkσ(n) . . . Φs2 −s1 Mkσ(1) u(0) (s1 ) + gkσ(1) ds1 . . . dsn α! σ∈Pn 0 0 0 where Pn is the group of permutations of {1, . . . , n}. Then, we obtain  ~ α |α|! C kuα kL2 H 1 ≤ CA M √ (kvkL2 + kf kL2 H −1 + kgkL2 H −1 ) T 0X α! X T X T X ~ = (C1 , C2 , . . . ), Ck = λk CA . where C Taking the weights to satisfy (2.14), it follows from Lemma 1.2 that kukRΩ L2 H 1 ≤ C(kvkL2 + kf kL2 H −1 + kgkL2 H −1 ) T 0X X T X T X where C depends only on R, A, M and T . Step 2. Fix an arbitrary α∗ ∈ J . Assume v = V ξα∗ , f = F ξα∗ , g = Gξα∗ ; in other words, the randomness of the data is localized to a single mode. Let u[α∗ ; V, F, G](t, x) be the solution. By linearity, the chaos expansion coefficients with indices of the form α∗ + α satisfy uα∗ +α [α∗ ; V, F, G] uα [(0); √Vα∗ ! , √Fα∗ ! , √G α∗ ! ] p = √ (α∗ + α)! α! and are zero otherwise. Then Z T ku[α∗ ; V, F, G](t)k2RΩ H 1 dt 0X 0 X q  α∗ +α  2 (α∗ + α)! V F G = |α∗ + α|! α! uα (0); √α∗ ! , √α∗ ! , √α∗ ! 2 1 α LT H0X 18 X q α∗ +α |α|!|α∗ |! (α∗ + α)! = ∗ ∗ ∗ kuα [(0); V, F, G]k2L2 H 1 α |α|!|α |! |α + α|! α!α ! T 0X ∗ qα ≤ ∗ ku[(0); V, F, G]k2RΩ L2 H 1 |α |! T 0X where the last inequality follows by Lemma 1.3. Step 3. ¯ Ω L2 and f, g ∈ R For the general case with random data, assume v ∈ R ¯ Ω L2 H −1 . The X T X solution can be written as X u= u[α∗ ; vα∗ , fα∗ , gα∗ ] α∗ Using the estimates from Step 2, X kukRΩ L2 H 1 ≤ ku[α∗ ; vα∗ , fα∗ , gα∗ ]kRΩ L2 H 1 T 0X T 0X α∗ !1/2 !1/2 X q α∗ |α∗ |! X ρ α∗  2 ≤C kv α ∗ k 2 + kfα∗ k 2 LX −1 + kg α ∗k 2 −1 α∗ |α∗ |! ρα∗ α∗ |α∗ |! LT HX LT HX   ≤ C kvkR¯ Ω L2 + kf kR¯ Ω L2 H −1 + kgkR¯ Ω L2 H −1 X T X T X where we have applied Cauchy-Schwartz inequality in the second inequality. The conver- P ∗  q α |α∗ |! gence of ∗ ∗ α |α |! ρ α ∗ follows from a sufficient condition such as (2.15). ¯ so u is a weak solution of (2.7) in the sense of Definition 3.1. Uniqueness Clearly, R ⊇ R, follows from the uniqueness of each equation in the propagator system.  µk Remark. The validity of the assumption that M := supk (1 ∨ λk ) < ∞ arises in some common examples. For example, taking Mk φ = uk ∆φ and gk = uk g, we have that µk , λk are both ∼ O(k). If M = ∞, then in the estimate for kuα kL2 H 1 in Step 1, we should T 0X ~ α ~ P ~ ) , and use the criterion k qk (λk CA ∨ µk )2 < 1 in replace the factor M λ by (λCA ∨ µ α place of (2.14). ¯ for Remark. If the input data is non-random, then it belongs to any weighted space R any ρ. In this case, condition (2.15) is automatically satisfied, and the condition for optimal solution weights R reduces to (2.14) alone. 3.3. Higher regularity of solutions. The weak solution of (2.7) is a generalized 1 , and we now investigate when it possesses better smoothness in the spatial process in H0X variable. Such a result will become useful in the analysis of the stochastic finite element 19 method in Chapter 3, because the spatial regularity of the solution is closely related to the convergence order of finite element schemes. In fact, within the limitations of our analysis in Chapter 3, we require at the least that the minimum spatial regularity of the solution u be H 2 (D)-smooth, and its time derivatives ut , utt be L2 (D) and H −1 (D) functions respectively. However, obtaining higher spatial regularity comes at the expense of worsening the weights R. In line with the strategy of the previous section, we derive higher regularity results from the analogous results in the deterministic theory. Thus, we will see, for example, that certain compatibility conditions at time t = 0 are necessary conditions for higher regularity to hold, except that the compatibility conditions in the stochastic case are more extensive than those in the deterministic case. We now recall a result from deterministic PDE. Theorem 3.4. (Evans [14], Theorems 5 and 6 in §7.1.32). Let A be a uniformly elliptic second order operator whose coefficients belong to HT1 WXk,∞ . Suppose u ∈ L2T H0X 1 −1 with ut ∈ L2T HX is the weak solution of ut + Au = f in D × (0, T ] u|∂D = 0 u(0) = v (i) Assume 1 v ∈ H0X , f ∈ L2T L2X . Then in fact u ∈ L2T HX 2 ∩ L∞ H 1 and u ∈ L2 L2 , and T 0X t T X   ess sup ku(t)kH 1 + kukL2 H 2 + kut kL2 L2 ≤ C0reg kvkH 1 + kf kL2 L2 0X T X T X 0X T X 0≤t≤T where the constant C0reg depends only on D, T and A. (ii) Fix m ≥ 1. Assume 2m+1 dk f 2m−2k v ∈ HX , ∈ L2T HX for k = 0, . . . , m dtk 2The statement of the results assumes that the operator A does not depend on time. A careful analysis of the proof shows that a similar result holds for time dependent operators under the assumptions on the coefficients described in this section. 20 and suppose the m-th order compatibility conditions hold:    1 , 1 ,..., v0 := V0 ∈ H0X V1 := f (0) − AV0 ∈ H0X   dm−1 f 1 . Vm := dtm−1 − AVm−1 ∈ H0X dk u 2m+2−2k Then dtk ∈ L2T HX for k = 0, . . . m + 1, and m k m k ! X d u X d f ≤ reg Cm kvkH 2m+1 + dtk 2m+2−2k X dtk 2m−2k k=0 L2T HX k=0 L2T ;HX reg where the constant Cm depends only on m, D, T and A. From Theorem 3.4(i), we can obtain the following higher regularity result for the sto- chastic equation (2.7), with deterministic input data. The case of random data can be shown in the same way as Steps 2 and 3 in the proof of Theorem 3.3. Corollary 3.5. Suppose u ∈ RΩ L2T H0X 1 is the weak solution of the SPDE (2.7). Also assume that v, f, g are deterministic with 1 v ∈ H0X , and f, g ∈ L2T L2X . ˜ satisfying Then for the weights R X (2) ρ˜k (λk C0reg )2 < 1, k ˜ Ω L2 H 2 and the weak solution u ∈ R T X   kukR˜ Ω L2 H 2 ≤ C kvkH 1 + kf kL2 L2 + kgkL2 L2 T X 0X T X T X Proof. The proof is similar to the proof of Theorem 3.3. The estimates for each uα are obtained by applying Theorem 3.4(i) to the propagator system.  No special compatibility conditions were necessary for Corollary 3.5, but it is unable to ensure boundedness of utt . Thus, we next show how to obtain a smoother solution and the boundedness of utt using the 1st order compatibility conditions. 21 Corollary 3.6. Suppose u ∈ RΩ L2T H0X 1 is the weak solution of the SPDE (2.7). Also assume that v, f, g are deterministic with 3 df dg v ∈ HX , and f, g ∈ L2T HX 2 , and , ∈ L2T L2X , dt dt and that the 1st order compatibility conditions hold for {v, f, gk }:    1 , 1 , v ∈ H0X f (0) − Av ∈ H0X (2.16)   1 Mk v + gk (0) ∈ H0X ∀k = 1, 2, . . . Then for the weights R′ satisfying X  2 (4) (2)  (2.17) ρ′k λk ∨ λk C1reg < 1, k the weak solution u ∈ R′Ω L2T HX 4 , u ∈ R′ L2 H 2 and u ∈ R′ L2 L2 and t Ω T X tt Ω T X kukR′ 2 4 + kut kR′ 2 2 + kutt kR′ 2 2 Ω LT H X Ω LT H X Ω LT LX  ≤ C kvkH 3 + kf kL2 H 2 + kgkL2 H 2 + kft kL2 L2 + kgt kL2 L2 X T X T X T X T X Proof. For α = (0), the (deterministic) compatibility conditions hold, and from The- orem 3.4(ii), ku(0) kL2 H 4 + ku(0),t kL2 H 2 + ku(0),tt kL2 L2 T X T X T X   reg ≤ C1 kvkH 3 + kf kL2 H 2 + kft kL2 L2 . X T X T X For α = εk , since we have assumed the coefficients of Mk to be sufficiently smooth (e.g., at least WX3,∞ ), so u(0) ∈ L2T HX4 implies that M u 2 2 k (0) + gk ∈ LT HX , and u(0),t ∈ LT HX 2 2 implies that (Mk u(0) + gk )t ∈ L2T L2X . The compatibility conditions for (Mk u(0) + gk ) t=0 = Mk v + gk (0) are also satisfied. Again applying Theorem 3.4(ii), kuεk kL2 H 4 + k(uεk )t kL2 H 2 + k(uεk )tt kL2 L2 T X T X T X   (4) (2) (2) (0) ≤ C1reg λk ku(0) kL2 H 4 + θk kgkL2 H 2 + λk k(u(0) )t kL2 H 2 + θk kgt kL2 L2 T X T X T X T X   (4) (2) ˜ ≤ (C1reg )2 (λk ∨ λk )M ˜ kvk 3 + kf k 2 2 + kft k 2 2 + kgk 2 2 + kgt k 2 2 H L H L L L H L L X T X T X T X T X ˜  (2) (µk ∨µk ) (0) ˜ = supk 1 ∨ where M (4) (2) . (The remark following Theorem 3.3 applies.) (λk ∨λk )C1reg 22 For |α| ≥ 2, we have Mk uα−εk ∈ L2T HX 2 and (M u 2 2 k α−εk )t ∈ LT LX . The compatibility conditions hold trivially, since uα−εk t=0 ≡ 0 whenever |α| ≥ 2. The usual computations give the estimates, kuα kL2 H 4 + kuα,t kL2 H 2 + kuεk ,tt kL2 L2 T X T X T X reg ~ ~(2) α (4) ˜ C1 (λ √∨ λ ) |α|! ˜ ≤ C1reg M α!   × kvkH 3 + kf kL2 H 2 + kft kL2 L2 + kgkL2 H 2 + kgt kL2 L2 . X T X T X T X T X The weighted norm kukR′ 2 4 < ∞ provided (2.17) holds.  Ω LT HX Due to the lower triangular property of the propagator system, the first order compat- ibility conditions for the stochastic parabolic equation call for additional conditions on the input data compared to the deterministic case. If the input data is smoother than what is assumed in Corollary 3.6, additional compatibility conditions are required on the derivatives {Dγ v, Dγ f, Dγ g} in order to further increase the spatial regularity of u, ut and utt , even if the boundedness of time derivatives beyond utt are not needed. Additionally, if the input data is random, similar arguments as Steps 2 and 3 in Theorem 3.3 extend Corollary 3.6 to the random input data case, this time with additional compatibility conditions on the modes {vα , fα , gα }. These results are summarized in the following theorem. Theorem 3.7. Suppose u ∈ RΩ L2T H0X 1 is the weak solution of the SPDE (2.7). For fixed m ≥ 2, also assume that ¯ Ω H m+1 , ¯ Ω L2T HX m df dg ¯ Ω L2T H m−2 , v∈R X and f, g ∈ R , and , ∈R X dt dt and that the compatibility conditions (2.16) hold for {Dγ vα , Dγ fα , Dγ gk,α }, for all α ∈ J , and all indices γ = (γ1 , . . . , γd ) with |γ| ≤ m − 2. Then for the weights R′ satisfying X  2 X q′ (m+2) (m)  (2.18) qk′ λk ∨ λk reg Cm <1 and k < 1, ρk k k we have for the weak solution m+2 m−2 (2.19) u ∈ R′Ω L2T HX , ut ∈ R′Ω L2T HX m , utt ∈ R′Ω L2T HX , 23 and kukR′ 2 m+2 + kut kR′ 2 m + kutt kR′ 2 m−2 Ω LT H X Ω LT H X Ω LT H X  ≤ C kvkR¯ Ω H m+1 + kf kR¯ Ω L2 H m + kgkR¯ Ω L2 H m + kft kR¯ Ω L2 H m−2 + kgt kR¯ Ω L2 H m−2 . X T X T X T X T X 24 CHAPTER 3 Error Analysis for the Stochastic Finite Element Method In this chapter, we study numerical solutions obtained with the stochastic finite element ˙ , method applied to a parabolic SPDE driven by a multiplicative abstract noise W du dt + Au + δ W˙ (Mu) = f on D × (0, T ] (3.1) u|∂D = 0, u|t=0 = v and derive a priori error estimates for the numerical solution. Here, A is a uniformly elliptic operator and A, M take the form P Au = − i,j Di (aij (x, t)Dj u) (3.2) P Mk u = i,j Di (σkij (x, t)Dj u) with aij , σkij measurable and uniformly bounded on D. ¯ So, equation 3.1 is a special case of equation 2.7. The stochastic finite element method combines discretization procedures from the clas- sical finite element theory in numerical analysis with stochastic analysis in order to obtain computable solutions of SPDEs. The variable of randomness is discretized by a Galerkin approximation of the Wiener chaos expansion. This reduces the propagator system to a finite system of deterministic PDE that is then solved using the finite element discretiza- tion. Additionally, thanks to the lower triangular property of the propagator system, the stochastic finite element method becomes an iterative procedure of applying the finite ele- ment method to each equation in the propagator system recursively. This formulation of the stochastic finite element method for the corresponding stochastic elliptic equation has been described in [65], while the formulation for the parabolic case is essentially the same [44]. An important question to address upon formulating the numerical algorithm is to quan- tify, a priori, the error of the numerical solution. As a consequence of the discretization 25 procedures, the numerical error estimates for the elliptic and parabolic problems are com- prised of two terms. One term represents the error from the stochastic discretization, while the other term represents the numerical error from the application of the deterministic fi- nite element method to each equation in the truncated propagator system. A feature of the error estimates that carries over from the deterministic theory to the stochastic case is the connection between the spatial regularity of the solution and the order of convergence of the finite element schemes—a smoother solution yields a higher order of convergence for the part of the numerical error coming from the finite element discretization. Moreover, the Malliavin calculus approach turns out to be an indispensable framework for the error analysis of the parabolic equations, because it avails us of the two stochastic adjoint opera- tors, the Malliavin derivative and the Malliavin divergence operators, satisfying the adjoint property (2.5). This provides a tool to investigate the stochastic finite element method in a completely analogous way to the deterministic theory. The main idea brought over from the deterministic theory is the definition of the so-called Ritz projection that comes from the finite element method applied to the corresponding elliptic problem; additionally, where the notion of invoking an adjoint problem is required, the Malliavin calculus provides exactly this tool of a stochastic adjoint problem. In this sense, one may construe this error analysis to be a direct generalization of the deterministic theory to the stochastic case. However, it should be noted that extensive research on variants of finite element meth- ods, such as hp-element methods, has produced highly efficient deterministic solvers. As such, it is frequently the case that the errors incurred by the stochastic finite element method are largely dominated by the error due to the stochastic discretization, rather than the spa- tial discretization. Nevertheless, it is our hope that the techniques described in this chapter will elucidate a way of using the Malliavin calculus as a framework for direct generalization of the numerical analysis. Before proceeding, we remark on the wealth of techniques in the literature that has been developed both for the stochastic analysis and numerical analysis of SPDEs. Though the basic conception of the discretization procedure is based on the basic protocols already familiar in the algorithms for deterministic PDE—finite differences, Galerkin approxima- tions, collocation methods, finite element methods, etc.—these methods differ essentially depending on the way stochasticity is modelled in the equations. SPDEs that are essentially 26 infinite dimensional Itˆo equations (for example, equations driven by cylindrical Brownian motion) are often treated by transforming the SPDE into an infinite system of SDE and applying the techniques of the usual Itˆo calculus. Numerical simulation of this type of SPDEs often involves discretizing the SDE system by finite differences in the time com- ponent of the Brownian motion increments [13, 25, 26, 33]. Stochastic Taylor expansions are used in [32, 36] for developing and analyzing high order methods. In problems of Uncertainty Quantification, equations that depend on finite dimensional noise (i.e. pertur- bation by finite number of random variables) have enjoyed the development of polynomial chaos, generalized polynomial chaos and stochastic collocation methods over the the past decade [22, 23, 67–69]. These methods make use of the Karhunen–Lo`eve expansion, an orthogonal stochastic expansion akin to the Wiener chaos expansion, likewise treating the stochasticity in the equation as independent variables. Within the realm of finite element methods for stochastic PDE, there has been much literature on the algorithms and analysis for both elliptic and parabolic SPDE. We describe a few studies that bear some connection to our present analysis. Convergence rates of the Wiener-Itˆ o expansions of white noise and the errors from the Galerkin approximations, sans spatial discretization, have been studied in [7, 10, 63]. In [38, 39], the finite element discretization for semilinear parabolic SPDE and linear stochastic wave equation, both with additive noise in the framework of Ito calculus, was studied. For general elliptic SPDEs, the stochastic finite element or stochastic collocation methods have been studied by [1–3, 21]. The white noise functional approach using primarily the Wick product model has been studied by [8, 29, 30, 52], appealing to a similar technique of transforming the SPDE into a deterministic system of PDEs. An analysis in which the existing deterministic finite element theory is extended to the stochastic setting has been studied by [61, 70]. This idea of extending the deterministic theory turns out to be similar in spirit to our present work. 1. Review of the error estimates for the finite element method for deterministic PDE We first describe the usual finite element set up for solving deterministic PDEs, and briefly review how the error estimates for elliptic and parabolic PDE are derived. The finite element set up will be used directly as the protocol for spatial discretization in the stochastic 27 finite element method, but beyond that, we will elucidate the principles governing how the deterministic theory is developed, that will become the conception for the subsequent error analysis. 1.1. The finite element approximation. Let D be a domain in Rd with smooth boundary and let Th be a family of quasi-uniform triangulations on D. Let (Kref , P, N ) −1 be a reference finite element. For K ∈ Th , let ShK = {z : z ◦ FK ∈ P(Kref )} where FK : Kref → K is affine. The finite element space is Sh = {z ∈ H01 (D) : z|K ∈ ShK , K ∈ Th } A property of Sh we assume is that there exists r ≥ 2 such that for h small,  (3.3) inf kv − zh kL2 + hk∇(v − zh )kL2 ≤ Chs kvkH s , for 1 ≤ s ≤ r zh ∈Sh whenever v ∈ H s ∩ H01 [62]. We also assume that, in particular, Sh consists of piecewise polynomials of degree at most r − 1, so that the inverse inequality holds, k∇zh kL2 ≤ Ch−1 kzh kL2 , ∀zh ∈ Sh . We denote the finite element basis of Sh by {Φl }l=1,...,dim Sh . 1.2. The finite element method for deterministic PDE. For illustration’s sake, we consider a simple parabolic equation, the heat equation on D ut − ∆u = f, on D u|∂D = 0 u(0, ·) = w and the corresponding elliptic problem −∆U = F, on D U |∂D = 0 We will give an overview of the derivation of the error estimates for the parabolic equation, a la Thom´ee, which utilizes error estimates for the corresponding elliptic problem as well ` as utilizes properties of the adjoint problem (which in this case coincides with the elliptic problem, since ∆ is self-adjoint.) 28 The finite element formulation for the elliptic problem is to find Uh ∈ Sh such that (∇Uh , ∇χ) = (F, χ), ∀χ ∈ Sh Then the following error estimate obtains, the proof of which is well documented in the literature and is not needed for our purposes. Theorem 1.1. Assume the solution U of the elliptic problem belongs to H s for some 1 ≤ s ≤ r. Then |Uh − U | ≤ Chs kU ks and |∇Uh − ∇U | ≤ Chs−1 kU ks . Similarly, the finite element formulation for the parabolic equation is to find uh (t) ∈ Sh , t ≥ 0, such that (uh,t , χ) + (∇uh , ∇χ) = (f, χ), ∀χ ∈ Sh , t > 0, with uh (0) = wh , where wh ∈ Sh is some approximation of w. This is a semi-discrete formulation where the time variable has not been discretized. We have the error estimate as follows. Theorem 1.2. Assume the initial condition w ∈ H r , and for simplicity take wh = Rh w. For the solution u of the parabolic problem, assume that u, ut ∈ H r . Then  Z t 1/2  r 2 ′ |uh (t) − u(t)| ≤ Ch kukr + kut kr dt , ∀t ≥ 0 0 We will highlight the key ideas of the proof to illustrate how the elliptic error estimates are being used to show the parabolic error estimates. We define the elliptic or Ritz projection Rh : H01 → Sh by (∇Rh v, ∇χ) = (∇v, ∇χ), ∀χ ∈ Sh In order words, Rh is the finite element approximation operator for the corresponding elliptic problem, which maps an exact solution v of the elliptic problem to the finite element approximation vh = Rh v. Then, the approximation error can be decomposed into the sum 29 of two terms,   uh (t) − u(t) = uh (t) − Rh u(t) + Rh u(t) − u(t) = θ(t) + π(t) In an obvious way, π(t) can be directly estimated using the elliptic error estimates, and in a less direct way, so too can θ(t) be estimated. Using the definitions of the weak and FE formulations, we can compute (θt , χ) + (∇θ, ∇χ) = −(πt , χ), ∀χ ∈ Sh , t > 0, and choosing χ = θ 1 d 2 |θ| + |∇θ|2 ≤ |πt ||θ| ≤ C|πt |2 + C ′ |θ|2 2 dx and the error estimates follow from applying Gronwall’s inequality.  The last equation in the sketch of the proof uses the L2 duality estimates, |(πt , θ)| ≤ |πt ||θ|, which is a sensible choice since elliptic error estimates provide knowledge of |πt |, and which leads to an order of convergence hr matching the norms of both kut kr and kut kr . However, this manner of estimates does not exploit the structure of the regularity properties of u, ut , . . . that solutions of parabolic problems possess. An alternative estimate is to use instead the duality pairing between H −1 and H 1 , |(πt , θ)| ≤ kπt k−1 kθk1 . This requires a different set of estimates in the negative norm, but also turns out to yield a higher order of convergence. 1.3. Using the adjoint problem for negative norm estimates. The adjoint prob- lem of the elliptic problem, which in the simple model problem happens to coincide with the elliptic problem itself, is used in a duality argument to yield error estimates in negative order norms. The feature of these estimates is the give-and-take between spatial regularity and order of convergence—one can estimate the error in a lower order Sobolev space in exchange for a higher order of convergence—though such trade-off is quite typical in finite element theory. Subsequently, we will use this to improve the order of convergence of the error estimates for the parabolic problem. 30 For a nonnegative integer q, we define the spaces H −q (D) to be the dual of H q (D) with respect to the inner product in L2 (D), with duality pairing h·, ·i. The norm is hv, φi kvk−q = sup φ∈H q kφkq We have the following analogue of the error estimates for the elliptic problem. Theorem 1.3. Let U ∈ H s for some 1 ≤ s ≤ r. Then kUh − U k−q ≤ Chq+s kU ks , for 0 ≤ q ≤ r − 2. To illustrate the duality argument, we give the highlights of the proof. The negative norm in the sense of the sup norm is to be estimated; to this end, for any φ ∈ H q , consider hUh − U, φi = (Uh − U, −∆ψ) = (∇(Uh − U ), ∇ψ) The existence of ψ is granted by the solution of the adjoint problem −∆ψ = φ with ψ|∂D = 0, and moreover has the property that kψkq+2 ≤ Ckφkq for any q ≥ 0. Consequently, by orthogonality of the error to Sh , the approximation property in Sh , and the elliptic error estimates, |hUh − U, φi| = |(∇(Uh − U ), ∇(ψ − χ))| for any χ ∈ Sh ≤ CkUh − U k1 inf kψ − χk1 χ∈Sh ≤ Chs−1 kU ks · hq+1 kψkq+2 ≤ Chq+s kU ks kφkq The result follows.  As noted above, the application of the negative norm estimates is to raise the order of convergence for the parabolic problem. Theorem 1.4. Let r ≥ 3. Assume w ∈ H r and for simplicity, take wh = Rh w. Also assume the compatibility conditions that yield u(t) ∈ H r+1 , ut (t) ∈ H r−1 for a.e. t ≥ 0. Then  Z t 1/2  r 2 ′ |uh (t) − u(t)| ≤ Ch ku(t)kr + kut kr−1 dt 0 31 To prove the theorem, we decompose the error into two terms similar to the previous proof of the parabolic estimates. The difference comes in estimating the term θ(t): 1 d 2 |θ| + (∇θ, ∇θ) ≤ kπt k−1 kθk1 ≤ Ckπt k2−1 + |∇θ|2 2 dt Integrating yields the desired estimate for θ. Together with the estimate for |π(t)| ≤ Chr kukr , the result follows.  2. The stochastic finite element method formulation We will formulate the stochastic finite element method for the linear parabolic SPDE (3.1). The stochastic finite element method adopts the same idea as in the deterministic case, by elucidating a finite dimensional stochastic finite element space and casting the weak formulation of the problem into that finite dimensional setting. As with many numerical schemes for SPDE, the stochastic finite element method considered here forms the stochastic finite element space as a tensor product space of the spatial and stochastic variables, to which well-developed discretization techniques for each variable are applied separately: a finite element approximation in the spatial variable and the Galerkin approximation in the stochastic variable. In our analysis, we consider only the semi-discrete case, in that the time variable is kept continuous, thus yielding a system of ODE. The fully discrete scheme can be created by applying a suitable time stepping algorithm to the system of ODE. Finite element approximation in space. We use the usual finite element set up described in Section 1.1; that is, Sh is a finite element space on a family Th of quasi-uniform triangu- lations, with the assumption that Sh is spanned by the FE basis {Φl }l=1,...,dim Sh consisting of piecewise polynomials of degree at most r − 1. Galerkin approximation in randomness. Letting JM,n := {γ ∈ J : |γ| ≤ n, dim(γ) ≤ M }, we define the truncated Wiener chaos space n X o S M,n = f = fγ ξγ : fγ ∈ R . γ∈JM,n 32 M represents the truncation of the white noise to a finite dimension, while n represents the highest polynomial degree of the Hermite polynomials that make up the Cameron-Martin basis. SFEM formulation. The stochastic finite element method for the parabolic problem is Find uM,n h ∈ Sh ⊗ S M,n such that duM,n M,n PM h dt , zh R∓1 2 + Auh + k=1 δ ξk (Mk uM,n h ), zh R∓1 H ∓1 (3.4) Ω LX Ω X = hhf, zh iiR∓1 H ∓1 Ω X for all zh ∈ S M,n ⊗ Sh , and for every t ∈ [0, T ]. Denote the numerical solution X X dim XSh uM,n h (x, t) = u ˆγ (x, t)ξγ = u ˆγ,l (t)Φl (x)ξγ γ∈JM,n γ∈JM,n l=1 Due to 3.2, solving (3.1) via the SFEM is equivalent to solving each equation in the truncated propagator system via FEM: for α ∈ JM,n ,  dˆ u  (0) (3.5) u(0) , zh ] = hf(0) , zh i, , zh + A[ˆ dt X√ M dˆ uα   (3.6) , zh + A[ˆ u α , zh ] + uα−εk , zh ] = hfα , zh i, αk M k [ˆ dt k=1 ˆα |t=0 = (vhM,n )α . The bilinear forms A, M k are the for all zh ∈ Sh , with initial conditions u bilinear forms associated with A, Mk . Note that by our assumptions, A is coercive and M k is bounded. The algorithm. Next, we write out the SFEM algorithm explicitly to show the resulting system of ODE. We define the mass and stiffness matrices identically to the usual FEM case, and also a noise matrix arising from the stochastic term: Mmass l′ l = (Φl , Φl′ ), Mstif l′ l f = A[Φl , Φl′ ], Mnoise k;l′ l = M k [Φl , Φl′ ]. The lower triangular discrete propagator system is solved iteratively. Let the vector of coefficients of the solution vector be ~u ˆγ = (ˆ ˆγ,dim Sh )T . Then, for γ = (0), uγ,1 , . . . , u Mmass (~u ˆ(0) )t + Mstif f ~u ˆ(0) = f~(0) 33 and for |γ| ≥ 1, X√   Mmass (~u ˆγ )t + Mstif f ~u ˆγ + γk Mnoise ~u ˆγ−εk + ~gk,γ−εk = f~γ k where f~γ = (hfγ , Φ1 i, . . . , hfγ , Φdim Sh i)T , and ~gk,γ = (hgk,γ , Φ1 i, . . . , hgk,γ , Φdim Sh i)T . Remark. We remark that the stochastic finite element formulation for the parabolic problem is identical to the formulation for the corresponding elliptic problem (3.8). For the elliptic problem, Find UhM,n ∈ Sh ⊗ S M,n such that M X (3.7) AUhM,n + δ ξk (Mk UhM,n ), zh R∓1 ∓1 = hhf, zh iiR∓1 H ∓1 Ω HX Ω X k=1 for all zh ∈ S M,n ⊗ Sh , and for every t ∈ [0, T ]. In this case, the implementation of the algorithm involves defining the stiffness and noise matrices, but not the mass matrix. 3. Error analysis for SPDE with time independent operators We first study (3.1) under the assumption that A, Mk do not depend on time. In ˙ (x) is restricted to become a purely spatial noise. The main particular, the white noise W goal of this section is to show the main error estimates for the parabolic problem (3.1) (see Theorem 3.9). This will be achieved by an analogous analysis as the deterministic theory, of going through the elliptic error estimates and adjoint problem. In view of this, we will begin by studying the formal stochastic adjoint problem as well as the negative norm error estimates for the elliptic SPDE, and finally stating and deriving the parabolic estimates. 3.1. The corresponding elliptic SPDE and the formal stochastic adjoint problem. The corresponding stochastic elliptic problem is AU + δ W˙ (MU ) = F in D (3.8) U |∂D = 0 34 ¯ Ω H −1 . For non-random F , [49] has shown the unique existence of the weak where F ∈ R X 1 . For arbitrary random F , an argument by induction yields the solution U in some RΩ H0X following result. ¯ Ω H −1 . Then there exists a unique weak solution U of (3.8) Theorem 3.1. Let F ∈ R X 1 , provided the weights r 2 = qα belonging to RΩ H0X α |α|! satisfy X ellip 2 X qk (3.9) qk (λk CA ) < 1, and < 1, ρ¯k k k Moreover, we have the bounds s p X ∞ Y ellip ellip βk |β|! (3.10) kUα kH 1 ≤ CA |α|! kFα−β kH −1 (λk CA ) . 0X X β!(α − β)! β≤α k=1 We first state a result on the boundedness of the stochastic operator in the LHS of equation (3.8) that will come in handy subsequently. r ∩R H 1 , r ≥ 1, where the weights satisfy P r 2 Lemma 3.2. Let χ ∈ RΩ HX Ω 0X k qk (λk ) < ∞. Then there exists C depending only on R, A, M such that kAχ + δ W˙ (Mχ)kRΩ H r−2 ≤ CkχkRΩ HXr . X Proof. We show the lemma for r = 1, for ease of notation; the proof for r > 1 is identical. By direct computation, X ∞ X √ kAχ + δ W˙ (Mχ)k2R −1 = rα2 kAχα + αk Mk χα−εk k2H −1 Ω HX X α k=1 ∞ !2 X X √ ≤ rα2 b CA kχα kH01 + αk λk kχα−εk kH01 α k=1 ∞ !2 X X √ b 2 ≤ 2(CA ) kχk2RΩ H 1 +2 rα2 αk λk kχα−εk kH01 0X α k=1 | {z } (∗) 35 b is the constant in kAφk b 1 where CA H −1 ≤ CA kφkH 1 , for all φ ∈ H0X . To estimate (∗), we X apply Jensen’s inequality to obtain  2 X X∞  αk |α|  (∗) = rα2  √ λk kχα−εk kH01  α k=1 |α| αk αk 6=0 X qα X ∞ αk |α|2 2 ≤ λk kχα−εk k2H 1 α |α|! k=1 |α| α k 0 αk 6=0 XX q α−εk = 1{αk 6=0} qk λ2k kχα−εk k2H 1 α (|α| − 1)! 0 k ! X X X = qk λ2k rα−εk kχα−εk k2H 1 = qk λ2k kχk2RΩ H 1 0 0X k α k αk 6=0 Hence, ! X kAχ + δ W˙ (Mχ)k2R H −1 ≤2 b 2 (CA ) + qk λ2k kχk2RΩ H 1 . Ω X 0X k  Let the operators A∗ , M∗k be the formal adjoints of A, Mk , respectively. The formal stochastic adjoint problem of (3.8) is A∗ ψ + M∗ · D W˙ ψ = φ on D (3.11) ψ|∂D = 0 for φ ∈ R−1 −1 Ω HX . Although for our error estimates, we consider only self-adjoint operators A, Mk of the form (3.2), the results in this section apply to nonself-adjoint operators as well. By definition, the term M∗ · D W˙ ψ can be formally written as ∞ X  √ M∗ · D W˙ ψ α = αk + 1M∗k ψα+εk , for α ∈ J k=1 where the infinite sum is interpreted as convergent in an appropriate space. Due to the adjoint property (2.5) between D W˙ and δ W˙ , we have the adjoint property between the operators hhχ, A∗ ψ + M∗ · D W˙ ψiiRΩ H 1 −1 −1 = hhAχ + δ W˙ (Mχ), ψiiRΩ H −1 ,R−1 H 1 . 0X ,RΩ HX X Ω 0X 36 Definition 3.3. A weak solution of (3.11), with φ ∈ R−1 −1 Ω HX , is a process ψ ∈ R−1 1 Ω H0X such that hhχ, A∗ ψ + M∗ · D W˙ ψiiRΩ H 1 −1 −1 = hhχ, φiiRΩ H 1 −1 −1 0X ,RΩ HX 0X ,RΩ HX 1 . for all χ ∈ RΩ H0X ellip Note that kU ∗ kH01 ≤ CA kF kH −1 for the solution of A∗ U ∗ = F , and kM∗k φkH −1 ≤ λk kφkH01 . Proposition 3.4. Suppose there exists {ψα , α ∈ J } belonging to H01 such that for all α, ∞ X √ −1 (i) αk + 1M∗k ψα+εk ∈ HX ; k=1 ∞ X √ (ii) A∗ ψα + αk + 1M∗k ψα+εk = φα in the weak sense. k=1 Let the weights R satisfy X ellip 2 1 (3.12) qk (λk CA ) < . 2 k Then there exists C depending on R, A∗ , M∗ , such that (3.13) kψkR−1 H 1 ≤ CkφkR−1 H −1 . Ω 0X Ω X Proof. From the deterministic elliptic estimates, ! ellip X√ kψα kH 1 ≤ CA kφα kH −1 + αk + 1kM∗k ψα+εk kH −1 0X k So X ellip 2 X rα−2 kψα k2H 1 ≤ 2(CA ) rα−2 kφα k2H −1 0X α α !2 X ellip 2 X √ −1 +2 (CA ) rα αk + 1λk kψα+εk kH 1 0X α k 37 In the second term, !2 ellip 2 X √ −1 (CA ) rα αk + 1λk kψα+εk kH 1 0X k p p s !2 X |α|! |α| + 1 αk + 1 1/2 ellip = q λk CA kψα+εk kH01 q α/2 1/2 qk |α| + 1 k k ! ! X X −2 2 αk + 1 ellip 2 ≤ rα+εk kψα+εk kH 1 qk (λk CA ) 0 |α| + 1 k k and XX αk + 1 −2 rα+ε kψα+εk k2H 1 α k 0 |α| + 1 k X X βk X X βk = rβ−2 kψβ k2H 1 = r−2 kψβ k2H 1 = kψk2R−1 H 1 0 |β| |β| β 0 Ω 0X k β:βk 6=0 β k Hence, ! X  ellip 2 ellip 2 1−2 qk (λk CA ) kψk2R−1 H 1 ≤ 2(CA ) kφk2R−1 H −1 . Ω 0X Ω X k The estimate follows from the condition (3.12).  Theorem 3.5. There exists a weak solution ψ ∈ R−1 1 Ω H0X to the adjoint problem (3.11) satisfying (3.13), provided (3.12) holds. Proof. The weak solution is constructed via the usual Galerkin approach. Fix an P integer p, and let φp := |α|≤p φα ξα . We will first construct the weak solution ψ p of (3.14) A∗ ψ p + M∗ · D W˙ ψ p = φp . Let ψαp = 0 if |α| > p. For |α| = p, define ψαp by the solution of A∗ ψαp = φα . For |α| < p, ∞ X √ p A∗ ψαp = φα − αk + 1M∗k ψα+ε k . k=1 The solvability of the equation for |α| = p follows from the usual deterministic theory, and kψαp kH01 ≤ CA kφα kH −1 . P √ p The solvability of the equation for |α| < p requires that k αk + 1M∗k ψα+ε k belongs to −1 HX , which we now verify. 38 (i) Denote by Φα the quantity ∞ X i Y −2 (α + εk1 + · · · + εkj−1 )kj + 1 (Φ(i) 2 α ) = rα+ε +···+εki kφα+εk1 +···+εki k2H −1 k 1 |α| + j k1 ,...,ki =1 j=1 (i) Clearly, Φα < ∞. If |α| = p − l, for l = 1, . . . , p, it is easy to show by induction on l that l ! ellip p X kψαp kH01 ≤ CA kφα kH −1 + rα−1 (l − 1)! 2i/2 qˆi/2 Φ(i) α i=1 P ellip 2 where qˆ = k qk (λk CA ) , and hence X√ l X p rα−2 k αk + 1M∗k ψα+ε k2 −1 k H ≤ (l − 1)! 2i qˆi Φ(i) α < ∞. k i=1 P √ p P p This verifies that k αk + 1M∗k ψα+ε k ∈ H −1 , and hence ψ p := α ψ α ξα is well-defined. By construction, ψ p solves equation (3.14). Moreover, by similar calculations as Propo- sition 3.4, ! X  ellip 2 ellip 2 1−2 qk (λk CA ) kψ p k2R−1 H 1 ≤ 2(CA ) kφp k2R−1 H −1 Ω 0X Ω X k ellip 2 ≤ 2(CA ) kφk2R−1 H −1 Ω X and by (3.12), the sequence ψ p is uniformly bounded in R−1 1 Ω H0X . Thus, there exists a weakly converging subsequence, say, with abuse of notation, ψ p ⇀ ψ weakly in R−1 1 Ω H0X . 1 . From Lemma 3.2, F := Aχ + δ (Mχ) belongs to Fix an arbitrary χ ∈ RΩ H0X ˙ W −1 RΩ H X . Then hhA∗ ψ + M∗ · D W˙ ψ, χii = hhψ, Aχ + δ W˙ (Mχ)ii = lim hhψ p , F ii p→∞ = lim hhA∗ ψ p + M∗ · D W˙ ψ p , χii = lim hhφp , χii = hhφ, χii. p→∞ p→∞ By definition, the solution ψ satisfies the hypothesis of Proposition 3.4, hence the esti- mate (3.13) holds.  Remark. Higher spatial regularity results follow as usual from the corresponding de- terministic results for each equation in the propagator. In a similar fashion to the proof of 39 Theorem 3.5, one can obtain higher regularity estimates such as kψkR−1 H r ≤ CkφkR−1 H r−2 Ω X Ω X for r ≥ 1, if φ ∈ R−1 r−2 Ω HX , and if the boundary ∂D and the coefficients of A, Mk are sufficiently smooth. 3.2. Error estimates for the corresponding elliptic SPDE. An extension of [65] to random forcing terms yields the following result for the approximation error of the SFEM approximation UhM,n of equation (3.8). 1 ∩ R H m+1 , where the weights satisfy Theorem 3.6. Suppose U ∈ RΩ H0X Ω X X ellip 2 1 X qk 1 (3.15) qk (λk CA ) < , and < . 2 ρ¯k 2 k k Then the error of approximation of the stochastic finite element method is given by (3.16) kU − UhM,n kRΩ H 1 ≤ CM,n hm kU kRΩ H m+1 + CkF kR¯ Ω H −1 QM,n (R, R) ¯ 0X X X Here, CM,n can be taken as   M +n ′ CM,n =C M and the constants C, C ′ are independent of h, M, n. The term s ˆW Q Qˆ n+1 ¯ = QM,n (R, R) + ˆ 2 1−Q (1 − Q) ˆ where X ellip 2 qk X ellip 2 qk ˆ= Q qk (λk CA ) + < 1, and ˆW = Q qk (λk CA ) + . ρ¯k ρ¯k k≥1 k>M Proof. The first part of this proof closely follows the proof in [65]. Denote the numer- P ical solution by UhM,n = α∈JM,p U ˆα ξα . We decompose the approximation error into two components, X X kU − UhM,n k2RΩ H 1 = ˆα k2 1 rα2 + kUα − U H kUα k2H 1 rα2 0X X X α∈JM,p α∈J \JM,p =: I1 + I2 40 For Term I1 , we use the definitions of the weak and numerical solution for each equation in the propagator system, M X M X √ √ ˆα + AU ˆα−ǫ vh = hfα , vh i = AUα + αk M k U αk Mk Uα−ǫk vh k k=1 k=1 for all vh ∈ Sh . Note that we are assuming complete knowledge of the forcing term F . An application of the approximation techniques in the classical finite element theory yields, (see the Online Supplementary Material of [65] for details), M X ˆα k 1 ≤ CˆA inf kUα − vh k 1 + √ ˆα−ε k 1 (3.17) kUα − U H H αk Ck kUα−εk − U k H X vh ∈Sh X X k=1 b C ellip ) and C := λ C ellip . By induction, where CˆA = (1 + CA A k k A X (3.18) ˆα k 1 ≤ CˆA kUα − U cα,β inf kUβ − vh kH 1 H X X vh ∈Sh β≤α where cα,β are constants depending on α, β. The following Lemma gives a possible choice for cα,β . ~ = (C1 , C2 , . . . ). Then the constants cα,β in (3.18) may be taken Lemma 3.7. Denote C as s  |α − β|! α ~ α−β cα,β =p C . (α − β)! β Proof. This is done by induction. Suppose X ˆγ k 1 ≤ CˆA kUγ − U cγ,β inf kUβ − vh kH 1 H X X vh ∈Sh β≤γ for all |γ| ≤ n − 1, dim γ ≤ M . Let |α| = n. Then the second term on the RHS of (3.17) is M X √ ˆα−ε k 1 αk Ck kUα−εk − U k H X k=1 M X X √ = CˆA αk Ck cα−εk ,β inf kUβ − vh kH 1 vh ∈Sh X k=1 β≤α−εk s  M X X √ |α − 1 − β|! α − εk ~ α−β = CˆA αk p C inf kUβ − vh kH 1 (α − εk − β)! β vh ∈Sh X k=1 β≤α−εk 41 s  M X X |α − 1 − β|! α = CˆA p ~ α−β inf kUβ − vh k 1 (αk − βk )C HX (α − β)! β vh ∈Sh k=1 β≤α−εk αk 6=0 s  XM X |α − 1 − β|! α ≤ CˆA p ~ α−β inf kUβ − vh k 1 (αk − βk )C HX k=1 β<α (α − β)! β vh ∈Sh αk 6=0 s  M X X |α − 1 − β|! α ~ α−β = CˆA (αk − βk ) p C inf kUβ − vh kH 1 (α − β)! β vh ∈Sh X β<α k=1 αk 6=0 s  X |α − β|! α ~ α−β = CˆA p C inf kUβ − vh kH 1 (α − β)! β vh ∈Sh X β<α X = CˆA cα,β inf kUβ − vh kH 1 vh ∈Sh X β<α Hence, M X ˆα k 1 ≤ CˆA inf kUα − vh k 1 + √ ˆα−ε k 1 kUα − U H H αk Ck kUα−εk − U k H X vh ∈Sh X X k=1 X ≤ CˆA cα,β inf kUβ − vh kH 1 . vh ∈Sh X β≤α  From Lemma 3.7 and denoting the constant in (3.3) by CF E , we obtain  2 X ˆα k2 1 rα2 ≤ h2m CF2 E CˆA kUα − U 2  cα,β kUβ kH m+1 rα  H X β≤α    X r2 X ≤ h2m CF2 E CˆA 2  α 2  c rβ2 kUβ k2H m+1  rβ2 α,β X β≤α β≤α   X |α|−1 ≤ h2m CF2 E CˆA 2  2 rα−β c2α,β  kU k2R H m+1 |β| Ω X β≤α So   X X X |α|−1 ˆα k2 1 rα2 ≤ h2m CF2 E CˆA kUα − U 2 kU k2R H m+1  2 rα−β c2α,β  HX Ω X |β| α∈JM,p α∈JM,p β≤α | {z } (∗) 42 |α|−1 α To estimate (∗), since |β| β < 1 due to Lemma 1.3, X X |α|−1   |α − β|!2 α ~ 2(α−β) 2 (∗) = rα−β C |β| (α − β)! β α∈JM,p β≤α X X |α − β|! ≤ ~ 2 )α−β (q 2 C (α − β)! α∈JM,p β≤α X X |β|! = ~ 2 )β (q 2 C β! β∈JM,p α≥β α∈JM,p X |β|! = ~ 2 )β (q 2 C × (#{α ∈ JM,p : α ≥ β}) β! β∈JM,p p X X    2 ~ 2 β n!M +p 2n = (q C ) × − β! M β! n=0 |β|=n dim β≤M   p   M +p X n M +p 1 ≤ [q]≤M ≤ M M 1 − qˆ n=0 PM 2 2 where [q]≤M := k=1 qk Ck = qˆ − qˆW . This gives the first term in the RHS of (3.16). For term I2 , we recall the estimates (3.10). We decompose the sum in Term I2 into X p n−1 X X X X n ∞ X X = + .     α∈J \JM,p n=0 i=0 |α(1) |=i n=p+1 i=0 |α(1) |=i α: α: |α(2) |=n−i |α(2) |=n−i Consider the innermost sum  s 2 X X ellip 2 α  X |β|! kUα k2H 1 rα2 ≤ (CA ) q ~β kFα−β kH −1 C  X X β!(α − β)! |α(1) |=i |α(1) |=i β≤α |α(2) |=n−i |α(2) |=n−i    ellip 2 X X X |β|! −2 ~ 2β ≤ (CA ) qα  2 r¯α−β kFα−β k2H −1   r¯α−β C  X β!(α − β)! |α(1) |=i β≤α β≤α |α(2) |=n−i X X  α−β ellip 2 ~ 2 β q |β|!|α − β|! ≤ (CA ) kF k2R¯ H −1 (q C ) Ω X ρ¯ β!(α − β)! |α(1) |=i β≤α |α(2) |=n−i i X X n−i X X  α−β ellip 2 ~ 2 β q |β|!|α − β|! = (CA ) kF k2R¯ H −1 (q C ) Ω X ρ¯ β!(α − β)! k=0 l=0 |β (1) |=k |γ (1) |=i−k |β (2) |=l |γ (2) |=n−i−l 43 We introduce the notation, for ρ = (ρ1 , ρ2 , . . . ), M X ∞ X [ρ]≤M = ρk , [ρ]>M = ρk . k=1 k=M +1 Then X kUα k2H 1 rα2 X |α(1) |=i |α(2) |=n−i     k + l h q ii−k h q in−i−l n − k − l i X X n−i ellip 2 ~ 2 k ~ 2 l ≤ (CA ) kF k2R¯ H −1 [q C ]≤M [q C ]>M X Ω k ρ¯ ≤M ρ¯ >M i−k k=0 l=0   h q i i  h q i n−i ellip 2 2 n ~ 2 ~ 2 ≤ (CA ) kF kR¯ H −1 [q C ]≤M + [q C ]>M + Ω X i ρ¯ ≤M ρ¯ >M The rest of the proof proceeds identically to the proof in [65], and we obtain the second term in the RHS of (3.16).  Remark. Define the (Ritz) projection ΠM,n h ¯ Ω H 1 → Sh ⊗ S M,n as the stochas- : R 0X tic finite element approximation operator for the stochastic elliptic problem (3.8). More ¯ Ω H 1 , the projection ΠM,n U is the stochastic finite element method’s precisely, for U ∈ R 0X h solution of the elliptic SPDE (3.8), satisfying M X M X (3.19) AU + δ ξk (Mk U ), z = A(ΠM,n h U ) + δ ξk (Mk (ΠM,n h U )), z k=1 k=1 ¯ Ω H 1 , and in view for all z ∈ Sh ⊗ S M,n . Due to Lemma 3.2, F := AU + δ W˙ (MU ) ∈ R 0X of 3.6, the estimates (3.16) hold with UhM,n = ΠM,n M,n h U . This also implies that Πh is a 1 into itself. continuous linear map from RΩ H0X We will also need error estimates in the L2 (D) and H −1 (D) norms. Proposition 3.8. Under the same assumptions as Theorem 3.6, the error of approxi- mation of the SFEM has the bounds (3.20) kU − UhM,n kRΩ H 1−k ≤ CM,n hm+k kU kRΩ H m+1 + CkF kR¯ Ω H −1 QM,n (R, R) ¯ X X X for k = 1, 2. 44 Proof. As in the proof of Theorem 3.6, X X U − UhM,p = ˆα )ξα + (Uα − U Uα ξα =: e1 + e2 , α∈JM,p α∈J \JM,p with ke1 kRΩ H 1 ≤ CM,n hm kU kRΩ H m+1 , and 0X X ¯ ke2 kRΩ H 1 ≤ CkF kR¯ Ω H −1 QM,n (R, R) 0X X We leave the estimate for e2 untouched. For e1 , we consider the two cases. Case: k = 1. Let ψ ∈ R−1 2 Ω HX be the solution of Aψ + M · D W 2 ˙ ψ = R e1 , with kψkR−1 H 2 ≤ CkR2 e1 kR−1 L2 = ke1 kRΩ L2 . Note that, in fact, ψ ∈ S M,n ⊗ HX 3 also. Ω X Ω X X Then, ke1 k2RΩ L2 = hhe1 , R2 e1 iiRΩ L2 −1 2 = hhe1 , R2 e1 iiRΩ H −1 , R−1 H 1 X X , RΩ LX X Ω X = hhe1 , Aψ + M · D W˙ ψiiRΩ H −1 , R−1 H 1 X Ω X = hhAe1 + δ W˙ (Me1 ), ψ − χiiRΩ H −1 , R−1 H 1 X Ω X for all χ ∈ S M,n ⊗ Sh . So ke1 k2RΩ L2 ≤ kAe1 + δ W˙ (Me1 )kRΩ H −1 inf kψ − χkR−1 H 1 X X χ∈S M,n ⊗Sh Ω X To estimate the first term, Lemma 3.2 implies that kAe1 + δ W˙ (Me1 )kRΩ H −1 ≤ Cke1 kRΩ H 1 X 0X To estimate the second term, we make use of the FE estimate (3.3), in particular 2 1 inf kΦ − χh kH 1 ≤ ChkΦkH 2 , ∀Φ ∈ HX ∩ H0X . χh ∈Sh 0X X This FE estimate is usually obtained by finding a projection operator Ih for which kΦ − Ih ΦkH 1 ≤ Ch2 kΦkH 3 , from which the desired estimate follows immediately. But here, we 0X X will show the estimate by constructing a near-infimizing χ. Fix ǫ > 0. For each α ∈ JM,n , 45 there exists χα ∈ Sh such that kψα − χα kH 1 ≤ inf kψα − χh kH 1 + κα (ǫ) ≤ Chkψα kH 2 + κα (ǫ) 0X χh ∈Sh 0X X P P where we choose κα (ǫ) = ǫ1/2 rα κ ¯ α , with ¯ 2α ακ = 12 . Set χ = α∈JM,n χα ξα ∈ S M,n ⊗ Sh . Then X  2 kψ − χk2R−1 H 1 ≤ rα−2 Chkψα kH 2 + κα (ǫ) ≤ Ch2 kψk2R−1 H 2 + ǫ Ω 0X X Ω X α∈JM,n and inf kψ − χkR−1 H 1 ≤ ChkψkR−1 H 2 ≤ Chke1 kRΩ H 1 χ∈S M,n ⊗Sh Ω 0X Ω X 0X Hence, ke1 k2RΩ L2 ≤ ke1 kRΩ H 1 ChkψkR−1 H 2 X 0X Ω X ≤ CM,n hm+1 kU kRΩ H m+1 ke1 kRΩ L2 . X X Case: k = 2. Since e1 ∈ S M,n ⊗ H0X 1 , we compute the norm |hhe1 , φii| |hhe1 , φii| ke1 kRΩ H −1 = sup = sup X φ∈R−1 1 kφkR−1 H 1 φ∈S M,n ⊗H 1 kφkR−1 H 1 Ω H0X Ω 0X 0X Ω 0X 1 , let ψ ∈ R−1 H 3 be the solution of Aψ + M · D ψ = φ, with For any φ ∈ S M,n ⊗ H0X Ω X ˙ W kψkR−1 H 3 ≤ CkφkR−1 H 1 . Note that, in fact, ψ ∈ S M,n ⊗ HX 3 also. Then, Ω X Ω 0X hhe1 , φii = hhe1 , Aψ + M · D W˙ ψii = hhAe1 + δ W˙ (Me1 ), ψ − χii for all χ ∈ S M,n ⊗ Sh , and by a similar argument in the previous case, we have that |hhe1 , φii| ≤ kAe1 + δ W˙ (Me1 )kRΩ H −1 inf kψ − χkR−1 H 1 X χ∈S M,n ⊗Sh Ω 0X ≤ Ch2 ke1 kRΩ H 1 kφkR−1 H 1 0X Ω 0X ≤ CM,n hm+2 kU kRΩ H m+1 kφkR−1 H 1 . Ω 0X The result follows.  3.3. The parabolic error estimates. 46 Theorem 3.9. Let m ≥ 2 be an even integer. Assume for the input data ¯ Ω H m+1 , v∈R ¯ Ω L2T HX f ∈R m , ¯ Ω L2T H m−2 , ft ∈ R X X ρ¯α with weights r¯α2 = |α|! , and assume that the appropriate compatibility conditions hold, so that m+2 −1 u ∈ R′Ω L2T H0X 1 ∩ R′Ω L2T HX , ut ∈ R′Ω L2T HX ∩ R′Ω L2T HX m , m−2 utt ∈ R′Ω L2T HX , ρ′ α where the weights ρ′α 2 = |α|! are chosen using the conditions (2.15) and (2.17). Also assume, for simplicity, that the discretized initial condition is vh = ΠM,n h v. Then, for every t ∈ (0, T ], we have the error estimate for the stochastic finite element solution uM,n h (t),   (3.21) keh (t)kRΩ L2 ≤ CM,n hm+1 kut kRΩ L2 H m + ku(t)kRΩ H m+1 X T X X   + CQM,n (R, R′ ) kft − utt kR′ L2 H −1 + kf (t) − ut (t)kR′ −1 Ω T X Ω HX qα where the weights R, rα 2 = |α|! , satisfy X  2 1 X qk 1 ellip (3.22) qk λ2k CA < , and < . 2 ρ′k 2 k k Proof. Let ΠM,n h denote the stochastic finite element approximation operator for the stochastic elliptic problem (3.8). In particular, M X M X AU + δ ξk (Mk U ), z = A(ΠM,n h U ) + δ ξk (Mk (ΠM,n h U )), z k=1 k=1 for all z ∈ S M,n ⊗ Sh . The error estimates (3.16) also imply that ΠM,n h is a continuous linear 1 into itself. map from RΩ H0X Decompose the error into     eh (t) := uM,n h (t) − u(t) = u M,n h (t) − ΠM,n h u(t) + ΠM,n h u(t) − u(t) = θ(t) + π(t). 47 Analysis for π. For every t ∈ (0, T ], we have that Au(t) + δ W˙ (Mu(t)) = f (t) − ut (t) ∈ m−1 R′Ω HX . Hence the elliptic estimates (3.16) and lower norm estimates (3.20) imply kπ(t)kRΩ L2 = kΠM,n h u(t) − u(t)kRΩ L2 X X ≤ CM,n hm+1 ku(t)kRΩ H m+1 + Ckf (t) − ut (t)kR′ −1 QM,n (R, R′ ) X Ω HX provided (3.22) holds. Analysis for θ. From the definitions of the numerical and weak solutions, M X hhθt , zii + Aθ + δ ξk (Mk θ), z k=1 M X = hhf, zii − (ΠM,n h u)t , z − AΠM,n h u+ δ ξk (Mk (ΠM,n h u)), z k=1 M X = hhf, zii − (ΠM,n h u)t , z − Au + δ ξk (Mk u), z ± hhut , zii k=1 = − (ΠM,n h u − u)t , z for all z ∈ S M,n ⊗ Sh . Choosing z = R2 θ, 1d X kθk2RΩ L2 + rα2 A[θα , θα ] 2 dt X α∈JM,n M X X √ ≤ k(ΠM,n h u − u)t kRΩ H −1 kθkRΩ H 1 + αk λk rα2 kθα−ǫk kH 1 kθα kH 1 X 0X 0X 0X α∈JM,n k=1 = (I) + (II) where λk are the constants in (2.11). For (II), M X X √ (II) = αk λk rα kθα−εk kH 1 rα kθα kH 1 X X α∈JM,n k=1  !2 1/2  1/2 X M X X √ ≤ αk λk rα kθα−εk kH 1   rα2 kθα k2H 1  X X α∈JM,n k=1 α∈JM,n 48  s 2 1/2 X XM   αk |α| 1/2   ≤  λk qk rα−εk kθα−εk kH 1   kθkRΩ H 1 |α| α k X X α∈JM,n k=1 αk 6=0  s 2 1/2 M  X X αk  |α| 1/2  ≤ λk qk rα−εk kθα−εk kH 1   kθkRΩ H 1 |α| αk X X α∈JM,n k=1 αk 6=0 where we applied Jensen’s inequality in the last inequality. Continuing,  1/2 M X X   (II) ≤  2 λ2k qk rα−ε kθ α−ε k2  kθkRΩ H 1 k k H 1  X k=1 α∈JM,n αk 6=0 M !1/2 X 1/2 ≤ λ2k qk kθk2RΩ H 1 := [q~λ2 ]≤M kθk2RΩ H 1 X X k=1 PM where [q~λ2 ]≤M = 2 k=1 λk qk . Then 1d kθk2RΩ L2 + CA coerc kθk2RΩ H 1 2 dt X 0X   M,n 2 1 ~ 2 1/2 ≤ ǫ0 k(Πh u − u)t kR H −1 + + [q λ ]≤M kθk2RΩ H 1 Ω X 4ǫ0 0X coerc is the coercivity constant in (2.8). By the first condition in (3.22), we can find where CA 1/2 ǫ0 such that 1 4ǫ0 + [q~λ2 ]≤M = CA coerc . So d kθk2RΩ L2 ≤ 2ǫ0 k(ΠM,n 2 h u − u)t kRΩ HX −1 dt X and Z t kθ(t)k2RΩ L2 ≤ kθ(0)k2RΩ L2 + 2ǫ0 k(ΠM,n 2 h u − u)t (s)kR −1 ds. X X Ω HX 0 Due to our assumption on the initial condition, vh = ΠM,n h v, the term θ(0) vanishes. The estimate for the second term in the last inequality is similar to the analysis for π(t), but since the norm appears inside a time integral, it suffices to show a bound for a.e. t. Since ΠM,n h 1 into itself, it follows that (ΠM,n u) = ΠM,n u . is a continuous linear map from RΩ H0X h t h t For a.e. s ∈ (0, T ], we have that m−2 Aut (s) + δ W˙ (Mut (s)) = ft (s) − utt (s) ∈ R′Ω HX . 49 Then (3.23) k(ΠM,n M,n h u − u)t (s)kRΩ H −1 = kΠh ut − ut (s)kRΩ H −1 X X ≤ CM,n hm+1 kut (s)kRΩ HXm + Ckft (s) − utt (s)kR′ −1 QM,n (R, R′ ) Ω HX for a.e. s, and hence 2 kθ(t)k2RΩ L2 ≤ CM,n h2(m+1) kut k2RΩ L2 H m + Ckft − utt k2R′ 2 −1 QM,n (R, R′ )2 X T X Ω LT HX for all t ∈ (0, T ]. Putting together the estimates for θ(t) and π(t), we obtain   keh (t)k2RΩ L2 ≤ CM,n 2 h2(m+1) kut k2RΩ L2 H m + ku(t)k2R H m+2 X T X Ω X   + CQM,n (R, R′ )2 kft − utt k2R′ L2 H −1 + kf (t) − ut (t)k2R′ −1 Ω T X Ω HX The constant C depends only on R, A, M and the elliptic estimate constant in (3.16).  Remarks. If the discrete initial condition vh is not ΠM,n h v, additional terms will arise from approximating the initial error, but those can be subsumed into the two main terms of the error estimate. If the boundary is not smooth enough, the use of regularity estimates for the stochastic adjoint problem in the proof of Proposition 3.8 will no longer hold. In this case, the application of the lower norm estimate to the term k(ΠM,n h u − u)t (s)kRΩ H −1 is no longerX valid, but we can nonetheless obtain a convergence rate of O(hm−1 ) in the first term of (3.21). In analogy to the deterministic equation case, the finite element convergence rate of hm+1 for the solution u ∈ RΩ HT1 HX m is optimal. Without invoking the stochastic adjoint problem, it is easy to obtain a convergence rate of hm−1 for the solution u ∈ RΩ HT1 HX m, which is two orders worse than optimal. The gain of two orders is achieved by extracting some crucial information from the estimates of lower norms, through the application of the stochastic adjoint problem in the duality technique. The term QM,n (R, R′ ) in the estimate (3.21) is, as usual, the error from truncating the Wiener chaos expansion up to JM,n . It arises from invoking the error estimates for the corresponding elliptic problem, and depends on the choice of the weighted space R in 50 which to bound the error, as well as on the weights R′ of the forcing term in the sense of the elliptic problem. It also implicitly assumes that R, R′ are related by the condition (3.22). However, the second inequality in (3.22) is a somewhat strict condition. If we consider the optimal weights R′ to behave like ρ′k ∼ k −(1+ǫ) λ−2 k for any ǫ > 0, then the optimal weights R can behave like qk ∼ k −(2+ǫ) λ−2 k for any ǫ > 0. Thus, the error estimate holds in a weighted space that is generally worse than the optimal space that the solution u belongs to. Additionally, the validity of the first and third term in the RHS of (3.21) requires the −1 boundedness of utt in the HX norm. This marks the departure of the SFEM from the deterministic FEM. 4. The SFEM for SPDE with time-dependent operators In this section, we extend the results in [44] to allow the noise term, as well as the operator A, to depend on time. As discussed in the previous chapter, such time-dependent noise encompasses a variety of equations driven by an abstract noise, such as equations with space-time white noise, or equations having two independent noise terms, one purely spatial and the other purely temporal. Most steps of the analysis for the time-independent noise case carry over to the time- dependent case, except for the step estimating the norm k(ΠM,n M,n h u − u)t kRΩ H −1 , where Πh X is the elliptic SFEM approximation operator. When the noise is purely spatial (and the operators A, Mk are independent of time), the equality (ΠM,n M,n h u − u)t = Πh ut − ut allows an immediate application of the elliptic error estimates in (3.23). On the other hand, if the noise depends on time, the elliptic SFEM approximation operator ΠM,n h (t) also depends on the time parameter, and by the product rule, (ΠM,n ˙ M,n M,n h (t)u(t))t = Πh (t)u(t) + Πh (t)ut (t). Thus, to complete the error estimates, we need to derive estimates for the time derivative ˙ M,n (t). of the SFEM approximation operator, Π h We assume the operators A(t), Mk (t) are of the form (3.2), with aij , σkij ∈ HT1 WXm+2,∞ for some m ≥ 2. Then aij , σkij are Lipschitz continuous in time and their time derivatives aij ij t , (σk )t exist a.e. Because much of the analysis hinges on studying the corresponding elliptic problem, which obviously has no time evolution, we emphasize that the operators 51 A(t), Mk (t) will also be understood as being parameterized by time t. We define the time- ˙ parameterized operators A(t), M˙ k (t) by X ˙ A(t)u =− Di (aij t (x, t)Dj u), i,j X M˙ k (t)u = Di ((σkij )t (x, t)Dj u). i,j (r) (r) We recall the constants CA and λk in (2.10), (2.11), and also define C˙ A and λ˙ k to be the constants in kwkL2 (0,T ;H01 (D)) ≤ C˙ A (kw0 kL2 (D) + kf kL2 (0,T ;H −1 (D)) ) for the weak solution w of the zero Dirichlet problem dw ˙ + A(t)w = f with w(0) = w0 ; and dt in (r) ˙ k (t)wkH r−2 (D) ≤ λ˙ kwkH r (D) , kM k ∀w ∈ H r (D), t ∈ (0, T ] (1) For brevity, we write λ˙ k = λ˙ k . 4.1. The stochastic elliptic problem with time parameterized operators. For fixed t ∈ (0, T ], define the operator L(t)U := A(t)U + δ W˙ (M(t)U ). We will study the time-parameterized elliptic problem L(t)U = F in D (3.24) U |∂D = 0. The elliptic SFEM approximation (Ritz) operator ΠM,n h (t) satisfies M X M X A(t)U + δ ξk (Mk (t)U ), z = A(t)(ΠM,n h (t)U ) + δ ξk (Mk (t)(ΠM,n h (t)U )), z k=1 k=1 for all z ∈ S M,n ⊗ Sh , and all t ∈ (0, T ]. We have termed ΠM,n h (t) the “approximation” operator because ΠM,n h (t)U produces a finite approximation of U . An alternative take on the SFEM approximation operator ΠM,n h (t)U is to consider instead the SFEM solution operator ThM,n (t)F , which creates from F a finite solution of the elliptic problem. We now define ThM,n (t). 52 ¯ Ω H −1 , and for the weights R For F ∈ R e satisfying X X X q˜k (3.25) q˜k λ2k CA 2 < 1, and < 1, ρ¯k k k ¯ Ω H −1 → R define T (t) : R e Ω H 1 to be the solution operator for the equation (3.24). That X 0X is, for each t ∈ (0, T ], U = U (t) := T (t)F is the weak solution of (3.24). Define the SFEM solution operator ThM,n (t) : R ¯ Ω H −1 → S M,n ⊗ Sh by T M,n (t)F ≡ X h ΠM,n h (t)U . Then, e(t) := ΠM,n M,n h (t)U − U = Th (t)F − T (t)F, and also et (t) = (T˙hM,n (t) − T˙ (t))F. The dot ˙ stands for time differentiation, and the derivatives T˙hM,n (t), T˙ (t) are understood in the weak sense, Z T Z T T˙ (t)ϕ(t)dt = − T (t)ϕ(t)dt ˙ 0 0 for all smooth functions ϕ. ¯ Ω H −1 and U ∈ From the usual elliptic error estimates (3.16), (3.20), given F ∈ R X e Ω H r+1 , with R, R e R¯ satisfying (3.25) with 1 on the RHS of both inequalities, we have that X 2 k(ThM,n (t) − T (t))F kRe Ω H 1−k ≤ CM,n hr+k kU kRe Ω H r+1 + CkF kR¯ Ω H −1 QM,n (R, ¯ e R) X X X for all t ∈ (0, T ], and for k = 0, 1, 2. The following two propositions show that similar estimates hold for (T˙hM,n (t) − T˙ (t))F , with k = 0, 2. e ΩH 1 ∩ R Proposition 4.1. Assume U (t) ∈ R e Ω H r+1 , with the weights r˜α2 = q˜α /|α|! 0X X satisfying (3.25), (3.28), and X X (r+1) 2 (3.26) q˜k λ˙ 2k < ∞ and q˜k (λ˙ k ) < ∞, k k Let the weights R : rα2 = q α /|α|! satisfy X 1 X qk 1 (3.27) 2 qk λ2k CA < and < , 2 q˜k 2 k k 53 and X (r+1) (r+1) 2 (3.28) qk λk CA < 1, k (r+1) (r+1) where CA is the constant in kwkH r+1 ≤ CA kA−1 wkH r−1 . Then, we have the estimate k(T˙hM,n (t) − T˙ (t))F kRΩ H 1 ≤ CM,n hr kU kRe Ω H r+1 + CkF kR¯ Ω H −1 QM,n (R, R). e X X X Proof. From the definitions of ThM,n (t) and T (t), hhL(t)e, χii = 0, ∀χ ∈ S M,n ⊗ Sh Differentiating both sides w.r.t. t, ˙ hhL(t)e + L(t)et , χii = 0, ∀χ ∈ S M,n ⊗ Sh ˙ where L(t) ˙ = A(t) + δ W˙ (M(t)) is the elliptic operator obtained by differentiating the coefficients of L w.r.t. t. Consider X X√ hhLet , R2 et iiR±1 H ∓1 = rα2 Aet,α + αk Mk et,α−εk , et,α H ∓1 Ω X X α k ! X X√ ≥ rα2 coerc CA ket,α k2H 1 − αk λk ket,α−εk kH 1 ket,α kH 1 0X 0X 0X α k so XX √ coerc CA ket k2RΩ H 1 ≤ hhLet , R2 et iiR±1 H ∓1 + rα2 αk λk ket,α−ǫk kH 1 ket,α kH 1 0X Ω X 0X 0X α k = (I) + (II) For Term (I), for any χ ∈ S M,n ⊗ Sh , ˙ + Let , χii ± hhLe, (I) = hhLet , R2 et ii + hhLe ˙ R2 et ii ˙ R2 et + χii − hhLe, = hhLet , R2 et + χii + hhLe, ˙ R2 et ii   ≤ kLet kRΩ H −1 + kLek ˙ −1 inf kR2 et + χkR−1 H 1 RΩ H X X χ∈S M,n ⊗Sh Ω X ˙ + kLek 2 RΩ H −1 kR et kR−1 H 1 . X Ω X 54 For Term (II), by Cauchy-Schwartz and Jensen’s inequalities,  !2 1/2 !1/2 X X √ X (II) ≤  rα αk λk ket,α−εk kH 1  rα2 ket,α k2H 1 0X 0X α k α   1/2 X X αk |α|2 q α−εk qk ≤  λ2 ket,α−εk k2H 1  ket kRΩ H 1 α |α| αk (|α| − 1)!|α| k 0X 0X k:αk 6=0 !1/2 XX = 1{αk 6=0} qk λ2k rα−ε 2 k ket,α−εk k2H 1 ket kRΩ H 1 0X 0X α k !1/2 X = qk λ2k ket k2RΩ H 1 . 0X k Combining,  !1/2  X coerc CA − qk λ2k  ket k2 1 RΩ H0X k   ˙ ≤ kLet kRΩ H −1 + kLek inf kR2 et + χkR−1 H 1 RΩ H −1 X X χ∈S M,n ⊗Sh Ω X ˙ + kLek 2 RΩ H −1 kR et kR−1 H 1 . X Ω X ellip coerc )−1 , and from (3.27), the LHS of the last equation is strictly positive, Since CA = (CA and we obtain a valid bound for ket k2R 1 . From Lemma 3.2, (3.25) and (3.26), we have Ω H0X ˙ kLχkRΩ H −1 ≤ CkχkRΩ H 1 and kLχk 1 X RΩ H −1 ≤ CkχkRΩ H 1 , for any χ ∈ RΩ H0X . Then 0X X 0X   ket k2RΩ H 1 ≤ C ket kRΩ H 1 + kekRΩ H 1 inf ket − R−2 χkRΩ H 1 0X 0X 0X χ∈S M,n ⊗Sh X + CkekRΩ H 1 ket kRΩ H 1 0X X In the “inf” term, because any χ ∈ S M,n ⊗ Sh has only finite non-zero Wiener chaos modes, infimizing ket − R−2 χkRΩ H 1 over χ ∈ S M,n ⊗ Sh is equivalent, by a simple rescaling, to X infimizing ket − χk ˜ ∈ S M,n ⊗ Sh . Thus, upon dividing through by ket kRΩ H 1 , ˜ RΩ H 1 over χ X 0X ket kRΩ H 1 ≤ C inf ket − χkRΩ H 1 0X χ∈S M,n ⊗Sh X ! ket − χkRΩ H 1 X + CkekRΩ H 1 1+ inf 0X χ∈S M,n ⊗Sh ket kRΩ H 1 0X ≤C inf ket − χkRΩ H 1 + CkekRΩ H 1 χ∈S M,n ⊗Sh X 0X 55 since inf χ∈S M,n ⊗Sh ket − χkRΩ H 1 ≤ ket kRΩ H 1 . By translation, we have that X 0X   ket kRΩ H 1 ≤ C inf kT˙ (t)F − χkRΩ H 1 + kekRΩ H 1 . 0X χ∈S M,n ⊗Sh X 0X P Continuing, the estimation of the “inf” term is as follows. An element φ = α φ α ξα P will be decomposed into φ = φ♯ + φ⊥ , where φ♯ = α∈JM,n φα ξα . Then, inf kT˙ (t)F − χkRΩ H 1 χ∈S M,n ⊗Sh X ≤ inf k(T˙ (t)F )♯ − χkRΩ H 1 + k(T˙ (t)F )⊥ kRΩ H 1 χ∈S M,n ⊗Sh X X = (III) + (IV ) Note that U˙ (t) = T˙ (t)F solves the equation L(t)U˙ (t) = −L(t)U ˙ e ΩH 1 , (t). Since U (t) ∈ R 0X ˙ it follows from Lemma 3.2 that L(t)U e Ω H −1 . So (IV ) is the error from Wiener chaos (t) ∈ R X truncation of U˙ , and by the same proof for Term I2 in Theorem 3.6, we have that (IV ) = kU˙ (t)⊥ kRΩ H 1 ≤ CQM,n (R, R)k ˙ e L(t)U (t)kRe Ω H −1 X X e ≤ CQM,n (R, R)kU (t)kRe Ω H 1 0X e ≤ CQM,n (R, R)kF kR¯ Ω H −1 . X For (III), we use the same arguments as the proof of Proposition 3.8 to obtain inf kU˙ (t)♯ − χkRΩ H 1 ≤ Chr kU˙ (t)♯ kRΩ H r+1 χ∈S M,n ⊗Sh X X ˙ ≤ Chr kL(t)U (t)kRe Ω H r−1 X ≤ Chr kU (t)kRe Ω H r+1 . X e Ω H r−1 into RΩ H r+1 ensured The 2nd inequality follows from the boundedness of L−1 from R X X by (3.28); the 3rd inequality follows from the boundedness of L˙ from R e Ω H r+1 into R X e Ω H r−1 X ensured by (3.26). Combining (III), (IV ) with the known estimates for kekRΩ H 1 , 0X   e ket kRΩ H 1 ≤ Chr kU kRe Ω H r+1 + CQM,n (R, R)kF kR¯ Ω H −1 0X X X   r ¯ + CM,n h kU kRe Ω H r+1 + CQM,n (R, R)kF kR¯ Ω H −1 X X 56 We abuse notation again to write CM,n in place of C(1 + CM,n ). The result follows by ¯ ≤ QM,n (R, R). noting that QM,n (R, R) e  Proposition 4.2. Assume, in addition to the conditions in Proposition 4.1, that X (3) (3) 1 X ∗,(3) (3.29) qk (λk CA )2 < and qk (λ˙ k )2 < ∞. 2 k k Then k(T˙hM,n (t) − T˙ (t))F kRΩ H −1 ≤ CM,n hr+2 kU kRe Ω H r+1 + CkF kR¯ Ω H −1 QM,n (R, R). e X X X Proof. From the proof of the ket kRΩ H 1 estimate, it is clear that the first term with 0X hr is due to the norm from ke♯t kS M,n ⊗H 1 , whereas the second term is due to the norm from 0X ke⊥ t k(S M,n )⊥ ⊗H0X 1 . As usual, we leave the second term untouched, and consider only the first term. We want to estimate ke♯t kRΩ H −1 . For any φ ∈ S M,n ⊗ H0X 1 , let ψ = ψ(t) ∈ S M,n ⊗ H 3 X X be the solution of L∗ (t)ψ = φ. ˙ ♯ , χii + Since hhLe♯ , χii = 0 for all χ ∈ S M,n ⊗ Sh , differentiating w.r.t. t gives that hhLe hhLe♯t , χii = 0. Hence hhe♯t , φii = hhe♯t , L∗ ψii = hhLe♯t , ψii = hhLe♯t + Le ˙ ♯ , ψ + χii − hhLe ˙ ♯ , ψii =: (V ) − (V I) for all χ ∈ S M,n ⊗ Sh . For Term (V), |(V )| ≤ kLe♯t + Le ˙ ♯k RΩ H −1 inf kψ − χkR−1 H 1 X χ∈S M,n ⊗Sh Ω X   ≤ C ke♯t kRΩ H 1 + ke♯ kRΩ H 1 h2 kφkR−1 H 1 0X 0X Ω X ≤ CM,n hr+2 kU kRe Ω H r+1 kφkR−1 H 1 . X Ω X For Term (VI), notice that L˙ ∗ ψ ∈ R−1 1 Ω HX , so |(V I)| = |hhe♯ , L˙ ∗ ψiiR±1 H ±1 | = |hhe♯ , L˙ ∗ ψiiR±1 H ∓1 | Ω X Ω X 57 ≤ ke♯ kRΩ H −1 kL˙ ∗ ψkR−1 H 1 ≤ Cke♯ kRΩ H −1 kψkR−1 H 3 X Ω X X Ω X ≤ CM,n hr+2 kU kRe Ω H r+1 kφkR−1 H 1 . X Ω X The penultimate inequality holds due to the second inequality in (3.29). Combining, |hhe♯t , φii| ke♯t kRH −1 = sup ≤ CM,n hr+2 kU kRe Ω H r+1 . X 1 φ∈S M,n ⊗H0X kφkR−1 H 1 X Ω X  4.2. The parabolic problem with time-dependent operators. We are now in the position to prove the parabolic estimates. 1 ∩ R′ L2 H m+2 be the solution to the stochastic par- Theorem 4.3. Let u ∈ R′Ω L2T H0X Ω T X abolic equation, and assume that the conditions in Theorem 3.9 hold. For the weights R, e : r˜α2 = q˜α /|α|! satisfies (3.26) and assume that (3.27) and (3.29) hold, where R X q˜k 1 < . qk′ 2 k Then we have the estimates kuM,n h − ukRΩ L2 X   ≤ CM,n hm+1 kf kR¯ Ω L2 H m + kft kR¯ Ω L2 H m−2 + kvkR¯ Ω H m+1 T X T X X   + CQM,n (R, R′ ) kf − ut kR′ L2 H −1 + kft − utt kR′ L2 H −1 + kf (t) − ut (t)kR′ −1 . Ω T X Ω T X Ω HX Proof. We set eh (t) := uM,n h (t) − u(t)     = uM,n h (t) − T M,n h (t)(f (t) − u t (t)) + (T M,n h (t) − T (t))(f (t) − u t (t)) = θ(t) + π(t) 58 Up to the point of equation (3.23), the proof of Theorem 3.9 is followed identically to yield estimates for θ, π. Thus, kπ(t)kRΩ L2 = k(ThM,n (t) − T (t))(f − ut (t))kRΩ L2 X X ≤ CM,n hm+1 ku(t)kRΩ H m+1 + Ckf (t) − ut (t)kR′ −1 QM,n (R, R′ ) X Ω HX and Z t kθ(t)k2RΩ L2 ≤ kθ(0)k2RΩ L2 + 2ǫ0 k(ΠM,n 2 h u − u)t (s)kR −1 ds. X X Ω HX 0 To estimate k(ΠM,n h u − u)t (s)kRΩ H −1 , X d M,n (ΠM,n h u − u)t (t) = (T (t) − T (t))(f (t) − ut (t)) dt h = (T˙hM,n (t) − T˙ (t))(f (t) − ut (t)) + (ThM,n (t) − T (t))(ft (t) − utt (t)) So from Proposition 4.2, k(ΠM,n h u − u)t (t)kRΩ H −1 X e ≤ CM,n hm+3 kf (t) − ut (t)kR′Ω HXm + CQM,n (R, R)kf (t) − ut (t)kR′ −1 Ω HX + CM,n hm+1 kft (t) − utt (t)kR′ m−2 e t (t) − utt (t)k ′ + CQM,n (R, R)kf −1 Ω HX R Ω HX Putting together the estimates for θ(t) = uM,n M,n M,n h (t)−Πh u(t) and π(t) = Πh u(t)−u(t), kuM,n h − ukRΩ L2 X   m+1 2 ≤ CM,n h h kukR′ L2 H m+2 + kut kR′ L2 H m + ku(t)kR′ H m+1 Ω T X Ω T X Ω X   + CQM,n (R, R′ ) kf − ut kR′ L2 H −1 + kft − utt kR′ L2 H −1 + kf (t) − ut (t)kR′ −1 Ω T X Ω T X Ω HX as desired.  5. Numerical Simulations We perform numerical simulations in order to test the error estimates (3.21). There are several aspects of the error estimates that are worth investigation. (1) The spectral convergence in n, coming from the second term in QM,n . This spectral convergence has been shown in [65] for the corresponding elliptic equation. 59 M +n  (2) The factor CM,n . Theorem 3.6 gives an upper bound for CM,n = C ′ M . How- ever, a tighter upper bound has been conjectured1 that for CM,n that is independent of M, n. (3) The optimality of the relationship between the convergence order hm+1 and the m+1 weighted space RΩ HX . How critical is the stochastic weights in the order of convergence? 5.1. A (Simple) Model Problem. To investigate the above questions, we simulate the 1D equation ∂u ∂t = ∆u + δ W˙ (∆u) + f, for x ∈ (0, π), t ∈ (0, T ] (3.30) u(0, t) = u(π, t) = 0, u(x, 0) = v(x) ˙ (x) = P uk (x)ξk is a spatial Gaussian white noise on L2 (0, π). For the CONS in where W k L2 (0, π), we take r 1 u1 (x) = , π r 2 uk (x) = cos(k − 1)π. π In the notation of (3.1), A = −∆ and M = (M1 , M2 , . . . ), where Mk u = (uk (x)ux (x))x . Then the propagator system is ∂uα X√ = ∆uα + αk Mk uα−εk , for x ∈ (0, π), t ∈ (0, T ] ∂t k uα (0, t) = uα (π, t) = 0, uα (x, 0) = vα (x). ellip (s) The relevant constants are: CA = (1 + C poinc ) = (1 + π); the constants λk in (s) kMk ukH s−2 ≤ λk kukH s are q 1 (s) λ1 = X, λ1 ∼ k 0 q 2 (s) λk = X, λk ∼ k s−1 1The author is grateful to Zhongqiang Zhang for bringing this conjecture to her attention. 60 For the numerical simulations, we solve the SPDE on the interval (0, π) up to final time T = 0.1. We fabricate the solution of the SPDE, by fixing the solution to be u(0) (x, t) = 0 |α|! X 2 uα (x, t) = √ sin(kx)(e−k t − 1) α! k≥1 αk 6=0 This fabricated solution u is rather spatially rough, but this is all the better for demon- strating some important features of the error estimates. Based on the fabricated solution, we reverse engineer the input data of the equation; that is, the solution u is obtained with zero initial conditions and with forcing term |α|! X 2 |α|! X X αk −k2 t fα = − √ k sin(kx) − √ (e − 1)(uk (x)(sin(lx))x )x α! k≥1 α! k≥1 l≥1 |α| αk 6=0 αk 6=0 αl −δkl 6=0 For the weighted space RL2 (Ω; H s ) with weights rα2 = q α /|α|!, and for any s ∈ N0 , (3.31) kukRΩ HXs provided qk ∼ k −(2s+1+δ) . Indeed, X qα X∞ X X 2 |α|! −k2 t kuα kH s = 1{αk 6=0} q α k(e − 1) sin(kx)k2H s |α|! α! α∈J n=0 |α|=n k≥1   X∞ X X X 2 |α|! |α|! = (e−k t − 1)2 k sin(kx)k2H s  qα − 1{αk =0} q α  α! α! n=0 k≥1 |α|=n |α|=n ∞ X X 2 = q n − (¯ (e−k t − 1)2 k sin(kx)k2H s (¯ q − qk ) n ) n=0 k≥1 P where q¯ = k≥1 qk . Since q¯n − (¯ q − qk )n ≤ n¯ q n−1 qk by the mean value theorem, ∞ X X 2 C 2s kuk2RΩ H s ≤ q n−1 n¯ (e−k t − 1)2 k qk X π n=0 k≥1 s−2 s−2 Similarly, we find that ut ∈ RΩ HX , and it follows by Lemma 3.2 that f ∈ RΩ HX . The weights. We consider u(t) belonging to R′Ω L2T HX 3 . For concreteness’ sake, we P take the weights (rα′ )2 = ρα /|α|!, with ρk = ρˆk −8 and ρˆ = 0.5 < k k −8 . Then according 61 to Theorem 3.9, we measure the error in the norm RΩ L2X with weights rα2 = q α /|α|!, and qk = qˆk −10 satisfying the two conditions (3.22), X ellip 2 1 X qˆk −10 1 qˆk −10 (λk CA ) < , and < . 2 ρˆk −9 2 k≥1 k≥1 So, we take ρˆ X −2 −1 π X  −10 −1 qˆ = 0.01 < min{ k , ellip 2 k }. 2 4(CA ) k k Code Specifications. For the finite element discretization, we implemented the cG(2) method. The domain (0, π) is partitioned into M + 1 uniform subintervals of length h := π/(M + 1). The subintervals are labelled Ii = [xi−1 , xi ] for i = 1, . . . M + 1, where the grid points are xi = ih for i = 0, . . . M + 1. The nodal basis Φl (x) = φl/2 (x) is given by φi−1/2 (x) = 1Ii (x)4h−2 (xi − x)(x − xi−1 ) φi (x) = 1Ii (x)2h−2 (x − xi−1/2 )(x − xi−1 ) + 1Ii+1 (x)2h−2 (x − xi+1/2 )(x − xi+1 ) for l = 1, . . . , 2M +1. The mass, stiffness and noise matrices are symmetric 5-banded sparse matrices given by Mmass l,l′ = (Φl , Φl′ ), and Mstif l,l′ f = (Φ′l , Φ′l′ ), and Mnoise ′ ′ k;l,l′ = (uk Φl , Φl′ ). Time stepping is implemented via a 2nd order Runge Kutta method u(1) = un + ∆tL(un ) 1 n  un+1 = u + u(1) + ∆tL(u(1) ) . 2 We took ∆t/h2 = 0.03. Then the maximal order of convergence in each deterministic equation is h3 . Test of Error Estimate. The basic error estimate has 3 terms2   N + p m+1 ˆ 1/2 + C Q ˆ (p+1)/2 . C h kukRΩ H m+1 + C Q N p X ˆN = P with m = 2. In the second term, Q 2 qk ∼ N −1 , by our choice of k≥N +1 qk (λk CA ) + q˜k qk ∼ k −10 . To investigate the relationship between N and h, in the first term we couple h and N by  setting h = N −p for p′ chosen in the following way. We estimate N p+p = (N +p)...(N +1) ′ p! ∼ 2The notation of some parameters has been reshuffled. The finite stochastic subspace is now J N,p . 62 Error (and order of convergence) from finite element approximation N=4 8 12 16 20 p = 2 1.1609e-11 6.0927e-13 8.5121e-14 2.6567e-14 9.9758e-15 (4.25) (4.85) (4.05) (4.39) p = 3 7.0086e-13 9.4364e-15 6.715e-16 1.0689e-16 2.3028e-17 (6.21) (6.52) (6.39) (6.88) p = 4 3.8629e-14 1.4847e-16 4.8386e-18 4.114e-19 5.7714e-20 (8.02) (8.44) (8.57) (8.80) Table 1. Absolute errors (and convergence orders) from the finite element part under the weights qk ∼ k −10 . (Values are squared of the error norm.) Error from white noise truncation (fixed p) N=4 8 12 16 20 p = 2 1.1903e-10 2.6821e-12 1.7075e-13 1.7931e-14 2.7177e-15 p = 3 1.1906e-10 2.6828e-12 1.708e-13 1.7936e-14 2.7185e-15 p = 4 1.1906e-10 2.6829e-12 1.708e-13 1.7936e-14 2.7185e-15 Order (5.47) (6.79) (7.83) (8.46) Table 2. Truncation error (and convergence order) for each fixed p. (Values are squared of the error norm.) N +p  ′ N p , and match the powers of N in the first two terms, p hm ∼ N p−mp v.s. N −1/2 . So, 1+2p we choose p′ = 2(m+1) . Then, the error estimate reduces to 2 terms ˆ (p+1)/2 . CN −1/2 + C Q In this way, we expect the convergence to be order −1/2 in N , and spectral in p. ′ 2p+1 Error results. We set h = 2N −p , where p′ = 2(m+1) with m = 2. We vary p = 2, 3, 4 and N = 4, 8, 12, 16, 20. Clearly, the mesh size h varies with both N and p. For the fabricated solution, it is fairly straightforward to compute the exact truncation error from the stochastic truncation. Thus, we focus our results mostly on the finite element error from the JN,p part. The tables 1 and 2 show the errors from the finite element approximation of the modes in JN,p and from the truncation of the white noise to N dimensions, up to polynomial order p. P 2 2 / N,p rα kuα k 2 . The weights The truncation error for fixed p is the error from the terms α∈J L X |α|≤p were taken to be qk ∼ k −10 . Notice that in some cases, both errors are comparable, while other times the finite element error is insignificant compared to the truncation error. For the error from the JN,p part, the orders of convergence is much higher than the O(N −1/2 ) that we tried to match. The order is approximately 5/2, 7/2 and 9/2 for p = 2, 3, 4 respectively, 63 Approximate order of covergence in N 8 6 4 Order of convergence 2 p=3 p=2 0 −2 −4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 s (power decay of qk~k−s) Figure 1. Order of convergence (in N ) for the JN,p part, as the power decay of the weights qk ∼ k −s varies. (Orders are computed for the square of the error norm.) which is exactly N p higher than expected. This over-correction calls into question whether  the term N p+p ∼ N p is indeed present in the error estimate. Calculating backwards, we might conjecture the term N τ hm+1 with τ = 0 instead. This conjecture warrants further investigation. Additionally, we note that the order of convergence in the truncation of white noise is order N −9 instead of N −1 . This is due to the fact that we are estimating the error in the RΩ L2X norm, with the weights required for u to possess higher spatial regularity, 3 , and which are worse than the weights required for u to be merely L2 in space. RΩ H X X Next, we investigate the relationship between the stochastic weights and the order of convergence in the JN,p part. For the stochastic weights R(s) characterized by the decay ′ qk ∼ k −s , (3.31) relates the power s to the spatial regularity in H s of the solution in (s) ′ s , and hence also to the convergence order hm+1 . Figure 1 shows the approximate RΩ H X ′ order of convergence in N as the power s varies, also with the coupling h = 2N −p for p′ found above. Two features are observed. First, the order of convergence is higher for p = 3  than for p = 2, possibly an artifact of the over-correction by the N p+p term. Second, for s . 8, the order of convergence is almost linear with slope ≈ 1, whereas for s & 10, the order of convergence plateaus off. The plateau is due to the convergence order being limited by the h3 from the cG(2) implementation. In view of (3.31), the linear portion of the graph corroborates qualitatively the trade-off between the stochastic weights on the one hand, and the spatial smoothness and order of convergence on the other. 64 In conclusion, we have seen qualitatively that the finite element order of convergence does depend on the spatial smoothness of the stochastically weighted solution, but further  work remains to be done to investigate (the absence of) the N p+p factor. 65 CHAPTER 4 Unbiased perturbations of the Navier-Stokes equations Stochastic perturbations of the Navier-Stokes equation have received much attention over the past few decades. Among the early studies of the stochastic Navier-Stokes equa- tions are those by Bensoussan and Temam [6], Foias et al. [17–19], Flandoli [15, 16], etc. Traditionally, the types of perturbations that were proposed includes stochastic forcing by a noise term such as a Gaussian random field or a cylindrical Wiener process, and are broadly accepted as a natural way to incorporate stochastic effects into the system. The stochastic Navier-Stokes equation is underpinned by a familiar physical basis, because it can be derived from Newton’s Second Law via the the fluid flow map, using a particular assumption on the stochasticity of the governing SODE of the flow map, known as the Kraichnan turbulence. (See [53,54] and the references therein.) However, due to the nonlinearity in the equations, the stochastic Navier-Stokes equation leads to a biased perturbation; that is, the mean solution of the stochastic equation does not coincide with the solution of the unperturbed equation. In fact, the mean solution and unperturbed solution can differ quite drastically, an observation that is also true for other nonlinear equations such as the stochastic Burgers equation. To derive a model for an unbiased perturbation of the Navier-Stokes equation, the nonlinear term must be modified. In [55], a quantized stochastic Navier-Stokes equation has been proposed as an unbiased perturbation. The quantized equation replaces the usual product in the nonlinear term with the Wick product, thereby turning the nonlinear term into a stochastic convolution. The replacement with the Wick product preserves the mean because of the identity E[u ⋄ v] = Eu Ev. Thus, this perturbation is unbiased in the sense that the mean Eu of the solution satisfies the unperturbed Navier-Stokes equation. 66 Apart from the interpretation of being an unbiased perturbation, the quantized sto- chastic Navier-Stokes equation also has a physical derivation based on Newton’s Second Law. This derivation likewise representation of the fluid flow map, but differs from the aforementioned derivation by the stochasticity assumption in the governing equation of the fluid flow map. In the quantized case, it can be shown that the flow map Φ(t, x) satisfies dΦ(t) = u⋄ (t, Φ(t)) dt for smooth functions u, where the function u⋄ is the Wick version of u. However, we will not delve into the study of the fluid flow map here. Thus, we will consider the quantized stochastic Navier-Stokes equation on an open bounded domain D ∈ Rd , d = 2, 3, driven by purely spatial noise,  ˙ (x), ut + ui ⋄ uxi + ∇P f = ν∆u + f (t, x) + σ i (t, x)uxi + ∇P g + g(t, x) ⋄ W (4.1) div u ≡ 0, u(0, x) = w(x), u|∂D = 0. where the diffusivity constant is ν > 0, and the functions f, g, σ are given deterministic Rd - valued functions. Here, the driving noise W˙ (x) = P ul (x)ξl is a stationary Gaussian white k noise on L2 (D), and we assume that supl kul kL∞ < ∞. Since we restrict the forcing term to be a stationary noise, it is natural to study the related steady solution of the stationary Navier-Stokes equation,  ¯i ⋄ u u u + f¯(x) + σ ¯xi + ∇P¯ f = ν∆¯ uxi + ∇P¯ g (x) + g¯(x) ⋄ W ¯ i (x)¯ ˙ (x) (4.2) ¯≡0 div u ¯|∂D = 0. u where f¯(x), g¯(x), σ ¯ (x) are given deterministic Rd -valued functions. The analysis of the quantized Navier-Stokes equation relies on studying the propagator system. The propagator system is a lower triangular system, thus the analysis is amenable to the same induction procedure used in the earlier chapters. As a side note, we remark that in comparison, the usual stochastic Navier-Stokes equation has a propagator system that is a full system of equations which, comparatively, are a much tougher beast to tackle using the Wiener chaos expansion. Additionally, apart from the zero-th chaos mode which, 67 being the mean, solves the deterministic Navier-Stokes equation, all higher modes in the propagator system solves a linearized Stokes equation. Thus, where a result is known for the deterministic Navier-Stokes equation, it is sometimes the case that an analogous result may be shown for the quantized equation. For instance, the existence of a unique stationary solution of (4.2) requires the same condition on the largeness of the viscosity ν as does the existence of a unique steady solution of the deterministic equation (4.8a). There is substantial theory on the steady solutions of the deterministic Stokes and Navier-Stokes equations, the long time convergence of a time-dependent solution to the steady solution, as well as other dynamical behaviour of the solution. In the subsequent sections, we begin to study some of these same questions for the quantized Navier-Stokes equation, focusing on the large viscosity case where the uniqueness of steady solutions and long time convergence has been established in the deterministic setting. We will study the existence of a unique stationary solution of (4.2) as well as the existence of a unique time- dependent solution of (4.1) on a finite time interval. The Wiener chaos expansion and the propagator system will be the central tool in obtaining a generalized solution, but to place the solution in a Kondratiev space involves a useful result invoking the Catalan numbers. The Catalan numbers arises naturally from the convolution of the Wiener chaos modes in the nonlinear term. It was used to study the Wick version of the Burgers equation [34]. The long time convergence of a time-dependent solution to a steady state solution is also presented, though the theory here is not complete—the convergence in the generalized sense and in a Kondratiev space are shown separately with differing sets of assumptions. Continuing work is done to reconcile the disparity. 1. Functional analysis framework To study equations (4.1) and (4.2), we adopt the variational/weak formulation in [59, 60]. Denote the following spaces V := {v ∈ C0∞ (D)n : div v = 0} V := closure of V in the H01 (D) norm ≡ {u ∈ H01 (D) : div u = 0} H := closure of V in the L2 (D) norm | · | V ′ := dual space of V w.r.t. inner product in H 68 Also denote the norms in V and V ′ by kwkV = |∇w| and kf kV ′ , respectively. The operator1 −∆ on H, defined on the domain dom(−∆), is symmetric positive definite and thus defines a norm via |∆w|, which is equivalent to the norm kwkH 2 . For m ∈ R, define the norms |w|m = |(−∆)m/2 w| on the closed subspace Vm := dom((−∆)m/2 ) = {v ∈ H m : div v = 0} The norms |w|m and kwkH m are equivalent. We have a constant c1 so that 1 c1 k · k2H 1 ≤ |w|21 ≤ k · k2H 1 . c1 Note that |w|1 = kwkV . Denote λ1 > 0 to be the smallest eigenvalue of −∆, then we have a Poincare inequality, (4.3) λ1 |v|2 ≤ kvk2V , for v ∈ V. Define the trilinear continuous form b on V × V × V by Z b(u, v, w) = uk ∂xk v j wj dx, D and the mapping B : V × V → V ′ by hB(u, v), wi = b(u, v, w). It is easy to check that b(u, v, w) = −b(u, w, v), and b(u, v, v) = 0 for all u, v, w ∈ V . B and b have many useful properties that follow from the following lemma. Lemma 1.1 (Lemma 2.1 in [59]). The form b is defined and is trilinear continuous on H m1 × H m2 +1 × H m3 , where mi ≥ 0 and d m1 + m2 + m3 ≥ 2 if mi 6= d2 , i = 1, 2, 3, (4.4) d m1 + m2 + m3 > 2 if mi = d2 , some i. 1Technically, the correct operator is Au := −P ∆u, where P is the orthogonal projection onto H. We abuse notation here and continue writing −∆. 69 In view of Lemma 1.1, let cb be the constant in |b(u, v, w)| ≤ cb |u|m1 |v|m2 +1 |w|m3 where mi satisfies (4.4). Also let cd , d = 2, 3, be the constants in 1/2 1/2 |b(u, v, w)| ≤ c2 |u|1/2 kukV kvkV |∆v|1/2 |w| if d = 2 1/2 |b(u, v, w)| ≤ c3 kukV kvkV |∆v|1/2 |w| if d = 3 for all u ∈ V , v ∈ dom(−∆), and w ∈ H (equations (2.31-32) in [59]). Other useful consequences of Lemma 1.1 is that B(·, ·) is a bilinear continuous operator from V × H 2 → L2 , and also from H 2 × V → L2 . In order to define the weak solution of (4.1), we recall that for a smooth function p, (∇p, v) = 0 for all v ∈ V . This leads us to define the weak solution by taking the test function space V , so that the pressure term drops out. Definition 1.2. A generalized weak solution of (4.1) is a generalized random element u ∈ D′ (L2 (0, T ; H01 (D))) such that  (4.5) ˙ (x), φii hhut + ui ⋄ uxi , φii = hhν∆u + f (t, x) + σ i (t, x)uxi + ∇P g + g(t, x) ⋄ W for all test functions φ ∈ D(V ). The pressure terms. The pressure terms P f , P g are determined from the velocity field u. The term P g is defined so that the white noise multiplies a divergence-free term. Specifically, define P g to be the solution of (4.6) −∆P g = ujxk σxkj + div g, and set β := ∇P g + uxk σ k + g. Clearly, div β = 0 as desired. Then, from the equation, P f solves ˙ (x)). −∆P f = ujxk ⋄ ukxj − div f − β ⋄ (∇W Using the Wiener chaos expansion, we will study equations (4.1) and (4.2) through the analysis of the propagator system of the QsNS equations. The technique is similar to 70 Chapter 3. We recall from (2.6) that X q  α (ui ⋄ ∂xi u)α = γ (uγ , ∇)uα−γ . 0≤γ≤α Applying the above formula, we obtain the propagator system of (4.1), ∂t u0 + B(u0 , u0 ) = ν∆u0 + f (4.7a) div u(0) = 0 , u(0) (0, x) = w(x), uα |∂D = 0. P q  α ∂t uα + B(uα , u(0) ) + B(u(0) , uα ) + 0<γ<α γ B(uγ , uα−γ ) P √ g  = ν∆uα + l αl ul (x) σ i ∂xi uα−ǫl + ∇Pα−ǫ l + 1α=ǫl g (4.7b) div uα = 0 uα (0, x) = 0, uα |∂D = 0 with equality holding in V ′ . Similarly, the propagator system of (4.2) is B(¯ u0 , u u0 + f¯ ¯0 ) = ν∆¯ (4.8a) , div u ¯(0) = 0, ¯(0) |∂D = 0 u P q  α B(¯ uα , u ¯(0) ) + B(¯ u(0) , u ¯α ) + 0<γ<α γ B(¯ uγ , u ¯α−γ ) P √ g  (4.8b) = ν∆uα + l αl ul (x) σ ¯α−ǫl + ∇P¯α−ǫ ¯ i ∂ xi u l + g¯α−ǫl div u ¯α = 0, ¯α |∂D = 0 u with equality holding in V ′ . The zeroth mode u(0) = Eu is the mean of (4.1) and solves the the unperturbed Navier- Stokes equations (4.7a). 2. Stationary QSNS Given deterministic functions f¯, g¯, σ ¯ ∈ L2 (D), we seek a weak/variational solution u ¯∈ D′ (V ) (generalized random element with values in V ) satisfying  (4.9) − νhh∆¯ u, ϕii + hh¯ ¯, ϕii = hhf¯, ϕii + hh σ u i ⋄ ∂ xi u ¯ + ∇P¯ g + g¯ ⋄ W ¯ i ∂ xi u ˙ (x), ϕii for all test random elements ϕ ∈ D(V ). 71 We will first show the existence and uniqueness of a generalized strong solution. Proposition 2.1. Assume the dimension d = 2, 3. Assume f¯, g¯, σ ¯ are deterministic functions satisfying (A0) f¯, g¯, σ ¯ ∈ H, (A1) ν 2 > cb kf¯kV ′ , (A2) g¯ ∈ H 1 (D), ¯ ∈ W 1,∞ (D). σ Then there exists a unique generalized strong solution u ∈ D′ (H 2 (D) ∩ V ) of (4.2). Remark. It is interesting to note that condition (A1) in Proposition 2.1, that ensures the existence of a generalized strong solution, is the same condition that ensure the unique- ness of the strong solution of the deterministic Navier-Stokes equation. Thus, Proposition 2.1 generalizes the analogous result in the deterministic Navier-Stokes theory, which is the special subcase when g¯ = σ ¯ = 0. Proof. Solution for α = (0). The equation for u ¯0 is the deterministic stationary Navier-Stokes equation, for which the existence and uniqueness of weak solutions is well-known [59, 60]. ¯0 ∈ V of (4.8a) satisfying From (A1), there exists a unique weak solution u 1 ¯ ν (4.10) k¯ u0 kV ≤ kf kV ′ < . ν cb Moreover, since f¯ ∈ L2 (D), then u ¯0 ∈ dom(−∆), with 2 ¯ c2d |∆¯ u0 | ≤ |f | + 3/2 |f¯|3 . ν 5 ν λ1 ¯0 (·, ·). Define the bilinear continuous form a The bilinear form a ¯0 on V × V by (4.11) ¯0 (u, v) = ν(∇u, ∇v) + b(u, u a ¯0 , v) + b(¯ u0 , u, v) where u ¯0 (x) is the solution of the stationary (deterministic) Navier-Stokes equation (4.8a) just found. Also define the mapping A¯0 : V → V ′ , by hA¯0 (u), vi = a ¯0 (u, v), for all v ∈ V. 72 Then (4.8b) can be written as X q  X√ g  A¯0 (¯ uα ) = − α γ B(¯ uγ , u ¯α−γ ) + αl ul (x) σ ¯α−ǫl + ∇P¯α−ǫ ¯ i ∂ xi u l + 1α=ǫl g¯ 0<γ<α l for |α| ≥ 1. To obtain the existence and uniqueness of uα , we intend to apply the Lax-Milgram ¯0 (·, ·). To do this, we first check the coercivity of a lemma to the bilinear form a ¯0 (·, ·) on V . Lemma 2.2. Assume (A1), and assume u0 solves (4.8a) with f ∈ V ′ . Then a ¯0 (·, ·) defined in (4.11) is coercive and bounded on V . Proof. Indeed, for any v ∈ V , ¯0 (v, v) = ν|∇v|2 + b(v, u a ¯0 , v) + b(¯ u0 , v, v) ≥ ν|∇v|2 − cb k¯ u0 kV kvk2V  ¯ 2, u0 kV kvk2V = βkvk = ν − cb k¯ V where β¯ := ν − cb k¯ u0 kV > 0 by (4.10). Next, a ¯0 (·, ·) is bounded, because |¯ a0 (v, w)| ≤ νkvkV kwkV + |b(v, u ¯0 , w)| + |b(¯ u0 , v, w)|  ≤ ν + cb k¯ u0 kV kvkV kwkV for any v, w ∈ V .  We continue with the proof of Proposition 2.1. ¯ǫl ∈ V such Solutions for α = ǫl . Equation (4.8b) in variational form reduces to finding u that  a uǫl , v) = hul σ ¯0 (¯ ¯0 + ∇P¯0g + g¯ , vi =: hGǫl , vi ¯ i ∂ xi u for all v ∈ V . To apply the Lax-Milgram lemma to (4.8b), we check that the forcing term  ¯0 + ∇P¯0g + g¯ ¯ i ∂ xi u Gǫl := ul σ belongs to V ′ . In fact, we have that Gǫl belongs to L2 (D). Indeed, due to assumption (A2), σ i ∂ xi u |¯ ¯0 | ≤ k¯ u0 kV , and from (4.6) we have kP¯0g kH 2 ≤ C(k¯ σ kL∞ k¯ σ kW 1,∞ k¯ u0 kV + k¯ g kH 1 ). 73 So, from (4.10),   |Gǫl | ≤ Ckul kL∞ k¯ σ kW 1,∞ k¯ u0 kV + k¯g kH 1 ν  ≤ Ckul kL∞ k¯ σ kW 1,∞ + k¯g kH 1 cb ¯ǫl ∈ V with the By the Lax-Milgram lemma, there exists a unique variational solution u estimate 1 ν  k¯ uǫl kV ≤ ¯ Ckul kL∞ k¯ σ kW 1,∞ + k¯ g kH 1 . β cb Additionally, by a standard technique in [60], there exists P¯ǫfl ∈ L2 (D) such that (4.8b) holds in V ′ . Next, observe that by the continuity property B : V × H 2 → L2 , uǫl = Gǫl − B(¯ −ν∆¯ ¯0 ) − B(¯ u ǫl , u ¯ǫl ) ∈ L2 (D) u0 , u ¯ǫl ∈ dom(−∆), and we have the estimate Hence, u 1  |∆¯ u ǫl | ≤ |Gǫl | + |B(¯ ¯0 )| + |B(¯ u ǫl , u u0 , u ¯ǫl )| ν 1  ≤ |Gǫl | + 2cb |∆¯ u0 | k¯ uǫl kV ν C supl kul kL∞  ν  2cb  ≤ k¯σ kW 1,∞ + k¯ g kH 1 1 + |∆¯ u 0 | ν β¯ cb β¯ ¯ =K and K ¯ f¯, g¯, σ ¯ = K(ν, ¯ ) does not depend on l. Solutions for |α| ≥ 2. Denote X√ g  Gα := α l ul σ ¯α−ǫl + ∇P¯α−ǫ ¯ i ∂ xi u l , l X q  α Fα := − γ B(¯ uγ , u ¯α−γ ) 0<γ<α ¯α ∈ V such that We first find u a uǫl , v) = hFα + Gα , vi ¯0 (¯ for all v ∈ V . 74 We prove by induction. Assume we have shown the existence of a unique solution ¯γ ∈ dom(−∆) for all |γ| ≤ n − 1. By a similar argument as above, we have Gα ∈ L2 (D) u with X√ |Gα | ≤ C αl kul kL∞ k¯ σ kW 1,∞ k¯ uα−ǫl kV < ∞. l Also, since B(·, ·) is a bilinear continuous from H 2 × H 2 → L2 , we deduce that Fα ∈ L2 (D) with X q  α |Fα | ≤ cb γ |∆¯ uγ | |∆¯ uα−γ | < ∞ 0<γ<α ¯α ∈ V with the Applying the Lax-Milgram lemma, there exists a unique solution u estimates 1 uα kV ≤ ¯ (|Gα | + |Fα |). k¯ β Finally, since −ν∆¯ uα = Fα + Gα − B(¯ ¯0 ) − B(¯ uα , u ¯α ) ∈ L2 (D), u0 , u we deduce that uα ∈ dom(−∆), with 1 |∆¯ uα | ≤ (|Fα | + |Gα | + |B(¯ ¯0 )| + |B(¯ uα , u u0 , u ¯α )|) ν 1 ≤ (|Fα | + |Gα | + 2cb k¯ uα kV |∆¯ u0 |) ν 1  2cb  ≤ (|Fα | + |Gα |) 1 + ¯ |∆¯ u0 | < ∞ ν β Hence, we have found a solution u ∈ D′ (H 2 (D) ∩ V ).  Next, we find the appropriate Kondratiev space to which the solution u belongs. As described previously, the estimation of the Kondratiev norm makes use of the recursion properties of the Catalan numbers. The details of how the Catalan number rescaling is used in our estimates is put off until Section 6, though in the proofs presented before that section, we will apply the results from there. Proposition 2.3. Assume (A0-2) hold. Then there exists q0 > 2, depending on ν, f¯, g¯, σ ¯ belongs to the Kondratiev space S−1,−q (H 2 (D) ∩ V ), for q > q0 . ¯ such that u Proof. For |α| ≥ 1, we have found in the proof of Proposition 2.1 estimates for |∆¯ uα |, |∆¯ ¯ u ǫl | ≤ K 75   1 X |∆¯ u γ | |∆¯ u α−γ | X 1 √ |∆¯ ¯0  uα | ≤ B √ p + 1σ6=0 1αl 6=0 p uα−ǫl kV  . k¯ α! 0<γ<α γ! (α − γ)! l (α − ǫ l )! ¯0 depends on ν, f¯, σ where B ¯ . Let Lǫl = 1 + |∆¯ uǫl |, and Lα = √1 |∆¯ uα | for |α| ≥ 2. Then α! the rest of the proof follows in a similar way to the proof of Lemma 3 in [55]: X ¯0 Lα ≤ B Lα−γ Lγ 0<γ<α and by the Catalan numbers method in Appendix 6,   2 2 |α| ¯ 2(|α|−1) K ¯ 2|α| (4.12) |∆¯ uα | ≤ α!C|α|−1 (2N)α B 0 α for |α| ≥ 1, and the result holds with q0 satisfying ∞ X (4.13) ¯02 K B ¯ 2 25−q0 i1−q0 = 1. i=1  3. The time-dependent QSNS (4.1) In this section, we consider for simplicity equation (4.1) with σ(t, x) = 0 and, wlog, P g = 0 and div g = 0. We will consider the time-dependent solution u(t) of (4.1) on a finite time interval [0, T ] if d = 2, 3, and also study its uniform boundedness on [0, ∞) for d = 2. The former result allows an arbitrarily large time interval, thereby ensuring a global-in-time solution. On the other hand, the latter result will become useful for showing the long-time convergence of the solution to a steady state solution. For any T < ∞, it is known that a strong solution u0 (t) of the deterministic Navier- Stokes equation (4.7a) exists on the finite interval [0, T ] if d = 2, and exists on [0, (T ∧ T1 )] for a specific T1 = T1 (u0 (0)) depending on u0 (0) if d = 3. Without further conditions, we have the following result for a generalized strong solution of the quantized Navier-Stokes equation. Lemma 3.1. For d = 2, 3, let T < ∞ if d = 2, or T ≤ T1 if d = 3. Assume the forcing terms f¯, g¯ and initial condition u(0) are deterministic functions satisfying (A0′ ) f, g ∈ L2 (0, T ; H), u(0) ∈ V. 76 Then there exists a unique generalized strong solution u(t) ∈ D′ (H 2 (D) ∩ V ) for a.e. t ∈ [0, T ]. Moreover, uα ∈ C([0, T ], V ) for all α. Proof. For α = (0), it is well-known that (4.7a) has a unique solution u0 , and u0 ∈ L2 ([0, T ]; dom(−∆)), u0 ∈ C([0, T ]; V ). The bilinear form a0 (t). For t ∈ [0, T ], define the bilinear continuous form a0 (t) on V × V by a0 (u, v; t) = ν(∇u, ∇v) + b(u, u0 (t), v) + b(u0 (t), u, v) where u0 (t, x) is the solution of the time-dependent (deterministic) Navier-Stokes equations given in (4.7a) just found. Also define the mapping A0 (t) : V → V ′ , for t ∈ [0, T ], by hA0 (t)u, vi = a0 (u, v; t), for all v ∈ V. Then (4.7b) can be written as X q  X√ g  α ∂t uα + A0 (t)uα + γ B(uγ , uα−γ ) = αl ul (x) σ i ∂xi uα−ǫl + ∇Pα−ǫ l + 1α=ǫl g 0<γ<α l This is a linearized Stokes equation of the form ∂t U + A0 (t)U = F . U |∂D = 0, U (0) = w Since u0 ∈ L2 ([0, T ]; dom(−∆)), it can be shown by standard compactness techniques that if F ∈ L2 (0, T ; H) and w ∈ V , then there exists a unique strong solution U ∈ L2 (0, T ; dom(−∆)) with U ∈ L2 (0, T ; dom(−∆)), Ut ∈ L2 (0, T ; H), and U ∈ C(0, T ; V ) We prove the lemma by induction. For |α| ≥ 1, assume that uγ ∈ L2 ([0, T ]; dom(−∆)) for all γ < α. We check for the RHS of (4.7b), X q  α − γ B(uγ , uα−γ ) + 1α=ǫl ul g ∈ L2 ([0, T ]; H) 0<γ<α 77 This follows from (A0′ ) and the fact that |B(uγ , uα−γ )| ≤ cb |∆uγ | |∆uα−γ |. It follows from linear theory that there exists a unique solution uα of (4.7b) with uα ∈ L2 ([0, T ]; dom(−∆)), ∂t uα ∈ L2 ([0, T ]; H), and uα ∈ C([0, T ]; V ).  Remark. If σ 6= 0, then in addition to (A0′ ), we must require g ∈ L2 (0, T ; H 1 (D)) and σ ∈ L2 (0, T ; W 1,∞ (D)). (Compare with (A2).) Next, we study ku(t)k−1,−q;H 2 on a finite interval [0, T ] as well as the uniform bound- edness of ku(t)k−1,−q;V for all time t ∈ [0, ∞). We recall the following established result on the uniform bounds of u0 in the V and H 2 (D) norms. Lemma 3.2. (Lemma 11.1 in [59]; see also [24]) Assume for the initial condition that u0 (0, ·) ∈ V , and assume f is continuous and bounded from [0, ∞) into H f ′ is continuous and bounded from [0, ∞) into V ′ Let u0 (t) be the strong solution of the deterministic Navier-Stokes equations (4.7a), defined on [0, ∞) if d = 2, or on [0, T1 ] if d = 3. Then (4.14a) sup ku0 (t)kV ≤ c′ (ku(0) (0, ·)kV , ν, f¯, D). t≥0 (4.14b) sup |∆u0 (t)| ≤ c′′ (τ, ku(0) (0, ·)kV , ν, f¯, D). t≥τ for any τ > 0. Proposition 3.3. (i) For d = 2, 3, let [0, T ] be a finite subinterval of [0, ∞) if d = 2, or of [0, T1 (u0 (0))] if d = 3. Assume (A0′ ) and assume (A1′′ ) ν > 4cb c′ . where c′ = c′ (ku(0, ·)kV , ν, f¯, D) in (4.14a). Then there exists some q1 > 2 depending on ν, c′ , cb and T , such that for q > q1 , u ∈ S−1,−q ( L2 (0, T ; dom(−∆)) ) ∩ S−1,−q (L∞ (0, T ; V )). 78 (ii) For d = 2, assume the hypothesis of Lemma 3.2, and assume g is bounded from [0, ∞) into H. Also assume 27 c4b c′4 (A1′ ) ν4 > λ1 where c′ = c′ (ku(0, ·)kV , ν, f¯, D) in (4.14a). Then there exists q2 > 2 depending on ν, c′ and cb , such that for q > q2 , sup ku(t)k−1,−q;V < ∞. t≥0 Proof. (i) For α = (0), (4.14a) and the usual deterministic theory implies that u0 ∈ L2 (0, T ; dom(−∆)) ∩ L∞ (0, T ; V ). For |α| = 1, α = ǫl , choose in (4.7b) the test function v = (−∆)uα , 1 d kuǫ k2 + ν|∆uǫl |2 ≤ |b(uǫl , u0 , ∆uǫl )| + |b(u0 , uǫl , ∆uǫl )| + |hul g, ∆uǫl i| 2 dt l V ≤ 2cb ku0 kV |∆uǫl |2 + |ul g| |∆uǫl | ν 1 ≤ 2cb c′ + |∆uǫl |2 + |ul g|2 2 2ν So Z T Z T 1 sup kuǫl (t)k2V + (ν − 4cb c ) ′ 2 |∆uǫl | dt ≤ |ul g|2 dt 0≤t≤T 0 ν 0 By (A1′′ ), Z T 1/2 2 Lǫl := sup kuǫl (t)kV + |∆uǫl | dt 0≤t≤T 0   Z T 1/2 1 1 2 ≤√ 1+ √ |ul g| dt ≤ K1 ν ν − 4cb c′ 0 where K1 does not depend on l. For |α| ≥ 2, 1 d kuα k2V + ν|∆uα |2 2 dt X q  α ≤ |b(uα , u0 , ∆uα )| + |b(u0 , uα , ∆uα )| + γ |b(uγ , uα−γ , ∆uα )| 0<γ<α X q  α ≤ 2cb ku0 kV |∆uα |2 + γ cb kuγ kV |∆uα−γ | |∆uα | 0<γ<α 79 So Z T 1 2 ′ sup kuα (t)kV + (ν − 2cb c ) |∆uα |2 dt 2 0≤t≤T 0 X q Z T 1/2  Z T 1/2 α 2 2 ≤ cb γ kuγ kV |∆uα−γ | dt |∆uα |2 dt 0<γ<α 0 0  2 c2b X q Z T 1/2 Z T ≤  α kuγ k2V |∆uα−γ |2 dt  +ν |∆uα |2 dt 2ν γ 2 0<γ<α 0 0 RT 1/2 ˜ γ := sup0≤t≤T kuγ (t)kV + and for L |∆uγ |2 dt , 0 Z T sup kuα (t)k2V + (ν − 4cb c′ ) |∆uα |2 dt 0≤t≤T 0  2 c2b  X q α   Z T 1/2 ≤ γ sup kuγ (t)kV |∆uα−γ |2 dt  ν 0≤t≤T 0 0<γ<α  2 c2b X q  ≤  α ˜γ L L ˜ α−γ  ν γ 0<γ<α Hence,   X q cb 1  ˜α ≤ √ L 1+ √ α ˜γ L L ˜ α−γ γ ν ν − 4cb c′ 0<γ<α Let Lα = ˜ . √1 L Then α! α X L α ≤ B1 Lγ Lα−γ 0<γ<α where B1 depends on ν and c′ . By the Catalan numbers method as discussion in Appendix 6,   √ |α| |α|−1 |α| kuα kL∞ (0,T ;V ) + k∆uα kL2 (0,T ;H) ≤ α!C|α|−1 B1 K1 α and the statement of the Proposition holds with q1 satisfying ∞ X B12 K12 25−q1 i1−q1 = 1. i=1 (ii) We now show the uniform boundedness of each mode uα for all t ≥ 0. For α = (0), this is shown in the estimates of (4.14a). For |α| = 1, α = ǫl , choose in (4.7b) the test 80 function v = (−∆)uα , 1d kuǫ k2 + ν|∆uǫl |2 ≤ |b(uǫl , u0 , ∆uǫl )| + |b(u0 , uǫl , ∆uǫl )| + |hul g, ∆uǫl i| 2 dt l V 1/2 ≤ 2cb ku0 kV kuǫl kV |∆uǫl |3/2 + |ul g| |∆uǫl | ε 1 1/2 2 ≤ |∆uǫl |2 + 2cb ku0 kV kuǫl kV |∆uǫl |1/2 + |ul g| 2 2ε ε 2c2 ku0 k2V 1 ≤ |∆uǫl |2 + b kuǫl kV |∆uǫl | + |ul g|2 2 2ε ε 3 4 2 cb ku0 kV4 1 ≤ (ε)|∆uǫl |2 + 3 kuǫl k2V + |ul g|2 ε ε Taking ε = ν2 , d 27 c4 4 kuǫl k2V + ν|∆uǫl |2 ≤ 3b ku0 k2V kuǫl k2V + |ul g|2 dt ν ν and from (4.3) and (4.14a), d  27 c4 c′4  4 kuǫl k2V ≤ b 3 − νλ 1 kuǫl k2V + |ul g|2 dt ν ν 4 ≤ −βkuǫl k2V + |ul g|2 ν 27 c4b c′4  where β := − ν3 − νλ1 > 0 by (A1′ ). By Gronwall’s inequality, Z T   4 4 kuǫl (T )k2V ≤ |ul g|2 e−β(T −s) ds ≤ kul k2L∞ kgk2L∞ (0,∞;H) 1 − e−βT 0 ν νβ for any T > 0. Also, 27 c4b c′2 4 |∆uǫl (t)|2 ≤ 4 kuǫl (t)k2V + 2 |ul g(t)|2 . ν ν It follows that Lǫl := sup (kuǫl (t)kV + |∆uǫl (t)|) ≤ K2 , t≥0 for all l, where the constant K2 is independent of l and t. For |α| ≥ 2, let Lα := √1 supt≥0 (kuα (t)kV + |∆uα (t)|. Then α! 1 d kuα k2V + ν|∆uα |2 2 dt X q  α ≤ |b(uα , u0 , ∆uα )| + |b(u0 , uα , ∆uα )| + γ |b(uγ , uα−γ , ∆uα )| 0<γ<α 81 1/2 X q  α ≤ 2cb ku0 kV kuα kV |∆uα |3/2 + γ cb kuγ kV |∆uα−γ | |∆uα |. 0<γ<α By similar computations, 1 d kuα k2V + ν|∆uα |2 2 dt  2 27 c4b 4c2b X q  α ≤ ku0 k4V kuα k2V +  γ kuγ kV |∆uα−γ | ν3 ν 0<γ<α   X q    2 27 c4b 4c 2 α ≤ 3 ku0 k4V kuα k2V + b  γ sup kuγ (s)kV sup |∆uα−γ (s)|  ν ν s≥0 s≥0 0<γ<α and so  2 d 4c2b X √ kuα k2V ≤ −βkuα k2V +  α!Lγ Lα−γ  . dt ν 0<γ<α By Gronwall’s inequality and triangle inequality,  2 Z T X √ 4c2 kuα (T )k2V ≤ b  α!Lγ Lα−γ e−β(T −s)/2  ds ν 0 0<γ<α  2 4c2b  X √ Z T 1/2 ≤ α!Lγ Lα−γ e−β(T −s) ds  ν 0 0<γ<α so 1 2c2 X √ sup kuα (T )kV ≤ √ b Lγ Lα−γ α! T ≥0 νβ 0<γ<α We have also,  2 27 c4b c′4 4c2b X √ |∆uα (t)|2 ≤ kuα (t)k2V +  α!Lγ Lα−γ  ν4 ν2 0<γ<α for any t ≥ 0. Hence, it follows that X L α ≤ B2 Lγ Lα−γ 0<γ<α 82 where B2 depends on ν, c′ and cb , but is independent of t. By the Catalan method in Appendix 6,   √ |α| |α|−1 |α| sup (kuα (t)kV + |∆uα (t)|) ≤ α!C|α|−1 B2 K2 t≥0 α for |α| ≥ 1, and the statement of the Proposition holds with q2 satisfying ∞ X B22 K22 25−q2 i1−q2 = 1. i=1  4. Long time convergence to the stationary solution In this section, we study the solutions u(t, x) of (4.1) and u ¯(x) of (4.2) with σ(t, x) = ¯ (x) = 0, and for simplicity consider the case with f (t, x) = f¯(x) and g(t, x) = g¯(x). σ ¯(x) as t → ∞, first in a We study the convergence of u(t, x) to the stationary solution u weak sense (in a generalized space D′ (H)) with some exponential rate of convergence in each mode, then in a strong sense (in some Kondratiev space S−1,−q (H)) using a compact embedding argument. The latter proof, unfortunately, is does not provide a rate of conver- gence. For time-dependent f, g, similar results can be obtained under suitable assumptions, but the exponential convergence of each mode is not guaranteed. Let z(t) := u(t) − u ¯. The propagator system for z is (4.15a) z0,t + B(u0 , u0 ) − B(¯ u0 , u ¯0 ) = ν∆z0 X q   (4.15b) zα,t + A0 (t; uα ) − A¯0 (¯ uα ) = − α γ B(uγ , uα−γ ) − B(¯ uγ , u ¯α−γ ) 0<γ<α with zα (0, x) = uα (0, x) − u ¯α (x), z|∂D = 0 and div zα ≡ 0, for all α. Proposition 4.1. Let d = 2. Assume (A0), (A0′ ), (A1), and assume  λ 3/4 2 ¯ c22 1 (A3) ν ′ > |f | + 3/2 |f¯|3 c2 ν ν 5 λ1 where c2 , c′2 are specific constants depending only on D. 83 Then the solution u(t) of (4.1) converges in D′ (H) to the solution u ¯ of (4.2), D ′ (H) u(t) −→ u ¯, as t → ∞. Remark. In the following proof, all computations follow through even when d = 3. So, a similar statement to Proposition 4.1 can be made for d = 3, provided a strong solution u(t) exists in D′ (H 2 ∩ V ) for all t > 0, and the zero-th mode u0 (t) satisfies the energy inequality [59] 1d |u0 (t)|2 + νku0 (t)k2V ≤ hf¯, u0 (t)i. 2 dt Remark. If f (t, x) and g(t, x) depend on time, then an additional condition for the proposition to hold is that f (t), g(t) converge to f¯, g¯ in H. Proof. For α = (0), the convergence for the deterministic Navier-Stokes equation is well-known: if u0 (t) is any weak solution of (4.7a) with initial condition u0 (0) ∈ H, then u0 (t) −→ u ¯(0) in H as t → ∞, provided (A3) holds. Moreover, |z0 (t)| decays exponentially, (4.16) |z0 (t)| ≤ |z0 (0)| e−¯ν t , c′2 where ν¯ := νλ1 − ν 1/3 u0 |4/3 |∆¯ > 0. (See e.g., Theorem 10.2 in [59]; the positivity of ν¯ follows from the fact that |∆¯ u0 | can be majorized by the RHS of (A3).) For α = ǫl , choosing the test function v = zǫl in the weak formulation of (4.15b), 1d |zǫ |2 + νkzǫl k2V 2 dt l ≤ |b(zǫl , u ¯0 , zǫl )| + |b(zǫl , z0 , zǫl )| + |b(z0 , u ¯ǫl , zǫl )| + |b(¯ uǫl , z0 , zǫl )| u0 kV kzǫl k2V + cb kz0 kL∞ |zǫl | kzǫl kV + 2cb |∆¯ ≤ cb k¯ uǫl | |z0 | kzǫl kV c2b 2c2 u0 kV kzǫl k2V + ≤ cb k¯ kz0 k2L∞ |zǫl |2 + εkzǫl k2V + b |∆¯ uǫl |2 |z0 |2 2ε ε ¯ So, where we have used the ε-inequality in the last line with any 0 < ε < β. 1 d c2 2c2 (4.17) |zǫl |2 + (β¯ − ε)kzǫl k2V ≤ b kz0 k2L∞ |zǫl |2 + b |∆¯ uǫl |2 |z0 |2 . 2 dt 2ε ε 84 ¯ Using the Poincare inequality (4.3) and taking ε = β2 , 2 2 d ¯ 1 |zǫ |2 ≤ 2cb kz0 k2 ∞ |zǫ |2 + 8cb |∆¯ |zǫl |2 + βλ uǫl |2 |z0 |2 L dt l β¯ l β¯ For some appropriately chosen t0 ∈ (0, ∞) to be discussed next, we apply Gronwall’s in- equality, RT Z T RT ϕ(t)dt |zǫl (T )|2 ≤ e t0 |zǫl (t0 )|2 + ψl (s)e s ϕ(t)dt ds t0 where 4c2 ¯ 1, ϕ(t) = ¯b kz0 (t)k2L∞ − βλ β 8c2 uǫl |2 |z0 (t)|2 . ψl (t) = ¯b |∆¯ β β¯2 λ1 The t0 is chosen large enough so that kz0 (t)k2L∞ < 4c2b whenever t ≥ t0 . Such t0 exists, because by (4.14b) and the Sobolev embedding z0 (t) ∈ C 1/2 is H¨older continuous with exponent γ < 1 and supt≥τ kz0 (t)kC γ ≤ c′′ is uniformly in t. Then due to (4.16), we deduce that in fact z0 (t, ·) −→ 0 uniformly on D as t → ∞. Consequently, we have that supt≥t0 ϕ(t) < 0. Set ϕ¯ > 0 satisfying   2ϕ¯ < min − sup ϕ(t), 2¯ ν . t≥t0 RT  Obviously, exp t0 ϕ(t)dt ≤ exp − 2ϕ(T ¯ − t0 ) . Moreover, from (4.16), 8c2 uǫl |2 |z0 (t0 )|2 e−2¯ν (t−t0 ) =: Cψl e−2¯ν (t−t0 ) −→ 0 |ψl (t)| ≤ ¯b |∆¯ β decays exponentially as t → ∞. Combining these results, Z T 2 −2ϕ(T ¯ −t0 ) 2 |zǫl (T )| ≤ e |zǫl (t0 )| + Cψl e−2¯ν (s−t0 ) e−2ϕ(T ¯ −s) ds t0 Cψl  −2φ(T ¯ −t0 ) −2¯  ≤ e−2ϕ(T ¯ −t0 ) |zǫl (t0 )|2 + e e ν (T −t0 ) −→ 0 ν − ϕ) 2(¯ ¯ as T → ∞. (In the first term, |zǫl (t0 )|2 has been shown to be finite for any finite t0 .) Since ϕ¯ < ν¯,   Cψl (4.18) |zǫl (T )|2 ≤ |zǫl (t0 )|2 + e−2ϕ(T ¯ −t0 ) =: Kǫ2l e−2ϕ(T ¯ −t0 ) ν − ϕ) 2(¯ ¯ 85 for T ≥ t0 . Kǫl does not depend on T . For |α| ≥ 2, we prove by induction. Fix α, and assume the induction hypothesis that: For each 0 < γ < α, for T ≥ t0 , 1−|γ| ϕ(T (4.19) |zγ (T )| ≤ Kγ e−2 ¯ −t 0) −→ 0 as T → ∞, where Kγ does not depend on T . We want to show that (4.19) also holds for α. From (4.15b) with test function v = zα , 1d |zα |2 + ν|∇zα |2 2 dt ≤ |b(zα , u ¯0 , zα )| + |b(zα , z0 , zα )| + |b(z0 , u ¯α , zα )| + |b(¯ uα , z0 , zα )| X q   α + γ |b(zγ , zα−γ , zα )| + |b(zγ , u¯α−γ , zα )| + |b(¯ uγ , zα−γ , zα )| 0<γ<α ¯ Similar to (4.17), using the ε-inequality with any 0 < ε < β/2, 1d c2 2c2 |zα |2 + (β¯ − 2ε)kzα k2V ≤ b kz0 k2L∞ |zα |2 + b |∆¯ uα |2 |z0 |2 2 dt 2ε ε c 2  X q α  2 + b γ kz α−γ kV + 2k¯ u α−γ kV |z γ | 1/2 4ε 0<γ<α ¯ Using the Poincare inequality and taking ε = β/4, d  4c2  16c2b |zα (t)|2 ≤ b kz k2 0 L ∞ − λ 1 ¯ β |z α | 2 + |∆¯uα |2 |z0 |2 dt β¯ β¯ 2c2  X q α  1/2  1/2 2 + ¯b γ kz α−γ V k + 2k¯u α−γ Vk |z γ | kz k γ V β 0<γ<α ≤ ϕ(t)|zα (t)|2 + ψα (t) where now 16c2b ψα (t) = ¯ uα |2 |z0 (t)|2 |∆¯ β    2c2b  X q α 2 X q  α + ¯ γ uα−γ kV kzγ (t)kV   kzα−γ (t)kV + 2k¯ γ |zγ (t)|  β 0<γ<α 0<γ<α 86 From the hypothesis (4.19),  X q  −|γ| 2ϕ(t−t  |ψα (t)| ≤ Cψα e−2¯ν (t−t0 ) + C˜ψα α γ Kγ e−2 ¯ 0) 0<γ<α where 16c2b Cψα = uα k2H 2 |z0 (t0 )|2 , k¯ β¯   2c2b  X q α  2 C˜ψα = ¯ γ uα−γ kV sup kzγ (s)kV  , sup kzα−γ (s)kV + 2k¯ β 0<γ<α s≥0 s≥0 and Cφα , C˜φα do not depend on t. By Gronwall’s inequality, Z T |zα (T )|2 ≤ e−ϕ(T ¯ −t0 ) |zα (t0 )|2 + ψα (s)e−ϕ(T ¯ −s) ds t0 Cψα X q  e−2 1−|γ| ϕ(T ¯ −t0 ) ≤ e−ϕ(T ¯ −t0 ) |zα (t0 )|2 + e−2ϕ(T ¯ −t0 ) + C˜ψα α γ Kγ ν − ϕ) 2(¯ ¯ 1 − 2−|γ| 0<γ<α −21−(|α|−1) ϕ(T ≤ Kα2 e ¯ −t 0) where Kα does not depend on T . Hence, 1−|α| ϕ(T (4.20) |zα (T )| ≤ Kα e−2 ¯ −t 0) for all T ≥ t0 . It follows that (4.19) holds also for α, and the result follows.  We proceed to deduce the long time convergence of u(t) in some Kondratiev space S−1,−q (H). The manner of estimates in Proposition 4.1 is not directly suited for applying the Catalan numbers method. Instead, we will use a compact embedding type argument in the following lemma to show the result. Lemma 4.2. Let uk ∈ S−1,−q (V ) be a sequence satisfying X rα   sup kukα k2V < ∞, α α! k that is, satisfying {uk } ∈ S−1,−q (ℓ∞ (V )). ˜ Then there exists a subsequence k˜N such that ukN converges in D′ (H) to some u ¯ ∈ D′ (H). Furthermore, if u ¯ ∈ S−1,−q (H), then the convergence is in S−1,−q (H). 87 Proof. The convergence in D′ (H) will follow easily from the fact that V is compactly embedded in H. Let JN = {α ∈ J : |α| ≤ N, and αi = 0 for i > N }. Since supk kuk0 kV < ∞, there exists a subsequence {kj0 }∞ k j=1 such that ku0 − u ¯0 kH → 0 for some u ¯0 ∈ H. N −1 ∞ Iteratively, for each N , there exists further subsequences {kjN }∞ j=1 ⊂ {kj }j=1 such that for every α ∈ JN , kukα − u ¯α kH → 0 ¯α ∈ H. In particular, for each N , we can find jN such that for some u kjN ¯α kH ≤ N −1 , kuα N − u for all α ∈ JN . P Consequently, choose the subsequence k˜N := kjNN and we have found the limit u ¯= αu ¯ α ξα . ˜ It follows that ukN → u ¯ in D′ (H). ¯ ∈ S−1,−q (H). Let ε > 0 be arbitrary. For any N , Now suppose u ˜ X rα ˜ X rα ˜ ¯k2−1,−q;H = kukN − u kukN − u ¯k2H + kukN − u ¯k2H = (I) + (II) α! α! α∈JN α∈J / N By our special choice of k˜N , there exists NI such that X rα ε (I) ≤ N −2 < whenever N > NI . α! 2 α∈JN From the hypothesis of the lemma, there exists NII such that X rα   X rα ε (II) ≤ 2 sup kuk k2V + 2 uk2H < k¯ whenever N > NII . α! k α! 2 α∈J / N α∈J / N ˜ Thus, kukN − u ¯k2−1,−q;H < ε whenever N > max{NI , NII }.  The hypothesis in Lemma 4.2 is stronger than requiring uk ∈ l∞ (S−1,−q (V )), thus it is a weaker statement of what might be construed as a compact embedding result for Kondratiev spaces. It is not shown whether S−1,−q (V ) is compactly embedded in S−1,−q (H). Nonetheless, it is sufficient for our purposes. Corollary 4.3. Let d = 2. Assume the hypotheses of Propositions 2.3 and 3.3(ii). Then, for the solutions u(t) and u ¯ of (4.1), (4.2), we have that u(t) −→ u ¯ in S−1,−q (H), as t → ∞, 88 for q > max{q0 , q2 }, where q0 , q2 are the numbers from Propositions 2.3, 3.3. Proof. In the proof of Proposition 3.3, we have in fact shown that u(t) belongs to the space S−1,−q (L∞ ([0, ∞); V )). Taking any sequence of times, tk → ∞, the sequence {u(tk )} satisfies the hypothesis of Lemma 4.2. So, there exists a subsequence of u(tk ) converging in S−1,−q (H) to u ¯. This is true for any sequence {tk }, hence u(t) −→ u ¯ in S−1,−q (H) as t → ∞.  5. Finite Approximation by Wiener Chaos Expansions In this section, we study the accuracy of the Galerkin approximation of the solutions of the quantized stochastic Navier-Stokes equations. The goal is to quantify the conver- gence rate of approximate solutions obtained from a finite truncation of the Wiener chaos expansion, where the convergence is in a suitable Kondratiev space. In relation to being a numerical approximation, quantifying the truncation error is the first step towards un- derstanding the error from the full discretization of the quantized stochastic Navier-Stokes equation. In what follows, we will consider the truncation error estimates for the steady solution (2N)−qα u|: for rα2 = ¯. Recall the estimate (4.12) for |∆¯ u α! , with q > q0 , we have   |α| rα2 |∆¯ uα | 2 ≤ 2 C|α|−1 (2N)(1−q)α B0−2 (B0 K)2|α| . α This estimate arose from the method of rescaling via Catalan numbers, and will be the estimate we use for the convergence analysis. For the time-dependent equation, similar analysis can be performed using the analogous Catalan rescaled estimate, and will not be shown. Let JM,P = {α : |α| ≤ P, dim(α) ≤ M }, where M, P may take value ∞. The projection P of u ¯M,P = α∈JM,P u ¯ into span{ξα , α ∈ JM,P } is u ¯ α ξα . ¯M,P can be written as ¯−u Then the error e = u X |∆e|2 = rα2 |∆¯ uα | 2 α∈J \JM,P ∞ X X = rα2 |∆¯ uα | 2 + rα2 |∆¯ uα | 2 |α|=P +1 {|α|≤P, |α≤M |<|α|} 89 ∞ X P |α|−1 X X X = rα2 |∆¯ uα | 2 + rα2 |∆¯ uα | 2 |α|=P +1 |α|=1 i=0 |α≤M |=i | {z } | {z } (IV ) (I) | {z } (II) | {z } (III) We define the following values ∞ X ˆ := 21−q B02 K 2 Q i1−q , i=1 M X ∞ X ˆ ≤M := 21−q B0 K 2 Q i1−q , ˆ >M := 21−q B0 K 2 Q i1−q . i=1 i=M +1 ˆ >M decays on the order of M 2−q . In particular, the term Q We proceed to estimate the terms (I)-(IV), by similar computations to Wan et al. For fixed 1 ≤ p ≤ P , |α| = p, and fixed i < p, X   |α| (I) ≤ 2 Cp−1 B0−2 (2N)(1−q)α (B0 K)2p α |α≤M |=i, |α>M |=p−i   2 −2 p ˆ i ˆ p−i = Cp−1 B0 Q≤M Q >M i Then for fixed 1 ≤ p ≤ P , |α| = p, p−1 X p−1   X p ˆ i ˆ p−i (II) = (I) ≤ 2 Cp−1 B0−2 Q≤M Q>M i i=0 i=0 2 = Cp−1 B0−2 (Q ˆp ) ˆp − Q ≤M And finally, P X P X (III) = (II) ≤ 2 Cp−1 ˆp ) ˆp − Q B0−2 (Q ≤M |α|=1 p=1 P 1 ˆ ˆ ≤M ) + 1 X 24p ˆp ) ˆp − Q ≤ ( Q − Q (Q ≤M B02 16πB02 p=2 (p − 1)3 Since Q ˆ p ≤ pQ ˆp − Q ˆ p−1 (Q ˆ−Q ˆ ≤M ) by the mean value theorem for x 7→ xp , ≤M XP ˆ p−1 1 ˆ 1 ˆ p24p Q (III) ≤ Q >M + Q >M B02 16πB02 p=2 (p − 1)3 90 XP ˆ p−1 1 ˆ 1 ˆ p(24 Q) ≤ 2 Q>M + Q >M B0 πB02 p=2 (p − 1)3 P X −1 1 ˆ ˆ p ≤ Q >M (24 Q) B02 p=0 To estimate Term (IV ), ∞ X X   |α| (IV ) ≤ 2 Cp−1 B0−2 (21−q B02 K 2 )p (N)(1−q)α α p=P +1 |α|=p ∞ X X p = B0−2 2 Cp−1 (21−q B02 K 2 )p i1−q p=P +1 i≥1 ∞ X ˆ P +1 24(p−1) ˆ p 1 (24 Q) ≤ B0−2 Q ≤ π(p − 1)3 16πB02 1 − 24 Q ˆ p=P +1 Putting the estimates together,  ˆ P +1 + M 2−q |∆e|2 ≤ C (24 Q) ˆ < 1 in (4.13), which ensured summability of the weighted norm Notice the condition 24 Q of the solution, is of course a required assumption for the convergence of the error estimate. 6. The Catalan numbers method The Catalan numbers method was used in the preceding sections to derive estimates for the norms in Kondratiev spaces. This method was previously described in [34, 55], but we restate it here just for the record. Lemma 6.1. Suppose Lα are a collection of positive real numbers indexed by α ∈ J , satisfying X Lα ≤ B Lγ Lα−γ . 0<γ<α Then  Y |α|−1 |α| Lα ≤ C|α|−1 B Lαǫii α i for all α, where Cn are the Catalan numbers. 91 Proof. The result is clearly true for α = ǫi . By induction, let |α| ≥ 2, and suppose the result is true for all γ < α. Then    X |γ| |α − γ|  Y αi  Lα ≤ C|γ|−1 C|α−γ|−1 B |α|−1 L ǫi γ α−γ 0<γ<α i |α|−1 X X n! (|α| − n)! |α|−1  Y αi  = Cn−1 C|α|−n−1 B L ǫi γ! (α − γ)! n=1 0<γ<α i |γ|=n |α|−1 X X |α|−1 α |α|! Y  = Cn−1 C|α|−n−1 B |α|−1 Lαǫii n γ α! n=1 0<γ<α i |γ|=n | {z } (∗) We claim that (∗) = 1, for any α and any n < |α|. Indeed, let Kα = (k1 , . . . , k|α| ) be the characteristic set of α. Each summand in (∗) is  −1 |α|! n! (|α| − n)! α! γ! (α − γ)! |α|! n! (|α|−n)! The term α! is the number of distinct permutations of Kα , whereas the term γ! (α−γ)! is the number of distinct permutations of Kα where only Kγ , Kα−γ has been permuted within themselves. On the other hand, the latter term is the number of distinct permutations of Kα corresponding to a particular γ, where the correspondence of a permutation of Kα to a γ ∈ {γ : 0 < γ < α, |γ| = n} can be made by taking Kγ to be the first n entries of that permutation of Kα . Thus, each summand in (∗) is the relative frequency of γ over all distinct permutations of Kα , and hence their sum must equal 1. To complete the proof, using the recursion property of the Catalan numbers, |α|−1 X   Y |α|! Lα ≤ Cn−1 C|α|−n−1 B |α|−1 Lαǫii α! n=1 i   Y |α|! = C|α|−1 B |α|−1 Lαǫii . α! i  If Lα satisfies the hypothesis of Lemma 6.1, and if Lǫi ≤ K for all i, then for r = (2N)−q , X X   |α| rα L2α ≤ 2 Cn−1 B 2(|α|−1) K 2|α| (2N)(1−q)α α |α|=n |α|=n 92   n X |α| (1−q)α = B −2 Cn−1 2 B 2 K 2 21−q N α |α|=n  X ∞ n 2 1−q n = B −2 Cn−1 2 2 B K 2 i(1−q) i=1 2 2n For large n, the Catalan numbers behave asymptotically like Cn ∼ √πn 3/2 . Hence, the sum P∞ P α 2 n=0 |α|=n r Lα converges for any q > max{q0 , 2}, where q0 satisfies ∞ X 2 2 5−q0 B K 2 i(1−q0 ) = 1. i=1 93 CHAPTER 5 Randomization of Incoherent Forcing for Improvement of Energy Approximations 1. Introduction In this chapter, we consider a linear SPDE ∂ ˙ Q (x), (5.1) v = Av + W x ∈ U, t > 0, ∂t and a system of deterministic PDEs ∂ (5.2) vi (x, t) = Avi (x, t) + ρi ei (x), x ∈ U, t > 0, for i = 1, 2, . . . , ∞, ∂t where U ⊂ Rd is an open bounded domain, A is a linear partial differential operator, ˙ (x) is a weighted spatial noise, given {ei , i ≥ 1} is an orthonormal basis in L2 (U ), and W by X ˙ (x) = W σi ei (x)ξi i≥1 with {ξi , i ≥ 1} being a set of independent Gaussian random variables and {σi , i ≥ 1} ˙ (x) is a standard spatial being a set of nonnegative weights (see (2.3)). If all σi = 1, W white noise; this case is presented in [43]. We assume that the initial conditions in (5.1) and (5.2) are zero. In fact, we recall from Definition 1.5, and the discussion therein, that (5.1) and (5.2) are equivalent in that X   vi (x, t) = E[v(x, t)ξi ] and v(x, t) = vi (x, t) W˙ , ei . L2 (U ) i≥1 The equivalence of (5.1) and (5.2) is a very simple implication of the Wiener chaos ex- pansion for SPDEs. System (5.2) is the propagator system for (5.1). Under very general assumptions, a solution of one of the two equations exists and is unique if and only if the P other has a unique solution (see [49] and Theorem 3.2). Moreover, if i≥1 σi2 < ∞, then 94 P 2 the solution is in L2 , and if i≥1 σi = ∞, then the solution is found in a Sobolev space with a negative index. The energy of a solution u of (5.1) is defined by X (5.3) E[v(t)] := Ekv(·, t)k2L2 (U ) = kvi (·, t)k2L2 (U ) . i≥1 Clearly, it is independent of the choice of the basis {ei , i ≥ 1} . Our main goal is to identify suitable bases {ei , i ≥ 1} as well as estimators vˆ(n) (x, t) = Pn i=1 vi (x, t)σi ξi such that the energy of vˆ(n) (x, t) efficiently approximates E[v(t)]. For a ˙ N (x), we want to study the behavior of the estimators as finite N -dimensional noise W N → ∞. Getting a little bit ahead of the story, we remark that, while the energy E[v(t)] does not depend on the choice of the basis, the rate of convergence of the approximate energy Pn 2 i=1 kvi (·, t)kL2 (U ) does and, sometimes, does so quite substantially. Approximating the energy kv(·, t)k2L2 (U ) for system (5.2), and similar systems, requires solving a large number of PDEs that differ only by the forcing terms. For example, the problem of efficient approximation of the energy comes up in the modeling of wave propaga- tion with incoherent sources [40], which appear in a wide range of problems in optics, such as those related to diffuse light [71]. Some popular examples include the Raman photonic crystal spectrometer [51], which is used to measure spatially incoherent light in environ- mental and biological sensing, as well as fluorescent or bioluminescent tomography [66], which has been used successfully to achieve in-vivo functional imaging in cancer research and drug monitoring. In modeling the performance of new designs for photonic crystal spectrometers, one has to compute the solutions of Maxwell equations, which govern the light propagation in the spectrometer, with spatially incoherent sources f (x). Similarly, current models in fluorescent tomography are based on solving a diffusion approximation of the well-known radiative transport equation, and due to the random phase value it is again natural to model the incoherent fluorescent light source by point sources. Therefore, engineers routinely model incoherence by solving very large systems of equations, each of them excited by a point mass function fi (x) = fi δxi (x), i = 1, . . . , N . On one hand, the incoherence property is well modelled by point sources in (5.2). On the other hand, the sheer number of required point sources sets a computational roadblock. 95 To mitigate the aforementioned numerical complications, it was proposed in [4, 5] to cir- cumvent the local scale problem by replacing the localized forcing terms with a new global scale forcing that efficiently consolidates most of the energy into just a few terms. This was implemented by replacing multiple Maxwell equations with point sources by a single ˙ N (x) = PN ξi ni (x), where {ξi , i = 1, 2, . . . , N } Maxwell equation driven by white noise W i=1 were independent standard Gaussian random variables and {ni , i = 1, 2, . . . , N } was a sub- set of a trigonometric basis. The numerical simulations presented in [4, 5] demonstrate a dramatic reduction in computational complexity in evaluating the energy kv(·, t)k2L2 while maintaining a similar level of accuracy of energy approximation. However, papers [4, 5] were not concerned with rigorous theoretical explanations of the validity of the proposed algorithm and the potential scope of its applicability. Thus, we present here a rigorous approach to the problem of efficient approximation of the energy (5.3) for systems of fairly general evolution equations (5.2). In section 3, we compare the efficiency of the small scale (point forcing) basis and the “large scale” A-eigenfunction basis, and we deduce our main result—the approximation of the energy using the latter basis yields a 1st order improvement over the former (see Theorem 3.1). In fact, we will show that the number of expansion terms under the eigen- function basis is O (1) in N , whereas under the point forcing basis it is O (N ). In section 4, we show numerical results for the one-dimensional heat equation under the point forcing and cosine bases that corroborate the theoretical results, and we also show results for the convection-diffusion equations that suggest the applicability of this method to a broader class of parabolic equations. We remark that the change of basis method is not the only way to tackle the determin- istic system. The key point is the randomization of the system (5.2) to the SPDE (5.1), which can then be handled by various methods, such as WCE or Monte Carlo simulation. While there are numerous works in the literature studying such equations with additive noise, most of these works use a single choice of basis, which is usually a generic basis in the case of white noise, or the basis derived from the Karhunen–Lo`eve expansion (e.g., [13,20]). We point out that, at least in the case of a self-adjoint operator A, the choice of new basis should be related to the eigenfunctions of A (see section 3), rather than to the basis arising from the Karhunen–Lo`eve expansion of the noise. Interestingly, [11, 22] specifically chose 96 to use a basis similar to the point forcing basis, but this was only to expedite the use of the finite element method. In [11], the stochastic term was handled by Monte Carlo simulation, and L2 -convergence properties of the solutions were studied. To the best of our knowledge, direct comparison of two bases has not received as much attention. 2. Change of Wiener chaos basis We introduce the framework that will lead up to the proposed change of Wiener chaos basis idea. Let U ⊂ Rd be an open bounded domain. Let −A be a positive definite self-adjoint elliptic operator of order 2m, equipped with either periodic or zero Dirichlet boundary conditions. (In the case of periodic boundary conditions, the domain U will be a torus Td .) We assume the dimensionality condition (5.4) 2m/d > 1/2. It is well known that −A has eigenfunctions {mi } that form an orthonormal basis in L2 (U ), and the corresponding eigenvalues {λi } behave asymptotically as [58] (5.5) λi ∼ i2m/d . We will refer to {mi } as the A-eigenfunction basis in L2 (U ). As an unbounded positive definite self-adjoint operator on L2 (U ), −A has a well-defined √ m or H m . Then Λ induces a Hilbert square root Λ = −A, which has domain D(Λ) = Hper 0 γ scale which we denote by HA , γ ∈ R, with norms ∞  X  1/2m 2γ (5.6) kφk2H γ = λj φ2j A j=1 PJ γ for φ of the form φ = j=1 φj mj , for some J ∈ N [41]. HA is the closure of the set of such γ φ in the norm k · kH γ . It can be shown that HA is equivalent to the usual Sobolev scale. A In particular, the norm k · kH −2m is equivalent to the Sobolev norm A |hφ, ψiH −2m ,H 2m | kφkH −2m := sup , ψ∈H·2m kψkH 2m where we denoted H·2m = Hper 2m or H 2m in the case of periodic or zero Dirichlet boundary 0 conditions, respectively. 97 In order to define a localized basis, we consider a partition of the domain U . Let N < ∞ be arbitrary. Let I = {Ii , i = 1, . . . , N } be a partition of U into (small) nonoverlapping subsets with Lebesgue measure |Ii | ∼ 1/N . We assume the family I is quasi-uniform in N . That is, there exist constants ρ1 , ρ2 such that max ri ≤ (ρ1 |U |)1/d N −1/d , i min εi ≥ (ρ2 |U |)1/d N −1/d , i where ri = diam(Ii ) and εi is the radius of the largest sphere Bi contained in Ii . The quasi- (N ) uniform assumption implies nondegeneracy, i.e., that there exists ρ3 such that 2εi ≥ ρ3 ri for all i, N . It then follows that ρ− |U |N −1 ≤ min |Ii | ≤ max |Ii | ≤ ρ+ |U |N −1 i i and ρ˜− (ri )d ≤ |Ii | ≤ ρ˜+ (ri )d , and hence εi ∼ ri ∼ N −1/d . We are now ready to introduce the two bases {ni } and {mi } that will be the focus of our comparative analysis. (1) Point forcing basis: 1 (5.7) ni (x) = p 1Ii (x) for i = 1, . . . , N, |Ii | and {ni }∞ ⊥ i=N +1 is any basis in SN , where SN = span{ni , i = 1, . . . , N }. (2) (Discrete) eigenfunction basis in SN : m1 = m1 , i−1 ! 1 X (5.8) mi = PN mi − (PN mi , mj )mj , i = 2, . . . , N, Zi j=1 where PN is the L2 projection onto SN and Zi is the normalization constant. In other words, {mi } is the Gram–Schmidt orthonormalization of the L2 projections of the first N eigenfunction basis elements onto SN . 98 Many quantities considered in this chapter, such as the definitions of the two bases, depend on the parameter N . The limit as N → ∞ is an object of study. However, in the rest of chapter, we will suppress explicitly writing this dependence on N if no ambiguity arises. ˙ Q (x) on L2 (U ) by the Wiener chaos expansion Define the Gaussian noise W X (5.9) ˙ Q (x) := W σi ni (x)ηi i≥1 where ηi ∼ i.i.d. N (0, 1), σi ≥ 0, and the covariance operator Q2 is defined by Qni = σi ni for i = 1, 2, . . . (see (2.3)). We do not restrict Q2 to being a nuclear operator, but in the case ˙ Q to be a finite N -dimensional noise, we will assume that Range Q ⊂ SN . where we desire W (N ) In this case, σi = σi are nonzero only for i = 1, . . . , N . We consider the equation ∂v ˙ Q (x)G(t) (5.10) = Av + W ∂t with zero initial conditions and either periodic or zero Dirichlet boundary conditions.1 Here, G(t) is a bounded function on [0, T ] satisfying Z t CG1 CG2 (5.11) ≤ e−λj (t−s) G(s)ds ≤ , for t ∈ (0, T ], λj 0 λj for j = 1, 2, . . . , where the constants CG1 , CG2 are independent of j and N , and CG2 is independent of t. At this point, we introduce the related equation driven by an infinite dimensional Gauss- ian noise, which will be used for the error analysis in section 3.1. We assume for the sequence (N ) (N ) {σi , i = 1, . . . , N } that supN supi≤N σi < ∞, and (N ) N →∞ σi −→ σi∗ uniformly for i ≤ N. (N ) That is, ∀ǫ > 0, ∃N0 such that if N > N0 , then |σi − σi∗ | < ǫ, ∀i ≤ N . For simplicity, we assume σi∗ = 1 for all i = 1, 2, . . . , but the results can be extended to any bounded sequence σi∗ . The uniform convergence of {σi } makes it possible to study the asymptotic behaviour of (5.10) through the related limiting SPDE driven by an infinite dimensional 1For simplicity, we will always assume zero initial conditions and periodic or zero Dirichlet boundary condi- tions, even when not explicitly stated. We also always take x ∈ U and t ∈ (0, T ] for arbitrary T < ∞. 99 noise. To this end, we consider the SPDE with Gaussian white noise on L2 (U ), with Q = I and W ˙ (x) = P ξi mi (x), i≥1 ∂u∗ ˙ (x)G(t). (5.12) = Au∗ + W ∂t Equation (5.12) is solved in the triple H −m ֒→ L2 ֒→ H m and should be understood in the weak sense. Its propagator system is ∂uˆ∗i (5.13) u∗i + mi (x)G(t). = Aˆ ∂t The equivalence of the propagator system to the weak solution can be shown. Moreover, there exists a solution u∗ such that u∗ (t) ∈ L2 (Ω; L2 (U )) for each t ∈ (0, T ], and the energy E[u∗ ] := ku∗ (t)k2L2 (Ω;L2 (U )) at any fixed t ∈ (0, T ] is finite. (See section 3.1.) The framework to allow us to change the basis of the Wiener Chaos expansion is ele- mentary. Direct computation gives that N X N X Qmj = Σjk mk , where Σjk = σi (ni , mj )(ni , mk ). k=1 i=1 The uniform convergence of {σi } implies that Σjk −→ δjk as N → ∞. We will write Σj in ˙ Q, place of Σjj . Then there are two equivalent WCEs for W N X N X ˙ Q (x) = W ni (x)σi ηi = mi (x)Σi ξi , i=1 i=1 where N X ξi = Σ−1 i σk (nk , mi )ηk . k=1 The ξi s are identically distributed standard Gaussian random variables, but in general they are not independent. The covariance matrix ρ = (ρij )N i,j=1 is PN k=1 Σik Σkj ρij := E[ξi ξj ] = . Σi Σj Clearly, ρ is symmetric positive definite for each N , and we have ρij −→ δij as N → ∞. In the case that σi ≡ 1 for all i = 1, . . . , N , the ξi s are i.i.d. standard Gaussian random variables, and the relationship between {ξi } and {ηi } reduces to the usual change of basis 100 formula, N X (5.14) ξi = (mi , ni )ηi . i=1 We remark that the second expansion in (2) is, strictly speaking, not a Wiener chaos expansion, because the ξi s are not orthogonal in L2 (Ω). It is a standard exercise to transform the expansion into an orthogonal expansion by a linear transformation of the ξi s. However, we will not do that here, but instead just work directly with the linearly independent expansion (2). By the change of basis formula, we write the solution of (5.10) in two expansions N X N X (5.15) u(x, t) = vˆi (x, t)ηi = u ˆi (x, t)ξi . i=1 i=1 Multiplying both sides of (5.10) by ηi or ξi , and taking expectation yields two equivalent propagator systems ∂ (5.16a) vˆi = Aˆ vi + σi ni (x)G(t) ∂t N X X N ∂ (5.16b) ρji u ˆj = ρji (Aˆ uj + Σj mj (x)G(t)) ∂t j=1 j=1 for i = 1, . . . , N . Since ρ is invertible, (5.16b) reduces to a simpler system ∂ (5.16b′) ˆi = Aˆ u ui + Σi mi (x)G(t) ∂t Then the energy of (5.10), E[u] := Ekuk2L2 , can be computed from the solutions of either system (5.16a) or (5.16b′) by a simple algebraic formula N X N X (5.17) E[u] = vi k2L2 kˆ = (ˆ ui , u ˆj )ρij i=1 i=1 In order to reduce the computational cost of computing the solutions of all N equations in the system (5.16a) or (5.16b′), we approximate the energy of the N -system by the energy of a truncated system. Truncating the systems (5.16) to n < N equations means to consider 101 the systems ∂ (5.18a) vˆi = Aˆ vi + σi ni (x)G(t), for i = 1, . . . , n ∂t ∂ (5.18b) ˆi = Aˆ u ui + Σi mi (x)G(t), for i = 1, . . . , n. ∂t System (5.18a) is the propagator system of ∂ (n) ˙ Pn Q (x)G(t) (5.19) v = Av (n) + W ∂t where Pn the projection into span{ni , i = 1, . . . , n}. System (5.18b) is the related system to ∂ (n) (5.20) u = Au(n) + Z˙ n (x)G(t) ∂t Pn where Z˙ n (x) = i=1 mi (x)Σi ξi . Obviously, (5.19) and (5.20) are different SPDEs with different energies, n X n X (n) E[v ]= vi k2L2 kˆ 6= E[u (n) ]= (ˆ ui , u ˆj )ρij . i=1 i,j=1 The energies E[v (n) ] and E[u(n) ] will be taken as an approximation to the true energy E[u]. The absolute and relative errors of the approximations will be denoted as N X (n) R[v (n) ] := E[u] − E[v (n) ] = vi k2L2 , kˆ and ¯ (n) ] = R[v ] R[v E[u] i=n+1 N X (n) R[u (n) ] := E[u] − E[u (n) ]= ui k2L2 , kˆ and ¯ (n) ] = R[u ] R[u E[u] i=n+1 for n ≤ N . We will compare the performance of the two bases using the relative error of the approximate energy. Given an allowable relative error r, let (5.21) ¯ (n) ] < r} and nP := inf{n : R[v ¯ (n) ] < r}. nE := inf{n : R[u be the minimum truncation sizes that achieves the relative error r. Define the improvement of the eigenfunction basis over the point forcing basis as nP /nC . 102 The improvement is an indication of the computational savings of using the eigenfunction basis for the relative error r. 3. Comparative error analysis and 1st order improvement In the foregoing section, all the quantities depend on the number N of subdivisions of U . In this section, we study the asymptotic behavior as N → ∞. We will formulate precise bounds on the relative error and compare the asymptotic behavior of the two bases as N → ∞. The main goal of this section is to show the 1st order improvement of the change of basis method, in the sense of the following theorem. Theorem 3.1. Given a relative error r ∈ (0, 1), we have, at worst, 1st order improve- ment as N → ∞. More precisely, there exist constants 0 < C0,min < C0,max ≤ 1, depending on r but independent of N , such that for every C0 ∈ [C0,min , C0,max ) there exists N0 = N0 (C0 ) > 0 such that nP ≥ C0 N nE whenever N > N0 . Moreover, N0 → ∞ as C0 ↑ C0,max . Obviously, 1st order improvement is the best one can hope for, simply because nP ≤ N and nE ≥ 1, so that nP /nE ≤ N . The result of Theorem 3.1 states that the constant in front of the 1st order improvement can vary in an interval, with a larger constant holding for larger N . A big part of the proof of Theorem 3.1 involves studying the decay in n of the relative ¯ (n) ] and R[u errors R[v ¯ (n) ]. Theorem 3.1 follows easily from two key facts: first, that the ¯ (n) ] for the point forcing basis decays no faster than linearly; second, that relative error R[v ¯ (n) ] for the eigenfunction basis decays no slower than superlinearly, on the relative error R[u the order of n−α with α > 0. (See Figures 1 and 2 for illustration and motivation.) We make these two statements more precise in the following two propositions. Proposition 3.2. For the solution of (5.12), define the relative error by P∞ i=n+1 kˆu∗i k2L2 ¯ ∗,(n) ] := R[u E[u∗ ] 103 1 0.9 Eigenfunction basis Point forcing basis 0.8 0.7 Relative Error 0.6 0.5 0.4 0.3 0.2 0.1 0 0 20 40 60 80 100 120 Number of coefficients, n Figure 1. Relative errors incurred R[u¯ (n) ] when the system is truncated to n coefficients, under the point forcing basis (dotted line) and the eigenfunc- tion (cosine) basis (solid line). The convection-diffusion equation was used to produce this data. (a) 0 (b) 0 10 10 ε = 0.01 −1 −1 ε = 0.1 10 10 −2 −2 10 10 −3 −3 10 10 N = 30 Relative error (log scale) Relative error (log scale) −4 −4 10 10 −5 −5 10 N = 60 10 −6 −6 10 10 N = 120 −7 −7 10 10 N = 240 −8 −8 10 10 N = 480 −9 −9 10 10 N = 960 −10 −10 10 10 0 1 2 3 0 1 2 3 10 10 10 10 10 10 10 10 Number of coefficients, n (log scale) Number of coefficients, n (log scale) Figure 2. (a) Relative errors on log-log axes for increasing values of N , under the cosine basis for the heat equation. (b) Relative errors for two values of diffusion coefficients ǫ = 0.1, 0.01. The graph for ǫ = 0.01 lies above the graph for ǫ = 0.1. for n = 1, 2, . . . . Then (5.22) ¯ ∗,(n) ] ∼ n−4m/d +1 . R[u Given a relative error r, n o d (5.23) ¯ ∗,(n) ] < r ∼ r d−4m n0 := inf n : R[u 104 as r ↓ 0. Proposition 3.3. There exists a constant C independent of n and N such that (5.24) ¯ (n) ] ≥ L(n) = 1 − nCN −1 , R[v where L(n) is a straight line passing through the point (0, 1) and with slope −CN −1 , which tends to 0 as N → ∞. To show the decay behavior of the relative errors, we will focus on finding bounds on the L2 norms of the solution modes u ˆi and vˆi . In order to be useful for explaining this contrasting behavior of the two bases, the bounds need to be sensitive to the localness or globalness of the basis and should provide accurate bounds on the solution modes. Error bounds involving kni k2L2 and kmi k2L2 are clearly insensitive to the choice of basis, since both norms equal 1. Standard methods for estimating the time evolution of ku(t)k2L2 , such as those involving Gronwall’s inequality, may also be inadequate. A case in point is the following. Suppose u(t) solves the heat equation ∂u = ∆u + f (x), x ∈ [0, X], ∂t with zero initial conditions and periodic boundary conditions (cf. section 4.1). Also assume that u(t) has periodic derivative. Then Z 1d ku(t)k2L2 = u(x, t)ut (x, t) dx 2 dt U ≤ −kux (t)k2L2 + ku(t)kH 1 kf kH −1  1 ≤ −kux (t)k2L2 + ku(t)k2L2 + kux (t)k2L2 + kf k2H −1 4 1 = ku(t)k2L2 + kf k2H −1 . 4 Thus, by Gronwall’s inequality, te2t ku(t)k2L2 ≤ kf k2H −1 = C(t)kf k2H −1 . 4 105 Using ni or the cosine basis mi in place of f , the energy of each mode is bounded by  2 X ui (t)k2L2 kˆ ≤ C(t)kmi k2H −1 ≤ C(t) , (i − 1)π vi (t)k2L2 ≤ C(t)kni k2H −1 ≤ C(t). kˆ Then the absolute error of the energy estimate of an n-equation truncated system (for fixed N ) decays on the order of     (n) 1 1 (n) N −n R[u ]∼O − , R[v ]∼O . n N N The estimate for R[v (n) ] is consistent with numerical results, but we will establish this result in more generality. But the estimate for R[u(n) ] is merely an upper bound and does not  predict the actual O n−3 decay (see Figure 2). 3.1. Fourier techniques for the eigenfunction basis. The Fourier expansion is an effective technique for obtaining exact error estimates. We begin by considering the limiting infinite dimensional equation (5.12) with covariance operator Q = I, and its propagator system (5.13). The solution of (5.12) has the Wiener chaos expansion: ∞ X ∞ X ∗ u (t) = ˆ∗i (t)ξi u = ˆ ˆ∗ij (t)mj ξi , u i=1 i,j=1 ˆ∗ = (ˆ where u ˆ u∗i , mj ) are the Fourier coefficients of u ˆ∗i with respect to the A-eigenfunction ij basis in L2 (U ). Note that only the low order modes {ξi mj }i,j≥1 are nonzero because the noise appears additively. From the propagator system (5.13), the Fourier coefficients solve the decoupled system of ODEs   du ˆ∗ ˆ ˆ∗ij + δij G(t), dt ˆij = −λj u (5.25)  uˆ ˆ∗ (0) = 0 ij for i, j = 1, 2, . . . . The solution satisfies Z t CG1 ˆ CG2 δij ˆ∗ij (t) = δij ≤u e−λj (t−s) G(s) ds ≤ δij . λj 0 λj 106 In other words, all the energy is concentrated on the modes {ξi mi }, and X∞ X∞ 2 1 ∗ 2 2 1 CG1 2 ≤ Eku (t)kL2 (U ) ≤ C G2 2, i=1 λ i λ i=1 i where the summations converge because of the asymptotic behavior of the eigenvalues (5.5) and the dimensionality condition (5.4). It follows that the dimensionality condition (5.4) is necessary and sufficient for u∗ (t) ∈ L2 (Ω; L2 (U )) to be square integrable. The consequence of this computation is that we have precise asymptotic estimates for the truncation error: ∞ X ∞ X (5.26) R[u∗,(n) ] := u∗i k2L2 = kˆ ˆ ˆ∗ii )2 ∼ n−4m/d +1 . (u i=n+1 i=n+1 ¯ ∗,(n) ] = R[u∗,(n) ]/E[u∗ ] as well, and we The asymptotics in (5.26) obviously hold for R[u obtain (5.22) in Proposition 3.2. Equation (5.23) then follows by taking the function inverse of the asymptotic bounds in (5.22). Thus, d d (5.27) Cr d−4m ≤ n0 ≤ C ′ r d−4m for r sufficiently small. ˙ Q defined We now look at the finite dimensional model, N < ∞, with nuclear Q and W ˙ Q can be viewed as a finite dimensional approximation of W in (5.9). W ˙ , and we have that E[u] → E[u∗ ]. Instead of (5.25), the relevant system of ODEs for (5.10) corresponding to the discrete eigenfunction basis in SN comes from (5.16b′):,   du ˆ ˆ dt ˆij = −λj u ˆij + Σi (mi , mj )G(t), (5.28)  uˆ ˆij (0) = 0 for i = 1, . . . , N , j = 1, 2, . . . . The solution is Z t ˆ u ˆij (t) = Σi (mi , mj ) e−λj (t−s) G(s) ds. 0 The behavior of the truncation error of system (5.28) can be expected to be close to that of (5.25). Indeed, since the Gram–Schmidt orthonormalized elements mi in (5.8) are finite sums of projections PN mk , and since PN mk −→ mk in L2 as N → ∞, it follows that 107 mi −→ mi in L2 as N → ∞ for each i = 1, 2, . . . . For any fixed j, (mi , mj ) −→ δij also. Thus, by (5.11) and the dominated convergence theorem, X Z t 2 −λj (t−s) (ˆ ui , u ˆi′ )ρii′ = ρii′ Σi Σi′ (mi , mj )(mi′ , mj ) e G(s)ds j≥1 0 X Z t 2 −λj (t−s) −→ δii′ δij δi′ j e G(s)ds u∗i k2L2 (U ) , = δii′ kˆ as N → ∞ j≥1 0 for all i, i′ = 1, 2, . . . . Since E[u] −→ E[u∗ ], it follows that N →∞ (5.29) R[u(n) ] = E[u] − E[u(n) ] −→ R[u∗,(n) ] = E[u∗ ] − E[u∗,(n) ] for every n. In particular, we deduce that for fixed truncation size n, the relative error incurred by truncating the large size N system must tend to the relative error incurred by truncating the infinite system.  ¯ (n) ] = O n− 4m We do not assert that R[u d +1 . Nonetheless, similar to n0 in (5.23), an asymptotic result for the minimum truncation size to achieve relative error r for the discrete eigenfunction basis, easily follows. Proposition 3.4. Let nE = nE (N, r) be defined as in (5.21). Then there is some r∗ ∈ (0, 1) such that for any relative error r < r∗ , there exists N (r) such that d d Cr d−4m ≤ nE ≤ C ′ r d−4m whenever N ≥ N (r). The constants C, C ′ are independent of r, N . Proof. From (5.27), there exists r∗ such that d d Cr d−4m ≤ n0 (r) ≤ C ′ r d−4m whenever r < r∗ . Fix δ ∈ (0, 1). For any r < r∗ /(1 + δ), choose N (r) such that ¯ (n) ¯ ∗,(n) ] < δr R[u ] − R[u holds for both n = n0 ((1 + δ)r) and n = n0 ((1 − δ)r). Then ¯ (n) ]|n=n ((1+δ)r) > r R[u and ¯ (n) ]|n=n ((1−δ)r) < r R[u 0 0 108 and n0 ((1 + δ)r) < nE ≤ n0 ((1 − δ)r). Hence d d d d C(1 + δ) d−4m r d−4m ≤ nE ≤ C ′ (1 − δ) d−4m r d−4m , and the result follows with r∗ /(1 + δ) in place of r∗ .  3.2. H −2m norm estimates for the point forcing basis. For the error analysis for the point forcing basis, similar computations for the system of ODEs for the point forcing basis coming from (5.16b′) show that vˆ ˆij := (ˆ vi , mj ) satisfies Z t vˆ ˆij (t) = σi (ni , mj ) e−λj (t−s) G(s) ds. 0 Clearly, the energy of the system is not concentrated on {vˆ ˆii , i = 1, . . . , N }, and ∞ X Z t 2 (5.30) vi k2L2 kˆ = σi2 (ni , mj )2 e −λj (t−s) G(s) ds . j=1 0 We have the following lemma. Lemma 3.5. Let ni be defined in (5.7) for i = 1, . . . , N . Then we have the bounds C1 N −2m/d ≤ kni kH −2m ≤ C2 N −1/2 , where C1 , C2 are independent of i and N . Proof. For the lower bound, consider the mollifier ζε with support in B(0, ε), and let αi be the center of the largest sphere Bi contained in Ii with radius εi . Then, denoting H·2m = Hper 2m or H 2m , 0 |hni , ψi| kni kH −2m = sup ψ∈H·2m (U ) kψkH 2m (U ) Z −1 ≥ kζεi (· − αi )kH 2m (U ) ni (x)ζεi (x − αi ) dx U = kζεi k−1 (n ∗ ζεi )(αi ) ≥ CN −2m/d . H 2m (Rd ) i −(2m+d/2) The last inequality holds because it can be computed that kζεi kH 2m ∼ εi and −d/2 (ni ∗ ζεi )(αi ) = ni (αi ) ∼ N 1/2 ∼ εi . 109 For the upper bound, 1 R |hni , ψi| |Ii | Ii |ψ| dx kni kH −2m = sup ≤ |Ii |1/2 sup ψ∈H·2m kψk H 2m ψ∈H·2m kψkH 2m ≤ CN −1/2 , where the C is independent of i and N , since by the Sobolev embedding every ψ belonging to H02m (U ) or Hper ¯ ). 2m (U ) also belongs to C 0,1/2 (U  Corollary 3.6. For each i = 1, . . . , N , we have the bounds C3 N −4m/d ≤ σi−2 kˆ vi k2L2 ≤ C4 N −1 , where C3 , C4 are independent of i and N . γ Proof. From the definition of the HA norm, (5.6), ∞   X CG1 2 vi k2L2 kˆ ≥ σi (ni , mj ) = σi2 CG1 2 kni k2H −2m λj A j=1 and similarly ∞  X 2 CG2 vi k2L2 kˆ ≤ σi (ni , mj ) = σi2 CG2 2 kni k2H −2m . λj A j=1 γ The result follows by the equivalence of the HA norms and the Sobolev norms, and from Lemma 3.5.  The lower bound in Corollary 3.6 gives another way to see that the solution of the finite system will not converge to a square integrable of the infinite system if the dimensionality condition (5.4) is not met. This lower bound also gives a lower bound for the relative error of the point forcing basis. Interestingly, a more informative lower bound on the relative error can be derived from the upper bound in Corollary 3.6: P E[u] − ni=1 kˆ vi k2L2 ¯ (n) R[v ] = E[u]  E[u] − C4 N −1 n supN supi≤N σi2 (5.31) ≥ E[u] C5 =1−n =: L(n). N E[u] 110 Since the constant C4 is independent of n and N , the relative error is bounded from below by a straight line L(n) passing through the point (0, 1), and with slope −C5 /(N E[u]) which tends to 0 as N → ∞. We have just shown Proposition 3.3. We now prove Theorem 3.1. Proof of Theorem 3.1. From (5.31), the linear lower bound L(n) attains relative error r for n ≥ nL , where (1 − r)E[u] ˜ )N nL = N = C(N C5 ¯ (n) ] ≥ L(n), is the value such that L(nL ) = r. Then, since R[u ˜ )N. nP ≥ nL = C(N ˜ ) depends on N because E[u] depends on N . We next show a series of estimates for C(N ˜ ) to remove the dependence on N . First, E[u{N } ] ≥ E[u{N =1} ] for all N , so C(N ˜ ) ≥ C(1) C(N ˜ for all N . Now, since E[u] → E[u∗ ], for any ǫ ∈ (0, E[u∗ ] − E[u{N =1} ]), there exists N (ǫ) such that E[u] ≥ E[u∗ ] − ǫ. So ∗ ˜ ) ≥ (1 − r)(E[u ] − ǫ) C(N C4 ˜ (1−r)E[u∗ ] whenever N > N (ǫ). Denote C(∞) = C5 . As ǫ ranges from 0 to E[u∗ ] − E[u{N =1} ], ˜ the right-hand side of the last inequality ranges from C(∞) ˜ to C(1). Clearly, N (ǫ) increases ˜ to ∞ as ǫ ↓ 0. In other words, for any C ∈ [C(1), ˜ C(∞)), there exists N (C) such that ˜ )N ≥ CN nP ≥ C(N ˜ whenever N > N (C). Moreover, N (C) increases to ∞ as C ↑ C(∞). ¯ ∗,(n) ] < r}, (5.23). Let ǫ0 = r − R[u For nE , recall n0 = inf{n : R[u ¯ ∗,(n0 ) ] > 0. From ¯ (n0 ) ] → R[u (5.29), R[u ¯ ∗,(n0 ) ] as N → ∞. So there exists N (n0 ) > 0 such that ¯ (n0 ) ] < R[u R[u ¯ ∗,(n0 ) ] + ǫ0 = r whenever N > N (n0 ). Hence, nE ≤ n0 if N > N (n0 ). 111 Combining the two inequalities for nP , nE , nP CN ≥ = C0 N nE n0 whenever N > N0 (C0 ) := max{N (C), N (n0 )}. Hence, C0 ∈ [C(1)/n0 , C(∞)/n0 ) and N0 → ∞ as C0 ↑ C(∞)/n0 .  The next result gives upper and lower bounds on the improvement in terms of the relative error r. Corollary 3.7. There exist r∗ ∈ (0, 1) and constants 0 < C∗,min < C∗,max ≤ C∗ ≤ 1 such that, for every r < r∗ and every C0 ∈ [C∗,min , C∗,max ), there exists N0 = N0 (r, C0 ) > 0 such that d nP d C0 r− d−4m N ≤ ≤ C∗ r− d−4m N nE whenever N > N0 . Moreover, N0 → ∞ as C0 ↑ C∗,max or as r ↓ 0. ˜ ) with Proof. In the proof of Theorem 3.1, the inequalities hold if we replace C(N (1−r ∗ )E[u] C˜∗ (N ) := C5 so that, for any C ∈ [C˜∗ (1), C˜∗ (∞)), there exists N (C) such that nP ≥ C˜∗ (N )N ≥ CN whenever N > N (C). Also N (C) increases to ∞ as C ↑ C˜∗ (∞). From Proposition 3.4, nP CN d ≥ d = C0 N r− d−4m nE C ′ r d−4m whenever N > N0 (r, C0 ). Also from Proposition 3.4, and since nP ≤ N , nP N d ≤ d = C∗ N r− d−4m . nE Cr d−4m  d d − d−4m − d−4m If r1 < r2 , then r1 < r2 , so Corollary 3.7 indicates that one would expect a slower convergence to 1st order improvement for a smaller relative error. This observation is in accordance with Trend (T3) in the numerical simulations. We also note that the interval endpoints in Theorem 3.1 and Corollary 3.7 are inversely proportional to C4 from Corollary 3.6, which is in turn inversely proportional to the norm of A. This point is corroborated 112 by the numerical result that showed that the improvement is better for a larger diffusivity constant (cf. Trend (T1)). 3.3. The non–self-adjoint case. If A is not self-adjoint or not positive definite, pre- cise bounds on the error decay such as those obtained in the previous section may not be readily available. But, under additional assumptions, we can still deduce certain as- ymptotic results similar to the positive definite self-adjoint case, including the result of 1st order improvement. To see this, let us consider again the SPDE (5.12) in the triple H m ֒→ L2 ֒→ H −m , where we assume A is a 2mth order non–self-adjoint elliptic operator. Also assume a more stringent dimensionality condition: (5.32) m/d > 1/2. 1 We decompose A = A0 + A1 into the symmetric part A0 = 2 (A + A∗ ) and the skew- symmetric part A1 = 12 (A − A∗ ), and we assume that −A0 is positive definite. Then −A0 generates an eigenfunction basis {mi } with eigenfunctions {λi } satisfying (5.5). Similarly γ 2 P∞ 2 γ/m , to (5.6), A0 defines a scale of Hilbert spaces HA 0 , with norm kφkH γ = j=1 (φ, mj ) λj A0 that is equivalent to the Sobolev scale Hγ. In the infinite dimensional case with white noise (5.12), the existence and uniqueness of the solution u∗ is shown in [57] because the asymptotics of the eigenvalues (5.5) and the ˙ ∈ L2 (Ω; H −m (U )). Applying the usual new dimensionality condition (5.32) imply that W ˆ∗i (t) is deterministic parabolic estimates to the propagator system (5.13), we have that u continuous in t, and u∗i (t)k2L2 (U ) ≤ CkGk2L2 (0,T ) kmi k2H −m ≤ C ′ kmi k2H −m ≤ C ′ λ−1 kˆ i . A0 for all t ∈ (0, T ]. Then we have a result analogous to (but weaker than) Proposition 3.2. For the error ∞ X ∗,(n) (5.33) R[u ] := u∗i k2L2 ≤ Cn−2m/d +1 , kˆ i=n+1 ¯ ∗,(n) ] < r}, and for n0 := min{n : R[u d n0 ≤ Cr d−2m . 113 In the finite dimensional case (5.10), we again have E[u] → E[u∗ ]. For the discrete eigenfunction basis, kˆ u∗i k2L2 = |kˆ ui k2L2 − kˆ u∗i kL2 | (kˆ ui kL2 − kˆ u∗i kL2 ) ≤ Ckˆ ui kL2 + kˆ ˆ∗i kL2 ui − u N →∞ ≤ C ′ kΣi mi − mi kH −m −→ 0 for each i ≤ N , and so n X (n) E[u ] − E[u∗,(n) ] = ui k2L2 − kˆ kˆ u∗i k2L2 −→ 0. i=1 ¯ (n) ] −→ R[u Hence, R[u ¯ ∗,(n) ] as N → ∞ for each n. For the point forcing basis, vi (t)k2L2 (U ) ≤ Cσi2 kni k2H −m ≤ C ′ N −1 , kˆ where the last inequality follows by an argument similar to the upper bound in Lemma 3.5. The proof of Theorem 3.1 follows through identically, so the statement of 1st order improvement applies to the non–self-adjoint case as well, provided (5.32) holds. However, this argument by parabolic estimates works only when (5.32) holds; the be- havior when 1/4 < m/d ≤ 1/2, which was covered in the self-adjoint case, is not addressed here. This should not be a surprise because the parabolic estimates are essentially Gronwall- type estimates, which we have noted in the beginning of section 3 to give suboptimal error bounds. The main difference between the two analyses is the estimation of the forcing terms in the H −m norm in the parabolic estimate case, rather than the H −2m norm in the self-adjoint case. Hence, the parabolic estimates provide only upper bounds on R[u∗,(n) ]  that are O n−2m/d +1 , which is less favorable and less precise than the o(n−4m/d +1 ) decay found in the self-adjoint case. Nonetheless, we conjecture that the asymptotic behavior of R[u∗,(n) ] should in principle be dominated by the self-adjoint part A0 , even though this is not reflected with the parabolic estimates (see section 4.3). 4. Examples and simulations The change of basis strategy is applied to some simple equations to illustrate the effi- ciency of the point forcing and cosine bases, (5.34), (5.35), for approximating the energy of the systems (5.16a,b). One of the equations considered is the heat equation, for which 114 we will observe results that corroborate the analysis in section 3. Although the analysis is asymptotic in nature, the 1st order convergence is already clear even for not-too-large system sizes. We also present numerical results for convection-diffusion equations that share very similar comparative properties to the pure diffusion case and extend the discussion to the connection with the pure convection equation. For our numerical simulations, we take the interval U = [0, X], and we let IN = {Ii , i = 1, . . . , N } be a uniform partition of U into intervals of length X/N . We consider the operator with A = ǫ∆ with periodic boundary conditions, whose eigenfunctions are the usual cosine basis. ǫ is a small diffusivity coefficient. The two bases on SN := span{ni , i = 1, . . . , N } are the following: (1) Point forcing basis: r N (5.34) ni (x) = 1I (x) for i = 1, . . . , N. X i (2) Cosine basis in SN : The eigenfunction basis in L2 ([0, X]) is the usual cosine basis: r 1 m1 (x) = , X r   2 (i − 1)πx mi (x) = cos , i = 2, 3, . . . . X X Define the cosine basis in SN as the Gram–Schmidt orthonormalization of the L2 projections of the first N cosine basis elements onto SN : m1 = m1 , i−1 ! 1 X (5.35) mi = PN mi − (PN mi , mj )mj , Zi j=1 where PN is the L2 projection onto SN and Zi is the normalization constant. In this example, we take G(t) = 1, and take the covariance Q = PN , so that σi ≡ σi∗ ≡ 1 ˙ N are for all i = 1, 2, . . . . Then Σi = 1 for all i = 1, . . . , N and the two WCEs for W N X N X ˙ N (x) = W ni (x)ηi = mi (x)ξi i=1 i=1 where ηi and ξi are related by the usual change of basis formula (5.14). As a side note, ˙ N is a finite truncation of the white noise W since W ˙ , it is well-known that we can give ξi 115 and ηi precise expressions: Z Z ξi := mi (x) dW (x) and ηi := ni (x) dW (x), U U where W (x) is a Brownian motion on U and from which the change of basis formula can be checked by direct computation. We study the equation ∂u ˙ N (x) (5.36) = ǫ∆u + W ∂t with zero initial conditions and periodic boundary conditions. Equations (5.15), (5.16), and (5.17) hold. Note that, strictly speaking, the analysis of section 3 does not apply to (5.36) because −∆ with periodic boundary conditions has an eigenvalue λ1 = 0 and thus is not strictly positive definite. Nonetheless, we can still apply the ideas from section 3 to obtain analogous results for the error decay and 1st order improvement. Equation (5.11) holds for j = 2, 3, . . . , while for j = 1, λ1 = 0, Z t e−λ1 (t−s) ds = t, 0 so equations (5.26), (5.27) and (5.29) hold also, as do Propositions 3.2 and 3.4. For the point forcing basis, the analogous result to Corollary 3.6 is C3 N −1 ≤ kˆ vi k2L2 ≤ C4 N −1 . Indeed, the lower bound is vi (t)k2L2 ≥ |vˆ kˆ ˆi,1 (t)|2 = (ni , m1 )2 t2 = N −1 t2 . For the upper bound, we integrate by parts backwards twice to find (5.37) (ni , λ−1 j mj ) = (−1) j−1 −1 ′ λj fi (X) + (fi , mj ), j ≥ 2, 116 for some function fi (x) such that fi′′ = ni . (This step takes the place of invoking the H −2 norm of ni .) It can be directly computed that    (i−1)X   0, x≤ N ,    q    fi (x) = 1 N x − (i−1)X , (i−1)X iX N N N , p 17X 4 1 so fi′ (X) = X/N and kfi k2L2 = · · · ≤ 15 N . Squaring (5.37) and summing over j, X X kni k2H −2 = λ−2 (ni , mj )2 ≤ 2f ′ (X)2 λ−2 2 j + 2kfi kL2 ≤ CN −1 j≥1 j≥1 vi k2L2 ≤ Ckni k2H −2 ≤ CN −1 . Hence kˆ 4.1. Heat equation. For the heat equation (5.36), we show in Figure 2(a) the relative error of the truncated system under the cosine basis for different values of N . We observe that, for each n, the relative error increases pointwise to a limit as N → ∞. We assume that the N = 960 error plot is representative of the error in the limit as N → ∞, at least for n not near 960. When n > 10, the relative error decays linearly on the log-log axes, with  ¯ (n) ] ∼ O n−3 . This same order of decay is seen for ǫ = 0.01 a gradient of ≈ −3; i.e., R[u only when n > 40 (Figure 2(b)), and the actual relative error is larger than for ǫ = 0.1. Both these orders of decay are consistent with (5.26) when m = d = 1. In contrast, the relative error decays linearly on the linear axes for the point forcing basis (cf. Figure 1) and does not exhibit the same limiting behavior as the error plots for the cosine basis do. In fact, in this case of periodic boundary conditions, the relative error plot is simply a straight line of slope −N −1 joining the points (0, 1) and (N, 0) because the vi k2L2 is equal. For a given level of relative error and for large values of N , energy of each kˆ nP for the point basis scales on the order of O(N ), whereas nE for the cosine basis scales with O(1). As a result, this implies the 1st order convergence seen in Table 1. Table 1(b) shows the improvements of the cosine basis for 5% error. We highlight several trends. (T1) For fixed N , the improvement increases for larger ǫ. This increase is most signifi- cant for large N . 117 Table 1. Improvement nP /nE in the number of basis elements required to attain 5% error. (a) Convection-diffusion equation (b) Heat equation N ǫ = 0.1 ǫ = 0.01 ǫ=0 ǫ = 0.1 ǫ = 0.01 30 2.6364 2.4167 2.4167 1.8125 1.0741 60 4.2846 4.1429 3.8000 3.4118 1.3571 120 8.7692 7.6000 4.7917 6.3889 2.3 240 17.5385 12.6667 8.4815 12.7222 4.3019 480 35.0769 21.7619 15.7241 25.3889 8.4444 960 70.2308 43.4762 30.4333 50.6667 16.8889 (T2) We have 1st order improvement: doubling N increases the improvement by a factor that approaches double as N becomes large. (T3) 1st order improvement is seen for a smaller error of 1% (data not shown), but the convergence to 1st order improvement is slower. Numerical scheme. The discontinuous Galerkin dG(1) scheme with a 2nd order Runge– Kutta time stepping scheme [12] was used in this computation. For each number N of forcing terms, we took N spatial grid points and used X = 2π, T = 0.5. The simulations were also done using a fixed number of grid points (960 grid points) for all values of N , but little difference was found in the quantitative and qualitative behaviors of the estimates. 4.2. Convection-diffusion equations. We applied the same change of basis method for the stochastic convection-diffusion equation ∂ ˙N (5.38) u + bux = ǫuxx + W ∂t with zero initial conditions and periodic boundary conditions. We performed simulations with constant convection speed b0 = 1.47 and small diffusive coefficients ǫ = 0.01, 0.1. Figure 1 shows the behavior of the relative errors of the two bases on linear axes. Under the point forcing expansion, the relative error of the truncated system decays linearly in n, whereas the relative error under the cosine expansion decays superlinearly. The improve- ment is also found for varying sizes of the full system, N = 30, 60, . . . , 960 (Table 1(a)). 4.3. Further remarks. As noted in section 3.3, it is not straightforward to deduce precise error estimates for general equations where A does not provide an eigenfunction basis. If the equation is simple enough, the error decay rate can be found from the explicit 118 solution. In the case of (5.38), Z ∂ ui k2L2 = kˆ 2ˆ ui (−bˆ ui,x + ǫˆ ui,xx + mi ) dx ∂t U Z = ˆ2i + 2ǫˆ bx u ui u ˆi,xx + 2ˆ ui mi dx. U If ǫ = 0, Z tZ Z t ui (t)k2L2 = kˆ kˆ ui (0)k2L2 + 2 u ˆi mi dx dt = 2 (ˆ ui (·, τ ), mi ) dτ, 0 U 0 ˆ so the error of each mode depends only on the coefficients u ˆii (τ ) := (ˆ ui (τ ), mi ) up to time t. By explicitly solving the convection equation,   ˆ X (i − 1)πtb u ˆii (t) = sin (i − 1)πb X and hence N X N X Z t (n) ˆ (N ) R[u ]= ui k2L2 kˆ = 2 u ˆii dτ i=n+1 i=n+1 0 N  X 2   X (i − 1)πtc =2 1 − cos (i − 1)πb X i=n+1     1 1 N →∞ 1 ∼O − −→ O . n N n  An approximately O n−1 decay for the pure diffusion case is seen in Figure 3—this is the decay rate predicted by the parabolic estimate analysis in section 3.3. If ǫ > 0, the decay rate seems to be a hybrid between the convection and the diffusion parts—for small  n, the O n−1 decay from the convection part dominates, while for large n the decay shows  better agreement with the O n−3 decay from the diffusion part. Evidently, the analysis in section 3.3 is unable to capture the intermediate and asymptotic behaviors of the error decay for the convection-diffusion equation. 119 0 10 −1 10 −2 10 Relative Error (log scale) −3 10 −4 10 Pure diffusion −5 10 Convection−diffusion −6 10 Pure convection −7 10 −8 10 0 1 2 3 10 10 10 10 Number of coefficients, n (log scale) Figure 3. Log scale plots of the relative errors incurred by the truncation of the convection-diffusion system, as well as the pure diffusion and the pure convection systems. The convection and diffusion coefficients are b = 6b0 and ǫ = 0.1, respectively. 120 Bibliography [1] Ivo Babuˇska, Fabio Nobile, and Ra´ ul Tempone. A stochastic collocation method for elliptic partial differential equations with random input data. SIAM J. Numer. Anal., 45(3):1005–1034 (electronic), 2007. [2] Ivo Babuˇska, Ra´ ul Tempone, and Georgios E. Zouraris. Galerkin finite element approximations of stochastic elliptic partial differential equations. SIAM J. Numer. Anal., 42(2):800–825, 2004. [3] Ivo Babuˇska, Ra´ ul Tempone, and Georgios E. Zouraris. Solving elliptic boundary value problems with uncertain coefficients by the finite element method: the stochastic formulation. Comput. Methods Appl. Mech. Engrg., 194(12-16):1251–1294, 2005. [4] Majid Badieirostami, Ali Adibi, Hao-Min Zhou, and Shui-Nee Chow. Efficient modeling of spatially incoherence sources based on wiener chaos expansion method for the analysis of photonic crystal spec- trometers. Proc. SPIE Int. Soc. Op. Eng., pages 648018–1–8, 2007. [5] Majid Badieirostami, Ali Adibi, Hao-Min Zhou, and Shui-Nee Chow. Wiener chaos expansion and simulation of electromagnetic wave propagation excited by a spatially incoherent source. Multiscale Model. Simul., 8(2):591–604, 2009/10. ´ [6] A. Bensoussan and R. Temam. Equations stochastiques du type Navier-Stokes. J. Functional Analysis, 13:195–222, 1973. [7] Fred Espen Benth and Jon Gjerde. Convergence rates for finite element approximations of stochastic partial differential equations. Stochastics Stochastics Rep., 63(3-4):313–326, 1998. [8] Fred Espen Benth and Thomas Gorm Theting. Some regularity results for the stochastic pressure equation of Wick-type. Stochastic Anal. Appl., 20(6):1191–1223, 2002. [9] R. H. Cameron and W. T. Martin. The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals. Ann. of Math. (2), 48:385–392, 1947. [10] Yanzhao Cao. On convergence rate of Wiener-Ito expansion for generalized random variables. Stochas- tics, 78(3):179–187, 2006. [11] Yanzhao Cao, Hongtao Yang, and Li Yin. Finite element methods for semilinear elliptic stochastic partial differential equations. Numer. Math., 106(2):181–198, 2007. [12] Bernardo Cockburn, George E. Karniadakis, and Chi-Wang Shu. The development of discontinuous Galerkin methods. In Discontinuous Galerkin methods (Newport, RI, 1999), volume 11 of Lect. Notes Comput. Sci. Eng., pages 3–50. Springer, Berlin, 2000. 121 [13] Arnaud Debussche and Jacques Printems. Weak order for the discretization of the stochastic heat equation. Math. Comp., 78(266):845–863, 2009. [14] Lawrence C. Evans. Partial differential equations, volume 19 of Graduate Studies in Mathematics. Amer- ican Mathematical Society, Providence, RI, 1998. [15] Franco Flandoli. Dissipativity and invariant measures for stochastic Navier-Stokes equations. NoDEA Nonlinear Differential Equations Appl., 1(4):403–423, 1994. [16] Franco Flandoli and Dariusz Gatarek. Martingale and stationary solutions for stochastic Navier-Stokes equations. Probab. Theory Related Fields, 102(3):367–391, 1995. [17] C. Foia¸s. Statistical study of Navier-Stokes equations. I, II. Rend. Sem. Mat. Univ. Padova, 48:219–348 (1973); ibid. 49 (1973), 9–123, 1972. [18] C. Foia¸s and R. Temam. Homogeneous statistical solutions of Navier-Stokes equations. Indiana Univ. Math. J., 29(6):913–957, 1980. [19] Ciprian Foias, Ricardo M. S. Rosa, and Roger Temam. A note on statistical solutions of the three- dimensional Navier-Stokes equations: the stationary case. C. R. Math. Acad. Sci. Paris, 348(5-6):347– 353, 2010. [20] Philipp Frauenfelder, Christoph Schwab, and Radu Alexandru Todor. Finite elements for elliptic prob- lems with stochastic coefficients. Comput. Methods Appl. Mech. Engrg., 194(2-5):205–228, 2005. [21] J. Galvis and M. Sarkis. Approximating infinity-dimensional stochastic Darcy’s equations without uni- form ellipticity. SIAM J. Numer. Anal., 47(5):3624–3651, 2009. [22] R. G. Ghanem and P. D. Spanos. Polynomial chaos in stochastic finite elements. J. Appl. Mech., 57:197– 202, 1990. [23] Roger G. Ghanem and Pol D. Spanos. Stochastic finite elements: a spectral approach. Springer-Verlag, New York, 1991. [24] Colette Guillop´e. Comportement ` a l’infini des solutions des ´equations de Navier-Stokes et propri´et´e des ensembles fonctionnels invariants (ou attracteurs). Ann. Inst. Fourier (Grenoble), 32(3):ix, 1–37, 1982. [25] Istv´ an Gy¨ ongy and David Nualart. Implicit scheme for quasi-linear parabolic partial differential equa- tions perturbed by space-time white noise. Stochastic Process. Appl., 58(1):57–72, 1995. [26] Istv´ an Gy¨ ongy and David Nualart. Implicit scheme for stochastic parabolic partial differential equations driven by space-time white noise. Potential Anal., 7(4):725–757, 1997. [27] Takeyuki Hida, Hui-Hsiung Kuo, J¨ urgen Potthoff, and Ludwig Streit. White noise, volume 253 of Mathematics and its Applications. Kluwer Academic Publishers Group, Dordrecht, 1993. An infinite- dimensional calculus. [28] Takeyuki Hida, Hui-Hsiung Kuo, J¨ urgen Potthoff, and Ludwig Streit. White noise, volume 253 of Mathematics and its Applications. Kluwer Academic Publishers Group, Dordrecht, 1993. An infinite- dimensional calculus. 122 [29] H. Holden, T. Lindstrøm, B. Øksendal, J. Ubøe, and T.-S. Zhang. The pressure equation for fluid flow in a stochastic medium. Potential Anal., 4(6):655–674, 1995. [30] Helge Holden, Tom Lindstrøm, Bernt Øksendal, Jan Ubøe, and Tu Sheng Zhang. Stochastic boundary value problems: a white noise functional approach. Probab. Theory Related Fields, 95(3):391–419, 1993. [31] Helge Holden, Bernt Øksendal, Jan Ubøe, and Tusheng Zhang. Stochastic partial differential equations. Probability and its Applications. Birkh¨ auser Boston Inc., Boston, MA, 1996. A modeling, white noise functional approach. [32] Arnulf Jentzen and Peter Kloeden. Taylor expansions of solutions of stochastic partial differential equa- tions with additive noise. Ann. Probab., 38(2):532–569, 2010. [33] Arnulf Jentzen and Peter E. Kloeden. Overcoming the order barrier in the numerical approximation of stochastic partial differential equations with additive space-time noise. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 465(2102):649–667, 2009. [34] S. Kaligotla and S. V. Lototsky. Wick product in stochastic burgers equation: a curse or a cure? Preprint, 2010. [35] Gopinath Kallianpur and Jie Xiong. Stochastic differential equations in infinite-dimensional spaces. Institute of Mathematical Statistics Lecture Notes—Monograph Series, 26. Institute of Mathematical Statistics, Hayward, CA, 1995. Expanded version of the lectures delivered as part of the 1993 Barrett Lectures at the University of Tennessee, Knoxville, TN, March 25–27, 1993, With a foreword by Balram S. Rajput and Jan Rosinski. [36] Peter E. Kloeden and Eckhard Platen. Numerical solution of stochastic differential equations, volume 23 of Applications of Mathematics (New York). Springer-Verlag, Berlin, 1992. [37] Yuri G. Kondratiev, Peter Leukert, and Ludwig Streit. Wick calculus in Gaussian analysis. Acta Appl. Math., 44(3):269–294, 1996. [38] Mih´ aly Kov´ acs, Stig Larsson, and Fredrik Lindgren. Strong convergence of the finite element method with truncated noise for semilinear parabolic stochastic equations with additive noise. Numer. Algo- rithms, 53(2-3):309–320, 2010. [39] Mih´ aly Kov´ acs, Stig Larsson, and Fardin Saedpanah. Finite element approximation of the linear sto- chastic wave equation with additive noise. SIAM J. Numer. Anal., 48(2):408–427, 2010. [40] John D. Kraus and Daniel Fleisch. Electromagnetics. McGraw-Hill Series in Electrical and Computer Engineering. McGraw-Hill, New York, 1991. [41] S. G. Kre˘ın, Yu. ¯I. Petun¯ın, and E. M. Sem¨enov. Interpolation of linear operators, volume 54 of Transla- tions of Mathematical Monographs. American Mathematical Society, Providence, R.I., 1982. Translated from the Russian by J. Sz˝ ucs. [42] Hui-Hsiung Kuo. White noise distribution theory. Probability and Stochastics Series. CRC Press, Boca Raton, FL, 1996. 123 [43] C. Y. Lee, B. L. Rozovskii, and H. M. Zhou. Randomization of forcing in large systems of PDEs for improvement of energy estimates. Multiscale Model. Simul., 8(4):1419–1438, 2010. [44] Chia Ying Lee and Boris Rozovskii. A stochastic finite element method for stochastic parabolic equations driven by purely spatial noise. Commun. Stoch. Anal., 4(2):271–297, 2010. [45] S. V. Lototsky and B. L. Rozovskii. A unified approach to stochastic evolution equations using the Skorokhod integral. Teor. Veroyatn. Primen., 54(2):288–303, 2009. [46] S. V. Lototsky, B. L. Rozovskii, and D. Seleˇsi. A note on generalized malliavin calculus. Preprint, 2010. [47] Sergey Lototsky and Boris Rozovskii. Stochastic differential equations: a Wiener chaos approach. In From stochastic calculus to mathematical finance, pages 433–506. Springer, Berlin, 2006. [48] Sergey V. Lototsky and Boris L. Rozovskii. Stochastic parabolic equations of full second order. In Topics in stochastic analysis and nonparametric estimation, volume 145 of IMA Vol. Math. Appl., pages 199– 210. Springer, New York, 2008. [49] Sergey V. Lototsky and Boris L. Rozovskii. Stochastic partial differential equations driven by purely spatial noise. SIAM J. Math. Anal., 41(4):1295–1322, 2009. [50] Paul Malliavin. Stochastic analysis, volume 313 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]. Springer-Verlag, Berlin, 1997. [51] Leonard Mandel and Emil. Wolf. Optical Coherence and Quantum Optics. Cambridge University Press, Cambridge, UK, 1995. [52] H. Manouzi. A finite element approximation of linear stochastic PDEs driven by multiplicative white noise. Int. J. Comput. Math., 85(3-4):527–546, 2008. [53] R. Mikulevicius and B. L. Rozovskii. Stochastic Navier-Stokes equations for turbulent flows. SIAM J. Math. Anal., 35(5):1250–1310, 2004. [54] R. Mikulevicius and B. L. Rozovskii. Global L2 -solutions of stochastic Navier-Stokes equations. Ann. Probab., 33(1):137–176, 2005. [55] R. Mikulevicius and B. L. Rozovskii. On quantized stochastic navier-stokes equations. Preprint, 2010. [56] David Nualart. The Malliavin calculus and related topics. Probability and its Applications (New York). Springer-Verlag, Berlin, second edition, 2006. [57] B. L. Rozovski˘ı. Stochastic evolution systems, volume 35 of Mathematics and its Applications (Soviet Se- ries). Kluwer Academic Publishers Group, Dordrecht, 1990. Linear theory and applications to nonlinear filtering, Translated from the Russian by A. Yarkho. [58] M. A. Shubin. Pseudodifferential operators and spectral theory. Springer Series in Soviet Mathematics. Springer-Verlag, Berlin, 1987. Translated from the Russian by Stig I. Andersson. [59] Roger Temam. Navier-Stokes equations and nonlinear functional analysis, volume 66 of CBMS-NSF Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 1995. 124 [60] Roger Temam. Navier-Stokes equations. AMS Chelsea Publishing, Providence, RI, 2001. Theory and numerical analysis, Reprint of the 1984 edition. [61] Thomas Gorm Theting. Solving parabolic Wick-stochastic boundary value problems using a finite ele- ment method. Stoch. Stoch. Rep., 75(1-2):49–77, 2003. [62] Vidar Thom´ee. Galerkin finite element methods for parabolic problems, volume 25 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006. [63] Gjermund V˚ age. Variational methods for PDEs applied to stochastic partial differential equations. Math. Scand., 82(1):113–137, 1998. ´ [64] John B. Walsh. An introduction to stochastic partial differential equations. In Ecole d’´et´e de probabilit´es de Saint-Flour, XIV—1984, volume 1180 of Lecture Notes in Math., pages 265–439. Springer, Berlin, 1986. [65] Xiaoliang Wan, Boris Rozovskii, and George Em Karniadakis. A stochastic modeling methodology based on weighted Wiener chaos and Malliavin calculus. Proc. Natl. Acad. Sci. USA, 106(34):14189–14194, 2009. [66] R. Weissleder. Molecular imaging in cancer. Science, 312:1168–1171, 2006. [67] Dongbin Xiu and George Em Karniadakis. Modeling uncertainty in steady state diffusion problems via generalized polynomial chaos. Comput. Methods Appl. Mech. Engrg., 191(43):4927–4948, 2002. [68] Dongbin Xiu and George Em Karniadakis. The Wiener-Askey polynomial chaos for stochastic differen- tial equations. SIAM J. Sci. Comput., 24(2):619–644 (electronic), 2002. [69] Dongbin Xiu and George Em Karniadakis. Modeling uncertainty in flow simulations via generalized polynomial chaos. J. Comput. Phys., 187(1):137–167, 2003. [70] Yubin Yan. Galerkin finite element methods for stochastic parabolic partial differential equations. SIAM J. Numer. Anal., 43(4):1363–1384 (electronic), 2005. [71] A. Yodh and B. Chance. Spectropscopy and imaging with diffusing light. Phys. Today, 48:34–40, 1995. 125