Geiser J. Decomposition Methods for Differential Equations: Theory and Applications
Подождите немного. Документ загружается.
Numerical Experiments 177
TABLE 6.41: Maximal temperatures for
numerical experiments, depending on the anisotropy
coefficients (α
1
r
,α
1
z
) of the insulation (cf. Figures 6.20
α
1
r
α
1
z
Maximal temperature
[K]
1 1 1273.18
1 10 1232.15
1-10, mixed 1-10, mixed 1238.38
10 1 918.35
6.5.3 Elastic Wave Propagation
In these sections, we focus on decoupling the wave equations motivated by a
realistic problem involving seismic sources and waves. The complicated wave
propagation process is decomposed into simpler processes with respect to the
directions of the propagation. We discretize and solve the simpler equations
by more accurate higher-order methods, see [103].
With the test examples, we verify the decomposition methods for the linear
acoustic wave equations. Moreover, the delicate initialization of the right-
hand side, while using a Dirac function, is presented with respect to the
decomposition method. The benefit of spatial splitting methods, applied to
the multidimensional operators, is discussed in several examples.
6.5.3.1 Real-Life Application of the Elastic Wave Propagation
In this application, we refer to our underlying model equations presented
in Section 1.3.3. Based on these equations we discuss the discretization and
decomposition.
We focus on a higher-order splitting method for a real-life application prob-
lem involving the approximations of time and space; we also focus on how the
Dirac function is approximated.
During numerical testing we have observed a need to reduce the allowable
time-step in the case where the ratio of λ over μ becomes too large. This
is likely from the influence of the explicitly treated mixed derivative. For
really high ratios (> 20) a reduction of 35% was necessary to avoid numerical
instabilities.
6.5.3.2 Initial Values and Boundary Conditions
In order to start the time-stepping scheme we need to know the values at
two earlier time levels. Starting at time t = 0, we know the value at level
n =0asU
0
= g
0
. The value at level n = −1 can be obtained by Taylor
© 2009 by Taylor & Francis Group, LLC
- 6.21).
178Decomposition Methods for Differential Equations Theory and Applications
expansion as shown below:
U
−1
= U
0
− τ∂
t
U
0
+
τ
2
n
2
∂
tt
U
0
−
τ
3
6
∂
ttt
U
0
+
τ
4
24
∂
tttt
U
0
+ O(τ
5
), (6.99)
where we use
∂
t
U
0
j,k
= g
1
j,k
, (6.100)
∂
tt
U
0
j,k
≈
1
ρ
M
4
g
0
j,k
)+f
j,k
, (6.101)
∂
ttt
U
0
j,k
≈
1
ρ
M
4
g
1
j,k
)+∂
t
f
0
j,k
, (6.102)
∂
tttt
U
0
j,k
≈
1
ρ
M
2
2
g
0
j,k
)+M
4
f
0
j,k
+ ∂
tt
f
0
j,k
, (6.103)
and also for (6.102) and (6.103).
The approximation of our right-hand side is given as
∂
t
f
0
j,k
≈
f
1
j,k
− f
−1
j,k
2τ
, (6.104)
∂
tt
f
0
j,k
≈
f
1
j,k
− 2f
0
j,k
+ f
−1
j,k
τ
2
n
. (6.105)
We are not considering the boundary value problem and thus we will not be
concerned with constructing proper difference stencils at grid points near the
boundaries of the computational domain. We have simply added a two-point-
thick layer of extra grid points at the boundaries of the domain and assigned
the correct analytical solution at all points in the layer for every time-step.
REMARK 6.17 For the Dirichlet boundary conditions, the splitting
method, see (4.35)–(4.37), also conserves the conditions. For the three equa-
tions (i.e., for U
∗
, U
∗∗
and for U
n+1
), we can use the same conditions.
For the Neumann boundary conditions and other boundary conditions of
higher order, we also have to split the boundary conditions with respect to
the split operators, see [154].
6.5.3.3 Test Example of the Two-Dimensional Wave Equation
The first test example is a two-dimensional example. The splitting method
is presented in Section 4.2.3 as well as the model equations.
We apply in the first test case a forcing function given as
f =
sin(t − x)sin(y) − 2μ sin(t − x)sin(y)
−(λ + μ)(cos(x)cos(t − y)+sin(t − x)sin(y)),
sin(t − y)sin(x) − 2Vs
2
sin(x)sin(t − y)
−(λ + μ)(cos(t − x)cos(y)+sin(y)sin(t − y))
T
, (6.106)
© 2009 by Taylor & Francis Group, LLC
Numerical Experiments 179
giving the analytical solution
U
true
=
sin(x − t)sin(y), sin(y − t)sin(x)
T
. (6.107)
Using the splitting method, we solved (1.10) on a domain Ω = [−1, 1]×[−1, 1]
and time interval t ∈ [0, 2]. We used two sets of material parameters; for the
first case ρ, λ,andμ were all equal to 1, for the second case ρ and μ were
1andλ was set to 14. Solving on four different grids with a refinement
factor of two in each direction between the successive grids, we obtained the
results shown in Table 6.42. For all test examples, the equidistant time-step
is given as τ =0.0063. The errors are measured in the L
∞
-norm defined as
||U
j,k
|| = max (max
j,k
|u
j,k
|, max
j,k
|v
j,k
|).
TABLE 6.42: Errors in max-norm for decreasing h and
smooth analytical solution U
true
. Convergence rate indicates
fourth-order convergence for the split scheme.
Grid step Time t =2,e
h
=err
U,L
∞
= ||U
n
− U
true
||
∞
h case 1 log
2
(
e
2h
e
h
) case 2 log
2
(
e
2h
e
h
)
0.05 1.7683e-07 2.5403e-07
0.025 1.2220e-08 3.855 2.1104e-08 3.589
0.0125 7.9018e-10 3.951 1.4376e-09 3.876
0.00625 5.0013e-11 3.982 9.2727e-11 3.955
As can be seen, we get the expected fourth-order convergence for problems
with smooth solutions.
To check the influence of the splitting error N
4,θ
on the error, we solved
the same problems using the nonsplit scheme (4.34). The results are shown in
Table 6.43. The errors are only marginally smaller than for the split scheme.
6.5.3.4 Singular Forcing Terms
In seismology and acoustics, it is common to use spatial singular forcing
terms that can look like
f = Fδ(x)g(t), (6.108)
where F is a constant direction vector. A numeric method for Equation (1.10),
see Chapter 1, needs to approximate the Dirac function correctly in order to
achieve full convergence. Obviously, we cannot expect convergence close to
the source as the solution will be singular for two- and three-dimensional
domains.
The analyses in [188] and [195] demonstrate that it is possible to derive
regularized approximations of the Dirac function that result in a point-wise
© 2009 by Taylor & Francis Group, LLC
180Decomposition Methods for Differential Equations Theory and Applications
TABLE 6.43: Errors in max-norm for
decreasing h and smooth analytical solution
U
true
and using the non-split scheme.
Comparing with Table 6.42 we see that the
splitting error is very small for this case.
Time t =2,
Grid step e
h
=err
U,L
∞
= ||U
n
− U
true
||
∞
h case 1 case 2
0.05 1.6878e-07 2.4593e-07
0.025 1.1561e-08 2.0682e-08
0.0125 7.4757e-10 1.4205e-09
0.00625 4.8112e-11 9.2573e-11
convergence of the solution away from the sources. Based on these analy-
ses, we define one second-order (δ
h
2
) and one fourth-order (δ
h
4
) regularized
approximation of the one-dimensional Dirac function:
δ
h
2
(˜x)=
1
h
⎧
⎨
⎩
1+˜x, −h ≤ ˜x<0,
1 − ˜x, 0 ≤ ˜x<h,
0, elsewhere,
(6.109)
δ
h
4
(˜x)=
1
h
⎧
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎩
1+
11
6
˜x +
5
8
˜x
2
+
1
6
˜x
3
, −2h ≤ ˜x<−h,
1+
1
2
˜x − ˜x
2
−
1
2
˜x
3
, −h ≤ ˜x<0,
1 −
1
2
˜x − ˜x
2
+
1
2
˜x
3
, 0 ≤ ˜x<h,
1 −
11
6
˜x +˜x
2
−
1
6
˜x
3
,h≤ ˜x<2h,
0, elsewhere,
(6.110)
where ˜x = x/h. The two- and three-dimensional Dirac functions are then
approximated as δ
h
2,4
(˜x)δ
h
2,4
(˜y)andδ
h
2,
4
(˜x)δ
h
2,4
(˜y)δ
h
2,4
(˜z). The chosen time
dependence was a smooth function given by
g(t)=
$
exp(−1/(t(1 − t))), 0 ≤ t<1,
0, elsewhere,
(6.111)
which is C
∞
. Using this forcing function, we can compute the analytical
solution by integrating the Green’s function given in [61]. The integration
was done using numerical quadrature routines from MATLAB. Figures 6.22
and 6.23 show examples of what the errors look like on a radius passing
through the singular source at time t =0.8 for different grid sizes h and
the two approximations δ
h
2
and δ
h
4
. We see that the error is smooth and
converges a small distance away from the source. However, using δ
h
2
limits
the convergence to second-order, while using δ
h
4
gives the full fourth-order
convergence away from the singular source. When t>1, the forcing goes
to zero and the solution will be smooth everywhere. Table 6.44 shows the
convergence behavior at time t =1.1 for four different grids. Note that the full
© 2009 by Taylor & Francis Group, LLC
Numerical Experiments 181
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−40
−35
−30
−25
−20
−15
−10
−5
0
distance from source
2−logarithm of |error|
FIGURE 6.22: The 2-logarithm of the error along a line going through
the source point for a point force located at x =0,y = 0, and approximated
in space by (6.110). Note that the error decays with O(h
4
)awayfromthe
source but not near it. The grid sizes were
h =0.05 (−·), 0.025 (·), 0.0125 (−), 0.00625 (∗). The numerical quadrature
had an absolute error of approximately 10
−11
≈ 2
−36
, so the error cannot be
resolved beneath that limit.
© 2009 by Taylor & Francis Group, LLC
182Decomposition Methods for Differential Equations Theory and Applications
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−40
−35
−30
−25
−20
−15
−10
−5
0
distance from source
2−logarithm of |error|
FIGURE 6.23: The 2-logarithm of the error along a line going through
the source point for a point force located at x =0,y = 0, and approximated
in space by (6.109). Note that the error only decays with O (h
2
)awayfrom
the source. The grid sizes were
h =0.05 (−·), 0.025 (·), 0.0125 (−), 0.00625 (∗).
© 2009 by Taylor & Francis Group, LLC
Numerical Experiments 183
TABLE 6.44: Errors in max-norm for
decreasing h and analytical solution U
true
.
Convergence rate approaches fourth-order
after the singular forcing term goes to zero of
the two-dimensional split scheme.
Time t =1.1,
Grid step e
h
=err
U,L
∞
= ||U
n
− U
true
||
∞
h case 1 log
2
(
e
2h
e
h
)
0.05 1.1788e-04
0.025 1.4146e-05 3.0588
0.0125 1.3554e-06 3.3836
0.00625 1.0718e-07 3.6606
0.003125 7.1890e-09 3.8981
convergence is achieved even if the lower-order δ
h
2
is used as an approximation
for the Dirac function.
The convergence rate approaches four as we refine the grids, even though
the solution was singular up to time t =1.
We implement the two-dimensional case in MATLAB on a single CPU.
Further computations are also done on multiple CPUs with optimal speed.
The spatial grid steps are given as
h =Δx =Δy =Δy ∈{0.05, 0.025, 0.0125, 0.00625, 0.003125},andthetime-
step is given as τ =0.0063. The material parameters for the elastic wave
propagation equations are λ =1,μ =1,andρ = 1. The total number of grid
points are about 26,000.
A visualization of the elastic wave propagation is given in Figures 6.24 and
6.25.
To visualize the influence of the singular point force for the two-dimensional
computation we use the x-component of the solution at time t =1overthe
normed domain [−1, 1] × [−1, 1].
6.5.3.5 Computational Cost of the Splitting Method
For a two-dimensional problem, the fourth-order explicit method (4.33) can
be implemented using approximately 160 floating point operations (flops) per
grid point.
The splitting method requires approximately 120 flops (first step) plus two
times 68 flops (second and third step) for a total of 256 flops. This increase of
about 60% in the number of flops is somewhat offset by the larger time-steps
allowed by the splitting method, especially for “nice” material properties,
making the two methods roughly comparable in computational cost.
© 2009 by Taylor & Francis Group, LLC
184Decomposition Methods for Differential Equations Theory and Applications
FIGURE 6.24:Thex-component of the solution for a singular point force
at time t = 1 and spatial grid step h =0.0125.
6.5.3.6 A Three-Dimensional Splitting Method
As the splitting method discussed for two dimensions in Section 4.2.3, we
extend the splitting method to three dimensions in the applications. There-
fore, a first discussion of the discretization and decomposition of the operator
is important. In three dimensions, a fourth-order difference approximation of
the operator becomes
M
4
=
⎛
⎜
⎜
⎜
⎝
(λ +2μ)
1 −
h
2
12
D
x
2
D
x
2
+ μ
1 −
h
2
12
D
y
2
D
y
2
+1−
h
2
12
D
z
2
D
z
2
(λ + μ)
1 −
h
2
6
D
x
2
D
x
0
1 −
h
2
6
D
y
2
D
y
0
(λ + μ)
1 −
h
2
6
D
x
2
D
x
0
1 −
h
2
6
D
z
2
D
z
0
(λ + μ)
1 −
h
2
6
D
x
2
D
x
0
1 −
h
2
6
D
y
2
D
y
0
(λ +2μ)
1 −
h
2
12
D
y
2
D
y
2
+ μ
1 −
h
2
12
D
x
2
D
x
2
+ μ
1 −
h
2
12
D
z
2
D
z
2
(λ + μ)
1 −
h
2
6
D
z
2
D
z
0
1 −
h
2
6
D
y
2
D
y
0
(λ + μ)
1 −
h
2
6
D
x
2
D
x
0
1 −
h
2
6
D
z
2
D
z
0
(λ + μ)
1 −
h
2
6
D
y
2
D
y
0
1 −
h
2
6
D
z
2
D
z
0
(λ +2μ)
1 −
h
2
12
D
z
2
D
z
2
+ μ
1 −
h
2
12
D
x
2
D
x
2
+1−
h
2
12
D
y
2
D
y
2
⎞
⎟
⎟
⎟
⎠
,
© 2009 by Taylor & Francis Group, LLC
Numerical Experiments 185
operating on grid functions U
n
j,k,l
defined at grid points x
j
,y
k
,z
l
,t
n
similarly
to the two-dimensional case. We can split M
4
into six parts; M
xx
, M
yy
,
M
zz
containing the three second-order directional difference operators, and
M
xy
, M
yz
, M
xz
containing the mixed difference operators.
We could split this scheme in a number of different ways depending on
how we treat the mixed derivative terms. We have chosen to implement the
following split scheme in three dimensions:
1.ρ
U
∗
j,k,l
− 2U
n
j,k,l
+ U
n−1
j,k,l
τ
2
n
= M
4
U
n
j,k,l
+ θf
n+1
j,k,l
+(1− 2θ)f
n
j,k,l
+ θf
n−1
j,k,l
),
2.ρ
U
∗∗
j,k,l
− U
∗
j,k,l
τ
2
n
= θM
xx
U
∗∗
j,k,l
− 2U
n
j,k,l
+ U
n−1
j,k,l
+
θ
2
(M
xy
+ M
xz
)
U
∗
j,k,l
− 2U
n
j,k,l
+ U
n−1
j,k,l
,
3.ρ
U
∗∗∗
j,k,l
− U
∗∗
j,k,l
τ
2
n
= θM
xx
U
∗∗∗
j,k,l
− 2U
n
j,k,l
+ U
n−1
j,k,l
+
θ
2
(M
xy
+ M
yz
)
U
∗∗
j,k,l
− 2U
n
j,k,l
+ U
n−1
j,k,l
,
4.ρ
U
n+1
j,k,l
− U
∗∗∗
j,k,l
τ
2
n
= θM
xx
U
n+1
j,k,l
− 2U
n
j,k,l
+ U
n−1
j,k,l
+
θ
2
(M
xz
+ M
yz
)
U
∗∗∗
j,k,l
− 2U
n
j,k,l
+ U
n−1
j,k,l
.
The properties such as splitting error, accuracy, stability, and so forth, for
the three-dimensional case are similar to the two-dimensional case treated in
the earlier sections.
6.5.3.7 Test Example of the Three-Dimensional Wave Equation
We have performed some numerical experiments with the three-dimensional
scheme in order to test the convergence and stability. We apply our splitting
scheme, presented in Section 6.5.3.6, and use a forcing
f =
− (−1+λ +4μ)sin(t − x)sin(y)sin(z) −
(λ + μ)cos(x)(2 sin(t)sin(y)sin(z)+cos(t)sin(y + z)),
−(−1+λ +4μ)sin(x)sin(t − y)sin(z) −
(λ + μ)cos(y)(2 sin(t)sin(x)sin(z)+cos(t)sin(x + z)),
−(λ + μ)cos(t − y)cos(z)sin(x) − sin(y)((λ + μ)c
os(t − x)cos(z)+
(−1+λ +4μ)sin(x)sin(t − z))
T
, (6.112)
© 2009 by Taylor & Francis Group, LLC
186Decomposition Methods for Differential Equations Theory and Applications
giving the analytical solution
U
true
=
sin(x − t)sin(y)sin(z),
sin(y − t)sin(x)sin(z),
sin(z − t)sin(x)sin(y)
T
. (6.113)
As in the earlier examples, we tested this for a number of different grid sizes.
Using the same two sets of material parameters as for the two-dimensional
case, we computed up until t = 2 and checked the maximum error for all
components of the solution. The results are given in Table 6.45.
TABLE 6.45: Errors in max-norm for decreasing h
and smooth analytical solution U
true
. Convergence rate
indicates fourth-order convergence of the
three-dimensional split scheme.
t =2,e
h
=err
U,L
∞
= ||U
n
− U
true
||
∞
h case 1 log
2
(
e
2h
e
h
) case 2 log
2
(
e
2h
e
h
)
0.1 4.2986e-07 1.8542e-06
0.05 3.5215e-08 3.61 1.3605e-07 3.77
0.025 3.0489e-09 3.53 8.0969e-09 4.07
0.0125 2.0428e-10 3.90 4.7053e-10 4.10
We implement the three-dimensional case in MATLAB on a single CPU.
Further computations are also done on multiple CPUs with optimal speed.
The spatial grid steps are given as
h =Δx =Δy =Δy ∈{0.1, 0.05, 0.025, 0.0125}, and the time-step is given as
τ = τ =0.0046. The parameters for the elastic wave propagation are given as
λ =14,μ =1,andρ = 1. We obtain a total number of grid points of about
110
6
.
The initial solution is given as U =0.
In Figure 6.25, we visualize the influence of the singular point force for
the three-dimensional computation using the y-component on a plane. The
solution is given at the time t = 1 and at the normed domain [−1, 1]×[−1, 1]×
[−1, 1]. Further, the isosurfaces for the solution U are also visualized to show
the wave fronts.
REMARK 6.18 Our splitting scheme has been shown to work well
in practice for different types of material properties. It is comparable to the
fully explicit fourth-order scheme (4.33) in terms of computational cost but
should be easier to implement, as no difference approximations of higher-order
operators are needed.
© 2009 by Taylor & Francis Group, LLC