On a Monte Carlo scheme for some linear stochastic partial differential equations

Takuya Nakagawa and Akihiro Tanaka

Abstract

The aim of this paper is to study the simulation of the expectation for the solution of linear stochastic partial differential equation driven by the space-time white noise with the bounded measurable coefficient and different boundary conditions. We first propose a Monte Carlo type method for the expectation of the solution of a linear stochastic partial differential equation and prove an upper bound for its weak rate error. In addition, we prove the central limit theorem for the proposed method in order to obtain confidence intervals for it. As an application, the Monte Carlo scheme applies to the stochastic heat equation with various boundary conditions, and we provide the result of numerical experiments which confirm the theoretical results in this paper.

MSC 2010: 35R60; 60H15; 60H35

A Appendix

A.1 Proof of the probability representation of (C1)

We consider the solution of SPDE (3.1) with condition (C1). It is well known that the fundamental solution of PDE (2.1) with L=t-a2x2 under condition (C1) is given by

(A.1)Gt(x,y)=gat(y-x)+gat(y+x)

(see e.g. [19]), and the fundamental solution is known to be the probability density function of the Wiener process on R+ reflected at 0. Therefore, we obtain the representation

m(t,x)=E[u0(|atZ1-x|)],

where Z1 is an RV with the standard normal distribution. Moreover, the representation of σ(t,x)2 follows from standard computations.

Lemma A.1

For any (t,x)R++×R+, there exist

• an RV 𝛽 with Beta(1,12) distribution,

• an RV Γ with Gamma(12,x2) distribution,

• an RV Z1 with standard normal distribution

such that (β,Γ,Z1) are mutually independent and

σ(t,x)2=E[tπah2(tβ,|at(1-β)2Z1-x|)+1axΓ-11[1/at,)(Γ)h2(t-1aΓ,12ΓZ1)].

Proof

Since the fundamental solution Gt(x,y) is given by (A.1), we have

σ(t,x)2=0tS(Gt-s(x,ξ)h(s,ξ))2dξds=0tS(ga(t-s)(ξ-x)+ga(t-s)(ξ+x))2h2(s,ξ)dξds=0tS(ga(t-s)2(ξ-x)+ga(t-s)2(ξ+x)+2ga(t-s)(ξ-x)ga(t-s)(ξ+x))h2(s,ξ)dξds=E1(t,x)+E2(t,x),

where E1(t,x) and E2(t,x) are defined by

E1(t,x)=0tS12πa(t-s)(ga(t-s)/2(ξ-x)+ga(t-s)/2(ξ+x))h2(s,ξ)dξds,
E2(t,x)=0tS1πa(t-s)exp(-x2a(t-s))ga(t-s)/2(ξ)h2(s,ξ)dξds.
Let β=βα,β be a Beta(α,β) distributed random variable which admits the probability density function

β(x;α,β)=1B(α,β)xα-1(1-x)β-1forx[0,1],

and βα,β and the Wiener process 𝑊 are mutually independent. Then we have

E1(t,x)=0tS12πat(1-s/t)(gat(1-s/t)/2(ξ-x)+gat(1-s/t)/2(ξ+x))h2(s,ξ)dξds=01St2πat(1-s)(gat(1-s)/2(ξ-x)+gat(1-s)/2(ξ+x))h2(ts,ξ)dξds=01t2πa(1-s)S(gat(1-s)/2(ξ-x)+gat(1-s)/2(ξ+x))h2(ts,ξ)dξds=01t2πa(1-s)E[h2(ts,|at(1-β)2Z1-x|)]ds=tB(1,1/2)2πa01β(s;1,1/2)E[h2(ts,|at(1-s)2Z1-x|)]ds=tπaE[h2(tβ,|at(1-β)2Z1-x|)].

Let Γ=Γ(α,β) be a Gamma(α,β) distributed random variable which admits the probability density function

Γ(x;α,β)=βαxα-1e-βxΓ(α)forx>0,α>0andβ>0,

and such that Γ(α,β) and the Wiener process 𝑊 are mutually independent. We have

E2(t,x)=0tS1πa(t-s)exp(-x2a(t-s))ga(t-s)/2(0,ξ)h2(s,ξ)dξds
=1atSs-3/2aπexp(-x2s)g12s(0,ξ)h2(t-1as,ξ)dξds
=1ax0Ss-11[1/(at),)(s)Γ(s;12,x2)g12s(0,ξ)h2(t-1as,ξ)dξds
=1ax0s-11[1/(at),)(s)Γ(s;12,x2)E[h2(t-1as,12sZ1)]ds
=1axE[Γ-11[1/(at),)(Γ)h2(t-1aΓ,12ΓZ1)],
which implies the statement. ∎

A.2 Proof of the probability representation of (C2)

We consider the solution of SPDE (3.1) with condition (C2). It is well known that the fundamental solution of PDE (2.1) with L=t-a2x2 under condition (C2) is given by

Gt(x,y)=gat(y-x)-gat(y+x).

Then we have the following representation:

m(t,x)=E[u0(atZ1+x)1R+(atZ1+x)-u0(atZ1-x)1R+(atZ1-x)],

where Z1 is an RV with standard normal distribution. Moreover, from an analogue of Lemma A.1, we obtain

σ(t,x)2=E[tπah2(tβ,|at(1-β)2Z1-x|)-1axΓ-11[1/at,)(Γ)h2(t-1aΓ,12ΓZ1)].

A.3 Proof of the probability representation of (C3)

We consider the solution of SPDE (3.1) under condition (C3). It is well known that the fundamental solution of PDE (2.1) with L=t-a2x2 with condition (C3) is given by

(A.2)Gt(x,y)=nZ{gat(y-x-2nl)+gat(-y-x-2nl)}

for any y[0,l] (see e.g. [19, p. 53]), and it is known that Gt(x,) has a Gaussian upper bound (see [20, p. 216]); thus Gt(x,) is square-integrable.

We now define the integer-valued functions K(x,) and Q(x,),

K(x,ν)=-x-ν2l+12,Q(x,ν)=-x-νl.

Then we obtain the following representation for m(t,x).

Lemma A.2

For any (t,x)R++×S, it holds that

m(t,x)=E[u0((-1)Q(x,atZ1)+1(x+2K(x,atZ1)l+atZ1))],

where Z1 is an RV with standard normal distribution.

The proof of Lemma A.2 follows from standard computation, and therefore we omit it.

In order to study the representation of σ2(t,x), we first consider the square of the fundamental solution (A.2). Note that gt is an even function and the linear map T(z1,z2):=(z1-z2,z1+z2) is in one-to-one correspondence for each z1,z2N; thus we obtain

(A.3)Gt(x,ξ)2=kZqZv,w=01gat(ξ+(-1)v(x+2kl))gat(ξ+(-1)w(x+2ql))=kZqZv,w=0112gat2(ξ+(-1)v(kl+x2)+(-1)w(ql+x2))gat2((x2+kl)-(-1)w+v(x2+ql))=kZqZv=0v=w112gat2(kl+x+(-1)vξ)gat2(ql)+kZqZv=0vw112gat2(kl+(-1)vξ)gat2(ql+x)=v=01I1(t,x,ξ,v)+v=01I2(t,x,ξ,v),

where I1 and I2 are defined by

(A.4)I1(t,x,ξ,v)=R212gat2(ly1+x+(-1)vξ)gat2(ly2)dy1dy2,
(A.5)I2(t,x,ξ,v)=R212gat2(ly1+(-1)vξ)gat2(ly2+x)dy1dy2.
The last equality in (A.3) follows since nZan=nZ-an1[n,n+1)(x)dx=-axdx for any sequence (an)nN satisfying |nZan|<.

By using the importance sampling technique and (A.3), we obtain the following representation for σ(t,x)2.

Lemma A.3

For any (t,x)R++×S, there exist RVs UiU(0,1) and ZiN(0,1) for each i=1,2 such that (U1,U2,Z1,Z2) are mutually independent and

σ(t,x)2=E[j=12{v=01tl2h2(t(1-U1),lU2)gatU12((-1)vlU2+δ1,jx)gatU12(δ2,jx)+v=01t2h2(t(1-U1),lU2)gatU12((-1)vlU2+δ1,jx)g1(φj(t,x,U1,Z2))g1(Z2)×1[0,1)c(φ~j(t,x,U1,Z2))+v=01t2h2(t(1-U1),lU2)g1(θj(v,t,x,U1,U2,Z1))g1(Z1)gatU12(δ2,jx)×1[0,1)c(θ~j(v,t,x,U1,U2,Z1))+v=01t2lh2(t(1-U1),lU2)g1(θj(v,t,x,U1,U2,Z1))g1(Z1)g1(φj(t,x,U1,Z2))g1(Z2)×1[0,1)c×[0,1)c(θ~j(v,t,x,U1,U2,Z1),φ~j(t,x,U1,Z2))}],

where δ1,j and δ2,j are Kronecker deltas, and θj,θ~j,φj and φ~j are defined by

θj(v,t,x,U1,U2,Z1):=2atU1(lθ~j(v,t,x,U1,U2,Z1)+(-1)vlU2+δ1,jx),
φj(t,x,U1,Z1):=2atU1(lφ~j(t,x,U1,Z2)+δ2,jx),
θ~j(v,t,x,U1,U2,Z1):=atU12Z1-δ1,jx-(-1)vlU2l,
φ~j(t,x,U1,Z2):=atU12Z2-δ2,jxl.

Proof

Using (A.3), we obtain

(A.6)σ(t,x)2=v=010l0tI1(t-s,x,ξ,v)h2(s,ξ)dsdξ+v=010l0tI2(t-s,x,ξ,v)h2(s,ξ)dsdξ.

We now consider the first term of (A.6). In order to derive a probability representation that satisfies (2.6) for p=2, we split the first term of (A.6) as follows:

v=010l0tI1(t-s,x,ξ,v)h2(s,ξ)dsdξ
=v=01t2l0101R2h2(t(1-s),lξ)gats2(ly1l+x+(-1)vlξ)gats2(ly2l)dy1dy2dsdξ
=t2l0101[0,l)×[0,l)h2(t(1-s),lξ)gats2(ly1l+x+(-1)vlξ)gats2(ly2l)dy1dy2dsdξ
+t2l0101[0,l)×[0,l)ch2(t(1-s),lξ)gats2(ly1l+x+(-1)vlξ)gats2(ly2l)dy1dy2dsdξ
+t2l0101[0,l)c×[0,l)h2(t(1-s),lξ)gats2(ly1l+x+(-1)vlξ)gats2(ly2l)dy1dy2dsdξ
+t2l0101[0,l)c×[0,l)ch2(t(1-s),lξ)gats2(ly1l+x+(-1)vlξ)gats2(ly2l)dy1dy2dsdξ
=tl20101h2(t(1-s),lξ)gats2(x+(-1)vlξ)gats2(0)dsdξ
+t20101[0,l)ch2(t(1-s),lξ)gats2(x+(-1)vlξ)gats2(ly2l)dy2dsdξ
+t20101[0,l)ch2(t(1-s),lξ)gats2(ly1l+x+(-1)vlξ)gats2(0)dy1dsdξ
+t2l0101[0,l)c×[0,l)ch2(t(1-s),lξ)gats2(ly1l+x+(-1)vlξ)gats2(ly2l)dy1dy2dsdξ.
On the other hand, using the importance sampling technique and Fubini’s theorem, we have

v=010l0tI1(t-s,x,ξ,v)h2(s,ξ)dsdξ=v=01t2l0101R2h2(t(1-s),lξ)gats2(ly1l+x+(-1)vlξ)gats2(ly2l)×gats2(y1+x+(-1)vlξ)gats2(y1+x+(-1)vlξ)gats2(y2)gats2(y2)dy1dy2dsdξ=v=01t2l0101E[h2(t(1-s),lξ)gats2(lats2Z1-x-(-1)vlξl+x+(-1)vlξ)gats2(ats2Z1)gats2(lats2Z2l)gats2(ats2Z2)]dsdξ=v=01t2lE[h2(t(1-U1),lU2)gatU12(latU12Z1-x-(-1)vlU2l+x+(-1)vlU2)gatU12(atU12Z1)gatU12(latU12Z2l)gatU12(atU12Z2)].

Note that gt(x)=g1(x/t)/t for any xR and t>0; thus we have

(A.7)v=010l0tI1(t-s,x,ξ,v)h2(s,ξ)dsdξ=E[v=01t2lh2(t(1-U1),lU2)g1(θ1(v,t,x,U1,U2,Z1))g1(Z1)g1(φ1(t,U1,Z2))g1(Z2)],

where θ1 and φ1 are defined by

θ1(v,t,x,U1,U2,Z1):=2atU1(latU12Z1-x-(-1)vlU2l+x+(-1)vlU2),
φ1(t,U1,Z2):=2atU1(latU12Z2l).
From the result of splitting the first term of (A.6) and (A.7), we obtain

v=010l0tI1(t-s,x,ξ,v)h2(s,ξ)dsdξ=E[v=01tl2h2(t(1-U1),lU2)gatU12(x+(-1)vlU2)gatU12(0)]+E[v=01t2h2(t(1-U1),lU2)gatU12(x+(-1)vlU2)g1(φ1(t,U1,Z2))g1(Z2)1[0,1)c(φ~1(t,U1,Z2))]+E[v=01t2h2(t(1-U1),lU2)g1(θ1(v,t,x,U1,U2,Z1))g1(Z1)gatU12(0)1[0,1)c(θ~1(v,t,x,U1,U2,Z1))]+E[v=01t2lh2(t(1-U1),lU2)g1(θ1(v,t,x,U1,U2,Z1))g1(Z1)g1(φ1(t,U1,Z2))g1(Z2)×1[0,1)c×[0,1)c(θ~1(v,t,x,U1,U2,Z1),φ~1(t,U1,Z2))],

where θ~1 and φ~1 are defined by

θ~1(v,t,x,U1,U2,Z1):=atU12Z1-x-(-1)vlU2l,φ~1(t,U1,Z2):=atU12Z2l.

In the same way,

(A.8)v=010l0tI2(t-s,x,ξ,v)h2(s,ξ)dsdξ=E[v=01t2lh2(t(1-U1),lU2)g1(θ2(v,t,U1,U2,Z1))g1(Z1)g1(φ2(t,x,U1,Z1))g1(Z2)],

where θ2 and φ2 are defined by

θ2(v,t,U1,U2,Z1):=2atU1(latU12Z1-(-1)vlU2l+(-1)vlU2),
φ2(t,x,U1,Z1):=2atU1(latU12Z2-xl+x),
and from splitting the second term of (A.6) as for the first term, we obtain
v=010l0tI2(t-s,x,ξ,v)h2(s,ξ)dsdξ
=E[v=01tl2h2(t(1-U1),lU2)gatU12((-1)vlU2)gatU12(x)]
+E[v=01t2h2(t(1-U1),lU2)gatU12((-1)vlU2)g1(φ2(t,x,U1,Z2))g1(Z2)1[0,1)c(φ~2(t,x,U1,Z2))]
+E[v=01t2h2(t(1-U1),lU2)g1(θ2(v,t,U1,U2,Z1))g1(Z1)gatU12(x)1[0,1)c(θ~2(v,t,U1,U2,Z1))]
+E[v=01t2lh2(t(1-U1),lU2)g1(θ2(v,t,U1,U2,Z1))g1(Z1)g1(φ2(t,x,U1,Z2))g1(Z2)
×1[0,1)c×[0,1)c(θ~2(v,t,U1,U2,Z1),φ~2(t,x,U1,Z2))],
where θ~2 and φ~2 are defined by

θ~2(v,t,U1,U2,Z1):=atU12Z1-(-1)vlU2l,φ~2(t,x,U1,Z2):=atU12Z2-xl.

Therefore, we have concluded the proof. ∎

Remark A.4

For 1<p2, the 𝑝-order moments of the probability representation in (A.7) and in (A.8) are not finite. Therefore, we have applied the importance sampling technique after splitting the space R2 to

[0,l)×[0,l)[0,l)×[0,l)c[0,l)c×[0,l)[0,l)c×[0,l)c.

A.4 Proof of the probability representation of (C4)

We consider the solution of SPDE (3.1) with condition (C4). It is well known that the fundamental solution of PDE (2.1) with L=t-a2x2 with condition (C4) is given by

Gt(x,y)=nZ{gat(y-x-2nl)-gat(-y-x-2nl)}

(see e.g. [19, p. 52]). Using the same arguments as in Subsection A.3, we have the following representation for m(t,x):

m(t,x)=E[(-1)Q(Wat)+1u0(X(Wat))].

Moreover, by using I1 and I2 defined in (A.4) and (A.5), respectively, we have

Gt2(x,ξ)=v=01I1(t,x,ξ,v)-v=01I2(t,x,ξ,v).

Hence we obtain the following lemma using the same arguments as in the proof of Lemma A.3.

Lemma A.5

For any (t,x)R++×S, there exist RVs UiU(0,1) and ZiN(0,1) for each i=1,2 such that (U1,U2,Z1,Z2) are mutually independent, respectively, and

σ(t,x)2=E[j=12(-1)j-1{v=01tl2h2(t(1-U1),lU2)gatU12((-1)vlU2+δ1,jx)gatU12(δ2,jx)
+v=01t2h2(t(1-U1),lU2)gatU12((-1)vlU2+δ1,jx)g1(φj(t,x,U1,Z2))g1(Z2)1[0,1)c(φ~j(t,x,U1,Z2))
+v=01t2h2(t(1-U1),lU2)g1(θj(v,t,x,U1,U2,Z1))g1(Z1)gatU12(δ2,jx)1[0,1)c(θ~j(v,t,x,U1,U2,Z1))
+v=01t2lh2(t(1-U1),lU2)g1(θj(v,t,x,U1,U2,Z1))g1(Z1)g1(φj(t,x,U1,Z2))g1(Z2)
×1[0,1)c×[0,1)c(θ~j(v,t,x,U1,U2,Z1),φ~j(t,x,U1,Z2))}],
where θj,θ~j,φj and φ~j are defined in Lemma A.3.

Acknowledgements

The authors would like to thank Professor Arturo Kohatsu-Higa for his valuable suggestions and fruitful discussions. The authors would also like to thank Associate Professor Dai Taguchi for his careful reading and advice. The authors would also like to thank an anonymous referee for his/her careful readings.

References

[1] E. J. Allen, S. J. Novosel and Z. Zhang, Finite element and difference approximation of some linear stochastic partial differential equations, Stoch. Stoch. Rep. 64 (1998), no. 1–2, 117–142. 10.1080/17442509808834159Search in Google Scholar

[2] A. Andersson, M. Kovács and S. Larsson, Weak error analysis for semilinear stochastic Volterra equations with additive noise, J. Math. Anal. Appl. 437 (2016), no. 2, 1283–1304. 10.1016/j.jmaa.2015.09.016Search in Google Scholar

[3] A. Andersson and S. Larsson, Weak convergence for a spatial approximation of the nonlinear stochastic heat equation, Math. Comp. 85 (2016), no. 299, 1335–1358. 10.1090/mcom/3016Search in Google Scholar

[4] C.-E. Bréhier, Influence of the regularity of the test functions for weak convergence in numerical discretization of SPDEs, J. Complexity 56 (2020), Article ID 101424. 10.1016/j.jco.2019.101424Search in Google Scholar

[5] C.-E. Bréhier and A. Debussche, Kolmogorov equations and weak order analysis for SPDEs with nonlinear diffusion coefficient, J. Math. Pures Appl. (9) 119 (2018), 193–254. 10.1016/j.matpur.2018.08.010Search in Google Scholar

[6] C.-E. Bréhier and L. Goudenège, Weak convergence rates of splitting schemes for the stochastic Allen–Cahn equation, BIT 60 (2020), no. 3, 543–582. 10.1007/s10543-019-00788-xSearch in Google Scholar

[7] C.-E. Bréhier, M. Hairer and A. M. Stuart, Weak error estimates for trajectories of SPDEs under spectral Galerkin discretization, J. Comput. Math. 36 (2018), no. 2, 159–182. 10.4208/jcm.1607-m2016-0539Search in Google Scholar

[8] D. Conus, A. Jentzen and R. Kurniawan, Weak convergence rates of spectral Galerkin approximations for SPDEs with nonlinear diffusion coefficients, Ann. Appl. Probab. 29 (2019), no. 2, 653–716. 10.1214/17-AAP1352Search in Google Scholar

[9] G. Da Prato and J. Zabczyk, Stochastic Equations in Infinite Dimensions, Encyclopedia Math. Appl. 44, Cambridge University, Cambridge, 1992. 10.1017/CBO9780511666223Search in Google Scholar

[10] A. Debussche, Weak approximation of stochastic partial differential equations: the nonlinear case, Math. Comp. 80 (2011), no. 273, 89–117. 10.1090/S0025-5718-2010-02395-6Search in Google Scholar

[11] R. Durrett, Probability: Theory and Examples, Cambridge University, Cambridge, 2010. 10.1017/CBO9780511779398Search in Google Scholar

[12] E. Hausenblas, Weak approximation of the stochastic wave equation, J. Comput. Appl. Math. 235 (2010), no. 1, 33–58. 10.1016/j.cam.2010.03.026Search in Google Scholar

[13] M. Hefter, A. Jentzen and R. Kurniawan, Weak convergence rates for numerical approximations of stochastic partial differential equations with nonlinear diffusion coefficients in UMD Banach spaces, preprint (2016), https://arxiv.org/abs/1612.03209. Search in Google Scholar

[14] L. Jacobe de Naurois, A. Jentzen and T. Welti, Weak convergence rates for spatial spectral Galerkin approximations of semilinear stochastic wave equations with multiplicative noise, preprint (2015), https://arxiv.org/abs/1508.05168. 10.1007/s00245-020-09744-6Search in Google Scholar

[15] L. Jacobe de Naurois, A. Jentzen and T. Welti, Lower bounds for weak approximation errors for spatial spectral Galerkin approximations of stochastic wave equations, Stochastic Partial Differential Equations and Related Fields, Springer Proc. Math. Stat. 229, Springer, Cham (2018), 237–248. 10.1007/978-3-319-74929-7_13Search in Google Scholar

[16] A. Jentzen and R. Kurniawan, Weak convergence rates for Euler-type approximations of semilinear stochastic evolution equations with nonlinear diffusion coefficients, Found. Comput. Math. 21 (2021), no. 2, 445–536. 10.1007/s10208-020-09448-xSearch in Google Scholar

[17] A. Lang and A. Petersson, Monte Carlo versus multilevel Monte Carlo in weak error simulations of SPDE approximations, Math. Comput. Simulation 143 (2018), 99–113. 10.1016/j.matcom.2017.05.002Search in Google Scholar

[18] J. Ma and J. Yong, Forward-Backward Stochastic Differential Equations and Their Applications, Lecture Notes in Math. 1702, Springer, Berlin, 1999. Search in Google Scholar

[19] A. D. Polyanin and V. E. Nazaikinskii, Handbook of Linear Partial Differential Equations for Engineers and Scientists, CRC Press, Boca Raton, 2015. 10.1201/b19056Search in Google Scholar

[20] J. B. Walsh, An introduction to stochastic partial differential equations, École d’été de probabilités de Saint-Flour XIV—1984, Lecture Notes in Math. 1180, Springer, Berlin (1986), 265–439. 10.1007/BFb0074920Search in Google Scholar

[21] X. Wang and S. Gan, Weak convergence analysis of the linear implicit Euler method for semilinear stochastic partial differential equations with additive noise, J. Math. Anal. Appl. 398 (2013), no. 1, 151–169. 10.1016/j.jmaa.2012.08.038Search in Google Scholar