Geiser J. Decomposition Methods for Differential Equations: Theory and Applications
Подождите немного. Документ загружается.
76Decomposition Methods for Differential Equations Theory and Applications
4.2.1 Classical ADI Method
Our classical method is based on the splitting method of [67] and [149].
We discuss the method for a simple wave equation given as
∂
2
c(x, y, t)
∂t
2
=(A + B) c(x, y, t)+f(x, y, t), for (x, y, t) ∈ Ω × [0,T], (4.7)
c(x, y, 0) = c
0
(x, y),
∂c(x, y, 0)
∂t
= c
1
(x, y), for (x, y) ∈ Ω, (Initial conditions),
c(x, y, t)=0, for (x, y, t) ∈ ∂Ω × (0,T), (Boundary conditions),
where the operators A = D
1
∂
2
∂x
2
and B = D
2
∂
2
∂y
2
are given, with D
1
,D
2
∈ R
+
.
We assume that c exists and belongs to C
4
(Ω × [0,T]).
We denote a grid l for the discretization with the nodes (iΔx, jΔy), where
Δx, Δy>0andi, j are integers. Our discretized domain is given as
Ω
h
=
l ∩ Ω. For the time discretization we denote τ
n
= t
n+1
−t
n
and n =0, 1,...,N
with t
N+1
= T .
The second-order ADI method is given as, see [149],
˜c − 2c(t
n
)+c(t
n−1
)=τ
2
n
A(η˜c +(1− 2η)c(t
n
)+ηc(t
n−1
)) (4.8)
+τ
2
n
Bc(t
n
)
+τ
2
n
(ηf(t
n+1
)+(1− 2η)f(t
n
)+ηf(t
n−1
)),
in Ω
h
× [t
n−1
,t
n+1
],
c(t
n+1
) − 2c(t
n
)+c(t
n−1
)=τ
2
n
A(η˜c +(1− 2η)c(t
n
)+ηc(t
n−1
)) (4.9)
+τ
2
n
B(ηc(t
n+1
)+(1− 2η)c(t
n
)+ηc(t
n−1
))
+τ
2
n
(ηf(t
n+1
)+(1− 2η)f(t
n
)+ηf(t
n−1
)),
in Ω
h
× [t
n−1
,t
n+1
],
where the result is given as c(t
n+1
) with the initial conditions c(t
n
)=c
n
,
dc(t
n
)
dt
=
dc
n
dt
as the known split approximation at the time level t = t
n
.We
have η ∈ [0, 0.5], and the boundary condition is given as c(x, y, t)=0onthe
boundary ∂Ω
h
× [t
n−1
,t
n+1
]. The local time-step is given as τ
n
= t
n+1
− t
n
.
The fully coupled explicit method is for η = 0, and the decoupled implicit
method is for 0 <η≤ 0.5, which is a mixing of explicit and implicit Euler
methods.
We obtain the following consistency result.
THEOREM 4.1
We apply the second-order discretization in time and space for our underlying
ADI method. We then obtain the consistency error
|Lc − L
τ
n
,h
c| = O(τ
2
n
+ h
2
), (4.10)
© 2009 by Taylor & Francis Group, LLC
Decomposition Methods for Hyperbolic Equations 77
for sufficiently smooth functions c ∈ C
4
(Ω×[0,T]). We have h =max{Δx, Δy},
Lc =
∂
2
c
∂t
2
−(A + B)c is the wave operator (d’Alembert operator), and the dis-
crete operator with respect to the ADI method is given as
L
τ
n
,h
c = ∂
+
t
∂
−
t
c − (D
1
∂
+
x
∂
−
x
+ D
2
∂
+
y
∂
−
y
)c + R
h
(c)
= ∂
+
t
∂
−
t
c − (D
1
∂
+
x
∂
−
x
+ D
2
∂
+
y
∂
−
y
)c
−ητ
2
n
(D
1
∂
+
x
∂
−
x
+ D
2
∂
+
y
∂
−
y
)∂
+
t
∂
−
t
c + η
2
τ
4
n
(∂
+
t
∂
−
t
(D
1
∂
+
x
∂
−
x
+ D
2
∂
+
y
∂
−
y
)c),
where the residual of the ADI method is given as
R
h
(c)=−ητ
2
n
(D
1
∂
+
x
∂
−
x
+D
2
∂
+
y
∂
−
y
)∂
+
t
∂
−
t
c+η
2
τ
4
n
(∂
+
t
∂
−
t
(D
1
∂
+
x
∂
−
x
+D
2
∂
+
y
∂
−
y
)c),
see [149].Here∂
+
p
is a forward and ∂
−
p
is a backward difference quotient with
respect to the variable p.
PROOF
We calculate the consistency error Lc − L
τ
n
,h
c,whichisgivenas
|Lc − L
τ
n
,h
c| = |(∂
t
∂
t
c − (D
1
∂
x
∂
x
+ D
2
∂
y
∂
y
)c)
−(∂
+
t
∂
−
t
c − (D
1
∂
+
x
∂
−
x
+ D
2
∂
+
y
∂
−
y
)c + R
h
(c))|
= O(τ
2
n
+ h
2
), (4.11)
for 0 <η≤ 0.5. Time and space are discretized as second-order methods, and
therefore, we obtain a second-order consistency result, see also [149].
REMARK 4.1 Higher-order ADI methods are more complicated to
develop, see [67]. Therefore, we concentrate on the LOD method, in which it
is possible to derive a fourth-order method.
4.2.2 LOD Method
The local one-dimensional (LOD) method is discussed in [133].
The intent of this method is to combine explicit and implicit discretiza-
tion methods to obtain simpler equations and gain a higher-order splitting
method, while the splitting error maintains the same order as the discretiza-
tion method.
The wave equation is given as
∂
2
c(x, y, t)
∂t
2
=(A
1
+ A
2
)c(x, y, t)
+ f (x, y, t), for (x, y, t) ∈ Ω × (0,T), (4.12)
c(x, y, 0) = c
0
(x, y),
∂c(x, y, 0)
∂t
= c
1
(x, y), for (x, y) ∈ Ω,
c(x, y, t)=0, for (x, y, t) ∈ ∂Ω × (0,T),
© 2009 by Taylor & Francis Group, LLC
78Decomposition Methods for Differential Equations Theory and Applications
where the unbounded operators are given as A
1
= D
1
∂
2
∂x
2
and A
2
= D
2
∂
2
∂y
2
,
with D
1
,D
2
∈ R
+
.
Further, we assume that the order of the spatial discretization is given by
˜
A
1
c ≈ A
1
c + O(h
p
), (4.13)
˜
A
2
c ≈ A
2
c + O(h
p
), (4.14)
with the order O(h
p
), p =2, 4.
Furthermore, we apply the second-order finite difference method for the
discretization in time and we obtain
∂
+
t
∂
−
t
c :=
c
n+1
− 2c
n
+ c
n−1
τ
2
n
. (4.15)
The LOD method for our Equation (4.12) is then
c
n+1,0
− 2c
n
+ c
n−1
= ητ
2
n
(
˜
A
1
+
˜
A
2
)c
n
, (4.16)
c
n+1,1
− c
n+1,0
= ητ
2
n
˜
A
1
(c
n+1,1
− 2c
n
+ c
n−1
), (4.17)
c
n+1
− c
n+1,1
= ητ
2
n
˜
A
2
(c
n+1
− 2c
n
+ c
n−1
), (4.18)
for η ∈ [0, 0.5], where the result is given as c(t
n+1
) with the initial conditions
c(t
n
)=c
n
,
dc(t
n
)
dt
=
dc
n
dt
as the known split approximation at the time level
t = t
n
.
Equation (4.16) is the explicit method and gains all operators. The equa-
tions (4.17)–(4.18) are implicit and obtain only one operator.
REMARK 4.2 The LOD method obtains the local one-dimensional
directions of each operator. Therefore, we can construct unconditional stable
methods.
In the next subsection, we discuss the stability and consistency of a higher-
order LOD method based on wave equations.
4.2.2.1 Stability and Consistency of the Fourth-Order Splitting
Method
The consistency of the fourth-order splitting method is given in the following
theorem.
We assume the following discretization orders of O(h
p
), p =2, 4forthe
discretization in space, in which h = h
x
= h
y
is the spatial step size. Further,
we assume a second-order time discretization, where τ
n
is the equidistant
time-step and t
n
= nτ
n
.
Then we obtain the following consistency of our method (4.16)–(4.18).
© 2009 by Taylor & Francis Group, LLC
Decomposition Methods for Hyperbolic Equations 79
THEOREM 4.2
The consistency of our method is given by
(
d
2
c(t)
dt
2
− Ac) − (∂
+
t
∂
−
t
c −
˜
Ac)=O(τ
4
n
), (4.19)
where ∂
+
t
∂
−
t
is a second-order discretization in time and
˜
A is the discretized
fourth-order spatial operator.
PROOF
We add Equations (4.16)–(4.18) and obtain, see also [133],
∂
+
t
∂
−
t
c
n
−
˜
A(θc
n+1
+(1−2θ)c
n
+ θc
n−1
)
−
˜
B(c
n+1
− 2c
n
+ c
n−1
)=0, (4.20)
where it holds that
˜
B = θ
2
τ
2
n
˜
A
1
˜
A
2
. (4.21)
Therefore, we obtain a splitting error of
˜
B(c
n+1
− 2c
n
+ c
n−1
).
For sufficient smoothness we have (c
n+1
− 2c
n
+ c
n−1
)=O(τ
2
n
), and we
obtain
˜
B(c
n+1
− 2c
n
+ c
n−1
)=O(τ
4
n
).
We thus obtain a fourth-order method, if the spatial operators are also
discretized as fourth-order terms.
DEFINITION 4.1
We define the operator P
±
(c
j
,θ)as
P
±
(c
j
,θ):=θ(
˜
Ac
j
,c
j
)+θ(
˜
Ac
j±1
,c
j±1
)+(1− 2θ)(
˜
Ac
j
,c
j±1
), (4.22)
where c
j
= c(t
j
) is the solution at time t
j
and θ ∈ [0, 0.5].
THEOREM 4.3
The stability of the fourth-order splitting method is given by
||(1 − τ
2
n
˜
B)
1/2
∂
+
t
c
n
||
2
+ P
+
(c
n
,θ)
≤||(1 − τ
2
n
˜
B)
1/2
∂
+
t
c
0
||
2
+ P
+
(c
0
,θ), (4.23)
where θ ∈ [0.25, 0.5].
The proof is given in the following manner.
PROOF
We must proof the theorem for a test function ∂
0
t
c
n
,where∂
0
t
is the
central difference.
© 2009 by Taylor & Francis Group, LLC
80Decomposition Methods for Differential Equations Theory and Applications
We choose ∂
0
t
c
n
as a test function and for n ≥ 1wehave
((1 − τ
2
n
˜
B)∂
+
t
∂
−
t
c
n
,∂
0
t
c
n
)
+(
˜
A(θc
n+1
+(1− 2θ)c
n
+ θc
n−1
),∂
0
t
c
n
)=0. (4.24)
Multiplying by τ
n
and summarizing over j yields
n
j=1
((1 − τ
2
n
˜
B)∂
+
t
∂
−
t
c
j
,∂
0
t
c
j
) τ
n
+(
˜
A(θc
j+1
+(1−2θ)c
j
+ θc
j−1
),∂
0
t
c
j
) τ
n
=0. (4.25)
We can derive the identities
((1 − τ
2
n
˜
B)∂
+
t
∂
−
t
c
j
,∂
0
t
c
j
) τ
n
=1/2||(1 − τ
2
n
˜
B)
1/2
∂
+
t
c
j
||
2
− 1/2||(1 − τ
2
n
˜
B)
1/2
∂
−
t
c
j
||
2
, (4.26)
(
˜
A(θc
j+1
− (1 − 2θ)c
j
+ θc
j−1
),∂
0
t
c
j
) τ
n
=1/2(P
+
(c
j
,θ) −P
−
(c
j
,θ)), (4.27)
and obtain the result
||(1 − τ
2
n
˜
B)
1/2
∂
+
t
c
n
||
2
+ P
+
(c
n
,θ)
≤||(1 − τ
2
n
˜
B)
1/2
∂
+
t
c
0
||
2
+ P
+
(c
0
,θ), (4.28)
see also the methodology of [133].
REMARK 4.3 For θ =
1
12
we obtain a fourth-order method.
For the local splitting error, we must use the multiplier
˜
A
1
˜
A
2
; thus, for large
constants we have an unconditionally small time-step.
4.2.3 LOD Method Applied to a System of Wave Equations
Here we discuss the elastic wave propagation given as a system of linear
wave equations and introduced in Chapter 2.
4.2.3.1 Basic Numerical Methods
The standard explicit difference scheme with central differences used every-
where was introduced in [4]. To save space, we first give an example in two
dimensions.
We discretize uniformly in space and time on the unit square. We then
get a grid Ω
h
with grid points x
j
= jh,y
k
= kh,whereh>0isthespatial
grid size, and a time interval [0,T] with the time points t
n
= nτ
n
,whereτ
n
© 2009 by Taylor & Francis Group, LLC
Decomposition Methods for Hyperbolic Equations 81
is the equidistant time-step. The indices are denoted as j ∈{0, 1,...,J},
k ∈{0, 1,...,K} and n ∈{0, 1,...,N} with J, K, N ∈ N
+
. If we define the
grid function as U
n
j,k
= U (x
j
,y
k
,t
n
)withU =(u, v)
T
or U =(u, v, w), the
basic explicit scheme is given as
ρ
U
n+1
j,k
− 2U
n
j,k
+ U
n−1
j,k
τ
2
n
= M
2
U
n
j,k
+ f
n
j,k
, (4.29)
where M
2
is a difference operator:
M
2
=
(λ +2μ)∂
2
x
+ μ∂
2
y
(λ + μ)∂
0
x
∂
0
y
v
j,k
(λ + μ)∂
0
x
∂
0
y
v
j,k
(λ +2μ)∂
2
y
+ μ∂
2
x
, (4.30)
and we use the standard difference operator notation:
∂
+
x
v
j,k
=
1
h
(v
j+1,k
− v
j,k
),∂
−
x
v
j,k
= ∂
+
x
v
j−1,k
,∂
0
x
=
1
2
∂
+
x
+ ∂
−
x
,
and
∂
2
x
= ∂
+
x
∂
−
x
.
M
2
approximates μ∇
2
+(λ + μ)∇(∇·) to second order. This explicit scheme
is stable for time-steps satisfying
τ
n
<
h
λ +3μ
. (4.31)
Replacing M
2
by M
4
, we have a fourth-order difference operator given by
M
4
=
⎛
⎝
(λ +2μ)
1 −
h
2
12
∂
2
x
∂
2
x
+ μ
1 −
h
2
12
∂
2
y
∂
2
y
(λ + μ)
1 −
h
2
6
∂
2
x
∂
0
x
1 −
h
2
6
∂
2
y
∂
0
y
(λ + μ)
1 −
h
2
6
∂
2
x
∂
0
x
1 −
h
2
6
∂
2
y
∂
0
y
(λ +2μ)
1 −
h
2
12
∂
2
y
∂
2
y
+ μ
1 −
h
2
12
∂
2
x
∂
2
x
⎞
⎠
, (4.32)
and using the modified equation approach [50], we obtain the explicit fourth-
order scheme:
ρ
U
n+1
j,k
− 2U
n
j,k
+ U
n−1
j,k
τ
2
n
= M
4
U
n
j,k
+ f
n
j,k
+
τ
2
n
12
(M
2
U
n
j,k
+ M
2
f
n
i,j
+ ∂
tt
f
n
i,j
), (4.33)
where M
2
is a second-order approximation to the
μ∇
2
+(λ + μ)∇(∇·)
2
operator (i.e., the right-hand-side operator squared). As M
2
only needs to be
second-order accurate, it has the same extent in space as M
4
, and no more
© 2009 by Taylor & Francis Group, LLC
82Decomposition Methods for Differential Equations Theory and Applications
grid points are used. The extra terms on the right-hand side are used to elim-
inate the second-order error in the approximation of the time derivative, and
we obtain a fully fourth-order scheme.
In [133], the following implicit scheme for the scalar wave equation was intro-
duced:
ρ
U
n+1
j,k
− 2U
n
j,k
+ U
n−1
j,k
τ
2
n
= M
4
θU
n+1
j,k
+(1− 2θ)U
n
j,k
+ θU
n−1
j,k
(4.34)
+θf
n+1
j,k
+(1− 2θ)f
n
j,k
+ θf
n−1
j,k
.
When θ =1/12, the error of this scheme is fourth order in time and space.
For this θ-value, it is, however, only conditionally stable.
In order to make it competitive with the explicit scheme (4.33), we will
provide operator-split versions of the implicit scheme (4.34). The presence of
the mixed derivative terms that couple different coordinate directions makes
creating split versions complicated.
4.2.3.2 Fourth-Order Splitting Method
Here, we present a fourth-order splitting method based on the basic scheme (4.34).
We split the operator M
4
into three parts, M
xx
, M
yy
,andM
xy
,wherewe
have
M
xx
=
⎛
⎝
(λ +2μ)
1 −
h
2
12
∂
2
x
∂
2
x
0
0 μ
1 −
h
2
12
∂
+
x
∂
−
x
∂
2
x
⎞
⎠
,
M
yy
=
⎛
⎝
μ
1 −
h
2
12
∂
2
y
∂
2
y
0
0(λ +2μ)
1 −
h
2
12
∂
2
y
∂
2
y
⎞
⎠
,
M
xy
= M
4
−M
xx
−M
yy
.
Our proposed splitting method has the following steps:
1.ρ
U
∗
j,k
− 2U
n
j,k
+ U
n−1
j,k
τ
2
n
= M
4
U
n
j,k
+θf
n+1
j,k
+(1− 2θ)f
n
j,k
+ θf
n−1
j,k
, (4.35)
2.ρ
U
∗∗
j,k
− 2U
∗
j,k
τ
2
n
= θM
xx
U
∗∗
j,k
− 2U
n
j,k
+ U
n−1
j,k
+
θ
2
M
xy
U
∗
j,k
− 2U
n
j,k
+ U
n−1
j,k
, (4.36)
3.ρ
U
n+1
j,k
− 2U
∗∗
j,k
τ
2
n
= θM
xx
U
n+1
j,k
− 2U
n
j,k
+ U
n−1
j,k
+
θ
2
M
xy
U
∗∗
j,k
− 2U
n
j,k
+ U
n−1
j,k
. (4.37)
© 2009 by Taylor & Francis Group, LLC
Decomposition Methods for Hyperbolic Equations 83
Here the first step is explicit, while the second and third steps treat the
derivatives along the coordinate axes implicitly and the mixed derivatives
explicitly. This is similar to how the mixed case is handled for parabolic
problems, see [22].
REMARK 4.4 For the Dirichlet boundary conditions, the splitting
method, see (4.35)–(4.37), also conserves the conditions. We can use for the
three equations, for U
∗
, U
∗∗
,andforU
n+1
, the same conditions, see (4.35)–
(4.37).
For the Neumann boundary conditions and other boundary conditions of
higher order, we also have to split the boundary conditions with respect to
the splitted operators, see [154].
In the next subsection, we discuss the stability and consistency for the LOD
method.
4.2.3.3 Stability and Consistency of the Fourth-Order Splitting
Method for the System of Wave Equations
The consistency of the fourth-order splitting method is given in the follow-
ing theorem.
We have for all sufficiently smooth functions U (x, t) the following discretiza-
tion order:
M
4
U = μ∇
2
U +(λ + μ)∇(∇·U)+O(h
4
), (4.38)
where O(h
4
) is fulfilled for each component.
In the following, we assume the operator O is vectorial (i.e., the accuracy
is given for each component). Furthermore, the split operators are also dis-
cretized with the same order of accuracy.
We obtain the following consistency result for the splitting method (4.35)–
(4.37).
THEOREM 4.4
The splitting method has a splitting error, which for smooth solutions U is
O(τ
4
n
), where we assume τ
n
= O(h) as a stability constraint.
PROOF
We assume here that f =(0, 0)
T
. We add the equations (4.35)–(4.37) and
obtain, as in the scalar case (see [133]), the following result for the discretized
equations:
∂
+
t
∂
−
t
U
n
j,k
−M
4
(θU
n+1
j,k
+(1− 2θ)U
n
j,k
+ θU
n−1
j,k
)
−N
4,θ
(U
n+1
j,k
− 2U
n
j,k
+ U
n−1
j,k
)=(0, 0)
T
, (4.39)
© 2009 by Taylor & Francis Group, LLC
84Decomposition Methods for Differential Equations Theory and Applications
where N
4,θ
= θ
2
τ
2
n
(M
xx
M
yy
+ M
xx
M
xy
+ M
xy
M
yy
)+θ
3
τ
4
n
M
xx
M
yy
M
xy
.
We therefore obtain a splitting error of N
4,θ
(U
n+1
j,k
− 2U
n
j,k
+ U
n−1
j,k
).
For sufficient smoothness, we have (U
n+1
j,k
−2U
n
j,k
+ U
n−1
j,k
)=O(τ
2
n
), and we
obtain N
4,θ
(U
n+1
− 2U
n
+ U
n−1
)=O(τ
4
n
).
It is important that the influence of the mixed terms can also be discretized
with a fourth-order method; therefore, the terms are canceled in the proof.
For the stability, we must denote an appropriate norm, which is in our case
the L
2
(Ω)-norm.
In the following, we introduce the notation of the norms.
REMARK 4.5 For our system, we extend the L
2
-norm
||U||
2
L
2
=(U, U)
L
2
=
Ω
U
2
dx, (4.40)
where U
2
= u
2
+ v
2
or U
2
= u
2
+ v
2
+ w
2
in two and three dimensions.
REMARK 4.6 The matrix
N
4,θ
= θ
2
τ
2
n
(M
xx
M
yy
+ M
xx
M
xy
+ M
xy
M
yy
)
+θ
3
τ
4
n
M
xx
M
yy
M
xy
, (4.41)
where M
xx
, M
yy
and M
xy
are symmetrical and positive-definite matrices.
Therefore, the matrix N
4,θ
is also symmetrical and positive-definite.
Furthermore, we can estimate the norms and define a weighted norm, see
[131] and [164].
REMARK 4.7 The energy norm is given as
(N
4,θ
U, U)
L
2
=
Ω
(N
4,θ
U · U ) dx. (4.42)
Consequently, we can write
||N
4,θ
U|| ≤ ω||U || , ∀ U ∈ H
d
, (4.43)
where ω ∈ R
+
is the weight and N
4,θ
is bounded. The dimension is d and H
is a Sobolev space, so each component of U is in H,see[65].
The stability of the fourth-order splitting method is given in the following
theorem.
© 2009 by Taylor & Francis Group, LLC
Decomposition Methods for Hyperbolic Equations 85
THEOREM 4.5
Let θ ∈ [0.25, 0.5], then the implicit time-stepping algorithm, see (4.29),
and the split procedure, see (4.35)–(4.37), are unconditionally stable. We can
estimate the split procedure iteratively as
||(1 − τ
2
n
˜ω)
1/2
∂
+
t
U
n
j,k
||
2
+ P
+
(U
n
j,k
,θ) ≤||(1 − τ
2
n
˜ω)
1/2
∂
+
t
U
0
j,k
||
2
+ P
+
(U
0
j,k
,θ), (4.44)
where we have
P
±
(U
n
j,k
,θ):=θ(M
4
U
n
j,k
,U
n
j,k
)+θ(M
4
U
n±1
j,k
,U
n±1
j,k
)+(1−2θ)(M
4
U
n
j,k
,U
n±1
j,k
),
and
P
±
≥ 0 for θ ∈ [0.25, 0.5].Further,1−τ
2
n
˜ω ∈ R
+
is the factor for the weighted
norm (I−τ
2
n
N
4,θ
)U ≤ (1 −τ
2
n
˜ω) U, ∀ U ∈ H
d
,whereI is the identity matrix
of H
d
.
We must now prove the iterative estimate for the split procedure.
PROOF
To obtain an energy estimate for the scheme, we multiply by a test function
∂
0
t
U
n
j,k
and summarize over each time level.
The following result is given for the discretized equations; see also Equation
(4.39).
(I−τ
2
n
N
4,θ
)∂
+
t
∂
−
t
U
n
j,k
−M
4
(θU
n+1
j,k
+(1−2θ)U
n
j,k
+ θU
n−1
j,k
)=(0, 0)
T
. (4.45)
For n ≥ 1 we can now rewrite Equation (4.45), and we multiply by a test
function ∂
0
t
U
n
j,k
. For the stability proof, we obtain
((I−τ
2
n
N
4,θ
)∂
+
t
∂
−
t
U
n
j,k
,∂
0
t
U
n
j,k
)
−(M
4
(θU
n+1
j,k
+(1−2θ)U
n
j,k
+ θU
n−1
j,k
),∂
0
t
U
n
j,k
)=0. (4.46)
Multiplying by τ
n
and summarizing over the time levels, we obtain
N
n=1
((I−τ
2
n
N
4,θ
)∂
+
t
∂
−
t
U
n
j,k
,∂
0
t
U
n
j,k
)τ
n
−
N
n=1
(M
4
(θU
n+1
j,k
+(1− 2θ)U
n
j,k
+ θU
n−1
j,k
),∂
0
t
U
n
j,k
)τ
n
=0. (4.47)
For each term of the sum we can derive the identities. For example, for
© 2009 by Taylor & Francis Group, LLC