- Email: [email protected]
Domain decomposition methods for linear and semilinear elliptic stochastic partial differential equations
Available online at www.sciencedirect.com
Applied Mathematics and Computation 195 (2008) 630–640 www.elsevier.com/locate/amc
Domain decomposition methods for linear and semilinear elliptic stochastic partial differential equations Kai Zhang
a,1
, Ran Zhang
a,b,*,2
, Yunguang Yin c, Shi Yu
a
a
School of Mathematics, Jilin University, Changchun, Jilin 130023, PR China Department of Mathematics, Dalian University of Technology, Dalian 116024, PR China College of Computer Science and Technology, Jilin University, Changchun, Jilin 130023, PR China b
c
Abstract In this paper, we study several overlapping domain decomposition methods for the numerical solutions of some linear and semilinear elliptic stochastic partial differential equations discretized by the finite element methods. In particular, we show that the algorithms converge and the convergence rates are independent of the finite element mesh parameter, as well as the number of subdomains used in the domain decomposition. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Domain decomposition method; Stochastic partial differential equation; Finite element method
1. Introduction Many physical and engineering phenomena are modeled by partial differential equations which often contain some levels of uncertainty. The complexity for the solution of these so-called stochastic partial differential equations (SPDE) is that SPDEs are able to more fully capture the behavior of interesting phenomena; it also means that the corresponding numerical analysis of the model will require new tools to model the systems, produce the solutions, and analyze the information stored within the solutions. In the last decade, many researchers have studied the different SPDEs and various numerical methods and approximation schemes for SPDEs have also been developed, analyzed, and tested [2–6,8,9,12,14,16–20,24,27]. In [3,18], the analysis based on the traditional finite element method was successfully used on SPDEs with random coefficients, using the tensor product between the deterministic and random variable spaces. Numerical methods for SPDEs with white noise and Brownian motion added to the forcing term have also been developed, analyzed, and tested by several authors [3,8,12,13,24].
*
1 2
Corresponding author. Address: School of Mathematics, Jilin University, Changchun, Jilin 130023, PR China. E-mail address: [email protected] (R. Zhang). The first author is supported by the 985 program of Jilin University. The second author is supported by NSFC (10671082, 10626026 and 10471054), J0630104 of Jilin University.
0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.05.009
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
631
In this paper, we study several Schwarz type overlapping domain decomposition methods for the following boundary value problem of linear/semilinear stochastic elliptic partial differential equation driven by a color noise: ( DuðxÞ þ f ðuðxÞÞ ¼ gðxÞ þ W_ ðxÞ; x 2 X; ð1:1Þ uðxÞ ¼ 0; x 2 oX; where X = [0, 1]2, f is a continuous function on satisfying certain regularity conditions given in Section 2, and g is a given deterministic function and W_ denote a color noise. The additive noises may appear in various forms, ranging from the space-time white noise to colored noises generated by some infinite dimensional Brownian motion with a prescribed covariance operator [21,23]. The existence and uniqueness of the weak solution for (1.1) have been established in [12] by converting the problem into the Hammerstein integral equation using the Green’s function. The difficulty in the error analysis of numerical approximations for a SPDE is the lack of regularity of its solution. For instance, as shown in [3], the required regularity conditions are not satisfied for the problem (1.1) for the standard error estimates of finite element methods. For one-dimension case, if W_ corresponds to the Brownian white noise, Allen, Novosel, and Zhang have show that the regularity estimates are usually very weak, and lead to very low order error estimates [3]. On the other hand, if the noise is more regular, Du and Zhang have prove that it is possible to get higher order of error estimates for the numerical solution. For the semilinear elliptic SPDE in 2D case, Cao et al. also obtain some convergence results. To the best of our knowledge, there exists few work in the literature which studies domain decomposition method for the linear and semilinear elliptic stochastic partial differential equations with noise in 2D case, except Jin et al. have present Schwarz type domain decomposition methods for the numerical solution of stochastic elliptic problems, whose coefficients are assumed to be a random field with finite variance [16]. In the present paper, for the linear elliptic SPDE problems, we shall show that the additive Schwarz domain decomposition method yields a preconditioner for the preconditioned conjugate gradient method, and we also show that the approximate solutions from the multiplicative Schwarz domain decomposition method converge to the exact solution of the linear system geometrically; for the semilinear case, we locally linearize the nonlinear equation via a Newton-like algorithm and then to solve the resulting linearized problems at each nonlinear iteration by a domain decomposition method, and show under certain assumptions that the mesh parameters independent convergence can be obtained. The rest of this paper is organized as follows. We first summarize some notations and theorems about the linear and semilinear elliptic SPDEs in 2D case in Section 2. In Section 3, we discuss some convergence results for the additive and multiplicative Schwarz domain decomposition methods for linear elliptic SPDEs. Further, we present the convergence results about the Schwarz domain decomposition methods for the semilinear elliptic SPDEs case in Section 4. In Section 5, the numerical results are presented to support the theoretical analysis. For Section 6, some concluding remarks are given. 2. Linear and semilinear elliptic SPDEs In this section, we summarize some notations and theorems about the linear and semilinear elliptic SPDEs in 2D case [8,12]. For the simplicity of presentation, we focus on the following color noise: W_ ðxÞ ¼
1 X 1 X i¼1
rij gij wi ðx1 Þwj ðx2 Þ;
x ¼ ðx1 ; x2 Þ;
j¼1
where {gij} are random variables satisfying gij N(0, 1), and covðgij ; gi0 j0 Þ ¼ Eðgij gi0 j0 Þ :¼ qiji0 j0 with {rij} to be 1 chosen. Let frsij g1 i;j¼1 approaches to frij gi;j¼1 as s ! 1 in some appropriate sense, then an approximation of W_ ðxÞ is 1 X 1 X W_ s ðxÞ ¼ rsij gij wi ðx1 Þwj ðx2 Þ; x ¼ ðx1 ; x2 Þ: i¼1
j¼1
632
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
For the convenience, we introduce the following notation: rs ¼ ðrs11 ; rs12 ; . . . ; rs21 ; rs22 ; . . . ÞT ;
r ¼ ðr11 ; r12 ; . . . ; r21 ; r22 ; . . . ÞT
and Q be the infinite matrix with Q ¼ ðqiji0 j0 Þ. For an integer s, let Qs be the infinite matrix with Qs ¼ s ððiji0 j0 Þ qiji0 j0 Þ (note that Q0 = Q), and define 1 X s 2 rij ri0 j0 ðiji0 j0 Þ qij qi0 j0 ; krkQs ¼ hr; riQs : hr1 ; r2 iQs ¼ rT1 Qs r2 ¼ i;i0 ;j;j0 ¼1
In order pffiffiffito get the higher order of the convergence, we specialize our color noise to the choice of fwk ðxi Þ ¼ 2 sin kpxi ; i ¼ 1; 2; k ¼ 1; 2; . . .g. Replacing W_ by W_ s in (1.1), we have the following ‘‘simple’’ problem: ( Dus ðxÞ þ f ðus ðxÞÞ ¼ gðxÞ þ W_ s ðxÞ; x 2 X; ð2:1Þ us ðxÞ ¼ 0; x 2 oX: Its variational form is: Find us 2 H 10 ðXÞ such that aðus ; vÞ ¼ ðF s ; vÞ;
v 2 H 10 ðXÞ;
ð2:2Þ 2
where F ¼ g þ W_ s , (Æ, Æ) denotes the inner product on L (X), and s
aðu; vÞ ¼ ðru; rvÞ þ ðf ðuÞ; vÞ:
ð2:3Þ s
H 10 ðXÞ
2
\ H ðXÞ and then establish an estimate of the We want to show that (2.2) has a unique solution u in error u us. To these ends, we shall assume that f satisfies the following conditions: (A1) There is a constant a < c such that 2
ðf ðsÞ f ðtÞÞðs tÞ P ajs tj ;
8s; t 2 R;
where c is the positive constant in the Poincare’s inequality (see [1]). (A2) f is Lipschitz continuous, i.e., there is a positive constant Lq such that jf ðsÞ f ðtÞj 6 Lq js tj;
8s; t 2 R:
Then, we have following results [8,12]: Lemma 1. Under assumptions (A1) and (A2), the variational problem (2.2) has a unique solution in H 10 ðXÞ \ H 2 ðXÞ, and for any integer q P 0, we have " #1=2 1 X 1 X q 2 s s _ EðkW kH q Þ 6 C ðrij ðijÞ Þ ð2:4Þ i¼1
j¼1
provided that the right hand side is convergent, where E(Æ) is the expectation and C is a constant independent of h. Theorem 2. Let u and us be the solution of (1.1) and (2.1) respectively. If f satisfies (A1), then Eðku us k1 Þ 6 Ckr rs kQ1 ;
ð2:5Þ
where C is a positive constant independent of u and h. We can get the semi-norm estimate by consider the derivative of the Green function, and the rest proof of Theorem 2 is same to the Theorem 3.1 of [12]. Then our aim is convert to use numerical methods to solve the ‘‘simple’’ problem (2.1), for example the standard finite element method. Theorem 3 [8,12]. Assume Vh contain all piecewise polynomials of degree r in H 10 ðXÞ and us 2 H 10 ðXÞ \ H rþ1 ðXÞ. Let u be the solution of (1.1) and ush 2 V h be the finite element solution of (2.1) respectively. If f satisfies (A1) and (A2), then
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
0
"
Eðku ush k1 Þ 6 C @kr rs kQ1 þ hr kgkH r1 þ hr
1 X 1 X
#1=2 1 A; ðrs ðijÞr1 Þ2 ij
i¼1
633
ð2:6Þ
j¼1
where C is a positive constant independent of u and h. The finite element method reduces an elliptic equation (2.1) to a linear system of equations. It is well known that as the number of triangles in the underlying triangulation is increased, which is indispensable for increasing the accuracy of the approximate solution, the size and condition number of the linear system increases. The Schwarz domain decomposition methods will enable us to solve the problem (2.1) efficiently [10,15,25,26]. We shall study several overlapping domain decomposition methods for the numerical solutions of (2.1) in linear and semilinear cases discretized by the finite element methods in the following sections. 3. Domain decomposition methods for linear elliptic SPDEs In this section, we shall concentrate on Schwarz domain decomposition methods for linear elliptic SPDEs. Schwarz’s method has several remarkable properties such as its simplicity, its variational interpretation, and its geometric convergence for very different classes of differential equations [10,15,25,26]. Currently, two variations of Schwarz’s alternating method have been under investigation. One is the additive Schwarz’s domain decomposition method and the other is the multiplicative Schwarz’s domain decomposition method. We shall apply these two methods to solve linear elliptic stochastic partial differential equation, i.e., to the ‘‘simple’’ problem: ( L1 us Dus ðxÞ þ aus ðxÞ ¼ gðxÞ þ W_ s ðxÞ; x 2 X; ð3:1Þ us ðxÞ ¼ 0; x 2 oX: Letting a1 ð; Þ be the symmetric bilinear form induced by L1, we assume that the differential operator is elliptic in the following sense: 2
2
C 1 kus k1 6 a1 ðus ; us Þ 6 C 2 kus k1
8u 2 H 10 ðXÞ
ð3:2Þ 1
for some positive constants C1 and C2, where the norm is the H -norm. Let Vh be the usual triangular finite element subspace of H 10 ðXÞ (inner product a1 ð; Þ and norm k ka1 ¼ a1 ð; Þ1=2 ) consisting of continuous piecewise polynomials of degree r. Following the Dryja–Widlund construction of the overlapping decomposition of Vh [11], the triangulation of X is introduced as follows. The e k : k ¼ 1; . . . ; pg. This forms a mesh TH of size H. region is first divided into nonoverlapping substructures f X h e k is further partitioned into finer partition fs : sh 2 Th g. The fine partitioning gives the fine mesh Next, each X h T of mesh size h and X ¼ Xh [sh 2Th sh . Both TH and Th are assumed to be regular. Along with this partitioning, we assume that we are given another sequence of (overlapping) subdomains Xk (1 6 k 6 p), which is e k ð1 6 k 6 pÞ, such that oXk aligns with the h-level mesh and the distance between oXk and o X ek a large one of X k is bounded from below by dk > 0. We denote the minimum of all dk by d. We assume that oX does not cut through any element sh 2 Th . For the subdomains meeting the boundary o X we cut off the part of Xk which is outside X and Vk Vh is the usual finite element space over Xk ; k ¼ 1; . . . ; p. We further assume that the subdomains {Xk} satisfy a limited overlap property, i.e., each point of X is contained in at most q0 subdomains where q0 is independent of H and h. Let V0 P Vh be a triangular finite element subspace defined on the coarse p grid. It is clear that X ¼ [pk¼1 Xk and V h ¼ k¼0 V k . Base on the decomposition of Vh discussed above, the corresponding finite element problem is: Find ush 2 V h , such that a1 ðush ; vÞ ¼ ðf~ ; vÞ;
v 2 V h;
which can be rewritten as Au ¼ b:
ð3:3Þ
634
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
For each ush 2 V h we let uk denote the vector obtained from u by choosing the coefficients with respect to the basis for Vk. Let Rk : Vh ! Vk be the rectangular matrix defined by Rk ush ¼ uk
8ush 2 V h
and Ak be the restriction of A to the subspace Vk given by vT Ak u ¼ a1 ðush ; vsh Þ
8ush ; vsh 2 V k ; k ¼ 1; . . . ; p:
Thus Ak is symmetric and positive definite for k ¼ 1; . . . ; p. The following are the Schwarz’s domain decomposition methods. 3.1. Additive Schwarz’s domain decomposition method Solve (3.3) by the preconditioned conjugate gradient method with preconditioner M defined by M¼
p X
RTk A1 k Rk :
ð3:4Þ
k¼1
In other words, writing M = DTD with D being an upper triangular matrix solve for x DADT x ¼ Db
ð3:5Þ T
by the conjugate gradient method and recover u = D x. Let be fxðmÞ ; m ¼ 0; . . . ; g be the sequence generated by applying the conjugate gradient method to the preconditioned linear system DADT. Then it is known that there exists a positive constant C such that 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1m T B condðDAD Þ 1C ð0Þ kxðmÞ xk 6 C @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3:6Þ A kx xk; condðDADT Þ þ 1 where cond(A) is the condition number of matrix A. Lemma 4 [26]. Under the construction of the subspaces of Vh and decomposition of u, there exists a constant C3 and C4 independent of the number P p of subdomains such that: (1) for any u 2 Vh , there exist a partition uk 2 V k ; k ¼ 0; 1; . . . ; p with u ¼ pk¼0 uk such that p X
uTk Ak uk 6 C 3 uT Au;
ð3:7Þ
k¼0
(2) for any partition uk 2 V k ; k ¼ 0; 1; . . . ; p with u ¼ uT Au 6 C 4
p X
Pp
k¼0 uk
such that
uTk Ak uk :
ð3:8Þ
k¼0
We are now ready to state the following theorem about the additive Schwarz’s domain decomposition method for solving the ‘‘simple’’ problem (3.1): Theorem 5. The matrix M is symmetric and positive definite and, with M = DTD, the condition number satisfies condðDADT Þ 6 C 3 C 4 :
ð3:9Þ T
From the above theorem, we can easily see that the preconditioned linear system DAD has a much better condition number than A and the convergence rate is given by (3.6). 3.2. Multiplicative Schwarz’s domain decomposition method ð0Þ
Starting with any initial guess uh in H1(X) which satisfies the homogenous boundary condition, we conðmÞ struct a domain decomposition sequence fuh ; m ¼ 0; . . .g as follows: for m ¼ 0; 1; . . ., and for k ¼ 1; . . . ; p,
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
(
ðmþðkþ1Þ=pÞ
a1 ðuh
ðmþðkþ1Þ=pÞ
uh
; vÞ ¼
ðx; yÞ ¼
R
ðmþðkþ1Þ=pÞ
bv dx dy; ðuh
Xk ðmþk=pÞ uh ðx; yÞ;
ðmþk=pÞ
uh
ÞjoXk ¼ 0
8v 2 V k ;
ðx; yÞ 2 X n Xk :
635
ð3:10Þ
ðmÞ
We have following result that the sequence fuh ; m ¼ 0; . . .g from the multiplicative Schwarz’s domain decomposition method converges geometrically to the exact solution us of the ‘‘simple’’ problem. Theorem 6. Let l ¼ 1 ðmÞ 2
kus uh k1 6
1 1
C 3 ð1þðq0 C 4 Þ2 Þ2
, then
C2 m s ð0Þ 2 l ku uh k1 : C1
ð3:11Þ
Furthermore, we have ðmÞ
Eðku uh k1 Þ 6 Clm=2 þ Ckr rs kQ1 :
ð3:12Þ
The ideas behind the proofs of these results follow from an analysis of the proof of a result stated for Sobolev spaces for a second order elliptic problem, Lemma 1 and Theorem 2. To end this section, we point out that there exist many other domain decompositions methods for the elliptic problems can also work for the linear elliptic SPDEs, here we only pick up two kinds as examples. 4. Domain decomposition methods for semilinear elliptic SPDEs Base on the same decomposition of Vh discussed in Section 3, we introduce and analyze some algorithms for the finite element solution of the following semilinear elliptic problem with homogenous Dirichlet boundary condition: Dus ðxÞ þ f ðus ðxÞÞ ¼ gðxÞ þ W_ s ðxÞ; x 2 X; ð4:1Þ us ðxÞ ¼ 0; x 2 oX; where the nonlinear term f satisfies the assumption (A1) and (A2) given in Section 2. Its finite element problem is: Find ush 2 V h such that a2 ðush ; vÞ ¼ ðF s ; vÞ;
v 2 V h;
ð4:2Þ
where F s ¼ g þ W_ s , and a2 ðu; vÞ ¼ ðru; rvÞ þ ðf ðuÞ; vÞ:
ð4:3Þ
As a direct consequence of assumptions (A1) and (A2), we can prove the following lemma, which will be used extensively in the convergence analysis in the sequel. Lemma 7. Assume f satisfies the assumption (A1) and (A2), then the functional a2 ð; Þ : H 10 ðXÞ H 10 ðXÞ ! R has following properties: there exist constants C1,C2 > 0, 2
a2 ðu; u vÞ a2 ðv; u vÞ P C 1 ku vk1 ; ja2 ðu; wÞ a2 ðv; wÞj 6 C 2 ku vk1 kwk1 for any u; v; w 2 H 10 ðXÞ. Let V h ¼ Spanf/1 ; . . . ; /n g and the finite element solution ush ¼ ! n X gi ðu1 ; . . . ; un Þ ¼ a2 uj /j ; /i ; f i ¼ ðF s ; /i Þ;
Pn
i¼1 ui /i .
Define
j¼1 T
G ¼ ðg1 ; . . . ; gn Þ and f^ ¼ ðf1 ; . . . ; fn ÞT . The rest of this section is devoted to the solution of the following nonlinear algebraic equation: FðuÞ ¼ GðuÞ f^ ¼ 0 and we consider the Newton-like method.
ð4:4Þ
636
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
4.1. Poisson–Schwarz–Newton method Here, we discuss a simple algorithm that combines the Schwarz preconditioning technique with a Newton’s method. Similar to the linear case, for each subspace Vk, let us define an operator pk : Vh ! Vk by pk ðuÞ ¼ RTk A1 k Rk GðuÞ; where Rk and Ak are the restriction and local solver respectively as Section 3. To define the additive Schwarz method, let us define p = p0 + p1 + + pp. We known that the following nonlinear equation [7]: f FðuÞ pðuÞ g~ ¼ 0; ð4:5Þ Pp s where P g~ ¼ k¼0 g~k and g~k ¼ pk uh , has a unique solution and is equivalent to Eq. (4.4). Let us define 1 M ¼ ð pk¼0 RTk A1 k Rk Þ , we see that M is symmetric and positive definite and the norm generated by M(k Æ kM), which is equivalent to the energy norm k Æ k1 [11]. Algorithm 1. For a suitable parameter k, iterate for m ¼ 0; 1; . . . until convergence ðmþ1Þ
uh
ðmÞ
ðmÞ
¼ uh þ kvh ;
Fðuh Þ ¼ M 1 Fðuh Þ. where vh ¼ f ðmÞ
ðmÞ
ðmÞ
Lemma 8. There exists two constants b1 and b2 such that 2
ðpðu þ vÞ pðuÞ; vÞM P b1 kvkM ;
2
2
kpðu þ vÞ pðuÞkM 6 b2 kvkM
8u; v 2 V h :
Using Lemma 8 and Theorem 3, we can get following convergence of the Algorithm 1: Theorem 9. If we choose 0 < k < 2b1/b2, where b1 and b2 are both defined in Lemma 8, then Algorithm 1 converges in the sense that ðmÞ
ð0Þ
kush uh k1 6 Clm=2 kush uh k1 :
ð4:6Þ
Here l = 1 kb2(2b1/b2 k) < 1 and C are independent of the mesh parameters h and H. Furthermore, we have 0 " #1=2 1 1 X 1 X ðmÞ r1 2 A: ð4:7Þ ðrsij ðijÞ Þ Eðku uh k1 Þ 6 C @lm=2 þ kr rs kQ1 þ hr kgkH r1 þ hr i¼1
j¼1
4.2. A Newton–Krylov–Schwarz method In this section, we use classical Newton as the outer iterative method and a Schwarz preconditioned Krylov subspace method as the inner iterative method for solving (4.2). For each u 2 Vh, let us define M 1 Q ðuÞ ¼
p X
RTk Q1 k ðuÞRk ;
k¼0
as the additive Schwarz preconditioner corresponding to the Jacobi operator Q(u) of G(u). Here Q1 k ðuÞ is the inverse of Q(u) in the subspace Vk and Rk : Vh ! Vk is the restriction operator. To solve for the mth Newton correction, we use nm steps of a Schwarz-preconditioned Krylov subspace iterative method with initial guess v(0) = 0. Let Tm be the iteration operator, i.e., at the lth Krylov iteration, the error is given by vðlÞ v ¼ T m ðvð0Þ vÞ:
ð4:8Þ
For the simplicity of presentation, we consider following iteration as an example of Krylov iteration: ðmÞ
ðmÞ
ðmÞ
l T m ðuh Þ ¼ ðI km M 1 Q ðuh ÞQðuh ÞÞ ;
ð4:9Þ
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
637
where the km are relaxation parameters. We assume that the operator Tm is bounded, i.e., there exists a constant 0 < qm < 1, such that kT m k1 6 qm :
ð4:10Þ
The estimate (4.10) is satisfied for a number of Krylov space methods, such as GMRES [22]. ð0Þ
Algorithm 2. For any given uh 2 V h , iterate for m ¼ 0; 1; . . . until convergence ðmÞ
ðmþ1Þ
Qðuh ÞðI T m Þ1 ðuh
ðmÞ ðmÞ uh Þ ¼ Gðuh Þ þ f^ :
ð4:11Þ
In practice, a damping parameter can usually be used in each outer iteration to accelerate the convergence of the Newton method. The parameters can be selected by using either a line search or a trust region approach. Before giving the convergence results we present a few auxiliary lemmas. Let A ¼ fa2 ð/i ; /j Þg; i; j ¼ 1; . . . ; n. We assume that Q(u) satisfies the Lipschitz condition, i.e., kQðuÞ QðvÞkA1 6 aku vkA : ð0Þ
Theorem 10. There exist constants c1 and c2 both sufficiently small, such that if kush uh k1 6 c1 and qm 6 c2 for all m, then ðmÞ
ð0Þ
kush uh k1 6 lm kush uh k1 :
ð4:12Þ
Here 0 6 l < 1 is a constant independent of the mesh parameters. Furthermore, we have 0 " #1=2 1 1 X 1 X ðmÞ r1 2 A: Eðku uh k1 Þ 6 C @lm þ kr rs kQ1 þ hr kgkH r1 þ hr ðrsij ðijÞ Þ i¼1
ð4:13Þ
j¼1
Proof. Using the properties of a2 ð; Þ and Q, for any v, w 2 Vh, we have, a2 ðQ1 ðvÞw; Q1 ðvÞwÞ 6 a1 a2 ðA1 w; A1 wÞ; a2 ðA1 QðvÞw; A1 QðvÞwÞ 6 a2 a2 ðw; wÞ: By definition of A-norm and G, we can get kGðv þ zÞ GðvÞ QðvÞzkA1 6 C 3 kzk2A : Let em ¼ ush uðmÞ , using relation (4.9), (4.11) and the above two estimate, we have 2
1
ðmþ1Þ
kemþ1 k1 6 Cðkemþ1 k1 þ qm ð1 qm Þ kuh
ðmÞ
uh k1 Þ;
which implies estimate (4.6). We can easily obtain (4.13) by (4.6) and Theorem 3.
h
5. Numerical results for some model equations In this section, we present numerical examples to demonstrate our theoretical results in the previous section. We will consider both linear and nonlinear problems (c.f. [8,12]). The normal random variables for W_ s shall be simulated by using the random number generator gsl-rangaussian of femlab. Theoretically, the number of samples M should be chosen so that the error generated by the Monte Carlo method is in the same magnitude of the errors generated by the finite element approxiðmÞ mation. we shall evaluate Eðuh Þ by using the Monte Carlo method to examine ðmÞ
e1 ðhÞ ¼ kEðuÞ Eðuh Þk21 to see if we have used enough samples. We also employ the following two types of errors:
638
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640 2
ðmÞ 2
e2 ðhÞ ¼ jEðkuk1 Þ Eðkuh k1 Þj;
ðmÞ 2
ðmÞ 2
e3 ðhÞ ¼ jEðkuh k1 Þ Eðkuh=2 k1 Þj;
to check error estimates for linear and nonlinear problems, respectively. Example 1. We test the performance of the additive and multiplicative Schwarz preconditioned algorithm for solving the model problem (1.1) in linear case on a two-dimensional domain X = [0, 1]2. Let the exact solution be uðx; yÞ ¼ sinðpxÞ sinðpyÞ in the absence of the color noise. The unit square will be triangulated as type-I 1 1 1 partition (c.f [8]) with h ¼ 14 ; 18 ; 16 ; 32 ; 64 and nh = 1. Here, we take the color noise as follows: Let the random variables {gij} being iid, namely 0; i > 64; j > 64; 1 64 ; rij ¼ qiji0 j0 ¼ Eðgij gi0 j0 Þ ¼ dii0 djj0 ; rij ¼ 3=2 rij ; else: ðijÞ Then by definition, we have kr r64 kQ1 6 6412 6 h2 and we choose the piecewise linear polynomial as Vh 2 2 (r = 1). For linear case f(u) = 0, we have EðuÞ ¼ u; kEðuÞk1 ¼ 1þ2p . 4 We take the coarse grid to be the 32 triangles of size h ¼ 14, whose union is the closure of X. Overlapping subdomains are construct by adjointing just enough fine elements to the coarse elements so that the overlap index d ¼ 161 ; 18 ; 14 ; 12. For the additive and the multiplicative preconditioners, the CG method was used to solve the discrete problems on the coarse mesh and on the subdomains. The conditioner numbers of the preconditioned system were obtained in Table 5.1 by using the Lanczos technique [22]. The computational results are displayed in Tables 5.2 and 5.3. Example 2. In this example, we consider the nonlinear problem with f(u) = sin(u) on the unit disk and the exact solution is uðxÞ ¼ sin pðx2 þ y 2 Þ in the absence of the color noise. We use the Algorithms 1 and 2 (given in Section 4) for the above problem with the stop criterion 106(m is large enough). The iterative numbers of Algorithms 1 and 2 are given in Table 5.4, and the final results are displayed in Table 5.5.
Table 5.1 Condition number of additive and multiplicative Schwarz DDM d
1 16
1 8
1 4
1 2
Additive Multiplicative
14.28 14.12
12.45 12.32
9.81 9.47
7.41 7.13
Table 5.2 Test of linear problem on unit by additive Schwarz DDM h
e1
1/8 1/16 1/32 1/64
2.13e2 5.56e3 1.48e3 3.85e4
Rate
Eðkush k2 Þ
e2
Rate
1.91 1.93 1.94
5.14981 5.16856 5.18601 5.18825
1.92e2 5.64e3 1.68e3 5.21e4
1.76 1.74 1.69
Table 5.3 Test of linear problem on unit by multiplicative Schwarz DDM ðmÞ
h
e1
Rate
Eðkuh k2 Þ
e2
Rate
1/8 1/16 1/32 1/64
2.06e2 5.47e3 1.43e3 3.69e4
1.91 1.93 1.95
5.16214 5.17612 5.18766 5.19045
1.88e2 5.51e3 1.67e3 5.21e4
1.77 1.72 1.68
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
639
Table 5.4 Iterative number of Algorithms 1 and 2 d
1 16
1 8
1 4
1 2
Alg 1 Alg 2
28 22
27 21
27 20
26 20
Table 5.5 Test of nonlinear problem on unit by Algorithm 1(left) and Algorithm 2(right) method ðmÞ
Eðkuh k2 Þ
e3
Rate
Alg 1: h 1/8 1/16 1/32 1/64
3.55994 3.57392 3.58611 3.58835
4.56e2 1.33e2 3.75e3 1.02e3
1.78 1.82 1.87
Alg 2: h 1/8 1/16 1/32 1/64
3.57379 3.58467 3.58501 3.58593
4.31e2 1.21e3 3.28e3 8.54e4
1.83 1.88 1.94
6. Conclusion In this paper, we develop several overlapping domain decomposition methods for the numerical solutions of some stochastic linear and semilinear elliptic equations driven by color noises: the additive and multiplicative Schwarz domain decomposition methods for linear case; two Newton-like Schwarz methods for the semilinear case. Both our theoretical analysis and numerical experiments establish the rates of convergence of Schwarz domain decomposition methods for finite element approximate solutions. References [1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [2] E.J. Allen, S.J. Novosel, Z. Zhang, Difference approximation of some linear stochastic partial differential equations, Stochastics and Stochastics Reports 64 (1998) 117–142. [3] I. Babuska, R. Tempone, G. Zouraris, Galerkin finite element approximations of stochastic elliptic partial differential equations, SIAM J. Numer. Anal. 42 (2004) 800–825. [4] R.E. Bank, D.J. Rose, Analysis of a multilevel iterative method for nonlinear finite element equations, Math. Comp. 39 (1982) 453– 465. [5] S. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer, Berlin, Heidelberg, New York,, 1994. [6] R. Buckdahn, E. Pardoux, Monotonicity methods for white noise driven quasi-linear PDEs, Prog. Prob. 22 (1990) 219–233. [7] X.C. Cai, M. Dryja, Domain decomposition methods for monotone nonlinear elliptic problems, Amer. Math. Soc. Providence, RI 180 (1994) 21–27. [8] Y.Z. Cao, H.T. Yang, L. Yin, Finite element methods for semilinear elliptic stochastic partial differential equations, Preprint. [9] T.F. Chan, T. Mathew, Domain decomposition algorithm, Acta Numer. (1994) 61–143. [10] T.F. Chan, J. Zou, A convergence theory of multilevel additive Schwarz methods on unstructured meshes, Numer. Algor. 13 (1996) 365–398. [11] M. Dryja, O.B. Widlund, Towards a unified theory of domain decomposition algorithm for elliptic problems, in: Third International Symposium on Domian Decompostion Methods for Partial Differential Equations, SIAM, Philadelphia, 1990. [12] Q. Du, T.Y. Zhang, Numerical approximation of some linear stochastic partial differential equations driven by special additive noise, SIAM J. Numer. Anal. 4 (2002) 1421–1445. [13] I. Gyo¨ngy, Lattice approximations for stochastic quasi-linear parabolic partial differential equations driven by space-time white noise II, Poten. Anal. 11 (1999) 1–37. [14] E. Hausenblas, Error analysis for approximation of stochastic differential equations driven by Poisson random measures, SIAM J. Numer. Anal. 40 (2002) 87–113.
640
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
[15] R. Hiptmair, A. Toselli, Overlapping and multilevel Schwarz methods for vector valued elliptic problems in three dimensions, IMA Vol. Math. Appl. 120 (2000) 181–208. [16] C. Jin, X.C. Cai, C.M. Li, Parallel domain decomposition methods for stochastic elliptic equations, Preprint. [17] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Cambridge University Press, Cambridge, UK, 1994. [18] A. Keese, H. Matthies, Parallel computation of stochastic groundwater flow, NIC Symposium 2004, in: Proceedings, John von Neumann Institute for Computing, Julich, NIC series, vol. 20, 2003, pp. 399–408. [19] P. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlag, New York, 1992. [20] G.N. Milstein, E. Platen, H. Schurz, Balanced implicit methods for stochastic systems, SIAM J. Numer. Anal. 35 (1998) 1010–1019. [21] G.D. Prato, J. Zabczyk, Stochastic Equations in Infinite Dimensions, Cambridge University Press, Cambridge, UK, 1992. [22] Y. Saad, Iterative methods for sparse linear systems, Philadelphia, SIAM, 2003. [23] J. Sancho, J.G. Ojalvo, H. Guo, Non-equilibrium Ginzburg–Landau model driven by colored noise, Phys. D 113 (1998) 331–337. [24] T. Shardlow, Weak convergence of a numerical method for a stochastic heat equation, BIT Numer. Math. 43 (2003) 179–193. [25] B.F. Smith, P. Bjorstad, W. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, 1996. [26] A. Toselli, O. Widlund, Domain Decomposition Methods, Springer, Berlin, 2004. [27] J.B. Walsh, An introduction to stochastic partial differential equations, Lecture Notes in Mathematics 1180, Springer-Verlag, 1986, pp. 265–439.
Applied Mathematics and Computation 195 (2008) 630–640 www.elsevier.com/locate/amc
Domain decomposition methods for linear and semilinear elliptic stochastic partial differential equations Kai Zhang
a,1
, Ran Zhang
a,b,*,2
, Yunguang Yin c, Shi Yu
a
a
School of Mathematics, Jilin University, Changchun, Jilin 130023, PR China Department of Mathematics, Dalian University of Technology, Dalian 116024, PR China College of Computer Science and Technology, Jilin University, Changchun, Jilin 130023, PR China b
c
Abstract In this paper, we study several overlapping domain decomposition methods for the numerical solutions of some linear and semilinear elliptic stochastic partial differential equations discretized by the finite element methods. In particular, we show that the algorithms converge and the convergence rates are independent of the finite element mesh parameter, as well as the number of subdomains used in the domain decomposition. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Domain decomposition method; Stochastic partial differential equation; Finite element method
1. Introduction Many physical and engineering phenomena are modeled by partial differential equations which often contain some levels of uncertainty. The complexity for the solution of these so-called stochastic partial differential equations (SPDE) is that SPDEs are able to more fully capture the behavior of interesting phenomena; it also means that the corresponding numerical analysis of the model will require new tools to model the systems, produce the solutions, and analyze the information stored within the solutions. In the last decade, many researchers have studied the different SPDEs and various numerical methods and approximation schemes for SPDEs have also been developed, analyzed, and tested [2–6,8,9,12,14,16–20,24,27]. In [3,18], the analysis based on the traditional finite element method was successfully used on SPDEs with random coefficients, using the tensor product between the deterministic and random variable spaces. Numerical methods for SPDEs with white noise and Brownian motion added to the forcing term have also been developed, analyzed, and tested by several authors [3,8,12,13,24].
*
1 2
Corresponding author. Address: School of Mathematics, Jilin University, Changchun, Jilin 130023, PR China. E-mail address: [email protected] (R. Zhang). The first author is supported by the 985 program of Jilin University. The second author is supported by NSFC (10671082, 10626026 and 10471054), J0630104 of Jilin University.
0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.05.009
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
631
In this paper, we study several Schwarz type overlapping domain decomposition methods for the following boundary value problem of linear/semilinear stochastic elliptic partial differential equation driven by a color noise: ( DuðxÞ þ f ðuðxÞÞ ¼ gðxÞ þ W_ ðxÞ; x 2 X; ð1:1Þ uðxÞ ¼ 0; x 2 oX; where X = [0, 1]2, f is a continuous function on satisfying certain regularity conditions given in Section 2, and g is a given deterministic function and W_ denote a color noise. The additive noises may appear in various forms, ranging from the space-time white noise to colored noises generated by some infinite dimensional Brownian motion with a prescribed covariance operator [21,23]. The existence and uniqueness of the weak solution for (1.1) have been established in [12] by converting the problem into the Hammerstein integral equation using the Green’s function. The difficulty in the error analysis of numerical approximations for a SPDE is the lack of regularity of its solution. For instance, as shown in [3], the required regularity conditions are not satisfied for the problem (1.1) for the standard error estimates of finite element methods. For one-dimension case, if W_ corresponds to the Brownian white noise, Allen, Novosel, and Zhang have show that the regularity estimates are usually very weak, and lead to very low order error estimates [3]. On the other hand, if the noise is more regular, Du and Zhang have prove that it is possible to get higher order of error estimates for the numerical solution. For the semilinear elliptic SPDE in 2D case, Cao et al. also obtain some convergence results. To the best of our knowledge, there exists few work in the literature which studies domain decomposition method for the linear and semilinear elliptic stochastic partial differential equations with noise in 2D case, except Jin et al. have present Schwarz type domain decomposition methods for the numerical solution of stochastic elliptic problems, whose coefficients are assumed to be a random field with finite variance [16]. In the present paper, for the linear elliptic SPDE problems, we shall show that the additive Schwarz domain decomposition method yields a preconditioner for the preconditioned conjugate gradient method, and we also show that the approximate solutions from the multiplicative Schwarz domain decomposition method converge to the exact solution of the linear system geometrically; for the semilinear case, we locally linearize the nonlinear equation via a Newton-like algorithm and then to solve the resulting linearized problems at each nonlinear iteration by a domain decomposition method, and show under certain assumptions that the mesh parameters independent convergence can be obtained. The rest of this paper is organized as follows. We first summarize some notations and theorems about the linear and semilinear elliptic SPDEs in 2D case in Section 2. In Section 3, we discuss some convergence results for the additive and multiplicative Schwarz domain decomposition methods for linear elliptic SPDEs. Further, we present the convergence results about the Schwarz domain decomposition methods for the semilinear elliptic SPDEs case in Section 4. In Section 5, the numerical results are presented to support the theoretical analysis. For Section 6, some concluding remarks are given. 2. Linear and semilinear elliptic SPDEs In this section, we summarize some notations and theorems about the linear and semilinear elliptic SPDEs in 2D case [8,12]. For the simplicity of presentation, we focus on the following color noise: W_ ðxÞ ¼
1 X 1 X i¼1
rij gij wi ðx1 Þwj ðx2 Þ;
x ¼ ðx1 ; x2 Þ;
j¼1
where {gij} are random variables satisfying gij N(0, 1), and covðgij ; gi0 j0 Þ ¼ Eðgij gi0 j0 Þ :¼ qiji0 j0 with {rij} to be 1 chosen. Let frsij g1 i;j¼1 approaches to frij gi;j¼1 as s ! 1 in some appropriate sense, then an approximation of W_ ðxÞ is 1 X 1 X W_ s ðxÞ ¼ rsij gij wi ðx1 Þwj ðx2 Þ; x ¼ ðx1 ; x2 Þ: i¼1
j¼1
632
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
For the convenience, we introduce the following notation: rs ¼ ðrs11 ; rs12 ; . . . ; rs21 ; rs22 ; . . . ÞT ;
r ¼ ðr11 ; r12 ; . . . ; r21 ; r22 ; . . . ÞT
and Q be the infinite matrix with Q ¼ ðqiji0 j0 Þ. For an integer s, let Qs be the infinite matrix with Qs ¼ s ððiji0 j0 Þ qiji0 j0 Þ (note that Q0 = Q), and define 1 X s 2 rij ri0 j0 ðiji0 j0 Þ qij qi0 j0 ; krkQs ¼ hr; riQs : hr1 ; r2 iQs ¼ rT1 Qs r2 ¼ i;i0 ;j;j0 ¼1
In order pffiffiffito get the higher order of the convergence, we specialize our color noise to the choice of fwk ðxi Þ ¼ 2 sin kpxi ; i ¼ 1; 2; k ¼ 1; 2; . . .g. Replacing W_ by W_ s in (1.1), we have the following ‘‘simple’’ problem: ( Dus ðxÞ þ f ðus ðxÞÞ ¼ gðxÞ þ W_ s ðxÞ; x 2 X; ð2:1Þ us ðxÞ ¼ 0; x 2 oX: Its variational form is: Find us 2 H 10 ðXÞ such that aðus ; vÞ ¼ ðF s ; vÞ;
v 2 H 10 ðXÞ;
ð2:2Þ 2
where F ¼ g þ W_ s , (Æ, Æ) denotes the inner product on L (X), and s
aðu; vÞ ¼ ðru; rvÞ þ ðf ðuÞ; vÞ:
ð2:3Þ s
H 10 ðXÞ
2
\ H ðXÞ and then establish an estimate of the We want to show that (2.2) has a unique solution u in error u us. To these ends, we shall assume that f satisfies the following conditions: (A1) There is a constant a < c such that 2
ðf ðsÞ f ðtÞÞðs tÞ P ajs tj ;
8s; t 2 R;
where c is the positive constant in the Poincare’s inequality (see [1]). (A2) f is Lipschitz continuous, i.e., there is a positive constant Lq such that jf ðsÞ f ðtÞj 6 Lq js tj;
8s; t 2 R:
Then, we have following results [8,12]: Lemma 1. Under assumptions (A1) and (A2), the variational problem (2.2) has a unique solution in H 10 ðXÞ \ H 2 ðXÞ, and for any integer q P 0, we have " #1=2 1 X 1 X q 2 s s _ EðkW kH q Þ 6 C ðrij ðijÞ Þ ð2:4Þ i¼1
j¼1
provided that the right hand side is convergent, where E(Æ) is the expectation and C is a constant independent of h. Theorem 2. Let u and us be the solution of (1.1) and (2.1) respectively. If f satisfies (A1), then Eðku us k1 Þ 6 Ckr rs kQ1 ;
ð2:5Þ
where C is a positive constant independent of u and h. We can get the semi-norm estimate by consider the derivative of the Green function, and the rest proof of Theorem 2 is same to the Theorem 3.1 of [12]. Then our aim is convert to use numerical methods to solve the ‘‘simple’’ problem (2.1), for example the standard finite element method. Theorem 3 [8,12]. Assume Vh contain all piecewise polynomials of degree r in H 10 ðXÞ and us 2 H 10 ðXÞ \ H rþ1 ðXÞ. Let u be the solution of (1.1) and ush 2 V h be the finite element solution of (2.1) respectively. If f satisfies (A1) and (A2), then
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
0
"
Eðku ush k1 Þ 6 C @kr rs kQ1 þ hr kgkH r1 þ hr
1 X 1 X
#1=2 1 A; ðrs ðijÞr1 Þ2 ij
i¼1
633
ð2:6Þ
j¼1
where C is a positive constant independent of u and h. The finite element method reduces an elliptic equation (2.1) to a linear system of equations. It is well known that as the number of triangles in the underlying triangulation is increased, which is indispensable for increasing the accuracy of the approximate solution, the size and condition number of the linear system increases. The Schwarz domain decomposition methods will enable us to solve the problem (2.1) efficiently [10,15,25,26]. We shall study several overlapping domain decomposition methods for the numerical solutions of (2.1) in linear and semilinear cases discretized by the finite element methods in the following sections. 3. Domain decomposition methods for linear elliptic SPDEs In this section, we shall concentrate on Schwarz domain decomposition methods for linear elliptic SPDEs. Schwarz’s method has several remarkable properties such as its simplicity, its variational interpretation, and its geometric convergence for very different classes of differential equations [10,15,25,26]. Currently, two variations of Schwarz’s alternating method have been under investigation. One is the additive Schwarz’s domain decomposition method and the other is the multiplicative Schwarz’s domain decomposition method. We shall apply these two methods to solve linear elliptic stochastic partial differential equation, i.e., to the ‘‘simple’’ problem: ( L1 us Dus ðxÞ þ aus ðxÞ ¼ gðxÞ þ W_ s ðxÞ; x 2 X; ð3:1Þ us ðxÞ ¼ 0; x 2 oX: Letting a1 ð; Þ be the symmetric bilinear form induced by L1, we assume that the differential operator is elliptic in the following sense: 2
2
C 1 kus k1 6 a1 ðus ; us Þ 6 C 2 kus k1
8u 2 H 10 ðXÞ
ð3:2Þ 1
for some positive constants C1 and C2, where the norm is the H -norm. Let Vh be the usual triangular finite element subspace of H 10 ðXÞ (inner product a1 ð; Þ and norm k ka1 ¼ a1 ð; Þ1=2 ) consisting of continuous piecewise polynomials of degree r. Following the Dryja–Widlund construction of the overlapping decomposition of Vh [11], the triangulation of X is introduced as follows. The e k : k ¼ 1; . . . ; pg. This forms a mesh TH of size H. region is first divided into nonoverlapping substructures f X h e k is further partitioned into finer partition fs : sh 2 Th g. The fine partitioning gives the fine mesh Next, each X h T of mesh size h and X ¼ Xh [sh 2Th sh . Both TH and Th are assumed to be regular. Along with this partitioning, we assume that we are given another sequence of (overlapping) subdomains Xk (1 6 k 6 p), which is e k ð1 6 k 6 pÞ, such that oXk aligns with the h-level mesh and the distance between oXk and o X ek a large one of X k is bounded from below by dk > 0. We denote the minimum of all dk by d. We assume that oX does not cut through any element sh 2 Th . For the subdomains meeting the boundary o X we cut off the part of Xk which is outside X and Vk Vh is the usual finite element space over Xk ; k ¼ 1; . . . ; p. We further assume that the subdomains {Xk} satisfy a limited overlap property, i.e., each point of X is contained in at most q0 subdomains where q0 is independent of H and h. Let V0 P Vh be a triangular finite element subspace defined on the coarse p grid. It is clear that X ¼ [pk¼1 Xk and V h ¼ k¼0 V k . Base on the decomposition of Vh discussed above, the corresponding finite element problem is: Find ush 2 V h , such that a1 ðush ; vÞ ¼ ðf~ ; vÞ;
v 2 V h;
which can be rewritten as Au ¼ b:
ð3:3Þ
634
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
For each ush 2 V h we let uk denote the vector obtained from u by choosing the coefficients with respect to the basis for Vk. Let Rk : Vh ! Vk be the rectangular matrix defined by Rk ush ¼ uk
8ush 2 V h
and Ak be the restriction of A to the subspace Vk given by vT Ak u ¼ a1 ðush ; vsh Þ
8ush ; vsh 2 V k ; k ¼ 1; . . . ; p:
Thus Ak is symmetric and positive definite for k ¼ 1; . . . ; p. The following are the Schwarz’s domain decomposition methods. 3.1. Additive Schwarz’s domain decomposition method Solve (3.3) by the preconditioned conjugate gradient method with preconditioner M defined by M¼
p X
RTk A1 k Rk :
ð3:4Þ
k¼1
In other words, writing M = DTD with D being an upper triangular matrix solve for x DADT x ¼ Db
ð3:5Þ T
by the conjugate gradient method and recover u = D x. Let be fxðmÞ ; m ¼ 0; . . . ; g be the sequence generated by applying the conjugate gradient method to the preconditioned linear system DADT. Then it is known that there exists a positive constant C such that 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1m T B condðDAD Þ 1C ð0Þ kxðmÞ xk 6 C @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3:6Þ A kx xk; condðDADT Þ þ 1 where cond(A) is the condition number of matrix A. Lemma 4 [26]. Under the construction of the subspaces of Vh and decomposition of u, there exists a constant C3 and C4 independent of the number P p of subdomains such that: (1) for any u 2 Vh , there exist a partition uk 2 V k ; k ¼ 0; 1; . . . ; p with u ¼ pk¼0 uk such that p X
uTk Ak uk 6 C 3 uT Au;
ð3:7Þ
k¼0
(2) for any partition uk 2 V k ; k ¼ 0; 1; . . . ; p with u ¼ uT Au 6 C 4
p X
Pp
k¼0 uk
such that
uTk Ak uk :
ð3:8Þ
k¼0
We are now ready to state the following theorem about the additive Schwarz’s domain decomposition method for solving the ‘‘simple’’ problem (3.1): Theorem 5. The matrix M is symmetric and positive definite and, with M = DTD, the condition number satisfies condðDADT Þ 6 C 3 C 4 :
ð3:9Þ T
From the above theorem, we can easily see that the preconditioned linear system DAD has a much better condition number than A and the convergence rate is given by (3.6). 3.2. Multiplicative Schwarz’s domain decomposition method ð0Þ
Starting with any initial guess uh in H1(X) which satisfies the homogenous boundary condition, we conðmÞ struct a domain decomposition sequence fuh ; m ¼ 0; . . .g as follows: for m ¼ 0; 1; . . ., and for k ¼ 1; . . . ; p,
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
(
ðmþðkþ1Þ=pÞ
a1 ðuh
ðmþðkþ1Þ=pÞ
uh
; vÞ ¼
ðx; yÞ ¼
R
ðmþðkþ1Þ=pÞ
bv dx dy; ðuh
Xk ðmþk=pÞ uh ðx; yÞ;
ðmþk=pÞ
uh
ÞjoXk ¼ 0
8v 2 V k ;
ðx; yÞ 2 X n Xk :
635
ð3:10Þ
ðmÞ
We have following result that the sequence fuh ; m ¼ 0; . . .g from the multiplicative Schwarz’s domain decomposition method converges geometrically to the exact solution us of the ‘‘simple’’ problem. Theorem 6. Let l ¼ 1 ðmÞ 2
kus uh k1 6
1 1
C 3 ð1þðq0 C 4 Þ2 Þ2
, then
C2 m s ð0Þ 2 l ku uh k1 : C1
ð3:11Þ
Furthermore, we have ðmÞ
Eðku uh k1 Þ 6 Clm=2 þ Ckr rs kQ1 :
ð3:12Þ
The ideas behind the proofs of these results follow from an analysis of the proof of a result stated for Sobolev spaces for a second order elliptic problem, Lemma 1 and Theorem 2. To end this section, we point out that there exist many other domain decompositions methods for the elliptic problems can also work for the linear elliptic SPDEs, here we only pick up two kinds as examples. 4. Domain decomposition methods for semilinear elliptic SPDEs Base on the same decomposition of Vh discussed in Section 3, we introduce and analyze some algorithms for the finite element solution of the following semilinear elliptic problem with homogenous Dirichlet boundary condition: Dus ðxÞ þ f ðus ðxÞÞ ¼ gðxÞ þ W_ s ðxÞ; x 2 X; ð4:1Þ us ðxÞ ¼ 0; x 2 oX; where the nonlinear term f satisfies the assumption (A1) and (A2) given in Section 2. Its finite element problem is: Find ush 2 V h such that a2 ðush ; vÞ ¼ ðF s ; vÞ;
v 2 V h;
ð4:2Þ
where F s ¼ g þ W_ s , and a2 ðu; vÞ ¼ ðru; rvÞ þ ðf ðuÞ; vÞ:
ð4:3Þ
As a direct consequence of assumptions (A1) and (A2), we can prove the following lemma, which will be used extensively in the convergence analysis in the sequel. Lemma 7. Assume f satisfies the assumption (A1) and (A2), then the functional a2 ð; Þ : H 10 ðXÞ H 10 ðXÞ ! R has following properties: there exist constants C1,C2 > 0, 2
a2 ðu; u vÞ a2 ðv; u vÞ P C 1 ku vk1 ; ja2 ðu; wÞ a2 ðv; wÞj 6 C 2 ku vk1 kwk1 for any u; v; w 2 H 10 ðXÞ. Let V h ¼ Spanf/1 ; . . . ; /n g and the finite element solution ush ¼ ! n X gi ðu1 ; . . . ; un Þ ¼ a2 uj /j ; /i ; f i ¼ ðF s ; /i Þ;
Pn
i¼1 ui /i .
Define
j¼1 T
G ¼ ðg1 ; . . . ; gn Þ and f^ ¼ ðf1 ; . . . ; fn ÞT . The rest of this section is devoted to the solution of the following nonlinear algebraic equation: FðuÞ ¼ GðuÞ f^ ¼ 0 and we consider the Newton-like method.
ð4:4Þ
636
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
4.1. Poisson–Schwarz–Newton method Here, we discuss a simple algorithm that combines the Schwarz preconditioning technique with a Newton’s method. Similar to the linear case, for each subspace Vk, let us define an operator pk : Vh ! Vk by pk ðuÞ ¼ RTk A1 k Rk GðuÞ; where Rk and Ak are the restriction and local solver respectively as Section 3. To define the additive Schwarz method, let us define p = p0 + p1 + + pp. We known that the following nonlinear equation [7]: f FðuÞ pðuÞ g~ ¼ 0; ð4:5Þ Pp s where P g~ ¼ k¼0 g~k and g~k ¼ pk uh , has a unique solution and is equivalent to Eq. (4.4). Let us define 1 M ¼ ð pk¼0 RTk A1 k Rk Þ , we see that M is symmetric and positive definite and the norm generated by M(k Æ kM), which is equivalent to the energy norm k Æ k1 [11]. Algorithm 1. For a suitable parameter k, iterate for m ¼ 0; 1; . . . until convergence ðmþ1Þ
uh
ðmÞ
ðmÞ
¼ uh þ kvh ;
Fðuh Þ ¼ M 1 Fðuh Þ. where vh ¼ f ðmÞ
ðmÞ
ðmÞ
Lemma 8. There exists two constants b1 and b2 such that 2
ðpðu þ vÞ pðuÞ; vÞM P b1 kvkM ;
2
2
kpðu þ vÞ pðuÞkM 6 b2 kvkM
8u; v 2 V h :
Using Lemma 8 and Theorem 3, we can get following convergence of the Algorithm 1: Theorem 9. If we choose 0 < k < 2b1/b2, where b1 and b2 are both defined in Lemma 8, then Algorithm 1 converges in the sense that ðmÞ
ð0Þ
kush uh k1 6 Clm=2 kush uh k1 :
ð4:6Þ
Here l = 1 kb2(2b1/b2 k) < 1 and C are independent of the mesh parameters h and H. Furthermore, we have 0 " #1=2 1 1 X 1 X ðmÞ r1 2 A: ð4:7Þ ðrsij ðijÞ Þ Eðku uh k1 Þ 6 C @lm=2 þ kr rs kQ1 þ hr kgkH r1 þ hr i¼1
j¼1
4.2. A Newton–Krylov–Schwarz method In this section, we use classical Newton as the outer iterative method and a Schwarz preconditioned Krylov subspace method as the inner iterative method for solving (4.2). For each u 2 Vh, let us define M 1 Q ðuÞ ¼
p X
RTk Q1 k ðuÞRk ;
k¼0
as the additive Schwarz preconditioner corresponding to the Jacobi operator Q(u) of G(u). Here Q1 k ðuÞ is the inverse of Q(u) in the subspace Vk and Rk : Vh ! Vk is the restriction operator. To solve for the mth Newton correction, we use nm steps of a Schwarz-preconditioned Krylov subspace iterative method with initial guess v(0) = 0. Let Tm be the iteration operator, i.e., at the lth Krylov iteration, the error is given by vðlÞ v ¼ T m ðvð0Þ vÞ:
ð4:8Þ
For the simplicity of presentation, we consider following iteration as an example of Krylov iteration: ðmÞ
ðmÞ
ðmÞ
l T m ðuh Þ ¼ ðI km M 1 Q ðuh ÞQðuh ÞÞ ;
ð4:9Þ
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
637
where the km are relaxation parameters. We assume that the operator Tm is bounded, i.e., there exists a constant 0 < qm < 1, such that kT m k1 6 qm :
ð4:10Þ
The estimate (4.10) is satisfied for a number of Krylov space methods, such as GMRES [22]. ð0Þ
Algorithm 2. For any given uh 2 V h , iterate for m ¼ 0; 1; . . . until convergence ðmÞ
ðmþ1Þ
Qðuh ÞðI T m Þ1 ðuh
ðmÞ ðmÞ uh Þ ¼ Gðuh Þ þ f^ :
ð4:11Þ
In practice, a damping parameter can usually be used in each outer iteration to accelerate the convergence of the Newton method. The parameters can be selected by using either a line search or a trust region approach. Before giving the convergence results we present a few auxiliary lemmas. Let A ¼ fa2 ð/i ; /j Þg; i; j ¼ 1; . . . ; n. We assume that Q(u) satisfies the Lipschitz condition, i.e., kQðuÞ QðvÞkA1 6 aku vkA : ð0Þ
Theorem 10. There exist constants c1 and c2 both sufficiently small, such that if kush uh k1 6 c1 and qm 6 c2 for all m, then ðmÞ
ð0Þ
kush uh k1 6 lm kush uh k1 :
ð4:12Þ
Here 0 6 l < 1 is a constant independent of the mesh parameters. Furthermore, we have 0 " #1=2 1 1 X 1 X ðmÞ r1 2 A: Eðku uh k1 Þ 6 C @lm þ kr rs kQ1 þ hr kgkH r1 þ hr ðrsij ðijÞ Þ i¼1
ð4:13Þ
j¼1
Proof. Using the properties of a2 ð; Þ and Q, for any v, w 2 Vh, we have, a2 ðQ1 ðvÞw; Q1 ðvÞwÞ 6 a1 a2 ðA1 w; A1 wÞ; a2 ðA1 QðvÞw; A1 QðvÞwÞ 6 a2 a2 ðw; wÞ: By definition of A-norm and G, we can get kGðv þ zÞ GðvÞ QðvÞzkA1 6 C 3 kzk2A : Let em ¼ ush uðmÞ , using relation (4.9), (4.11) and the above two estimate, we have 2
1
ðmþ1Þ
kemþ1 k1 6 Cðkemþ1 k1 þ qm ð1 qm Þ kuh
ðmÞ
uh k1 Þ;
which implies estimate (4.6). We can easily obtain (4.13) by (4.6) and Theorem 3.
h
5. Numerical results for some model equations In this section, we present numerical examples to demonstrate our theoretical results in the previous section. We will consider both linear and nonlinear problems (c.f. [8,12]). The normal random variables for W_ s shall be simulated by using the random number generator gsl-rangaussian of femlab. Theoretically, the number of samples M should be chosen so that the error generated by the Monte Carlo method is in the same magnitude of the errors generated by the finite element approxiðmÞ mation. we shall evaluate Eðuh Þ by using the Monte Carlo method to examine ðmÞ
e1 ðhÞ ¼ kEðuÞ Eðuh Þk21 to see if we have used enough samples. We also employ the following two types of errors:
638
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640 2
ðmÞ 2
e2 ðhÞ ¼ jEðkuk1 Þ Eðkuh k1 Þj;
ðmÞ 2
ðmÞ 2
e3 ðhÞ ¼ jEðkuh k1 Þ Eðkuh=2 k1 Þj;
to check error estimates for linear and nonlinear problems, respectively. Example 1. We test the performance of the additive and multiplicative Schwarz preconditioned algorithm for solving the model problem (1.1) in linear case on a two-dimensional domain X = [0, 1]2. Let the exact solution be uðx; yÞ ¼ sinðpxÞ sinðpyÞ in the absence of the color noise. The unit square will be triangulated as type-I 1 1 1 partition (c.f [8]) with h ¼ 14 ; 18 ; 16 ; 32 ; 64 and nh = 1. Here, we take the color noise as follows: Let the random variables {gij} being iid, namely 0; i > 64; j > 64; 1 64 ; rij ¼ qiji0 j0 ¼ Eðgij gi0 j0 Þ ¼ dii0 djj0 ; rij ¼ 3=2 rij ; else: ðijÞ Then by definition, we have kr r64 kQ1 6 6412 6 h2 and we choose the piecewise linear polynomial as Vh 2 2 (r = 1). For linear case f(u) = 0, we have EðuÞ ¼ u; kEðuÞk1 ¼ 1þ2p . 4 We take the coarse grid to be the 32 triangles of size h ¼ 14, whose union is the closure of X. Overlapping subdomains are construct by adjointing just enough fine elements to the coarse elements so that the overlap index d ¼ 161 ; 18 ; 14 ; 12. For the additive and the multiplicative preconditioners, the CG method was used to solve the discrete problems on the coarse mesh and on the subdomains. The conditioner numbers of the preconditioned system were obtained in Table 5.1 by using the Lanczos technique [22]. The computational results are displayed in Tables 5.2 and 5.3. Example 2. In this example, we consider the nonlinear problem with f(u) = sin(u) on the unit disk and the exact solution is uðxÞ ¼ sin pðx2 þ y 2 Þ in the absence of the color noise. We use the Algorithms 1 and 2 (given in Section 4) for the above problem with the stop criterion 106(m is large enough). The iterative numbers of Algorithms 1 and 2 are given in Table 5.4, and the final results are displayed in Table 5.5.
Table 5.1 Condition number of additive and multiplicative Schwarz DDM d
1 16
1 8
1 4
1 2
Additive Multiplicative
14.28 14.12
12.45 12.32
9.81 9.47
7.41 7.13
Table 5.2 Test of linear problem on unit by additive Schwarz DDM h
e1
1/8 1/16 1/32 1/64
2.13e2 5.56e3 1.48e3 3.85e4
Rate
Eðkush k2 Þ
e2
Rate
1.91 1.93 1.94
5.14981 5.16856 5.18601 5.18825
1.92e2 5.64e3 1.68e3 5.21e4
1.76 1.74 1.69
Table 5.3 Test of linear problem on unit by multiplicative Schwarz DDM ðmÞ
h
e1
Rate
Eðkuh k2 Þ
e2
Rate
1/8 1/16 1/32 1/64
2.06e2 5.47e3 1.43e3 3.69e4
1.91 1.93 1.95
5.16214 5.17612 5.18766 5.19045
1.88e2 5.51e3 1.67e3 5.21e4
1.77 1.72 1.68
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
639
Table 5.4 Iterative number of Algorithms 1 and 2 d
1 16
1 8
1 4
1 2
Alg 1 Alg 2
28 22
27 21
27 20
26 20
Table 5.5 Test of nonlinear problem on unit by Algorithm 1(left) and Algorithm 2(right) method ðmÞ
Eðkuh k2 Þ
e3
Rate
Alg 1: h 1/8 1/16 1/32 1/64
3.55994 3.57392 3.58611 3.58835
4.56e2 1.33e2 3.75e3 1.02e3
1.78 1.82 1.87
Alg 2: h 1/8 1/16 1/32 1/64
3.57379 3.58467 3.58501 3.58593
4.31e2 1.21e3 3.28e3 8.54e4
1.83 1.88 1.94
6. Conclusion In this paper, we develop several overlapping domain decomposition methods for the numerical solutions of some stochastic linear and semilinear elliptic equations driven by color noises: the additive and multiplicative Schwarz domain decomposition methods for linear case; two Newton-like Schwarz methods for the semilinear case. Both our theoretical analysis and numerical experiments establish the rates of convergence of Schwarz domain decomposition methods for finite element approximate solutions. References [1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [2] E.J. Allen, S.J. Novosel, Z. Zhang, Difference approximation of some linear stochastic partial differential equations, Stochastics and Stochastics Reports 64 (1998) 117–142. [3] I. Babuska, R. Tempone, G. Zouraris, Galerkin finite element approximations of stochastic elliptic partial differential equations, SIAM J. Numer. Anal. 42 (2004) 800–825. [4] R.E. Bank, D.J. Rose, Analysis of a multilevel iterative method for nonlinear finite element equations, Math. Comp. 39 (1982) 453– 465. [5] S. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer, Berlin, Heidelberg, New York,, 1994. [6] R. Buckdahn, E. Pardoux, Monotonicity methods for white noise driven quasi-linear PDEs, Prog. Prob. 22 (1990) 219–233. [7] X.C. Cai, M. Dryja, Domain decomposition methods for monotone nonlinear elliptic problems, Amer. Math. Soc. Providence, RI 180 (1994) 21–27. [8] Y.Z. Cao, H.T. Yang, L. Yin, Finite element methods for semilinear elliptic stochastic partial differential equations, Preprint. [9] T.F. Chan, T. Mathew, Domain decomposition algorithm, Acta Numer. (1994) 61–143. [10] T.F. Chan, J. Zou, A convergence theory of multilevel additive Schwarz methods on unstructured meshes, Numer. Algor. 13 (1996) 365–398. [11] M. Dryja, O.B. Widlund, Towards a unified theory of domain decomposition algorithm for elliptic problems, in: Third International Symposium on Domian Decompostion Methods for Partial Differential Equations, SIAM, Philadelphia, 1990. [12] Q. Du, T.Y. Zhang, Numerical approximation of some linear stochastic partial differential equations driven by special additive noise, SIAM J. Numer. Anal. 4 (2002) 1421–1445. [13] I. Gyo¨ngy, Lattice approximations for stochastic quasi-linear parabolic partial differential equations driven by space-time white noise II, Poten. Anal. 11 (1999) 1–37. [14] E. Hausenblas, Error analysis for approximation of stochastic differential equations driven by Poisson random measures, SIAM J. Numer. Anal. 40 (2002) 87–113.
640
K. Zhang et al. / Applied Mathematics and Computation 195 (2008) 630–640
[15] R. Hiptmair, A. Toselli, Overlapping and multilevel Schwarz methods for vector valued elliptic problems in three dimensions, IMA Vol. Math. Appl. 120 (2000) 181–208. [16] C. Jin, X.C. Cai, C.M. Li, Parallel domain decomposition methods for stochastic elliptic equations, Preprint. [17] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Cambridge University Press, Cambridge, UK, 1994. [18] A. Keese, H. Matthies, Parallel computation of stochastic groundwater flow, NIC Symposium 2004, in: Proceedings, John von Neumann Institute for Computing, Julich, NIC series, vol. 20, 2003, pp. 399–408. [19] P. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlag, New York, 1992. [20] G.N. Milstein, E. Platen, H. Schurz, Balanced implicit methods for stochastic systems, SIAM J. Numer. Anal. 35 (1998) 1010–1019. [21] G.D. Prato, J. Zabczyk, Stochastic Equations in Infinite Dimensions, Cambridge University Press, Cambridge, UK, 1992. [22] Y. Saad, Iterative methods for sparse linear systems, Philadelphia, SIAM, 2003. [23] J. Sancho, J.G. Ojalvo, H. Guo, Non-equilibrium Ginzburg–Landau model driven by colored noise, Phys. D 113 (1998) 331–337. [24] T. Shardlow, Weak convergence of a numerical method for a stochastic heat equation, BIT Numer. Math. 43 (2003) 179–193. [25] B.F. Smith, P. Bjorstad, W. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, 1996. [26] A. Toselli, O. Widlund, Domain Decomposition Methods, Springer, Berlin, 2004. [27] J.B. Walsh, An introduction to stochastic partial differential equations, Lecture Notes in Mathematics 1180, Springer-Verlag, 1986, pp. 265–439.