On the basis of integral representations we propose fast numerical methods to solve the Cauchy problem for the stochastic wave equation without boundaries and with the Dirichlet boundary conditions. The algorithms are exact in a probabilistic sense

Andreas Martin

Gerhard Winkler‡

Martin A.

Orsingher E.

Sergei M. Prigarin†

Walsh J. B.

English

Crossref

Exact and fast numerical algorithms for the stochastic wave equation

On the basis of integral representations we propose fast numerical methods to solve the Cauchy problem for the stochastic wave equation without boundaries and with Dirichlet boundary conditions. The algorithms are exact in a probabilistic sense

Martin, A.

Prigarin, S.M.

Winkler, G.

PuSH

Exact and fast numerical algorithms for the stochastic wave equation.

On the basis of integral representations we propose fast numerical methods to solve the Cauchy problem for the stochastic wave equation without boundaries and with the Dirichlet boundary conditions. The algorithms are exact in a probabilistic sense. 1 Statement of the problem and integral representations Consider the stochastic wave equatio

Sergei M. Prigarin

Gerhard Winkler

CiteSeerX

Martin, Andreas

Prigarin, Sergej M.

Winkler, Gerhard

EconStor (ZBW Kiel)

Prigarin, Sergej

Universität München: Elektronischen Publikationen

Martin, Prigarin, Winkler:
Exact and Fast Numerical Algorithms for the
Stochastic Wave Equation
Sonderforschungsbereich 386, Paper 349 (2003)
Online unter: http://epub.ub.uni-muenchen.de/
Projektpartner
Exact and Fast Numerical Algorithms for the
Stochastic Wave Equation ∗
Andreas Martin†, Sergei M. Prigarin‡and Gerhard Winkler §
April 9, 2003
Abstract
On the basis of integral representations we propose fast numerical
methods to solve the Cauchy problem for the stochastic wave equation
without boundaries and with the Dirichlet boundary conditions. The
algorithms are exact in a probabilistic sense.
1 Statement of the problem and integral rep-
resentations
Consider the stochastic wave equation
∂2X
∂t2
(t, x)− a2∂
2X
∂x2
(t, x) = g(t, x) + f(t, x) dW (t, x) (1)
with random initial conditions, where W is a Gaussian white noise on the
plane. We give shortly the precise statement of the problem. The aim of
∗AMS Subject Classification. Primary, 60h15, 68U20.
Key words and phrases. stochastic wave equation, numerical simulation , isonormal pro-
cess.
This work was partly supported by a DAAD Doktorandenstipendium im Rahmen des
gemeinsamen Hochschulsonderprogramms III von Bund und Ländern.
†Institute of Biomathematics and Biometry, GSF-National Research Centre for Envi-
ronment and Health, Ingolstädter Landstraße 1, 85764 Neuherberg/München, Germany.
E-mail: martin@gsf.de
‡Institute of Computational Mathematics and Mathematical Geophysics, pr. Lavren-
tieva 6, Novosibirsk, 630090, Russia. E-mail: smp@osmf.sscc.ru
§Institute of Biomathematics and Biometry, GSF-National Research Centre for Envi-
ronment and Health, Ingolstädter Landstraße 1, 85764 Neuherberg/München, Germany.
E-mail: gwinkler@gsf.de
1
1 STATEMENT OF THE PROBLEM 2
the paper is to construct numerical algorithms in order to solve the Cauchy
problem for equation (1) without boundaries and with the Dirichlet boundary
conditions, i.e., to simulate realizations of the random field X.
First, we present the necessary formalism for the stochastic wave equation
and integral representations for the solutions that underlay our numerical
methods. The details can be found in [5], [1] and [2].
Let (Ω, F , P) be a complete probability space, W be an isonormal process
indexed by L2(R+×R,B(R+×R), λ) and {Ft; t ∈ R+} be a nullset augmented
filtration which satisfies
– W (1A) is Ft-measurable whenever A ∈ B([0, t]× R)
with λ(A) < ∞;
– W (1A) is independent of Ft whenever A ∈ B((t,∞)× R)
with λ(A) < ∞.
Moreover, let F = {F (x); x ∈ R} and H = {H(x); x ∈ R} be almost surely
continuous and F -measurable processes such that
∫
R E{|F (x)|
2}ϕ(x) dx and∫
R E{|H(x)|
2}ϕ(x) dx are finite for any smooth ϕ with compact support.
Finally, assume f and g to be real valued functions on R+ × R such that
fϕ is square integrable and gϕ is integrable for any smooth ϕ with compact
support, respectively.
Similar to the deterministic case, the solution of equation (1) with initial
conditions (F, H) on R+×R can be represented in terms of Green’s function
Ga (hereafter a > 0),
Ga(t, x, s, y) =
{
1
2a
if 0 ≤ s ≤ t and |x− y| ≤ a(t− s)
0 otherwise.
The solution is given by
X(t, x) = u0(t, x) +
∫
R+
∫
R
Ga(t, x, s, y)g(s, y) dy ds (2)
+
∫
R+
∫
R
Ga(t, x, s, y)f(s, y) dW (s, y),
where u0(t, x) = (F (x− at) + F (x + at))/2 +
∫ x+at
x−at H(y)/2a dy.
Remark 1.1 We require the integrability conditions on F and H to have
first and second moments for the random field X. An example for a homo-
geneous wave equation with random initial conditions is studied in [4].
1 STATEMENT OF THE PROBLEM 3
Remark 1.2 The stochastic integral in (2) is defined by∫
R+
∫
R
h(s, y) dW (s, y) := W (h)
for any h ∈ L2(R+ × R,B(R+ × R), λ) and we use the integral notation for
convenience.
Remark 1.3 It holds X(0, ·) = F almost surely, and for f = 0 we have
∂X/∂t(0, ·) = H. However, for f 6= 0 the last equation has only formal
meaning, because the random field X is not differentiable.
For the Cauchy problem on R+ × [0, 1] with Dirichlet boundary conditions
Green’s function ΣGa has the following form
ΣGa(t, x, s, y) =
∑
k
[
Ga(t, x + 2k; s, y)1[0,1](y)
−Ga(t,−x + 2k; s, y)1[0,1](y)
]
,
and instead of u0 we consider
u0,D(t, x) =
1
2
∑
k
[
F (x + 2k − at)1(0,1)(x + 2k − at)
+ F (x + 2k + at)1(0,1)(x + 2k + at)
− F (−x + 2k − at)1(0,1)(−x + 2k − at)
− F (−x + 2k + at)1(0,1)(−x + 2k + at)
]
+
1
2a
∑
k
[ ∫ (x+2k+at)∧1
(x+2k−at)∨0
H(y) dy −
∫ (−x+2k+at)∧1
(−x+2k−at)∨0
H(y) dy
]
.
In this case the solution of (1) with initial conditions (F, H) is
X(t, x) = u0,D(t, x) +
∫
R+
∫ 1
0
ΣGa(t, x, s, y)g(s, y) dy ds (3)
+
∫
R+
∫ 1
0
ΣGa(t, x, s, y)f(s, y) dW (s, y).
In both cases the solutions are Gaussian random fields and the integral rep-
resentations give the basis to construct “exact” numerical solutions.
2 NUMERICAL ALGORITHMS 4
2 Numerical algorithms
A realization of the solution will be simulated on a grid
G = {(tj, xk) = (jh, akh)}
with j ∈ {0, 1/2, 1, 3/2, . . .}, k ∈ {. . . ,−1,−1/2, 0, 1/2, 1, . . .}, where h is
a time step, ah is a spatial step and for a grid point the indices j, k are
simultaneously either integers or fractional numbers.
In addition to the values Xjk = X(jh, akh) on the grid G the algorithm
operates with the values µjk on an “adjacent” grid G∗ = {(tj, xk) = (jh, akh)}
with the same sets of indices, but one of the indices is always integer while the
other is fractional. Moreover, µjk is given by µ
j
k =
∫ xk+ah/2
xk−ah/2
Hj(y) dy, where
Hj(y) is the second part of the initial conditions for the Cauchy problem at
time tj (see Remark 1.3).
2.1 Cauchy problem without boundaries
The algorithm uses the Markov property of X with respect to time, i.e. the
values for the next time level depend only on the values obtained at the
previous time level. Therefore the complexity of the algorithm is linear.
Initialisation. Set j = 0. The values X0k are taken from the initial
condition, i.e. X0k = F (akh) = X(0, akh) = X(t0, xk). Moreover, set µ
0
k =∫ xk+ah/2
xk−ah/2
H(y) dy.
Step 1. Equation (3) gives an explicit rule to calculate Xj+1/2:
X
j+1/2
k+1/2 =
(
1
2
(
Xjk + X
j
k+1
)
+
1
2a
µjk+1/2
)
+ ξ
j+1/2
k+1/2, (4)
where
ξ
j+1/2
k+1/2 =
∫ tj+1/2
tj
∫ xk+1
xk
Ga(tj+1/2, xk+1/2, s, y)g(s, y) dy ds
+
∫ tj+1/2
tj
∫ xk+1
xk
Ga(tj+1/2, xk+1/2, s, y)f(s, y) dW (s, y)
is a Gaussian random variable with mean∫ tj+1/2
tj
∫ xk+1
xk
Ga(tj+1/2, xk+1/2, s, y)g(s, y) dy ds
and variance∫ tj+1/2
tj
∫ xk+1
xk
G2a(tj+1/2, xk+1/2, s, y)f
2(s, y) dy ds.
2 NUMERICAL ALGORITHMS 5
Step 2. The next step is to calculate µ
j+1/2
k+1 . To this end, one observes
that Xj+1k+1 can be expressed in two different ways,
Xj+1k+1 =
(
1
2
(
X
j+1/2
k+1/2 + X
j+1/2
k+3/2
)
+
1
2a
µ
j+1/2
k+1
)
+ ξj+1k+1 (5)
and
Xj+1k+1 =
(
1
2
(
Xjk + X
j
k+2
)
+
1
2a
(
µjk+1/2 + µ
j
k+3/2
))
(6)
+
(
ξ
j+1/2
k+1/2 + ξ
j+1/2
k+3/2 + δ
j
k+1 + ξ
j+1
k+1
)
,
where
δjk+1 =
1
2a
∫ tj+1/2
tj
∫ xk+1+a(s−tj)
xk+1−a(s−tj)
g(s, y) dy ds (7)
+
1
2a
∫ tj+1/2
tj
∫ xk+1+a(s−tj)
xk+1−a(s−tj)
f(s, y) dW (s, y)
is a Gaussian random variable with mean given by the first double integral
in (7) and variance(
1
2a
)2 ∫ tj+1/2
tj
∫ xk+1+a(s−tj)
xk+1−a(s−tj)
f 2(s, y) dy ds.
A combination of (5) and (6) yields
µ
j+1/2
k+1 = 2a
[
1
2
(
Xjk + X
j
k+2
)
+
1
2a
(
µjk+1/2 + µ
j
k+3/2
)
(8)
+
(
ξ
j+1/2
k+1/2 + ξ
j+1/2
k+3/2
)
+ δjk+1 −
1
2
(
X
j+1/2
k+1/2 + X
j+1/2
k+3/2
)]
.
Thus, at step two the values µ
j+1/2
k+1 are simulated according to expression (8).
Note, that all the components of the expression, except δjk+1, are defined at
the previous steps of the algorithm, and the random variable ξj+1k+1 appears
in (5) and (6) but cancels out in (8).
Cycling. Set j = j + 1/2 and go back to step 1.
Proposition 2.1 The sequence (Xj,k) on the grid G fulfils Xj,k = X(tj, xk)
with probability one.
Proof. On the grid G Green’s function Gα can be decomposed into triangles
with corners given by the grid points. Thus induction on the time level tj
gives the assertion. 2
2 NUMERICAL ALGORITHMS 6
2.2 Cauchy problem with boundary conditions
This algorithm is a modification of the previous one. As ah is a step size in
space, here we assume that N = 1/ah is an integer.
Initialisation. Set j = 0 and determine X0k for k ∈ {0, . . . , N} and
µ0k+1/2 for k ∈ {0, . . . , N −1} by the initial conditions F and H, respectively.
Step 1. The values X
j+1/2
k+1/2 for k ∈ {0, . . . , N−1} are calculated according
to equation (4).
Step 2. The values µ
j+1/2
k+1 for k ∈ {0, . . . , N−2} are calculated according
to equation (8).
Step 3. The values Xj+1k for k ∈ {1, . . . , N −1} are calculated according
to equation (4) while Xj+10 and X
j+1
N are defined by the boundary conditions,
i.e., Xj+10 = X
j+1
N = 0.
Step 4. Now, the values µj+1k+1 have to be calculated. For k ∈ {12 ,
3
2
, . . . , N−5
2
}
equation (8) is used. To find µj+11/2 and µ
j+1
N−1/2 the boundary conditions should
be applied. The value X
j+3/2
1/2 can be expressed in two different ways (note
that Xj+10 = 0),
X
j+3/2
1/2 =
(
1
2
(
Xj+10 + X
j+1
1
)
+
1
2a
µj+11/2
)
+ ξ
j+3/2
1/2
=
(
1
2
Xj+11 +
1
2a
µj+11/2
)
+ ξ
j+3/2
1/2 ,
X
j+3/2
1/2 =
(
1
2
(
−Xj+1/21/2 + X
j+1/2
3/2
)
+
1
2a
µ
j+1/2
1
)
+
(
ξj+11 + δ
j+1/2
1/2 + ξ
j+3/2
1/2
)
.
From these two representations we obtain the formula for simulation,
µj+11/2 = 2a
[
1
2
(
−Xj+1/21/2 + X
j+1/2
3/2
)
+
1
2a
(
µ
j+1/2
1
)
+
(
ξj+11
)
+ δ
j+1/2
1/2 −
1
2
(
Xj+11
)]
.
Similarly, for µj+1N−1/2 we have
µj+1N−1/2 = 2a
[
1
2
(
X
j+1/2
N−3/2 −X
j+1/2
N−1/2
)
+
1
2a
(
µ
j+1/2
N−1
)
+
(
ξj+1N−1
)
+ δ
j+1/2
N−1/2 −
1
2
(
Xj+1N−1
)]
.
Cycling. Set j = j + 1 and go back to Step 1.
The sequence (Xj,k) fulfils Xj,k = X(tj, xk) with probability one as before.
3 CONCLUSION 7
2.3 Results of simulation
Results of numerical simulations are presented in Figs. 1 and 2. In addition,
in Fig. 3 we present an example of the sample size dependence of the empir-
ical correlation between two points of the solution for the Cauchy problem
without boundaries.
Figure 1: Realization of the solution of the Cauchy problem without bound-
aries; zero initial conditions and harmonic initial conditions.
Figure 2: Realization of the solution of the Cauchy problem with Dirichlet
boundary conditions; zero initial conditions and harmonic initial conditions.
3 Conclusion
The two proposed methods are “exact” in the sense that the finite-dimensional
distributions of the field X on the grid G coincide with the finite-dimensional
distributions of the numerical solution. Moreover, the complexity of the al-
gorithms is linear and they are therefore attractive for numerical simulations.
Remark 3.1 “Exact” algorithms to solve the stochastic Klein-Gordon equa-
tion
∂2X/∂t2(t, x)− ∂2X/∂x2(t, x) = αX(t, x) + dW (t, x)
REFERENCES 8
Figure 3: Empirical correlations with confidence intervals of level 0.95 via
the sample size; exact correlation is equal to 0.25.
were discussed in [3]. But the algorithms proposed in that paper don’t have
the Markov property and, practically, result in a general simulation of Gaus-
sian random vectors with correlated components. Because of this fact, it can
turn out that the algorithms are not feasible for large grids.
References
[1] R. Carmona and D. Nualart. Random nonlinear wave equations: Smooth-
ness of the solutions. Probab. Theory Relat. Fields, 79(4):469–508, 1988.
[2] A. Martin. Hyperbolic stochastic partial differential equations: Small
balls and simulation; propagation of singularities. PhD thesis, Aachen:
Shaker Verlag. Neuherberg: Institut f. Biomathematik u. Biometrie, GSF
- Forschungszentrum f. Umwelt u. Gesundheit, 2002.
[3] A. Martin, S.M. Prigarin, and G. Winkler. Exact numerical algorithms for
linear stochastic wave equation and stochastic Klein-Gordon equation. In
G.A. Mikhailov (ed.) et al., Proceedings to the International Conference
on Computational Mathematics. ICM&MG Publisher, Prospect Lavren-
tieva 6, 630090 Novosibirks, Russia. 2002.
[4] E. Orsingher. Vibrations with random initial conditions. Boll. Unione
Mat. Ital., VI. Ser., B, 4:541–556, 1985.
[5] J.B. Walsh. An introduction to stochastic partial differential equations.
In École d’été de probabilités de Saint-Flour XIV - 1984, Lect. Notes
Math. 1180, 265-437. 1986.


Exact and Fast Numerical Algorithms for the Stochastic Wave Equation

Open Access LMU

Martin, Prigarin, Winkler:Exact and Fast Numerical Algorithms for theStochastic Wave EquationSonderforschungsbereich 386, Paper 349 (2003)Online unter: http://epub.ub.uni-muenchen.de/ProjektpartnerExact and Fast Numerical Algorithms for theStochastic Wave Equation ∗Andreas Martin†, Sergei M. Prigarin‡and Gerhard Winkler §April 9, 2003AbstractOn the basis of integral representations we propose fast numericalmethods to solve the Cauchy problem for the stochastic wave equationwithout boundaries and with the Dirichlet boundary conditions. Thealgorithms are exact in a probabilistic sense.1 Statement of the problem and integral rep-resentationsConsider the stochastic wave equation∂2X∂t2(t, x)− a2∂2X∂x2(t, x) = g(t, x) + f(t, x) dW (t, x) (1)with random initial conditions, where W is a Gaussian white noise on theplane. We give shortly the precise statement of the problem. The aim of∗AMS Subject Classification. Primary, 60h15, 68U20.Key words and phrases. stochastic wave equation, numerical simulation , isonormal pro-cess.This work was partly supported by a DAAD Doktorandenstipendium im Rahmen desgemeinsamen Hochschulsonderprogramms III von Bund und La¨ndern.†Institute of Biomathematics and Biometry, GSF-National Research Centre for Envi-ronment and Health, Ingolsta¨dter Landstraße 1, 85764 Neuherberg/Mu¨nchen, Germany.E-mail: martin@gsf.de‡Institute of Computational Mathematics and Mathematical Geophysics, pr. Lavren-tieva 6, Novosibirsk, 630090, Russia. E-mail: smp@osmf.sscc.ru§Institute of Biomathematics and Biometry, GSF-National Research Centre for Envi-ronment and Health, Ingolsta¨dter Landstraße 1, 85764 Neuherberg/Mu¨nchen, Germany.E-mail: gwinkler@gsf.de11 STATEMENT OF THE PROBLEM 2the paper is to construct numerical algorithms in order to solve the Cauchyproblem for equation (1) without boundaries and with the Dirichlet boundaryconditions, i.e., to simulate realizations of the random field X.First, we present the necessary formalism for the stochastic wave equationand integral representations for the solutions that underlay our numericalmethods. The details can be found in [5], [1] and [2].Let (Ω, F , P) be a complete probability space,W be an isonormal processindexed by L2(R+×R,B(R+×R), λ) and {Ft; t ∈ R+} be a nullset augmentedfiltration which satisfies– W (1A) is Ft-measurable whenever A ∈ B([0, t]× R)with λ(A) <∞;– W (1A) is independent of Ft whenever A ∈ B((t,∞)× R)with λ(A) <∞.Moreover, let F = {F (x); x ∈ R} and H = {H(x); x ∈ R} be almost surelycontinuous and F -measurable processes such that ∫R E{|F (x)|2}ϕ(x) dx and∫R E{|H(x)|2}ϕ(x) dx are finite for any smooth ϕ with compact support.Finally, assume f and g to be real valued functions on R+ × R such thatfϕ is square integrable and gϕ is integrable for any smooth ϕ with compactsupport, respectively.Similar to the deterministic case, the solution of equation (1) with initialconditions (F,H) on R+×R can be represented in terms of Green’s functionGa (hereafter a > 0),Ga(t, x, s, y) ={12aif 0 ≤ s ≤ t and |x− y| ≤ a(t− s)0 otherwise.The solution is given byX(t, x) = u0(t, x) +∫R+∫RGa(t, x, s, y)g(s, y) dy ds (2)+∫R+∫RGa(t, x, s, y)f(s, y) dW (s, y),where u0(t, x) = (F (x− at) + F (x+ at))/2 +∫ x+atx−at H(y)/2a dy.Remark 1.1 We require the integrability conditions on F and H to havefirst and second moments for the random field X. An example for a homo-geneous wave equation with random initial conditions is studied in [4].1 STATEMENT OF THE PROBLEM 3Remark 1.2 The stochastic integral in (2) is defined by∫R+∫Rh(s, y) dW (s, y) :=W (h)for any h ∈ L2(R+ × R,B(R+ × R), λ) and we use the integral notation forconvenience.Remark 1.3 It holds X(0, ·) = F almost surely, and for f = 0 we have∂X/∂t(0, ·) = H. However, for f 6= 0 the last equation has only formalmeaning, because the random field X is not differentiable.For the Cauchy problem on R+ × [0, 1] with Dirichlet boundary conditionsGreen’s function ΣGa has the following formΣGa(t, x, s, y) =∑k[Ga(t, x+ 2k; s, y)1[0,1](y)−Ga(t,−x+ 2k; s, y)1[0,1](y)],and instead of u0 we consideru0,D(t, x) =12∑k[F (x+ 2k − at)1(0,1)(x+ 2k − at)+ F (x+ 2k + at)1(0,1)(x+ 2k + at)− F (−x+ 2k − at)1(0,1)(−x+ 2k − at)− F (−x+ 2k + at)1(0,1)(−x+ 2k + at)]+12a∑k[ ∫ (x+2k+at)∧1(x+2k−at)∨0H(y) dy −∫ (−x+2k+at)∧1(−x+2k−at)∨0H(y) dy].In this case the solution of (1) with initial conditions (F,H) isX(t, x) = u0,D(t, x) +∫R+∫ 10ΣGa(t, x, s, y)g(s, y) dy ds (3)+∫R+∫ 10ΣGa(t, x, s, y)f(s, y) dW (s, y).In both cases the solutions are Gaussian random fields and the integral rep-resentations give the basis to construct “exact” numerical solutions.2 NUMERICAL ALGORITHMS 42 Numerical algorithmsA realization of the solution will be simulated on a gridG = {(tj, xk) = (jh, akh)}with j ∈ {0, 1/2, 1, 3/2, . . .}, k ∈ {. . . ,−1,−1/2, 0, 1/2, 1, . . .}, where h isa time step, ah is a spatial step and for a grid point the indices j, k aresimultaneously either integers or fractional numbers.In addition to the values Xjk = X(jh, akh) on the grid G the algorithmoperates with the values µjk on an “adjacent” grid G∗ = {(tj, xk) = (jh, akh)}with the same sets of indices, but one of the indices is always integer while theother is fractional. Moreover, µjk is given by µjk =∫ xk+ah/2xk−ah/2 Hj(y) dy, whereHj(y) is the second part of the initial conditions for the Cauchy problem attime tj (see Remark 1.3).2.1 Cauchy problem without boundariesThe algorithm uses the Markov property of X with respect to time, i.e. thevalues for the next time level depend only on the values obtained at theprevious time level. Therefore the complexity of the algorithm is linear.Initialisation. Set j = 0. The values X0k are taken from the initialcondition, i.e. X0k = F (akh) = X(0, akh) = X(t0, xk). Moreover, set µ0k =∫ xk+ah/2xk−ah/2 H(y) dy.Step 1. Equation (3) gives an explicit rule to calculate Xj+1/2:Xj+1/2k+1/2 =(12(Xjk +Xjk+1)+12aµjk+1/2)+ ξj+1/2k+1/2, (4)whereξj+1/2k+1/2 =∫ tj+1/2tj∫ xk+1xkGa(tj+1/2, xk+1/2, s, y)g(s, y) dy ds+∫ tj+1/2tj∫ xk+1xkGa(tj+1/2, xk+1/2, s, y)f(s, y) dW (s, y)is a Gaussian random variable with mean∫ tj+1/2tj∫ xk+1xkGa(tj+1/2, xk+1/2, s, y)g(s, y) dy dsand variance∫ tj+1/2tj∫ xk+1xkG2a(tj+1/2, xk+1/2, s, y)f2(s, y) dy ds.2 NUMERICAL ALGORITHMS 5Step 2. The next step is to calculate µj+1/2k+1 . To this end, one observesthat Xj+1k+1 can be expressed in two different ways,Xj+1k+1 =(12(Xj+1/2k+1/2 +Xj+1/2k+3/2)+12aµj+1/2k+1)+ ξj+1k+1 (5)andXj+1k+1 =(12(Xjk +Xjk+2)+12a(µjk+1/2 + µjk+3/2))(6)+(ξj+1/2k+1/2 + ξj+1/2k+3/2 + δjk+1 + ξj+1k+1),whereδjk+1 =12a∫ tj+1/2tj∫ xk+1+a(s−tj)xk+1−a(s−tj)g(s, y) dy ds (7)+12a∫ tj+1/2tj∫ xk+1+a(s−tj)xk+1−a(s−tj)f(s, y) dW (s, y)is a Gaussian random variable with mean given by the first double integralin (7) and variance(12a)2 ∫ tj+1/2tj∫ xk+1+a(s−tj)xk+1−a(s−tj)f 2(s, y) dy ds.A combination of (5) and (6) yieldsµj+1/2k+1 = 2a[12(Xjk +Xjk+2)+12a(µjk+1/2 + µjk+3/2)(8)+(ξj+1/2k+1/2 + ξj+1/2k+3/2)+ δjk+1 −12(Xj+1/2k+1/2 +Xj+1/2k+3/2)].Thus, at step two the values µj+1/2k+1 are simulated according to expression (8).Note, that all the components of the expression, except δjk+1, are defined atthe previous steps of the algorithm, and the random variable ξj+1k+1 appearsin (5) and (6) but cancels out in (8).Cycling. Set j = j + 1/2 and go back to step 1.Proposition 2.1 The sequence (Xj,k) on the grid G fulfils Xj,k = X(tj, xk)with probability one.Proof. On the grid G Green’s function Gα can be decomposed into triangleswith corners given by the grid points. Thus induction on the time level tjgives the assertion. 22 NUMERICAL ALGORITHMS 62.2 Cauchy problem with boundary conditionsThis algorithm is a modification of the previous one. As ah is a step size inspace, here we assume that N = 1/ah is an integer.Initialisation. Set j = 0 and determine X0k for k ∈ {0, . . . , N} andµ0k+1/2 for k ∈ {0, . . . , N −1} by the initial conditions F and H, respectively.Step 1. The valuesXj+1/2k+1/2 for k ∈ {0, . . . , N−1} are calculated accordingto equation (4).Step 2. The values µj+1/2k+1 for k ∈ {0, . . . , N−2} are calculated accordingto equation (8).Step 3. The values Xj+1k for k ∈ {1, . . . , N −1} are calculated accordingto equation (4) while Xj+10 and Xj+1N are defined by the boundary conditions,i.e., Xj+10 = Xj+1N = 0.Step 4. Now, the values µj+1k+1 have to be calculated. For k ∈ {12 , 32 , . . . , N−52 }equation (8) is used. To find µj+11/2 and µj+1N−1/2 the boundary conditions shouldbe applied. The value Xj+3/21/2 can be expressed in two different ways (notethat Xj+10 = 0),Xj+3/21/2 =(12(Xj+10 +Xj+11)+12aµj+11/2)+ ξj+3/21/2=(12Xj+11 +12aµj+11/2)+ ξj+3/21/2 ,Xj+3/21/2 =(12(−Xj+1/21/2 +Xj+1/23/2)+12aµj+1/21)+(ξj+11 + δj+1/21/2 + ξj+3/21/2).From these two representations we obtain the formula for simulation,µj+11/2 = 2a[12(−Xj+1/21/2 +Xj+1/23/2)+12a(µj+1/21)+(ξj+11)+ δj+1/21/2 −12(Xj+11)].Similarly, for µj+1N−1/2 we haveµj+1N−1/2 = 2a[12(Xj+1/2N−3/2 −Xj+1/2N−1/2)+12a(µj+1/2N−1)+(ξj+1N−1)+ δj+1/2N−1/2 −12(Xj+1N−1)].Cycling. Set j = j + 1 and go back to Step 1.The sequence (Xj,k) fulfilsXj,k = X(tj, xk) with probability one as before.3 CONCLUSION 72.3 Results of simulationResults of numerical simulations are presented in Figs. 1 and 2. In addition,in Fig. 3 we present an example of the sample size dependence of the empir-ical correlation between two points of the solution for the Cauchy problemwithout boundaries.Figure 1: Realization of the solution of the Cauchy problem without bound-aries; zero initial conditions and harmonic initial conditions.Figure 2: Realization of the solution of the Cauchy problem with Dirichletboundary conditions; zero initial conditions and harmonic initial conditions.3 ConclusionThe two proposed methods are “exact” in the sense that the finite-dimensionaldistributions of the field X on the grid G coincide with the finite-dimensionaldistributions of the numerical solution. Moreover, the complexity of the al-gorithms is linear and they are therefore attractive for numerical simulations.Remark 3.1 “Exact” algorithms to solve the stochastic Klein-Gordon equa-tion∂2X/∂t2(t, x)− ∂2X/∂x2(t, x) = αX(t, x) + dW (t, x)REFERENCES 8Figure 3: Empirical correlations with confidence intervals of level 0.95 viathe sample size; exact correlation is equal to 0.25.were discussed in [3]. But the algorithms proposed in that paper don’t havethe Markov property and, practically, result in a general simulation of Gaus-sian random vectors with correlated components. Because of this fact, it canturn out that the algorithms are not feasible for large grids.References[1] R. Carmona and D. Nualart. Random nonlinear wave equations: Smooth-ness of the solutions. Probab. Theory Relat. Fields, 79(4):469–508, 1988.[2] A. Martin. Hyperbolic stochastic partial differential equations: Smallballs and simulation; propagation of singularities. PhD thesis, Aachen:Shaker Verlag. Neuherberg: Institut f. Biomathematik u. Biometrie, GSF- Forschungszentrum f. Umwelt u. Gesundheit, 2002.[3] A. Martin, S.M. Prigarin, and G. Winkler. Exact numerical algorithms forlinear stochastic wave equation and stochastic Klein-Gordon equation. InG.A. Mikhailov (ed.) et al., Proceedings to the International Conferenceon Computational Mathematics. ICM&MG Publisher, Prospect Lavren-tieva 6, 630090 Novosibirks, Russia. 2002.[4] E. Orsingher. Vibrations with random initial conditions. Boll. UnioneMat. Ital., VI. Ser., B, 4:541–556, 1985.[5] J.B. Walsh. An introduction to stochastic partial differential equations.In E´cole d’e´te´ de probabilite´s de Saint-Flour XIV - 1984, Lect. NotesMath. 1180, 265-437. 1986.

An introduction to stochastic partial diﬀerential equations.

Exact numerical algorithms for linear stochastic wave equation and stochastic Klein-Gordon equation.

Hyperbolic stochastic partial diﬀerential equations: Small balls and simulation; propagation of singularities.

Random nonlinear wave equations: Smoothness of the solutions.

Vibrations with random initial conditions.

https://epub.ub.uni-muenchen.de/1725/1/paper_349.pdf

Exact and Fast Numerical Algorithms for the Stochastic Wave Equation

Abstract

Similar works

Full text

Available Versions