A new efficient operator splitting method for stochastic Maxwell equations
aa r X i v : . [ m a t h . NA ] F e b A NEW EFFICIENT OPERATOR SPLITTING METHOD FOR STOCHASTICMAXWELL EQUATIONS
CHUCHU CHEN, JIALIN HONG, AND LIHAI JIA
BSTRACT . This paper proposes and analyzes a new operator splitting method for stochasticMaxwell equations driven by additive noise, which not only decomposes the original multi-dimensional system into some local one-dimensional subsystems, but also separates the deter-ministic and stochastic parts. This method is numerically efficient, and preserves the symplec-ticity, the multi-symplecticity as well as the growth rate of the averaged energy. A detailed H -regularity analysis of stochastic Maxwell equations is obtained, which is a crucial prerequi-site of the error analysis. Under the regularity assumptions of the initial data and the noise, theconvergence order one in mean square sense of the operator splitting method is established.
1. I
NTRODUCTION
This paper proposes and analyzes a new operator splitting method for solving the followingthree-dimensional stochastic Maxwell equations with additive noise:d E ( t ) = curl H ( t ) d t + λ d W ( t ) , ( t , x ) ∈ ( , T ] × D , (1.1a)d H ( t ) = − curl E ( t ) d t + λ d W ( t ) , ( t , x ) ∈ ( , T ] × D , (1.1b) E ( , x ) = E ( x ) , H ( , x ) = H ( x ) , x ∈ D , (1.1c)where D = ( a − , a + ) × ( a − , a + ) × ( a − , a + ) ⊆ R is a cuboid, E and H represent the electric fieldand the magnetic field, respectively, and λ i ∈ R , i = , , describe the amplitude of the noise.Here, { W ( t ) } t ∈ [ , T ] is a standard Q -Wiener process with respect to a filtered probability space ( Ω , F , { F t } ≤ t ≤ T , P ) with Q being a symmetric, positive definite operator with finite trace on aHilbert space U .Stochastic Maxwell equations provide the foundation in stochastic electromagnetism and sta-tistical radiophysics etc., when considering thermodynamic fluctuations or wave propagation inrandom field (cf. [25]). It is of special importance to develop efficient numerical methods forsimulating stochastic Maxwell equations in large scale and long time computations. In the lastdecade, various numerical methods for stochastic Maxwell equations have been developed, an-alyzed, and tested in the literature. For instance, see [6, 26] for the finite element method anddiscontinuous Galerkin method, [1] for Wiener chaos expansion method, [7,8,12] for temporallysemi-discrete methods, [10, 21, 22] for stochastic multi-symplectic methods.Taking the multi-dimensional character of stochastic Maxwell equations into consideration,the operator splitting method (or the dimension splitting method) is one of the most promisingmethods in the reduction of computational costs. The splitting technique is introduced to the Key words and phrases. stochastic Maxwell equations, operator splitting method, mean square convergence order.The research of C. Chen and J. Hong were supported by the NNSFC (NOs. 11971470, 11871068, 12022118, and12031020), the research of L. Ji was supported by the NNSFC (NOs. 11601032, and 11971458). context of stochastic partial differential equations (SPDEs) in [2] for Zakai equation in filtering,and is extended to the general linear SPDEs in [18], in order to decompose the vector field intothe stochastic and deterministic parts. Afterwards, there are many papers using the splittingtechnique, and we refer to [13, 14, 24] for the case of stochastic Schr¨odinger equation, to [4, 5]for the case of stochastic Allen–Cahn equation, and to [3, 15] for the case of stochastic Navier–Stokes equations, etc.In this paper, we propose a new operator splitting method for stochastic Maxwell equations,which combines the superiority of the decomposability of the original multi-dimensional systeminto some local one-dimensional subsystems (see c.f. [11, 20, 23] for the case of deterministicMaxwell equations) and the feature of the separation of the deterministic and stochastic parts(see c.f. [2, 18]). To be specific, stochastic Maxwell equations (1.1) are split into the followingsubsystems d E [ ] ( t ) = λ d W ( t ) d H [ ] ( t ) = λ d W ( t ) , d E [ ] ( t ) = − ∂ x H [ ] ( t ) d t d H [ ] ( t ) = − ∂ x E [ ] ( t ) d t , d E [ ] ( t ) = ∂ x H [ ] ( t ) d t d H [ ] ( t ) = ∂ x E [ ] ( t ) d t , (1.2) d E [ ] ( t ) = λ d W ( t ) d H [ ] ( t ) = λ d W ( t ) , d E [ ] ( t ) = ∂ y H [ ] ( t ) d t d H [ ] ( t ) = ∂ y E [ ] ( t ) d t , d E [ ] ( t ) = − ∂ y H [ ] ( t ) d t d H [ ] ( t ) = − ∂ y E [ ] ( t ) d t , (1.3) d E [ ] ( t ) = λ d W ( t ) d H [ ] ( t ) = λ d W ( t ) , d E [ ] ( t ) = − ∂ z H [ ] ( t ) d t d H [ ] ( t ) = − ∂ z E [ ] ( t ) d t , d E [ ] ( t ) = ∂ z H [ ] ( t ) d t d H [ ] ( t ) = ∂ z E [ ] ( t ) d t . (1.4)Thus an approximation of (1.1) in the interval [ t n − , t n ] is split into three steps: solving thesubsystem (1.2) and taking its solution at time t n as the initial value at time t n − while solvingthe subsystem (1.3) again on [ t n − , t n ] , then repeating it to the subsystem (1.4). Namely, theapproximation { ( ˜ E ( t n ) , ˜ H ( t n )) } ≤ n ≤ N is defined recurrently by ( ˜ E ( t n ) , ˜ H ( t n )) = Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n ( ˜ E ( t n − ) , ˜ H ( t n − )) , n ≥ , (1.5)where ( ˜ E ( t ) , ˜ H ( t )) = ( E , H ) , and { Ψ [ j ] s , t : 0 ≤ s ≤ t } , j = , , λ i ( i = ,
2) are such thatthe components E [ ] and H [ ] in subsystem (1.2) only contain the stochastic part whose exactsolutions can be given explicitly, and other components E [ ] , E [ ] , H [ ] , H [ ] only contain the de-terministic part (i.e., two one-dimensional wave equations) whose spatial derivative is in onespace direction. The same holds for (1.3) and (1.4). This feature significantly reduces the com-plexity of the problem and improves the computational efficiency. Moreover, we show thateach subsystem is still a stochastic Hamiltonian PDE preserving both the infinite-dimensionalstochastic symplectic structure and stochastic multi-symplectic conservation law. If each sub-system is solved exactly, then the operator splitting method (1.5) preserves the averaged energyexactly.In order to study the mean square convergence of the operator splitting method (1.5), a rele-vant prerequisite is to provide the regularity theory of the solution of (1.1). Utilizing the analysisfor mixed inhomogeneous boundary value problems for the Laplacian on a cuboid ( [16, Lemma3.1]), the L ( Ω ; H ( D ) ) -regularity ( H -regularity for short) for the solution of (1.1) is proved PERATOR SPLITTING METHOD 3 with certain assumptions being made only on the initial data and the noise; see Proposition 2.2.Our main result states that for each T >
0, there exists a constant C independent of τ and N suchthat max ≤ n ≤ N h E (cid:16) k E ( t n ) − ˜ E ( t n ) k L ( D ) + k H ( t n ) − ˜ H ( t n ) k L ( D ) (cid:17)i ≤ C τ , (1.6)where τ is the uniform step size. Furthermore, it follows from a straightforward modification ofthe proof for this error estimate that the result still holds if we perturb the order of the subsystems,for example, the splitting method by using another order Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n and so on.The paper is organized as follows. In Section 2, we introduce the basic setting and give theregularity analysis of stochastic Maxwell equations. Section 3 presents the formulation of theoperator splitting method and analyzes the preservation of stochastic symplecticity and multi-symplecticity by the phase flow of each subsystem. In Section 4, we establish error estimates inmean square sense of the operator splitting method.2. P ROPERTIES OF STOCHASTIC M AXWELL EQUATIONS
This section presents the regularity analysis of stochastic Maxwell equations. The regularityin L ( Ω ; H k ( D ) ) with k = , C a generic positive constant which may be different from line to line, but independent of thestep size τ and the partition number N .2.1. Notation.
In this paper, we adopt the following conventions. We write Id for the identityoperator and v · w for the Euclidean inner product. E = ( E , E , E ) ⊤ and H = ( H , H , H ) ⊤ represent the electric field and the magnetic field, respectively, and λ i = ( λ i , λ i , λ i ) ⊤ ∈ R , i = , , describe the amplitude of the noise.We work on a finite time interval [ , T ] , T >
0, and on the cuboid D = ( a − , a + ) × ( a − , a + ) × ( a − , a + ) ⊆ R with Lipschitz boundary Γ : = ∂ D . Denote Γ ± j : = { x ∈ ¯ D | x j = a ± j } and Γ j : = Γ − j ∪ Γ + j for j ∈ { , , } . Then Γ = Γ ∪ Γ ∪ Γ . Let n be the unit outward normal to Γ .Let H s ( D ) , s ∈ R be the classical Sobolev space. The basic Hilbert space we work with is V : = L ( D ) with inner product h· , ·i V and norm k · k V . Let us define two spaces related to thecurl operator H ( curl , D ) = (cid:8) v ∈ L ( D ) : curl v ∈ L ( D ) (cid:9) , H ( curl , D ) = { v ∈ H ( curl , D ) : n × v | Γ = } endowed with the norm k v k H ( curl , D ) = (cid:16) k v k L ( D ) + k curl v k L ( D ) (cid:17) . We enforce a perfectly electric conducting (PEC) boundary condition on D , n × E = , on [ , T ] × Γ . (2.1)An important skew-adjoint operator which will be often used throughout the paper is theMaxwell operator M : = (cid:18) − curl 0 (cid:19) (2.2) CHUCHU CHEN, JIALIN HONG, AND LIHAI JI with domain D ( M ) = (cid:26)(cid:18) EH (cid:19) ∈ V : M (cid:18) EH (cid:19) = (cid:18) curl H − curl E (cid:19) ∈ V , n × E (cid:12)(cid:12)(cid:12) Γ = (cid:27) = H ( curl , D ) × H ( curl , D ) . (2.3)The corresponding graph norm is k v k D ( M ) : = ( k v k V + k Mv k V ) / . Recursively we could definethe domain D ( M k ) : = { u ∈ D ( M k − ) : M k − u ∈ D ( M ) } of the k -th power of the operator M , k ∈ N , with D ( M ) = V .Let U : = { f ∈ L ( D ) : f = Γ } . If we denote an orthonormal basis of U by { e k } k ∈ N , theKarhunen-Lo`eve expansion yields W ( t , x ) = ∑ k ∈ N Q e k ( x ) β k ( t ) , t ∈ [ , T ] , x ∈ D , with { β k } k ∈ N being a family of independent real-valued Brownian motions. Remark 2.1.
We point out that, under the assumption n · H | Γ = , one can verify that the PECboundary condition (2.1) is equivalent to n · H = , on [ , T ] × Γ . (2.4) In fact, since Q h | Γ = , for any h ∈ U , (2.1) together with (1.1b) leads to d ( n · H ( t )) = − n · curl E ( t ) d t + n · λ d W ( t )= [ ∇ · ( n × E ) − E · curl n ] d t + ∞ ∑ k = n · λ Q e k d β k ( t ) = , due to the identity ∇ · ( u × v ) = v · curl u − u · curl v for the second equation. The last equationholds because of (2.1) and the fact that n can be written as the gradient of a parametrizationof Γ and then curl ( ∇ f ) = for some scalar function f . Thus one has n · H = n · H = on [ , T ] × Γ . Well-posedness.
Let us give a short comment here on the well-posedness of stochasticMaxwell equations and the uniform boundedness of the solution in L ( Ω ; D ( M k )) -norm. Wenote that, with the operator M defined in (2.2), stochastic Maxwell equations (1.1) can be rewrit-ten as d Z ( t ) = M Z ( t ) d t + λ d W ( t ) , t ∈ ( , T ] ; Z ( ) = Z , (2.5)where Z = ( E ⊤ , H ⊤ ) ⊤ , Z = ( E ⊤ , H ⊤ ) ⊤ and λ = ( λ ⊤ , λ ⊤ ) ⊤ . Lemma 2.1. [17]
The Maxwell operator M defined in (2.2) with domain D ( M ) is closed, skew-adjoint on V , and generates a unitary C -semigroup S ( t ) = e tM on V for t ∈ R + . Similar to the well-posedness and uniform boundedness results from [8, Propositions 2.1 and2.2], we have the following theorem.
Theorem 2.1.
Let Z be an F -measurable V -valued random variable satisfying k Z k L ( Ω ; V ) < ∞ , and let Q be a symmetric, positive definite, trace class operator on U . Then stochasticMaxwell equations (2.5) have a unique mild solution Z ∈ L ( Ω ; C ([ , T ] ; V )) satisfying Z ( t ) = S ( t ) Z + Z t S ( t − s ) λ d W ( s ) , P -a . s ., (2.6) for each t ∈ [ , T ] . PERATOR SPLITTING METHOD 5
Furthermore, if for any k ∈ N , Q ∈ HS ( U , H k ( D )) and k Z k L ( Ω , D ( M k )) < ∞ , then the mildsolution Z ∈ L ( Ω ; C ([ , T ] ; D ( M k ))) and satisfies E h sup t ∈ [ , T ] k Z ( t ) k D ( M k ) i ≤ C (cid:16) + E k Z k D ( M k ) (cid:17) , (2.7) where the constant C may depend on T , | λ | and k Q k HS ( U , H k ( D )) . Here HS ( U , H k ( D )) denotesthe family of Hilbert-Schmidt operators from U to H k ( D ) . The averaged energy of stochastic Maxwell equations (1.1) evolutes linearly with the rate | λ | Tr ( Q ) (see [10, Theorem 2.1]): E k Z ( t ) k V = E k Z k V + t | λ | Tr ( Q ) . (2.8)From the proof of [10, Theorem 2.2], one gets the divergence evolution laws for the stochasticMaxwell equations. Let H ( div , D ) = { v ∈ ( L ( D )) : ∇ · v ∈ L ( D ) } . Lemma 2.2.
Assume that E , H ∈ H ( div , D ) , Q ∈ HS ( U , H ( D )) . The divergence of system (1.1) satisfies ∇ · E ( t ) = ∇ · E + λ · ( ∇ W ( t )) , ∇ · H ( t ) = ∇ · H + λ · ( ∇ W ( t )) . (2.9)2.3. H k -regularity results. In this part, we prove the regularities of the mild solution (2.6) in L ( Ω ; H k ( D ) ) -norm ( H k -regularity, for short) with k = ,
2, which are essential in the meansquare convergence analysis.
Proposition 2.1.
Assume that Q ∈ HS ( U , H ( D )) and Z ∈ L ( Ω ; H ( D ) ) . Then the solution (2.6) has the H -regularity, i.e., E k Z ( t ) k H ( D ) ≤ C (cid:16) + E k Z k H ( D ) (cid:17) , where the constant C may depend on T , | λ | and k Q k HS ( U , H ( D )) .Proof. The H -regularity is deduced by utilizing the fact that v ∈ H ( curl , D ) ∩ H ( div , D ) belongsto H ( D ) if v × n = v · n = Γ . Moreover, the H -norm of v is dominated by k v k H ( D ) ≤ C (cid:16) k v k L ( D ) + k curl v k L ( D ) + k ∇ · v k L ( D ) (cid:17) . (2.10)Meanwhile, from (2.9), it can be verified that there exists a constant C depending on T , | λ | , k Q k HS ( U , H ( D )) such that E k ∇ · E ( t ) k L ( D ) ≤ E k ∇ · E k L ( D ) + E k λ · ( ∇ W ( t )) k L ( D ) ≤ C ( + E k ∇ · E k L ( D ) ) . (2.11)Similarly, it holds E k ∇ · H ( t ) k L ( D ) ≤ C ( + E k ∇ · H k L ( D ) ) (2.12)with C depending on T , | λ | and k Q k HS ( U , H ( D )) . CHUCHU CHEN, JIALIN HONG, AND LIHAI JI
Since n × E = n · H = [ , T ] × Γ , from (2.10), we obtain E k Z ( t ) k H ( D ) ≤ C (cid:18) E k E ( t ) k L ( D ) + E k H ( t ) k L ( D ) + E k curl E ( t ) k L ( D ) + E k curl H ( t ) k L ( D ) + E k ∇ · E ( t ) k L ( D ) + E k ∇ · H ( t ) k L ( D ) (cid:19) ≤ C (cid:18) E k Z ( t ) k D ( M ) + E k ∇ · E ( t ) k L ( D ) + E k ∇ · H ( t ) k L ( D ) (cid:19) ≤ C ( + E k Z k H ( D ) ) , due to (2.7), (2.11) and (2.12). The proof is completed. (cid:3) In our error analysis we still need the H -regularity of the solution (2.6), whose proof re-lies on the following lemma about the mixed inhomogeneous boundary value problems for theLaplacian on a cuboid; see [16, Lemma 3.1] for the detailed proof. Lemma 2.3. [16, Lemma 3.1]
Let j ∈ { , , } and Γ ∗ = Γ \ Γ j . Take f ∈ L ( D ) and g ∈ H / ( Γ j ) : = (cid:0) L ( Γ j ) , H ( Γ j ) (cid:1) / , . If there is a unique function v ∈ H Γ ∗ ( D ) solving Z D v ϕ d x + Z D ∇ v · ∇ϕ d x = Z D f ϕ d x + Z Γ + j g ϕ d σ − Z Γ − j g ϕ d σ , (2.13) for all ϕ ∈ H Γ ∗ ( D ) , then the solution v belongs to H ( D ) ∩ H Γ ∗ ( D ) and satisfies v − ∆ v = f onD, ∂ n v = g on Γ j , and k v k H ( D ) ≤ C (cid:0) k f k L ( D ) + k g k H / ( Γ j ) (cid:1) with the constant C only dependingon D. Following the approach in [16, 20], we introduce a subspace H ( D ) of H ( D ) as H ( D ) = n f ∈ H ( D ) | tr Γ ′ f ∈ H / ( Γ ′ ) , for all faces Γ ′ of D o , and achieve higher regularity of the solution (2.6). Proposition 2.2.
Assume that Q ∈ HS (cid:0) U , H ( D ) ∩ H ( D ) (cid:1) , ∇ Q : U → H ( D ) , and Z ∈ L (cid:0) Ω ; H ( D ) (cid:1) , ∇ · E ∈ L ( Ω ; H ( D )) . Then the solution (2.6) has the H -regularity, i.e., E k Z ( t ) k H ( D ) ≤ C ( + E k Z k H ( D ) ) . For any t ∈ [ , T ] , the field ( E ( t ) , H ( t )) has the tracesE j = E k = , ∂ j E j = ∂ k E j = ∂ j E k = ∂ k E k = , on Γ i , H i = , ∂ j H i = ∂ k H i = , on Γ i , for all permutations ( i , j , k ) of ( , , ) .Proof. Step 1. The H -regularity follows from Proposition 2.1. Moreover, the asserted zero-order traces for E and H are a direct consequence of the boundary conditions (2.1) and (2.4),respectively. The first-order traces result from the established zero-order traces and the following H -regularity.From (2.9), one gets E k ∇ · E ( t ) k H ( D ) + E k ∇ · H ( t ) k H ( D ) ≤ C ( + k Z k H ( D ) ) , PERATOR SPLITTING METHOD 7 where the constant C may depend on T , | λ | and k Q k HS ( U , H ( D )) . Moreover from assumptionsof the initial data and noise, we have that for any t ∈ [ , T ] , ∇ · E ( t ) , ∇ · H ( t ) ∈ L ( Ω ; H ( D )) .We note that M Z = (cid:18) − curl ( curl E ) − curl ( curl H ) (cid:19) , and ∆ E = ∇ ( ∇ · E ) − curl ( curl E ) . It follows from (2.7) with k = ∆ E ∈ L ( Ω ; L ( D ) ) ,and E k ∆ E ( t ) k L ( D ) ≤ C ( E k ∇ · E ( t ) k H ( D ) + E k Z ( t ) k D ( M ) ) . The field H can be estimated similarly. Standard interior elliptic regularity then leads to E ( t ) , H ( t ) ∈ L ( Ω ; H ( D ) ) . Step 2.
We first consider the H -regularity of the first component E of E . Set Γ ∗ = Γ ∪ Γ .From Step 1 it implies that f : = ( I − ∆ ) E ∈ L ( D ) . By employing cut-off and mollificationin y , z directions, one can approximate a given ϕ ∈ H Γ ∗ ( D ) in H ( D ) by a smooth ψ havingsupport in [ a − , a + ] × [ a − + η , a + − η ] × [ a − + η , a + − η ] for some small η : = η ( ψ ) >
0. Foreach κ ∈ ( , η ) , define D κ : = ( a − + κ , a + − κ ) × ( a − + κ , a + − κ ) × ( a − + κ , a + − κ ) and denoteby Γ ± ( κ ) those open faces of D κ that contain the points of the form ( a ∓ ± κ , y , z ) . Integrationby parts and the support of ψ yield that Z D E ψ d x + Z D ∇ E · ∇ψ d x = lim κ → Z D κ E ψ d x + Z D κ ∇ E · ∇ψ d x = lim κ → (cid:20) Z D κ ψ ( I − ∆ ) E d x + Z ∂ D κ ψ∇ E · n d σ (cid:21) = Z D ψ f d x ± lim κ → Z Γ ± ( κ ) ψ∂ x E d σ = Z D ψ f d x ± lim κ → Z Γ ± ( κ ) ψ (cid:0) ∇ · E − ∂ y E − ∂ z E (cid:1) d σ = Z D ψ f d x + Z Γ + ψ∇ · E d σ − Z Γ − ψ∇ · E d σ , (2.14)where the last step is due to integration by parts once more and the fact that ψ vanishes on theboundary of Γ ± ( κ ) , as well as E , E on Γ . Lemma 2.3 and ∇ · E ∈ L ( Ω ; H ( D )) (see Step1 ) lead to E ∈ L ( Ω ; H ( D )) . In the same manner, one observes that E , E ∈ L ( Ω ; H ( D )) .Furthermore E k E j k H ( D ) ≤ C (cid:18) E k E j k L ( D ) + E k ∆ E j k L ( D ) + E k ∇ · E k H / ( Γ j ) (cid:19) ≤ C (cid:18) E k ∇ · E k H ( D ) + E k Z k D ( M ) + E k ∇ · E k H / ( Γ j ) (cid:19) . Step 3.
Now we consider H , and set Γ ∗ = Γ and e f : = ( I − ∆ ) H ∈ L ( D ) . Here we usethe homogeneous version of Lemma 2.3 (see also [20, Lemma 3.6]). As in Step 2 , we take asmooth function ψ ∈ H ( D ) having support in [ a − + η , a + − η ] × [ a − , a + ] × [ a − , a + ] for some CHUCHU CHEN, JIALIN HONG, AND LIHAI JI small η : = η ( ψ ) >
0. Choose κ ∈ ( , η ) so that ψ vanishes around Γ ± ( κ ) , then Z D H ψ d x + Z D ∇ H · ∇ψ d x = lim κ → (cid:20) Z D κ H ψ d x + Z D κ ∇ H · ∇ψ d x (cid:21) = lim κ → (cid:20) Z D κ ψ ( I − ∆ ) H d x + Z ∂ D κ ψ∇ H · n d σ (cid:21) = Z D ψ e f d x + lim κ → Z ∂ D κ h ψ∇ H · n − (cid:0) ( curl H ) × n (cid:1) · ( ψ , , ) i d σ = Z D ψ e f d x + lim κ → Z ∂ D κ ψ∂ x H · n d σ = Z D ψ e f d x ± lim κ → (cid:20) Z Γ ± ( κ ) ψ∂ x H d σ + Z Γ ± ( κ ) ψ∂ x H d σ (cid:21) = Z D ψ e f d x , (2.15)where in the last step, we use again integration by parts and the fact that H , H vanishes on Γ , Γ , respectively. Hence E k H k H ( D ) ≤ C (cid:0) E k H k L ( D ) + E k ∆ H k L ( D ) (cid:1) ≤ C (cid:0) E k ∇ · H k H ( D ) + E k Z k D ( M ) (cid:1) . Components H , H can be treated similarly.Combining Steps 1-3 and (2.7), we have E k Z ( t ) k H ( D ) ≤ C (cid:16) E k Z ( t ) k D ( M ) + E k ∇ · E ( t ) k H ( D ) + E k ∇ · H ( t ) k H ( D ) + ∑ Γ ′ ∈ Γ E k ∇ · E ( t ) k H / ( Γ ′ ) (cid:17) ≤ C ( + E k Z k H ( D ) ) . Thus the proof is finished. (cid:3)
3. O
PERATOR SPLITTING METHOD
In this section, we introduce the dimension splitting of the Maxwell operator. By employ-ing a corresponding decomposition of λ i ( i = , ) , we get three subsystems, which enjoy thesuperiority of the decomposability of the original multi-dimensional problem into some localone-dimensional subproblems and the feature of the separation of the deterministic and stochas-tic parts. Moreover, each subsystem is a stochastic Hamiltonian system.3.1. Operator splitting.
Let us now describe how to split three-dimensional stochastic Maxwellequations (1.1) into some one-dimensional subsystems. We split the ‘curl’ operator intocurl = curl x + curl y + curl z , (3.1)where curl x = − ∂ x ∂ x , curl y = ∂ y − ∂ y , curl z = − ∂ z ∂ z , PERATOR SPLITTING METHOD 9 are one-dimensional diffenrential operators (see e.g. [11,23] for the case of deterministic Maxwellequations). Define operators M x = (cid:20) x − curl x (cid:21) , M y = (cid:20) y − curl y (cid:21) , M z = (cid:20) z − curl z (cid:21) , (3.2)on V endowed with domains D ( M x ) = (cid:8) u ∈ V : M x u ∈ V , u = u = Γ (cid:9) , D ( M y ) = (cid:8) u ∈ V : M y u ∈ V , u = u = Γ (cid:9) , D ( M z ) = (cid:8) u ∈ V : M z u ∈ V , u = u = Γ (cid:9) , respectively. Note that Mu = M x u + M y u + M z u for u ∈ D ( M x ) ∩ D ( M y ) ∩ D ( M z ) ⊂ D ( M ) . Then, stochastic Maxwell equations (1.1) can bedecomposed into the following three subsystemsd Z [ ] ( t ) = M x Z [ ] ( t ) d t + λ [ ] d W ( t ) , (3.3)d Z [ ] ( t ) = M y Z [ ] ( t ) d t + λ [ ] d W ( t ) , (3.4)d Z [ ] ( t ) = M z Z [ ] ( t ) d t + λ [ ] d W ( t ) , (3.5)where λ [ ] = ( λ , , , λ , , ) ⊤ , λ [ ] = ( , λ , , , λ , ) ⊤ , λ [ ] = ( , , λ , , , λ ) ⊤ .3.2. Properties of subsystems.
The splitting of the original system (1.1) into the three subsys-tems (3.3)-(3.5) gives an approach to simulating stochastic Maxwell equations effectively andefficiently in large scale and long time computations.To show this clearly, we rewrite (3.3)-(3.5) into the following componentwise forms, respec-tively, d E [ ] ( t ) = λ d W ( t ) d H [ ] ( t ) = λ d W ( t ) , d E [ ] ( t ) = − ∂ x H [ ] ( t ) d t d H [ ] ( t ) = − ∂ x E [ ] ( t ) d t , d E [ ] ( t ) = ∂ x H [ ] ( t ) d t d H [ ] ( t ) = ∂ x E [ ] ( t ) d t , (3.6) d E [ ] ( t ) = λ d W ( t ) d H [ ] ( t ) = λ d W ( t ) , d E [ ] ( t ) = ∂ y H [ ] ( t ) d t d H [ ] ( t ) = ∂ y E [ ] ( t ) d t , d E [ ] ( t ) = − ∂ y H [ ] ( t ) d t d H [ ] ( t ) = − ∂ y E [ ] ( t ) d t , (3.7) d E [ ] ( t ) = λ d W ( t ) d H [ ] ( t ) = λ d W ( t ) , d E [ ] ( t ) = − ∂ z H [ ] ( t ) d t d H [ ] ( t ) = − ∂ z E [ ] ( t ) d t , d E [ ] ( t ) = ∂ z H [ ] ( t ) d t d H [ ] ( t ) = ∂ z E [ ] ( t ) d t . (3.8)Note that (3.6)-(3.8) possess similar characters and structures. Taking (3.6) for an example, forcomponents E [ ] and H [ ] , they only contain the stochastic part and admit the explicit formu-las of the exact solutions. And for other components E [ ] , H [ ] and E [ ] , H [ ] , they only con-tain the deterministic part with spatial derivative being in one space direction. They are twoone-dimensional wave equations, which can be solved independently. These characters of thesplitting will lead to a dramatic reduction of computational costs in solving stochastic Maxwellequations. Well-posedness.
The skew-adjointness of operators M x , M y and M z is shown in the fol-lowing lemma. Lemma 3.1.
Operators M α ( α = x , y , z ) with domain D ( M α ) are skew-adjoint on V .Proof. To prove the skew-adjointness of M α it is enough to show that M α is a skew-symmetricoperator and that Id ± M α has dense range. To show the skew-symmetry of M x , we take ψ =( u ⊤ , v ⊤ ) ⊤ and ˜ ψ = ( ˜ u ⊤ , ˜ v ⊤ ) ⊤ in D ( M x ) . The integration by parts formula then implies h M x ψ , ˜ ψ i V = Z D ( − ˜ u ∂ x v + ˜ u ∂ x v + ˜ v ∂ x u − ˜ v ∂ x u ) d x d y d z = Z D ( v ∂ x ˜ u − v ∂ x ˜ u − u ∂ x ˜ v + u ∂ x ˜ v ) d x d y d z = − h ψ , M x ˜ ψ i V , due to the boundary conditions in the definition of D ( M x ) . Thus M x is skew-symmetric andanalogously for M y , M z .To check the density of Id ± M x , we have to show thatran ( Id ± M x ) = V . (3.9)Because C ∞ ( D ) is dense in V , we infer that (3.9) is equivalent to show that for every f ∈ C ∞ ( D ) there is a g = ( E ⊤ , H ⊤ ) ⊤ ∈ D ( M x ) such that ( Id ± M x ) g = f , (3.10)or equivalently, E = f , E ∓ ∂ x H = f , E ± ∂ x H = f , H = f , H ∓ ∂ x E = f , H ± ∂ x E = f . It yields E − ∂ xx E = f ± ∂ x f = : e f ∈ L ( D ) , E − ∂ xx E = f ∓ ∂ x f = : e f ∈ L ( D ) . In order to solve these equations, we introduce the operator A x = ∂ xx with domain D ( A x ) = { u ∈ L ( D ) : ∂ x u , A x u ∈ L ( D ) , u = Γ } . The Lax-Milgram lemma thus yields the existence of E , E ∈ D ( A x ) such that E j − A x E j = e f j , j = ,
3. Defining H = f ± ∂ x E , H = f ∓ ∂ x E , we obtain a solution g = ( E ⊤ , H ⊤ ) ⊤ ∈ D ( M x ) of (3.10). Similarly, we can get the results for M y and M z . Thus the proof is finished. (cid:3) By applying Stone’s theorem, operators M x , M y and M z generate unitary C -semigroups S x ( t ) = e tM x , S y ( t ) = e tM y and S z ( t ) = e tM z on V for t ∈ R + , respectively. Therefore each subsys-tem has a unique mild solutions given by Z [ ] ( t ) = S x ( t ) Z [ ] + λ [ ] Z t S x ( t − s ) d W ( s ) , (3.11) PERATOR SPLITTING METHOD 11 Z [ ] ( t ) = S y ( t ) Z [ ] + λ [ ] Z t S y ( t − s ) d W ( s ) , (3.12) Z [ ] ( t ) = S z ( t ) Z [ ] + λ [ ] Z t S z ( t − s ) d W ( s ) , (3.13) P -a.s., respectively. Proposition 3.1.
Let Z [ j ] be F -measurable V -valued random variables satisfying k Z [ j ] k L ( Ω ; V ) < ∞ for j = , , , and let Q be a symmetric, positive definite, trace class operator on U . Then E k Z [ j ] ( t ) k V = E k Z [ j ] k V + t | λ [ j ] | Tr ( Q ) . Stochastic symplecticity and multi-symplecticity.
The phase flow of stochastic Maxwellequations (1.1) conserves the infinite-dimensional stochastic symplectic structure ( [8]) and sto-chastic multi-symplectic conservation law ( [10, 21]). This part investigates these geometricstructures for subsystems (3.3)-(3.5).To present the formulation of stochastic Hamiltonian system for each subsystem, we define H [ j ] = Z D (cid:16) | E [ j ] | + | H [ j ] | (cid:17) d x , H [ j ] = Z D (cid:16) λ j H [ j ] j − λ j E [ j ] j (cid:17) d x , where j = , ,
3. Then subsystems (3.3)-(3.5) are reformulated as " d E [ j ] ( t ) d H [ j ] ( t ) = " α − curl α δ H [ j ] δ E [ j ] δ H [ j ] δ H [ j ] d t + " Id − Id δ H [ j ] δ E [ j ] δ H [ j ] δ H [ j ] ◦ d W ( t ) , (3.14)for ( α , j ) ∈ { ( x , ) , ( y , ) , ( z , ) } . Lemma 3.2.
The phase flows of subsystems (3.3) , (3.4) and (3.5) preserve symplectic structures ω [ j ] ( t ) = ω [ j ] ( ) , P -a . s ., j = , , , respectively, where ω [ j ] ( t ) : = Z D d E [ j ] ( t , x ) ∧ curl α d H [ j ] ( t , x ) d x , are the integral over the space domain of the differential 2-forms with ‘d’ denoting the exteriorderivative for ( α , j ) ∈ { ( x , ) , ( y , ) , ( z , ) } , respectively.Proof. Taking the exterior differential on both sides of (3.14) and utilizing the skew-adjointnessof M α on D ( M α ) yield the conclusion. (cid:3) We remark that the canonical formulations of stochastic Hamiltonian system for subsystems(3.3)-(3.5) are " d E [ j ] ( t ) d H [ j ] ( t ) = " Id − Id δ H [ α ] δ E [ j ] δ H [ α ] δ H [ j ] d t + " Id − Id δ H [ j ] δ E [ j ] δ H [ j ] δ H [ j ] ◦ d W ( t ) , (3.15)with H [ α ] = Z D (cid:16) E [ j ] · curl α E [ j ] + H [ j ] · curl α H [ j ] (cid:17) d x , where ( α , j ) ∈ { ( x , ) , ( y , ) , ( z , ) } respectively. In this case, the symplectic structures ω [ j ] ( t ) : = Z D d E [ j ] ( t , x ) ∧ d H [ j ] ( t , x ) d x , j = , , u [ j ] =( H [ j ] , H [ j ] , H [ j ] , E [ j ] , E [ j ] , E [ j ] ) ⊤ , j = , ,
3, and denote skew-symmetric matrices F = " Id − Id , K j = " D j D j , with D = −
10 1 0 , D = − , D = − , and Hamiltonian S ( u [ j ] ) = λ j H [ j ] − λ j E [ j ] . Then, (3.3)-(3.5) can be written as F d u [ j ] + K j ∂ α u [ j ] d t = ∇ S ( u [ j ] ) ◦ d W , (3.16)with ( α , j ) ∈ { ( x , ) , ( y , ) , ( z , ) } , respectively. Similar to the proof of [10, Theorem 2.3], weget the following stochastic multi-symplecticity. Lemma 3.3.
Subsystems (3.3) , (3.4) and (3.5) possess stochastic multi-symplectic conservativelaws, respectively, i.e., for ( α , j ) ∈ { ( x , ) , ( y , ) , ( z , ) } , d ω [ j ] + ∂ α κ [ j ] d t = , P -a . s ., where ω [ j ] ( t , x ) = du [ j ] ∧ Fdu [ j ] , κ [ j ] ( t , x ) = du [ j ] ∧ K j du [ j ] , are differential 2-forms associated with skew-symmetric matrices F and K j .
4. E
RROR ESTIMATES OF THE OPERATOR SPLITTING METHOD
This section presents the error analysis of the operator splitting method. Moreover, the errorscaused by temporally semi-discretizing subsystems are also analyzed.Let t n = n τ ( n = , , · · · , N ) be a uniform partition of [ , T ] with τ = T / N , and let { Ψ [ j ] s , t :0 ≤ s ≤ t } ( j = , , ) be solution flows of subsystems (3.3)-(3.5), respectively. Let ˜ Z ( t n − ) =( ˜ E ( t n − ) , ˜ H ( t n − )) denote the approximated solution at time t n − , and let ˜ Z ( t ) = Z . The de-tailed algorithm for solving stochastic Maxwell equations (1.1) on the interval [ t n − , t n ] with n = , · · · , N is given by • Stage 1.
Take ˜ Z ( t n − ) as the initial data, solve (3.3) or (3.6), and obtain Ψ [ ] t n − , t n (cid:0) ˜ Z ( t n − ) (cid:1) ; • Stage 2.
Take Ψ [ ] t n − , t n ˜ Z ( t n − ) as the initial data, solve (3.4) or (3.7), and obtain Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n (cid:0) ˜ Z ( t n − ) (cid:1) ; PERATOR SPLITTING METHOD 13 • Stage 3.
Take Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n ˜ Z ( t n − ) as the initial data, solve (3.5) or (3.8), and obtain˜ Z ( t n ) = Ψ t n − , t n (cid:0) ˜ Z ( t n − ) (cid:1) with Ψ t n − , t n : = Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n .Thus, the composed solution at time t n = n τ reads˜ Z ( t n ) = Ψ t n − , t n ◦ Ψ t n − , t n − ◦ · · · ◦ Ψ t , t ◦ Ψ t , t ( Z ) . (4.1)Note that all the results in this paper still hold for other splittings by changing the order ofsubsystems (3.3)-(3.5), for example, by using another order Ψ t n − , t n : = Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n ◦ Ψ [ ] t n − , t n and so on.If Ψ [ j ] s , t denotes the exact solution flow of each subsystem, then based on Proposition 3.1, itcan be seen that the averaged energy is preserved exactly by the splitting method (4.1), i.e., E k ˜ Z ( t n ) k V = E k Z k V + t n | λ | Tr ( Q ) . Splitting error.
In this subsection we denote Ψ [ j ] s , t the exact solution flow of each subsystemand establish the mean square convergence analysis of the splitting algorithm for stochasticMaxwell equations (1.1). To present the formulation of Ψ t n , t n + , we note that the mild solutionof the subsystem (3.3) satisfies Ψ [ ] t n , t n + ( ˜ Z ( t n )) = S x ( τ ) ˜ Z ( t n ) + Z t n + t n S x ( t n + − r ) λ [ ] d W ( r ) . Similarly for (3.4) starting from Ψ [ ] t n , t n + ( ˜ Z ( t n )) at time t n , Ψ [ ] t n , t n + ◦ Ψ [ ] t n , t n + ( ˜ Z ( t n )) = S y ( τ ) Ψ [ ] t n , t n + ( ˜ Z ( t n )) + Z t n + t n S y ( t n + − r ) λ [ ] d W ( r )= S y ( τ ) S x ( τ ) ˜ Z ( t n ) + Z t n + t n S y ( τ ) S x ( t n + − r ) λ [ ] d W ( r )+ Z t n + t n S y ( t n + − r ) λ [ ] d W ( r ) , and then for (3.5) starting from Ψ [ ] t n , t n + ◦ Ψ [ ] t n , t n + ( ˜ Z ( t n )) at time t n ,˜ Z ( t n + ) = Ψ t n , t n + ( ˜ Z ( t n )) = Ψ [ ] t n , t n + ◦ Ψ [ ] t n , t n + ◦ Ψ [ ] t n , t n + ( ˜ Z ( t n ))= S z ( τ ) Ψ [ ] t n , t n + ◦ Ψ [ ] t n , t n + ( ˜ Z ( t n )) + Z t n + t n S z ( t n + − r ) λ [ ] d W ( r )= S z ( τ ) S y ( τ ) S x ( τ ) ˜ Z ( t n ) + Z t n + t n S z ( τ ) S y ( τ ) S x ( t n + − r ) λ [ ] d W ( r )+ Z t n + t n S z ( τ ) S y ( t n + − r ) λ [ ] d W ( r ) + Z t n + t n S z ( t n + − r ) λ [ ] d W ( r ) . (4.2)Denote the splitting error at time t n by e n : = Z ( t n ) − ˜ Z ( t n ) , n = , , · · · , N . We now give themean square convergence analysis of the error. Theorem 4.1.
Under the same assumption as in Proposition 2.2, for sufficiently small τ , thereexists a constant C independent of τ such that max ≤ n ≤ N E (cid:2) k e n k V (cid:3) ≤ C τ . (4.3) Proof.
Since λ = λ [ ] + λ [ ] + λ [ ] , the mild solution of (2.5) is Z ( t n + ) = S ( τ ) Z ( t n ) + Z t n + t n S ( t n + − r ) λ [ ] d W ( r )+ Z t n + t n S ( t n + − r ) λ [ ] d W ( r ) + Z t n + t n S ( t n + − r ) λ [ ] d W ( r ) . (4.4)Subtracting (4.2) from (4.4), we obtain e n + = S z ( τ ) S y ( τ ) S x ( τ ) e n + ( S ( τ ) − S z ( τ ) S y ( τ ) S x ( τ )) Z ( t n )+ Z t n + t n [ S ( t n + − r ) − S z ( τ ) S y ( τ ) S x ( t n + − r )] λ [ ] d W ( r )+ Z t n + t n [ S ( t n + − r ) − S z ( τ ) S y ( t n + − r )] λ [ ] d W ( r ) (4.5) + Z t n + t n [ S ( t n + − r ) − S z ( t n + − r )] λ [ ] d W ( r )= : S z ( τ ) S y ( τ ) S x ( τ ) e n + I + I + I + I . Next, we give the estimates of terms I j , j = , . . . , Estimate of term I Following [19], for the generator F of a C -semigroup and a real number τ ≥
0, we define thebounded operators α ( τ F ) = e τ F and α k ( τ F ) = Z e ( − ξ ) τ F ξ k − ( k − ) ! d ξ , for k ≥ . These operators satisfy the recurrence relation α k ( τ F ) = k ! Id + τ F α k + ( τ F ) , k ≥ . Notice that S ( τ ) Z ( t n ) is the solution of the problem dd t Z ( t ) = ( M z + M y + M x ) Z ( t ) at time t = t n . It can be rewritten by the variation-of-constants formula Z ( t n + ) = S ( τ ) Z ( t n )= S z ( τ ) Z ( t n ) + Z t n + t n S z ( s ) M y S ( t n + − s ) Z ( t n ) d s + Z t n + t n S z ( s ) M x S ( t n + − s ) Z ( t n ) d s = S z ( τ ) Z ( t n ) + τ S z ( τ ) M y Z ( t n ) − Z t n + t n sS z ( s ) (cid:0) M z M y − M y M (cid:1) S ( t n + − s ) Z ( t n ) d s + τ S z ( τ ) M x Z ( t n ) − Z t n + t n sS z ( s ) (cid:0) M z M x − M x M (cid:1) S ( t n + − s ) Z ( t n ) d s , (4.6) where the integration by parts formula is utilized in the last step. For the term S z ( τ ) Z ( t n ) , weuse the relations S y ( τ ) = α ( τ M y ) = Id + τ M y + τ α ( τ M y ) M y and S x ( τ ) = α ( τ M x ) = Id + τ M x + τ α ( τ M x ) M x to get S z ( τ ) Z ( t n ) = S z ( τ ) h S y ( τ ) − τ M y − τ α ( τ M y ) M y i Z ( t n ) PERATOR SPLITTING METHOD 15 = S z ( τ ) S y ( τ ) Z ( t n ) − τ S z ( τ ) M y Z ( t n ) − τ S z ( τ ) α ( τ M y ) M y Z ( t n )= S z ( τ ) S y ( τ ) h S x ( τ ) − τ M x − τ α ( τ M x ) M x i Z ( t n ) − τ S z ( τ ) M y Z ( t n ) − τ S z ( τ ) α ( τ M y ) M y Z ( t n )= S z ( τ ) S y ( τ ) S x ( τ ) Z ( t n ) − τ S z ( τ ) S y ( τ ) M x Z ( t n ) − τ S z ( τ ) M y Z ( t n ) − τ S z ( τ ) S y ( τ ) α ( τ M x ) M x Z ( t n ) − τ S z ( τ ) α ( τ M y ) M y Z ( t n ) . Substituting it into (4.6) and using the relation Id − S y ( τ ) = − τ M y α ( τ M y ) yield I = (cid:16) S ( τ ) − S z ( τ ) S y ( τ ) S x ( τ ) (cid:17) Z ( t n )= − τ S z ( τ ) α ( τ M y ) M y M x Z ( t n ) − τ S z ( τ ) S y ( τ ) α ( τ M x ) M x Z ( t n ) − τ S z ( τ ) α ( τ M y ) M y Z ( t n ) − Z t n + t n sS z ( s ) (cid:0) M z M y − M y M (cid:1) S ( t n + − s ) Z ( t n ) d s − Z t n + t n sS z ( s ) (cid:0) M z M x − M x M (cid:1) S ( t n + − s ) Z ( t n ) d s . It follows from the H -regularity of continuous solution Z ( t ) (see Proposition 2.2) that E k I k V ≤ C τ E k Z ( t n ) k H ( D ) ≤ C τ . (4.7)(ii) Estimates of terms I , I and I For the term I , Itˆo isometry yields E k I k V = Z t n + t n (cid:13)(cid:13)(cid:13)(cid:16) S ( t n + − s ) − S z ( τ ) S y ( τ ) S x ( t n + − s ) (cid:17) λ [ ] ◦ Q (cid:13)(cid:13)(cid:13) HS ( U , V ) d s (4.8) ≤ C Z t n + t n (cid:13)(cid:13)(cid:13)(cid:16) S z ( τ ) S y ( τ ) (cid:16) S ( t n + − s ) − S x ( t n + − s ) (cid:17) λ [ ] ◦ Q (cid:17) (cid:13)(cid:13)(cid:13) HS ( U , V ) d s + C Z t n + t n (cid:13)(cid:13)(cid:13)(cid:16) S z ( τ ) ( Id − S y ( τ )) S ( t n + − s ) (cid:17) λ [ ] ◦ Q (cid:13)(cid:13)(cid:13) HS ( U , V ) d s + C Z t n + t n (cid:13)(cid:13)(cid:13) ( Id − S z ( τ )) S ( t n + − s ) λ [ ] ◦ Q (cid:13)(cid:13)(cid:13) HS ( U , V ) d s ≤ C ( | λ [ ] | , k Q k HS ( U , H ( D )) ) τ , where we have used the fact that S ( t ) and S α ( t ) , α = x , y , z are unitary C -semigroups and theinequalities k S ( s ) − Id k L ( D ( M ) , V ) = s k α ( sM ) M k L ( D ( M ) , V ) ≤ C τ , ∀ s ≤ τ and k S α ( s ) − Id k L ( D ( M α ) , V ) = s k α ( sM α ) M α k L ( D ( M α ) , V ) ≤ C τ , ∀ s ≤ τ . Estimates of terms I and I are similar with E k I k V + E k I k V ≤ C ( | λ [ ] | , | λ [ ] | , k Q k HS ( U , H ( D )) ) τ . (4.9)Taking E k · k V on both sides of the error equation (4.5) and combining estimates (4.7)-(4.9),we obtain E k e n + k V ≤ E k S z ( τ ) S y ( τ ) S x ( τ ) e n k V + E h S z ( τ ) S y ( τ ) S x ( τ ) e n , I i V + C E (cid:0) k I k V + k I k V + k I k V + k I k V (cid:1) ≤ ( + C τ ) E k S z ( τ ) S y ( τ ) S x ( τ ) e n k V + C τ E k I k V + C E (cid:0) k I k V + k I k V + k I k V + k I k V (cid:1) ≤ ( + C τ ) E k e n k V + C τ , due to the unitarity of semigroups S α ( α = x , y , z ).Gronwall’s inequality completes the proof. (cid:3) Temporally semi-discretized error.
In this subsection, each subsystem in
Stages 1-3 isdiscretized temporally by using numerical methods, for example, the implicit Euler method,the midpoint method, the exponential Euler method, etc. Notice that the midpoint method andthe exponential Euler method preserve the stochastic symplectic structure. We refer readersto [7, 8, 12] for the analysis of these methods for stochastic Maxwell equations, and to [9] for theprobabilistic superiority of stochastic symplectic methods. Let ( E n − , H n − ) denote the approx-imated solution at time t n − , and let ( E , H ) = ( E , H ) and ∆ W n = W ( t n ) − W ( t n − ) . We takethe implicit Euler method for an example, and thus Stages 1-3 become: • Stage ′ . Use the implicit Euler method to temporally solve (3.6), ( E n − , ∗ = E n − + λ ∆ W n H n − , ∗ = H n − + λ ∆ W n , ( E n − , ∗ = E n − − τ∂ x H n − , ∗ H n − , ∗ = H n − − τ∂ x E n − , ∗ , ( E n − , ∗ = E n − + τ∂ x H n − , ∗ H n − , ∗ = H n − + τ∂ x E n − , ∗ . • Stage ′ . Use the implicit Euler method to temporally solve (3.7), ( E n − , ∗∗ = E n − , ∗ + λ ∆ W n H n − , ∗∗ = H n − , ∗ + λ ∆ W n , ( E n − , ∗∗ = E n − , ∗ + τ∂ y H n − , ∗∗ H n − , ∗∗ = H n − , ∗ + τ∂ y E n − , ∗∗ , ( E n − , ∗∗ = E n − , ∗ − τ∂ y H n − , ∗∗ H n − , ∗∗ = H n − , ∗ − τ∂ y E n − , ∗∗ . • Stage ′ . Use the implicit Euler method to temporally solve (3.8), ( E n = E n − , ∗∗ + λ ∆ W n H n = H n − , ∗∗ + λ ∆ W n , ( E n = E n − , ∗∗ − τ∂ z H n H n = H n − , ∗∗ − τ∂ z E n , PERATOR SPLITTING METHOD 17 ( E n = E n − , ∗∗ + τ∂ z H n H n = H n − , ∗∗ + τ∂ z E n . To avoid the confusion of notations, below we use Φ IEt , t + τ ( Z ) , Φ Mt , t + τ ( Z ) and Φ EEt , t + τ ( Z ) todenote the temporally semi-discrete solutions at t + τ starting from Z at time t of stochasticMaxwell equations (1.1) via the combination of the splitting and the implicit Euler method, themidpoint method, the exponential Euler method, respectively.Similar as (4.2), the one-step temporal approximations read: • Implicit Euler method Φ IEt , t + τ ( Z ) = S IE τ , z S IE τ , y S IE τ , x Z + S IE τ , z S IE τ , y S IE τ , x λ [ ] ∆ W + S IE τ , z S IE τ , y λ [ ] ∆ W + S IE τ , z λ [ ] ∆ W , where S IE τ , α = ( Id − τ M α ) − with α = x , y , z , and ∆ W = W ( t + τ ) − W ( t ) . • Midpoint method Φ Mt , t + τ ( Z ) = S M τ , z S M τ , y S M τ , x Z + S M τ , z S M τ , y T M τ , x λ [ ] ∆ W + S M τ , z T M τ , y λ [ ] ∆ W + T M τ , z λ [ ] ∆ W , where S M τ , α = ( Id − τ M α ) − ( Id + τ M α ) , T M τ , α = ( Id − τ M α ) − with α = x , y , z , and ∆ W = W ( t + τ ) − W ( t ) . • Exponential Euler method Φ EEt , t + τ ( Z ) = S z ( τ ) S y ( τ ) S x ( τ ) Z + S z ( τ ) S y ( τ ) S x ( τ ) λ [ ] ∆ W + S z ( τ ) S y ( τ ) λ [ ] ∆ W + S z ( τ ) λ [ ] ∆ W , where S α ( τ ) = e τ M α with α = x , y , z , and ∆ W = W ( t + τ ) − W ( t ) .For the analysis of the numerical error of the above three numerical methods, the followingtwo lemmas are introduced. Lemma 4.1.
It holds (cid:12)(cid:12) S IE τ , α (cid:12)(cid:12) L ( V , V ) ≤ , (cid:12)(cid:12) S M τ , α (cid:12)(cid:12) L ( V , V ) = , (cid:12)(cid:12) T M τ , α (cid:12)(cid:12) L ( V , V ) ≤ , (cid:12)(cid:12) T M τ , α − Id (cid:12)(cid:12) L ( V , V ) ≤ C τ . The proof of this lemma is standard, thus is omitted.
Lemma 4.2.
For the implicit Euler method and the midpoint method, there exists a constant Cindependent of τ such that (i) k S ( τ ) Z − S IE τ , z S IE τ , y S IE τ , x Z k V ≤ C τ k Z k H ( D ) . (ii) k S ( τ ) Z − S M τ , z S M τ , y S M τ , x Z k V ≤ C τ k Z k H ( D ) . Proof. (i)
It suffices to prove k S z ( τ ) S y ( τ ) S x ( τ ) Z − S IE τ , z S IE τ , y S IE τ , x Z k V ≤ C τ k Z k H ( D ) , based on the estimate of the term I in the proof of Theorem 4.1. For α = x , y , z , denote a = τ M α and ξ = ( Id − a ) − . Then ξ = Id + ξ a = Id + ( Id + ξ a ) a = Id + a + ξ a , i.e., S IE τ , α = Id + τ M α + τ ξ M α , where ξ is a bounded operator. Recall the relation S α ( τ ) = Id + τ M α + τ α ( τ M α ) M α , we have for any u ∈ D ( M α ) k ( S α ( τ ) − S IE τ , α ) u k V = k τ α ( τ M α ) M α u − τ ξ M α u k V ≤ C τ k u k D ( M α ) . Therefore, k S z ( τ ) S y ( τ ) S x ( τ ) Z − S IE τ , z S IE τ , y S IE τ , x Z k V ≤ k ( S z ( τ ) − S IE τ , z ) S y ( τ ) S x ( τ ) Z k V + k S IE τ , z ( S y ( τ ) − S IE τ , y ) S x ( τ ) Z k V + k S IE τ , z S IE τ , y ( S x ( τ ) − S IE τ , x ) Z k V ≤ C τ k Z k H ( D ) . (ii) Similarly, by denoting a = τ M α and ξ = ( Id − a ) − , we have S M τ , α =( Id + τ ξ M α )( Id + τ M α ) = Id + τ M α + τ ξ M α + τ ξ M α = Id + τ M α + τ ( Id + τ ξ M α ) M α + τ ξ M α = Id + τ M α + τ ξ M α . Thus k ( S α ( τ ) − S M τ , α ) u k V = k τ α ( τ M α ) M α u − τ ξ M α u k V ≤ C τ k u k D ( M α ) . The rest step is similar as in (i) . The proof is completed. (cid:3)
Denote the discrete composed solution at t n by Z n = Φ t n − , t n ◦ Φ t n − , t n − ◦ · · · ◦ Φ t , t ( Z ) with Φ t , t + τ ∈ { Φ IEt , t + τ , Φ Mt , t + τ , Φ EEt , t + τ } . Set ˆ e n : = Z ( t n ) − Z n . The mean square convergenceresult is stated below. Theorem 4.2.
Under the same assumption as in Proposition 2.2, for sufficiently small τ , thereexists a constant C independent of τ such that max ≤ n ≤ N (cid:16) E k ˆ e n k V (cid:17) / ≤ C τ . Proof.
We take the midpoint method as an example, since the proofs of the implicit Euler methodand the exponential Euler method are similar. Recall that Z n + = S M τ , z S M τ , y S M τ , x Z n + S M τ , z S M τ , y T M τ , x λ [ ] ∆ W n + + S M τ , z T M τ , y λ [ ] ∆ W n + + T M τ , z λ [ ] ∆ W n + , where ∆ W n + = W ( t n + ) − W ( t n ) . Subtracting it from (4.4) leads toˆ e n + = S M τ , z S M τ , y S M τ , x ˆ e n + ˆ I + ˆ I + ˆ I + ˆ I , PERATOR SPLITTING METHOD 19 with ˆ I = (cid:16) S ( τ ) − S M τ , z S M τ , y S M τ , x (cid:17) Z ( t n ) , ˆ I = Z t n + t n h S ( t n + − r ) − S M τ , z S M τ , y T M τ , x i λ [ ] d W ( r ) , ˆ I = Z t n + t n h S ( t n + − r ) − S M τ , z T M τ , y i λ [ ] d W ( r ) , ˆ I = Z t n + t n h S ( t n + − r ) − T M τ , z i λ [ ] d W ( r ) . Utilizing Lemma 4.2 (ii) , E k ˆ I k V ≤ C τ E k Z ( t n ) k H ( D ) ≤ C τ . For the term ˆ I , E (cid:13)(cid:13) ˆ I (cid:13)(cid:13) V = Z t n + t n (cid:13)(cid:13)(cid:13)(cid:16) S ( t n + − s ) − S M τ , z S M τ , y T M τ , x (cid:17) λ [ ] ◦ Q (cid:13)(cid:13)(cid:13) HS ( U , V ) d s ≤ C Z t n + t n (cid:13)(cid:13)(cid:13)(cid:16) S M τ , z S M τ , y (cid:16) S ( t n + − s ) − T M τ , x (cid:17) λ [ ] ◦ Q (cid:17) (cid:13)(cid:13)(cid:13) HS ( U , V ) d s + C Z t n + t n (cid:13)(cid:13)(cid:13)(cid:16) S M τ , z (cid:0) Id − S M τ , y (cid:1) S ( t n + − s ) (cid:17) λ [ ] ◦ Q (cid:13)(cid:13)(cid:13) HS ( U , V ) d s + C Z t n + t n (cid:13)(cid:13)(cid:13) (cid:0) Id − S M τ , z (cid:1) S ( t n + − s ) λ [ ] ◦ Q (cid:13)(cid:13)(cid:13) HS ( U , V ) d s ≤ C ( | λ [ ] | , k Q k HS ( U , H ( D )) ) τ , due to Lemma 4.1. Similarly, we can obtain E (cid:13)(cid:13) ˆ I (cid:13)(cid:13) V ≤ C ( | λ [ ] | , k Q k HS ( U , H ( D )) ) τ , E (cid:13)(cid:13) ˆ I (cid:13)(cid:13) V ≤ C ( | λ [ ] | , k Q k HS ( U , H ( D )) ) τ . Combing the above results, we get E k ˆ e n + k V ≤ E k S M τ , z S M τ , y S M τ , x ˆ e n k V + E h S M τ , z S M τ , y S M τ , x ˆ e n , ˆ I i V + C E (cid:0) k ˆ I k V + k ˆ I k V + k ˆ I k V + k ˆ I k V (cid:1) ≤ ( + C τ ) E k S M τ , z S M τ , y S M τ , x ˆ e n k V + C τ E k ˆ I k V + C E (cid:0) k ˆ I k V + k ˆ I k V + k ˆ I k V + k ˆ I k V (cid:1) ≤ ( + C τ ) E k e n k V + C τ . Thus the proof is completed by using Gronwall’s inequality. (cid:3) R EFERENCES [1] M. Badieirostami, A. Adibi, H. Zhou, and S. Chow. Wiener chaos expansion and simulation of electromagneticwave propagation excited by a spatially incoherent source.
Multiscale Model. Simul. , 8(2):591–604, 2009/10.[2] A. Bensoussan and R. Glowinski. Approximation of Zakai equation by the splitting up method. In
Stochasticsystems and optimization (Warsaw, 1988) , volume 136 of
Lect. Notes Control Inf. Sci. , pages 257–265. Springer,Berlin, 1989.[3] H. Bessaih and A. Millet. Strong L convergence of time numerical schemes for the stochastic two-dimensionalNavier-Stokes equations. IMA J. Numer. Anal. , 39(4):2135–2167, 2019.[4] C.-E. Br´ehier, J. Cui, and J. Hong. Strong convergence rates of semidiscrete splitting approximations for thestochastic Allen-Cahn equation.
IMA J. Numer. Anal. , 39(4):2096–2134, 2019.[5] C.-E. Br´ehier and L. Gouden`ege. Analysis of some splitting schemes for the stochastic Allen-Cahn equation.
Discrete Contin. Dyn. Syst. Ser. B , 24(8):4169–4190, 2019.[6] C. Chen. A symplectic discontinuous Galerkin full discretization for stochastic Maxwell equations. arXiv:2009.09880 , 2020.[7] C. Chen, J. Hong, and L. Ji. Mean-square convergence of a semidiscrete scheme for stochastic Maxwell equa-tions.
SIAM J. Numer. Anal. , 57(2):728–750, 2019.[8] C. Chen, J. Hong, and L. Ji. Runge-Kutta semidiscretizations for stochastic Maxwell equations with additivenoise.
SIAM J. Numer. Anal. , 57(2):702–727, 2019.[9] C. Chen, J. Hong, D. Jin, and L. Sun. Asymptotically-preserving large deviations principles by stochastic sym-plectic methods for a linear stochastic oscillator.
SIAM J. Numer. Anal. , 59(1):32–59, 2021.[10] C. Chen, J. Hong, and L. Zhang. Preservation of physical properties of stochastic Maxwell equations withadditive noise via stochastic multi-symplectic methods.
J. Comput. Phys. , 306:500–519, 2016.[11] W. Chen, X. Li, and D. Liang. Energy-conserved splitting finite-difference time-domain methods for Maxwell’sequations in three dimensions.
SIAM J. Numer. Anal. , 48(4):1530–1554, 2010.[12] D. Cohen, J. Cui, J. Hong, and L. Sun. Exponential integrators for stochastic Maxwell’s equations driven by Itˆonoise.
J. Comput. Phys. , 410:109382, 21, 2020.[13] J. Cui and J. Hong. Analysis of a splitting scheme for damped stochastic nonlinear Schr¨odinger equation withmultiplicative noise.
SIAM J. Numer. Anal. , 56(4):2045–2069, 2018.[14] J. Cui, J. Hong, Z. Liu, and W. Zhou. Strong convergence rate of splitting schemes for stochastic nonlinearSchr¨odinger equations.
J. Differential Equations , 266(9):5625–5663, 2019.[15] P. D¨orsek. Semigroup splitting and cubature approximations for the stochastic Navier-Stokes equations.
SIAMJ. Numer. Anal. , 50(2):729–746, 2012.[16] J. Eilinghoff and R. Schnaubelt. Error estimates in L ∼ schnaubelt/ media/adi- strong.pdf.[17] K. Engel and R. Nagel. One-Parameter Semigroups for Linear Evolution Equations . Springer, Berlin, 2000.[18] I. Gy¨ongy and N. Krylov. On the splitting-up method and stochastic partial differential equations.
Ann. Probab. ,31(2):564–591, 2003.[19] E. Hansen and A. Ostermann. Dimension splitting for evolution equations.
Numer. Math. , 108(4):557–570,2008.[20] M. Hochbruck, T. Jahnke, and R. Schnaubelt. Convergence of an ADI splitting for Maxwell’s equations.
Numer.Math. , 129(3):535–561, 2015.[21] J. Hong, L. Ji, and L. Zhang. A stochastic multi-symplectic scheme for stochastic Maxwell equations withadditive noise.
J. Comput. Phys. , 268:255–268, 2014.[22] J. Hong, L. Ji, L. Zhang, and J. Cai. An energy-conserving method for stochastic Maxwell equations withmultiplicative noise.
J. Comput. Phys. , 351:216–229, 2017.[23] L. Kong, J. Hong, and J. Zhang. Splitting multisymplectic integrators for Maxwell’s equations.
J. Comput.Phys. , 229(11):4259–4278, 2010.[24] J. Liu. Order of convergence of splitting schemes for both deterministic and stochastic nonlinear Schr¨odingerequations.
SIAM J. Numer. Anal. , 51(4):1911–1932, 2013.[25] S. M. Rytov, Yu. A. Kravtsov, and V. I. Tatarski˘ı.
Principles of statistical radiophysics. 3 . Springer-Verlag,Berlin, 1989. Elements of random fields, Translated from the second Russian edition by Alexander P. Repyev.[26] K. Zhang.
Numerical studies of some stochastic partial differential equations . ProQuest LLC, Ann Arbor, MI,2008. Thesis (Ph.D.)–The Chinese University of Hong Kong (Hong Kong).
PERATOR SPLITTING METHOD 21
LSEC, ICMSEC, A
CADEMY OF M ATHEMATICS AND S YSTEMS S CIENCE , C
HINESE A CADEMY OF S CI - ENCES , AND S CHOOL OF M ATHEMATICAL S CIENCES , U
NIVERSITY OF C HINESE A CADEMY OF S CIENCES ,B EIJING
HINA
Email address : [email protected] LSEC, ICMSEC, A
CADEMY OF M ATHEMATICS AND S YSTEMS S CIENCE , C
HINESE A CADEMY OF S CI - ENCES , AND S CHOOL OF M ATHEMATICAL S CIENCES , U
NIVERSITY OF C HINESE A CADEMY OF S CIENCES ,B EIJING
HINA
Email address : [email protected] I NSTITUTE OF A PPLIED P HYSICS AND C OMPUTATIONAL M ATHEMATICS , B
EIJING
EOPLE ’ S R E - PUBLIC OF C HINA . Email address ::