Numerical Methods for Partial Differential Equations with Random - - PowerPoint PPT Presentation
Numerical Methods for Partial Differential Equations with Random Data Howard Elman University of Maryland Outline I. Problem statement and discretization Example: diffusion equation with random diffusion coefficient Discretization by
1
coefficient
Forcing function f Boundary data Viscosity ν in Navier-Stokes equations
D N D D d
n u a g u R f u a D D D D D \
) ( ,
in ) (
a = a(x,ω) a random field For each fixed x, a(x,ω) a random variable Uncertainty / randomness: Other possibly uncertain quantities :
D
g
2
div grad ) grad (
2
u f p u u u
3
Random field a(x,ω) Mean μ(x) = E(a(x,·)) Variance Covariance function c(x,y) = E( (a(x,·)- μ(x)) (a(y,·)- μ(y)) ) is finite
: , D y x
2 2)
) , ( ( ) ( x a E x
D in ) (
u a
4
2 1
a
Problem is well-posed Assumptions:
Sample a(x,ω) at all x , solve in usual way Standard weak formulation: find such that
1 D E
1
0 D
E
D D
for all Multiple realizations (samples) of a(x,·) Multiple realizations of u Statistical properties of u Problem: convergence is slow, requires many solves
5
= identically distributed uncorrelated random variables with mean 0 and variance 1
1
r r r r
D in ) (
u a
Covariance function is finite random field (diffusion coefficient) has Karhunen-Loève expansion:
r
D
r r x
a ), ( mean )) , ( ( ) ( ) ( x a E x x a
= eigenfunctions/eigenvalues of covariance operator
6
Requires: m large enough so that the fluctuation of a is well-represented, i.e. is small
1 1 / m
~ Principal components analysis
m r r r r
1
D in ) (
u a
Truncated Karhunen-Loève expansion:
7
More precisely: error from truncation is
2 1 2
| | | | D D
m j j
Choose m to make this small
Introduce a weak formulation analogous to finite elements in space that handles the “stochastic” component of the problem
Devise a special strategy for sampling ξ that converges more quickly than Monte Carlo simulation; derived from interpolation
8
m r r r r
1
Ghanem, Spanos, Babuška, Deb, Oden, Matthies, Keese, Karniadakis, Xiu, Hesthaven, Tempone, Nobile, Webster, Schwab, Todor, Ernst, Powell, Furnival, E., Ullmann, Rosseel, Vandewalle
Probability space (Ω,F, P)
2 P
{square integrable functions wrt dP(ω)} Inner product on :
2 P
Use to concoct weak formulation on product space
2 1 P E
Find such that
2 1 P E
for all
2 1 P E
D
Solution u=u(x,ω) is itself a random field
9
Truncated Karhunen-Loève expansion
m r r r r
1
Stochastic weak formulation uses
) (
) ( ) , ( ) ( ) , (
D D
d dx v u x a dP dx v u a v u a
ξ plays the role of a Cartesian coordinate Bilinear form entails integral over image of random variables ξ Require joint density function associated with ξ
10
11
D D
d dx fv d dx v u x a ) ( ) ( ) , (
Find such that
2 1 P E
for all
2 1 P E
Like an ordinary Galerkin (or Petrov-Galerkin) problem on a (d+m)-dimensional “continuous” space d = dimension of spatial domain m = dimension of stochastic space
12
D D
Finite dimensional spaces:
for example: piecewise linear on triangles
for example: polynomial chaos = m-variate Hermite polynomials (orthogonal wrt Gaussian measure)
x
N j j h
H S
1 1
} { by spanned ), (D
N l l p
L T
1 2
} { by spanned ), (
p h hp hp hp hp
x
N j N l l j jl hp
1 1
Discrete weak formulation:
13
2 2
) (
2
L Tp
Underlying space:
2 2 1 1 M M
Let polynomial of degree j orthogonal wrt
) ( k k j
k
Examples: if ρk ~ Gaussian measure Hermite polynomials ρk ~ uniform distribution Legendre polynomials Any ρk can be handled computationally (Gautschi) Rys polynomials
)} ( ) ( ) ( {
) ( 2 ) 2 ( 1 ) 1 (
2 1
m m j j j
m
q q q
spanned by Orthogonality of basis functions sparsity of coefficient matrix
q k Γ kq q l r lq r q l lq k j r jk r k j jk r r
r m 1 r D D D
Properties of A:
} and { }
r
G
r
A
14
m r r r r
x a x a x a
1
) ( ) ( ) ( )) ( , (
15
) (
2
L Tp )} ( ) ( ) ( {
) ( 2 ) 2 ( 1 ) 1 (
2 1
m m j j j
m
q q q
spanned by Full tensor product basis:
,m , i p ji 1 ,
Dimension:
m
p ) 1 (
“Complete” polynomial basis:
p j j j
m
2 1
Dimension:
m+p p = (m+p)! m! p! Too large
2 1 N
More manageable Order these in a systematic way
16
) (
2
L Tp )} ( ) ( ) ( {
) ( 2 ) 2 ( 1 ) 1 (
2 1
m m j j j
m
q q q
spanned by
3 ) ( , 1 ) ( , ) ( , 1 ) (
3 3 2 2 1
H H H H
Orthogonal (Hermite) polynomials in 1D: “Complete” polynomial basis:
p j j j
m
2 1 2 5 1 3 1 4 2 1 3 1 2 1
) ( 3 ) ( 1 ) ( ) ( 1 ) (
2 3 2 10 1 2 2 9 2 2 8 2 2 1 7 2 1 6
3 ) ( ) 1 ( ) ( ) 1 ( ) ( ) 1 ( ) ( ) (
m=2, p=3
m+p p
5 2 = =10 Gives basis set:
17
Nξ= (m+p)! m!p! 10! 6!4! = = 210 For m-variate polynomials of total degree p:
using orthogonality of stochastic basis functions
18
N j j j l M l l hp
1 1 1 1
First moment of u (expected value): Free! Similarly for second moment / covariance
1 1 1
l N l l l N l N j l j jl hp
x
19
1 1 1
l N l l l N l N j l j jl hp
x
) ) , ( ( x u P
hp
E.g.:
at some point x Sample ξ
1 N l l l hp
Evaluate Precomputed Repeat Not free, but no solves required
Given as above
20
r m r r r
x a x a x a
1
) ( ) ( ) , (
Let ξ be a specified realization (~ Monte Carlo) Weak formulation:
D D
r m r r r
1
Discretize in space in usual way. Stochastic collocation: choose special set from considerations of interpolation
) ( ) 2 ( ) 1 (
N
Advantage: Spatial systems are decoupled
Given and v(ξ), consider an interpolant
21
) ( ) 2 ( ) 1 ( N
where Lagrange interpolating polynomial
) ( jk j k
D D
r m r r r
1
If solves the discrete (in space) version of
1 ) ( k N k k h hp
) (k h
u ), ( ) ( ) ( ) )( (
1 ) (
N k k k
) (k
with then the collocated solution is
22
k N k k h hp
1 ) (
Not free but can be precomputed
Obtained by sampling, cheap ) ( ) ( ) , (
1 ) ( k N k k h hp
L x u x u Solution
1D interpolating polynomial
p k j
k j
23
1 ) (
N k k k
One choice of
2 1
2 1
m k k k k k
m
Advantage: easy to construct Disadvantage: “curse of dimensionality,” dimension = (p+1)
m
Given: 1D interpolation rule
) ( 1 ) ( ) ( ) ( k j m j k j k k
k
Multidimensional rule above is induced by fully populated multidimensional grid . Derived from (1D) grid
) ( ) ( 1 ) ( k m k k
k
) ( ) 2 ( ) 1 ( m
Alternative: multidimensional sparse grid (Smolyak)
p i i m p i i i
m m
1 2 1
1 ) ( ) ( ) (
) (
k k
24
25
Example of sparse grid for m=3, p=16 For v of the form interpolating function takes the form
2 2 ) 1 ( ) ( 1 1 ) 1 ( ) (
2 2 1 1 1
i i p i i i i
m
) 1 ( ) ( m m i i
m m
2 2 1 1 m m
For sparse grid and a tensor product polynomial of total degree at most p,
Theorem (Novak, Ritter, Wasilkowski, Wozniakowski) That is: sparse grid interpolation evaluates the set of complete m-variate polynomials exactly Overhead: number of sparse grid points to achieve this (= # stochastic dof) is larger than for Galerkin
m+p p ≈ 2
p
m+p p vs.
p j j j q q q v
m m m j j j
m
2 1 ) ( 2 ) 2 ( 1 ) 1 (
), ( ) ( ) ( ) (
2 1
26
Monte-Carlo:
)) ( ) ( ( )) ( ) ( ( ) ( ) (
hp h h hp
u E u E u E u E u E u E
, r r c
p
1
2
2 1
2 1
h s h h h s
(Babuška, Tempone, Zouraris, Nobile, Webster) Convergence is slow wrt number of samples but independent of number of random variables m Stochastic Galerkin and Collocation: Rule of thumb: the same p gives the same error (for all versions of SG and collocation) More dof for collocation than SG Exponential in polynomial degree p Constants (c2 , r) depend on m
27
Monte-Carlo methods: Many samples needed for statistical quantities Many systems to solve Systems are independent Statistical quantities are free (once data is accumulated) Stochastic Galerkin methods: One large system to solve Statistical quantities are free or (relatively) cheap Stochastic collocation methods: Systems are independent Fewer systems than Monte Carlo More degrees of freedom than Galerkin Statistical quantities are (relatively) cheap
s r r h h s
x u s u E
1 ) (
) ( 1 ) ( With s realizations: Convergence is slow but independent of m Similar convergence behavior Faster than MC Depends on m
28
For both: compute a discrete solution, a random field
) ( ) ( ) ( ) ( ) , (
1 1 1 l N l l N l N j l j jl hp
L x u L x u x u
x
Stochastic Galerkin:
N l l l N l N j l j jl hp
x u x u x u
x
1 1 1
) ( ) ( ) ( ) ( ) , (
Stochastic Collocation: Postprocess to get statistics
29
Stochastic Galerkin: Solve one large system of order Nx x Nξ m+p p Nξ=( ) Frequently cited as a problem for this methodology Stochastic Collocation: Solve Nξ “ordinary” algebraic systems (of order Nx), one for each sparse grid point Here:
) ( ) (
Galerkin p n collocatio
Some savings possible
30
spatial discretization parameter h
, ,
) ( ) ( h h r r
A A A A
31
(E. & Furnival)
, ,
) 2 ( ) 2 ( h h r r
A A A A
spatial discretization parameter 2h Solving Au=f
Ω r q l lq r D k j r jk r r r
r m 1 r
Develop MG algorithm for spatial component of the problem
32
) ( 1 ) ( ) ( 1 ) (
) (
h h h h
f Q u A Q I u
for i=0,1,… for j=1:k k smoothing steps end Restriction Solve Coarse grid correction Prolongation end
) (
) ( ) ( ) ( ) 2 ( h h h h
u A f r R
) 2 ( ) 2 ( ) 2 ( h h h
r c A
) 2 ( ) ( ) ( h h h
c u u P
Prolongation and restriction:
T T
P R R I P I , , P R P induced by natural inclusion in spatial domain Let Q = smoothing operator
,
) (
N Q A h
) ( , ) ( ] ) ( [
) (
2 ) ( 1 ) ( increases k A k h h
k y y k y A Q I A
h
Establish approximation property and smoothing property
y y C y A A
h
A h h
] ) ) ( [
2 1 ) 2 ( 1 ) (
) (
R P(
33
) ( ) ( 1 ) ( 1 ) 2 ( 1 ) ( ) 1 (
] ) ( [ )] ) ) [(
i k h h h h i
e A Q I A A A e R P(
Error propagation matrix:
) ( ) ( ) (
|| || ) ( || ] ) ( [ || || ] ) ( [ )] ) ) ( || || ||
) ( 2 ) ( ) ( 1 ) ( ) ( ) ( 1 ) ( 1 ) 2 ( 1 ) ( ) 1 (
h h h
A i i k h h A i k h h h h A i
e k C e A Q I A C e A Q I A A A e R P( Analysis is:
34
2 ) ( ) ( 2 ) ( 2 2 2 2 / 1 2 2 2 ) 2 ( ) ( 1 ) 2 ( 1 ) (
2 2 2 ) ( ) (
L L L a h a h h h h h a h h A h h A h h
h h
D D D
Approximability Property of mass matrix Regularity “Standard” MG analysis for deterministic problem:
35
Introduce semi-discrete space
p
T H ) (
1 0 D
Weak formulation
p p p p p p
u T H v v v u a Solution ) ( all for ) ( ) , (
1 0 D
) ( ) ( 2
2 2
L L p a hp p
u D Ch u u
D
Similarly for other steps used for deterministic analysis Discrete stochastic space
a h p a p h a p h hp A h h
u u u u u u y A A
h
2 , 2 1 ) 2 ( 1 ) (
] ) ) ( [
) (
R P(
Approximation (in 2D): Established using best approximation property of and interpolant
hp
u ) , ( ) , ( ~
j p j p
x u x u
36
spatial discretization size h
derives from basis of multivariate polynomials of total degree p, orthogonal wrt probability measure ρ(ξ)dξ
r
G ,
m 1 r r r r
G a G a G Maximum eigenvalue η = max root of orthogonal polynomial, bounded for bounded measure CG iteration is an option ), ( ) ( ) (
x1 1 x1 1 1 x 1 x1 1
m 1 r m 1 r r r r r
a a G a a
37
Polynomial degree # terms (m) in KL-expansion m=1 m=2 m=3 m=4 p=1 8 8 8 8 p=2 8 8 8 8 p=3 9 9 9 9 p=4 9 10 10 10
h=1/16
Polynomial degree m=1 m=2 m=3 m=4 p=1 7 7 8 8 p=2 8 8 8 8 p=3 8 8 9 9 p=4 9 9 9 9
h=1/32
Ω r q l lq r k j r jk r r r
r m 1 r D
Solving Au=f Preconditioner for use with CG:
k j D
(Kruger, Pellissetti, Ghanem) Deterministic diffusion, from mean
38
Theorem : For μ constant, the Rayleigh quotient satisfies
1 r m r r
Consequence: dictates convergence of PCG
m r r r r
1
Recall
m 1 r r r
39
Sketch of Proof
1 r m r r
40
c(p) bounded by largest root of scalar orthogonal polynomial
m 1 r r r
D
r r r
D
r
r
In spatial domain: From stochastic component: as above
Replace action of with multigrid preconditioner
1
MG MG
,
(Le Maitre, et al.) Analysis:
MG MG
Spectral equivalence
to diffusion operator
2 1 1 2
41
f u a ) (
# Samples s Max SFEM 100 1000 10,000 40,000 Mean .06311 .06361 .06330 .06313 .06313 Variance 2.360(-5) 2.161(-5) 2.407(-5) 2.258(-5) 2.316(-5)
Compare Monte-Carlo simulation with SFEM, for N.B.: No negative samples of diffusion obtained in MC Solve one system
Solve s systems of size 225
42
Recall, for stochastic collocation
43
1 ) ( k N k k h hp
Discrete solution Obtained by solving
D D
r m r r r
1
For set of samples situated in a sparse grid
) (k
Advantage of this approach: simpler (decoupled) systems Disadvantage: larger stochastic space for comparable accuracy larger by factor approximately
p
m (#KL) p Galerkin Collocation Sparse Collocation Tensor 4 1 2 3 4 5 15 35 70 9 41 137 401 16 81 256 625 10 1 2 3 4 11 66 286 1001 21 221 1582 8,801 1024 59,049 1,048,576 9,765,625 30 1 2 3 4 31 496 5456 46,376 61 1861 37,941 582,801 1.07(9) 2.06(14) 1.15(18) 9.31(20) ~ size of coarse grid space for MG / Version 1 # systems for collocation MG / Version II
solution costs
Solution algorithm for both discretizations: Preconditioned conjugate gradient with mean-based preconditioning, using AMG for the approximate diffusion solve ¹Estimated using a high-degree (p=10) Galerkin solution. (E., Miller, Phipps, Tuminaro)
46
46
polynomial degree for SG “level” for collocation produces comparable errors
p=1 p=2 p=3 p=4 p=5
Accuracy: for fixed m=4: similar p=
p=3 p=5 p=4 p=2 p=1 p=4 p=5 p=6 p=3 p=2 p=1
Degrees of freedom
Performance:
M=3 Error Error
47
47 47
p m=5 m=10 m=15 m=5 m=10 m=15 1 .058 .147 .32- .069 .163 .286 2 .269 1.20 3.80 .532 2.13 5.08 3 1.20 13.14 51.45 2.41 16.99 57.98 4 3.50 53.79 168.11 8.31 102.60 493.04 5 6.51 117.73 24.56 515.75
CPU times for larger m = #KL terms:
Galerkin Collocation
Performed on a serial machine with C code and CG/AMG code from Trilinos Observation: Galerkin faster, more so as number of stochastic variables (KL terms) grows
) , ( min
) , (
x c
e a x a
Nonlinear
48
For the problem discussed, based on a KL expansion, has a linear dependence on the stochastic variable ξ Other models have nonlinear dependence. For example
m r r r r
x a x a x c
1
) ( ) ( ) , (
In particular: coercivity is guaranteed with this choice For Gaussian c, called a log-normal distribution
49
) ( ) ( ) ( ) , (
1 M r r r r
x a x a x a
For stochastic Galerkin, need a finite term expansion for a
Note: not ξr
M 1 r r r
matrix
j i r ij r
Less sparse More importantly: # terms M will be larger perhaps as large as 2Nξ mvp will be more expensive
comes from for each sparse grid point
50
Collocation is less dependent on this expansion
D
k )
) ( ) (k
) (k
Many matrices to assemble, but mvp is not a difficulty
51
coefficients
models that simulate sampling at potentially lower cost
stochastic collocation method, were presented, each with some advantages
continues in this direction