Docstoc

Finite element and wavelet approximation of a parabolic stochastic

Document Sample
Finite element and wavelet approximation  of a parabolic stochastic Powered By Docstoc
					 Finite element and wavelet approximation of a
parabolic stochastic partial differential equation
                       Stig Larsson
            Department of Mathematical Sciences
             Chalmers University of Technology
                       ¨
                  and Goteborg University
                     ¨
                    Goteborg, Sweden




                                                  EFEF, Marseille, 2007 – p. 1/2
   Outline

Parabolic stochastic partial differential equation with additive noise
 du − ∆u dt = dW,       x ∈ D, t > 0


  u = 0,                 x ∈ ∂D, t > 0


  u(0) = u0

     Abstract framework
         Mild solution
         Finite element approximation
         Error estimates
     Approximation of the noise
         Survey of the literature
         Wavelet multiresolution




                                                                     EFEF, Marseille, 2007 – p. 2/2
   Co-workers

Yubin Yan
Mihaly Kovacs
Li Bin
Asta Hall
Fredrik Lindgren




                   EFEF, Marseille, 2007 – p. 3/2
   Abstract framework

  du + Au dt = dW,            t>0
  u(0) = u0

     H = L2 (D), · , (·, ·), D ⊂ Rd , bounded domain
                               1
     A = −∆, D(A) = H 2 (D) ∩ H0 (D)
     u(t), H-valued random process on probability space (Ω, F, P)
     W (t), H-valued Wiener process
     E(t) = e−tA , analytic semigroup generated by −A

Mild solution (stochastic convolution):
                      t
u(t) = E(t)u0 +           E(t − s) dW (s),   t≥0
                  0




                                                                EFEF, Marseille, 2007 – p. 4/2
   Wiener process

                   ∞                                               ∞
                        1/2                                              1/2
    W (t) =            γl βl (t)el ,             W (t, x, ω) =          γl     βl (t, ω)el (x)
               l=1                                                l=1

    Qel = γl el ,        γl > 0,          {el } ON basis
    covariance operator Q : H → H, self-adjoint, positive definite,
    bounded, linear operator, cov(W (t)) = tQ
    βl (t), independent identically distributed, real-valued, Brownian
    motions

Two important cases:
    Tr(Q) < ∞: W (t) is an H-valued Wiener process
        ∞                         2       ∞                       ∞
               1/2
    E         γl       βl (t)el       =         γl Eβl (t)2 = t         γl = t Tr(Q) < ∞
        l=1                               l=1                     l=1

    Q = I: W (t) is not H-valued, since Tr(I) = ∞, “white noise”



                                                                                                 EFEF, Marseille, 2007 – p. 5/2
     Stochastic integral

                                t
u(t) = E(t)u0 +                     E(t − s) dW (s),          t≥0
                            0

                                                                        t
      We can define the stochastic integral                                  B(s) dW (s)
                                                                    0
      if B(s)Q1/2 is a Hilbert-Schmidt operator on H
      Isometry property
                  t                        2            t
      E               B(s) dW (s)              =E           B(s)Q1/2        2
                                                                            HS   ds
              0                                     0

Hilbert-Schmidt operator B, B                       HS      <∞
     2                ∞               2
 B   HS   =           l=1   Bϕl           , {ϕl } arbitrary orthonormal basis in H

Da Prato and Zabczyk, Stochastic Equations in Infinite Dimensions



                                                                                          EFEF, Marseille, 2007 – p. 6/2
     Regularity


|v|β = Aβ/2 v ,           ˙
                          H β = D(Aβ/2 ),               β∈R
     2
 v         ˙
     L2 (Ω,H β )
                   = E(|v|2 ) =
                          β                     |Aβ/2 v|2 dx dP(ω),       β∈R
                                       Ω    D



Theorem 1. If         A(β−1)/2 Q1/2         HS   < ∞ for some β ≥ 0, then
 u(t)         ˙
        L2 (Ω,H β )   ≤C        u0          ˙
                                      L2 (Ω,H β )   + A(β−1)/2 Q1/2       HS



Remark:                           ˙                   ˙
                   W (t) ∈ L2 (Ω, H −(1−β) ) ⊂ L2 (Ω, H −1 )

Two cases:

      If Q1/2        2
                     HS   = Tr(Q) < ∞, then β = 1
                                         ∂2
      If Q = I, d = 1, A =             − ∂x2 ,   then A(β−1)/2       HS   < ∞ for β < 1/2
                                        −(1−β)
        A(β−1)/2       2
                       HS   =        j λj           ≈    j   j −(1−β)2/d < ∞ iff d = 1, β < 1/2

                                                                                        EFEF, Marseille, 2007 – p. 7/2
   The finite element method


     triangulations {Th }0<h<1 , mesh size h
     finite element spaces {Sh }0<h<1
           1       ˙
     Sh ⊂ H0 (D) = H 1
     Sh continuous piecewise linear functions
     Ah : Sh → Sh , discrete Laplacian, (Ah ψ, χ) = (∇ψ, ∇χ), ∀χ ∈ Sh
     Ph : L2 → Sh , orthogonal projection, (Ph f, χ) = (f, χ), ∀χ ∈ Sh

  uh (t) ∈ Sh ,   uh (0) = Ph u0
  duh + Ah uh dt = Ph dW

More rigorously, with Eh (t) = e−tAh ,

 uh (t) ∈ Sh , uh (0) = Ph u0

                               t
 uh (t) = Eh (t)Ph u0 +
                                  Eh (t − s)Ph dW (s)
                           0


                                                                    EFEF, Marseille, 2007 – p. 8/2
   Error estimates for the deterministic problem


   ut + Au = 0,            t>0                uh,t + Ah uh = 0,     t>0
   u(0) = v                                   uh (0) = Ph v
u(t) = E(t)v                               uh (t) = Eh (t)Ph v

Denote
Fh (t)v = Eh (t)Ph v − E(t)v,                       |v|β = Aβ/2 v

We have, for 0 ≤ β ≤ 2,

      Fh (t)v ≤ Chβ |v|β , t ≥ 0
             t                      1/2
                           2
                 Fh (s)v       ds         ≤ Chβ |v|β−1 , t ≥ 0
         0

V. Thomée, Galerkin Finite Element Methods for Parabolic Problems



                                                                          EFEF, Marseille, 2007 – p. 9/2
   Strong convergence in L2 norm


Theorem 2. If    A(β−1)/2 Q1/2      HS   < ∞ for some β ∈ [0, 2], then
 uh (t) − u(t)   L2 (Ω,H)   ≤ Chβ    u0          ˙
                                           L2 (Ω,H β )   + A(β−1)/2 Q1/2   HS


                             2
Recall:     uh (t) − u(t)    L2 (Ω,H)    = E( uh (t) − u(t) 2 )

Two cases:

     If Q1/2     2
                 HS   = Tr(Q) < ∞, then the convergence rate is O(h).
                                 ∂2
     If Q = I, d = 1, A =      − ∂x2 ,   then the rate is almost O(h1/2 ).

No result for Q = I, d ≥ 2.




                                                                                EFEF, Marseille, 2007 – p. 10/2
   Strong convergence: proof

                         t
u(t) = E(t)u0 +              E(t − s) dW (s)
                     0
                                   t
uh (t) = Eh (t)Ph u0 +                 Eh (t − s) Ph dW (s)
                               0

Fh (t) = Eh (t)Ph − E(t)
                                            t
uh (t) − u(t) = Fh (t)u0 +                      Fh (t − s) dW (s) = e1 (t) + e2 (t)
                                        0


 Fh (t)u0 ≤ Chβ |u0 |β                 (deterministic error estimate)

=⇒      e1 (t)   L2 (Ω,H)     ≤ Chβ u0                    ˙
                                                    L2 (Ω,H β )




                                                                                      EFEF, Marseille, 2007 – p. 11/2
     Strong convergence: proof

               t                                             t
                                       2
E
                   B(s) dW (s)             =E                    B(s)Q1/2   2
                                                                                         ds (isometry)
                                                                            HS
            0                                             0
          t                           1/2

                             2
                   Fh (s)v       ds                ≤ Chβ |v|β−1 , with v = Q1/2 ϕl (deterministic)
       0
=⇒
                                                t                            2                t
                2
    e2 (t)      L2 (Ω,H)      =E                    Fh (t − s) dW (s)            =                Fh (t − s)Q1/2    2
                                                                                                                    HS     ds
                                            0                                             0
                                   ∞            t                                                   ∞
                              =                         Fh (t − s)Q1/2 ϕl        2
                                                                                     ds ≤ C              h2β |Q1/2 ϕl |2
                                                                                                                       β−1
                                  l=1       0                                                      l=1
                                                    ∞
                              = Ch2β                      A(β−1)/2 Q1/2 ϕl           2
                                                                                         = Ch2β A(β−1)/2 Q1/2                2
                                                                                                                             HS
                                                l=1


If Tr(Q) < ∞, we may choose β = 1, otherwise β < 1.


                                                                                                                   EFEF, Marseille, 2007 – p. 12/2
   Time discretization

  du + Au dt = dW ,                 t>0
  u(0) = u0

The implicit Euler method: k = ∆t, tn = nk, ∆W n = W (tn ) − W (tn−1 )

  U n ∈ Sh ,   U 0 = Ph u0
  U n − U n−1 + kAh U n = Ph ∆W n ,

U n = Ekh U n−1 + Ekh Ph ∆W n ,             Ekh = (I + kAh )−1

                      n
                               n−j+1
       n
U n = Ekh Ph u0 +             Ekh    Ph ∆W j
                    j=1


                              tn
u(tn ) = E(tn )u0 +                E(tn − s) dW (s)
                          0


                                                                 EFEF, Marseille, 2007 – p. 13/2
   Error estimates for deterministic parabolic problem

             n
Denote Fn = Ekh Ph − E(tn )
We have the following estimates for 0 ≤ β ≤ 2:

     Fn v ≤ C(k β/2 + hβ )|v|β
          n                1/2
                       2
      k         Fj v             ≤ C(k β/2 + hβ )|v|β−1
          j=1




                                                          EFEF, Marseille, 2007 – p. 14/2
   Strong convergence


Theorem 3. If A(β−1)/2 Q1/2         HS   < ∞ for some β ∈ [0, 2], then, with
en = U n − u(tn ),

 en   L2 (Ω,H)   ≤ C(k β/2 + hβ )    u0         ˙
                                          L2 (Ω,H β )   + A(β−1)/2 Q1/2   HS




J. Printems (2001) (only time-discretization)
Y. Yan, BIT (2004), SIAM J. Numer. Anal (2005)




                                                                               EFEF, Marseille, 2007 – p. 15/2
Approximating the noise


 No result for Q = I, d ≥ 2.
 Increments ∆W j are not directly computable: infinite series with
 unknown eigenfunctions el .




                                                             EFEF, Marseille, 2007 – p. 16/2
   Approximating the noise: white noise in 1-D


Q = I,   d = 1,        D = (0, 1)
          ∂u        ∂2W
SPDE:        + Au =
          ∂t        ∂x∂t
Piecewise constant approximation:
                       tj+1        xi+1                              tj+1    xi+1
∂2W      1                                ∂2W           1
     :=                                        dx dt =                              dW
∂x∂t    ∆x∆t       tj         xi          ∂x∂t         ∆x∆t        tj       xi

                   tj+1        xi+1
         1
ηij := √                              dW ∈ N (0, 1)
        ∆x∆t      tj          xi

                                              N   M
∂2W           ∂2W            1                             √
     (x, t) ≈      (x, t) =                             ηij ∆x∆t χ[xi ,xi+1 ] (x) χ[tj ,tj+1 ] (t)
∂x∂t          ∂x∂t          ∆x∆t              i=0 j=0

         ∂u
          ˆ        ∂2W
PDE:           u
            + Aˆ =
         ∂t        ∂x∂t

                                                                                         EFEF, Marseille, 2007 – p. 17/2
   Approximating the noise: white noise in 1-D


         ∂u
          ˆ        ∂2W
PDE:           u
            + Aˆ =
         ∂t        ∂x∂t
                                                 ˆ
Finite element or finite difference approximation U :
                            1/p
   ˆ
 E|Uij − u(xi , tj )|p            ≤ C ∆t1/4 + ∆x1/2

Gyöngy (1999)
Allen, Novosel, and Zhang (1998)
Davie and Gaines (2000) (also lower bounds)
Walsh (2005)
Kossioris and Zouraris

Proof technique:
                1                               t       1
                                                                           ∂2W
u(x, t) =           G(x, y, t)u0 (y) dy +                   G(x, y, t − s)      (y, s) dy ds
            0                               0       0                      ∂x∂t


                                                                                         EFEF, Marseille, 2007 – p. 18/2
   Approximating the noise: colored noise in rectangle


 A(β−1)/2 Q1/2   HS   < ∞,       D = (0, 1)d
Explicit eigenfunctions:
Ael = λl el
Qel = γl el
know: λl ≈ l2/d ,     assume: γl ≈ l−α
                          ∞                 ∞
 A(β−1)/2 Q1/2   2
                 HS   =         λβ−1 γl ≈
                                 l                l(β−1)d/2−α < ∞
                          l=1               l=1

if α > 1 − (β − 1)d/2




                                                                    EFEF, Marseille, 2007 – p. 19/2
   Approximating the noise: colored noise in rectangle


 A(β−1)/2 Q1/2      HS   < ∞,     D = (0, 1)d ,       explicit eigenfunctions

1. Spectral Galerkin approximation in x, difference method in t
            N                                     N
                   1/2
W N (t) =         γj βj (t)ej ,   uN (x, t) =           ˆ
                                                        uj (t)ej (x)
            j=1                                   j=1


Shardlow (1999)
Lord and Rougemont (2003)
Müller-Gronbach and Ritter (2004), (2005) (lower bounds)




                                                                                EFEF, Marseille, 2007 – p. 20/2
     Approximating the noise: colored noise in rectangle


 A(β−1)/2 Q1/2        HS   < ∞,        D = (0, 1)d ,   explicit eigenfunctions

2. Finite element approximation in x, difference method in t
              J
                     1/2
W J (t) =           γj βj (t)ej
              j=1

duJ + Ah uJ dt = Ph dW J
  h       h
                               1/2
 E   uJ (t)
      h       − uh (t)     2
                                     ≤ Chβ A(β−1)/2 Q1/2    HS

if J ≥ Nh = dim(Sh )

Du and Zhang (2002) (d = 1)
Yan (2004)




                                                                                 EFEF, Marseille, 2007 – p. 21/2
   Approx. noise: colored noise in general domain


Wavelet multiresolution may provide a solution.

Initial attempt:      d = 1,   Haar basis
φ = φ0,0 = χ[0,1] ,      ψ(x) = φ(2x) − φ(2x − 1),        ψj,k (x) = 2j/2 ψ(2j x − k)
              2j −1, ∞
{φ0,0 , ψj,k }k=0, j=0     is an ON basis in H = L2 ((0, 1)).

                                  ∞ 2j −1
           1/2                                1/2
W (t) = γ0,0 β0,0 (t)φ0,0 (x) +             γj,k βj,k (t)ψj,k (x)
                                  j=0 k=0


                                    J 2j −1
            1/2                                1/2
W J (t) = γ0,0 β0,0 (t)φ0,0 (x) +             γj,k βj,k (t)ψj,k (x)
                                    j=0 k=0




                                                                             EFEF, Marseille, 2007 – p. 22/2
   Wavelet multiresolution

                                    J 2j −1
            1/2                                1/2
W J (t) = γ0,0 β0,0 (t)φ0,0 (x) +             γj,k βj,k (t)ψj,k (x)
                                    j=0 k=0


Framework of Yan is directly applicable. For example:

When γj,k = 1 we have white noise and the same piecewise constant
approximation as Allen, Novosel and Zhang. The rate is O(h1/2 ). (Li Bin)

For general domains D with d > 1 we may have to give up orthogonality.




                                                                      EFEF, Marseille, 2007 – p. 23/2
   Biorthogonal expansion


Frame: a countable set {φj }j∈J such that

      2
a f       ≤         |(f, φj )|2 ≤ b f   2

              j∈J


Then there is a dual frame {φ∗ }j∈J such that
                             j


f=          (f, φ∗ )φj
                 j
      j∈J


We make a formal expansion:

                                             1/2
W (t) =           (W (t), φ∗ )φj =
                           j                γj βj (t)φj
            j∈J                      j∈J




                                                          EFEF, Marseille, 2007 – p. 24/2
   Biorthogonal expansion


Let J ⊂ J be a finite subset. We investigate conditions under which we
can prove

  J
 WA (t) − WA (t) L2 (Ω,H) → 0, as J → J
  J         J
 WAh (t) − WA (t) L2 (Ω,H) → 0, as h → 0

together with error estimates, for the stochastic convolutions

                 t                                           t
                                                 1/2
WA (t) =             E(t − s) dW (s) =          γj               E(t − s)φj dβj (t)
             0                            j∈J            0
                 t                                               t
 J                                  J             1/2
WA (t)   =           E(t − s) dW (s) =           γj                  E(t − s)φj dβj (t)
             0                             j∈J               0
                     t                                                     t
 J                                                        1/2
WAh (t) =                Eh (t − s)Ph dW J (s) =         γj                    E(t − s)Ph φj dβj (t)
                 0                                 j∈J                 0




                                                                                                 EFEF, Marseille, 2007 – p. 25/2
   Error estimates (colored noise)


For continuous piecewise linear finite elements:

Theorem 4. Let −A be a self-adjoint nonpositive operator with compact inverse
                                                                                             β−1
generating an analytic semigroup on H . If T r(Q)    < ∞ and {φj }j∈J ⊂ D(A                   2      )
(β ≤ 2) is a frame for H , then


      E( WA (t) − WAh (t) 2 ) ≤ C
                   J                                                       ˜ ˜
                                                           (A−1 φk , φl )(Qφk , φl )
                                       k∈J \J l∈J \J
                                                          β−1            β−1
                                  + Ch  2β
                                                     (A    2    φk , A    2          ˜ ˜
                                                                               φl )(Qφk , φl )
                                             k,l∈J



To control the above sums more is needed!




                                                                                          EFEF, Marseille, 2007 – p. 26/2
   Support, cancellation, smoothness


Suppose J ↔ (j, k), j = 0, 1, . . . and k = 0, 1, ..., K(j, d). Usually
K(j, d) ∼ 2dj .
     Locality
                                ˜
     diam supp φj,k ∼ diam supp φj,k ∼ 2−j , j ≥ j0
     Cancellation property
     |(f, φj,k )| ≤ C2−j(s+d/2) f       W s,∞ (supp φj,k ) ,   s ≤ m, j ≥ j0 ,
                                                                   ˜
     and
          ˜
     |(f, φj,k )| ≤ C2−j(s+d/2) f                              s ≤ m, j ≥ j0 ,
                                        W s,∞ (supp φj,k ) ,
                                                    ˜

     Smoothness (inverse inequality)
      φj,k H s (Ω) ≤ C2sj φj,k L2 (Ω) , s ≤ γ,
     and
      ˜
      φj,k   H r (Ω)
                              ˜
                       ≤ C2rj φj,k   L2 (Ω) ,   r ≤ γ.
                                                    ˜



                                                                                 EFEF, Marseille, 2007 – p. 27/2
   Results in 1D (colored noise)

                                                              1
Let J ↔ {(j, k) : j = 0, 1, . . . , N } and (Qf )(x) :=           q(x, y)f (y) dy.
                                                          0

     Haar Basis (no smoothness)
     Theorem 5. Let    H := L2 [0, 1] and Au := −(au′ )′ + cu with
     a ≥ a0 > 0, c ≥ 0 smooth functions. Let q ∈ W 1,∞ ([0, 1]2 ) and {φj,k } be the
     Haar basis. Then, with 2−N := h,
                  J
     E( WA (t) − WAh (t) 2 ) ≤ Ch2 .
     Smooth wavelet basis
                   ∈ W 3,∞ ([0, 1]2 ) and assume locality, cancelation for
     Theorem 6. Let q
     m ≥ 1, m ≥ 2 and smoothness with γ = γ = 1. Then, with 2−N := h,
            ˜                              ˜
                  J
     E( WA (t) − WAh (t) 2 ) ≤ Ch4 | ln h|2 .

Good news: such wavelet basis exist!



                                                                               EFEF, Marseille, 2007 – p. 28/2
Future work


 Higher dimensions
 Add time stepping
 Implementation
                                       J
 Weak convergence: |E[f (WA (t)) − f (WAh (t))]|




                                                   EFEF, Marseille, 2007 – p. 29/2

				
DOCUMENT INFO