Multi-stage waveform Relaxation and Multisplitting Methods for Differential Algebraic Systems

0.14 MB·

by  Juergen Geiser

Multi-stage waveform Relaxation and
Multisplitting Methods for Differential
Algebraic Systems
6 Ju¨rgen Geiser1
1
0
RuhrUniversityof Bochum, Department of Electrical Engineering and Information
2
Technology, Universit¨atsstraße 150, D-44801 Bochum, Germany,

n [email protected]
a
J

4

Abstract. We are motivated to solve differential algebraic equations
]
A withnewmulti-stageandmultisplittingmethods.Themulti-stagestrat-
egy ofthewaveform relaxation (WR)methodsaregiven withouterand
N
inner iterations. While the outer iterations decouple the initial value
h. problemofdifferentialalgebraicequations(DAEs)intheformofAdy(t)+
dt
at By(t) = f(t) to MAdykd+t1(t) +M1yk+1(t) = N1yk(t)+NAdy+k(t)f(t),
m where A = MA −NA, B = M1 −N1. The inner iterations decouple
further M = M −N and M = M −N with additional iterative
[ 1 2 2 2 3 3
processes, such that we result to invert simpler matrices and accelerate

1 the solver process. The multisplitting method use additional a decom-
v position oftheouteriterativeprocesswithparallelalgorithms, basedon
5 thepartitionofunity,suchthatwecouldimprovethesolvermethod.We
9
discuss thedifferent algorithms and present a first experiment based on
4
a DAEsystem.
0
0
.
1 Keywords. Numerical Analysis, Multi-level Waveform relaxation methods,
0
Multi-splitting methods, Differential-Algebraic equations.
6
AMSsubjectclassifications35K28,35K27,35K20,47D60,47H20,65M12,
1
: 65M55.
v
i
X
1 Introduction
r
a
We are motivated to accelerate solver methods for differential-algebraic equa-
tions (DAEs). Many mathematical methods are based on such combination of
differential and algebraicequations,e.g., simulations of the power systems, con-
strained mechanical systems, singular perturbations, see [6] and [2].
We start with respect to assume, that the considered partial differential-
algebraicequations (PDAEs) can be semi-discretizedand written into an initial
value problem of differential algebraic equations (DAEs) in the form:
dy(t)
A +By(t)=f(t), t∈[t0,T],y(t0)=y0, (1)
dt
2
whereA∈Cm×m isasingularandB ∈Cm×m isanon-singularcomplexmatrix
with rank m and f(t):[t0,T]→Cm is a sufficient smooth right hand side.
For solving such problem initial value problems, waveform-relaxation(WR)
methods are developed and investigated by many authors, see [3], [5] and [7].
Themainideaistodecomposethepartitionsoflargesystemsintoiteratively
coupled smaller subsystems and solve such subsystems independently over the
integration intervals (also called time-windows), see [7]. For the performance
of the algorithms, since recent years two-stage strategy is introduced in WR
methods, means a first splitting in blocks for pure parallel splitting. For each
processor, we apply additional a splitting , i.e., an inner splitting instead of a
direct method, see [1].
We propose additional a multistage splitting, means, we assume to split ad-
ditional the inner splitting such that we could also preform the inner splitting
with parallel splitting, while only the last inner splitting is done serial.
Thenewclassofmulti-stagewaveform-relaxation(MSWR) methodsaredis-
cussedinthefollowing,withrespecttothethree-stagewaveform-relaxation(TH-
SWS) method.
For simplification, we deal with
dy(t)
+By(t)=f(t), t∈[t0,T],y(t0)=y0, (2)
dt
where B =M1−N1 and the outer iteration is obtained as
dyk+1(t)
+M1yk+1(t)=N1yk(t)+f(t), (3)
dt
yk+1(t0)=y0, k =1,2,..., (4)
where B =M1−N1, then we apply the first inner iteration:
dzν+1(t)
+M2zν+1(t)=N2zν(t)+N1yk(t)+f(t), (5)
dt
zν+1(t0)=y0k, k=1,2,...,ν =0,...,νk−1, (6)
where M1 =M2−N2 and we obtain yk+1 =zνk.
The last or second inner iteration is given as:
dz˜µ+1(t)
+M3z˜µ+1(t)=N3z˜µ(t)+N2zν(t)+N1yk(t)+f(t), (7)
dt
z˜µ+1(t0)=zν+1(t0)=y0k, (8)
k =1,2,...,ν =0,...,ν −1,µ=0,...,µ −1,
k νk
where M2 = M3−N3 and we obtain yk+1 = zνk =z˜µνk . We have at least µνk
inner first iterations within the inner second iteration ν and within the k+1
k
outer iterations. Means we deal with three-stage iterative methods.
For more flexibility in the approaches, we apply a multisplitting method,
which is based on the partition of the unit, see [4].
3
We decompose the solution into several units, means we have:
L
yk+1 = E yp,k+1, (9)
p
p=1
X
L
E =I, (10)
p
p=1
X
whereI ∈Cm×m istheunitmatrixandE arediagonalandthediagonalentries
p
are given as E ≥0 for i=1,...,m.
p,ii
We can solve L independent waveform relaxation schemes in parallel, while
the synchronization or update is done with the Equation (9).
The paper is outlined as following. In Section 2, we discuss the different
hierarchyofsolvermethods.ThenumericalexperimentsarepresentedinSection
3. The conclusions are done in Section 4.
2 Hierarchy of Solver schemes
In the following, we deal with the following hierarchy of solver schemes:
1. One-stage WR schemes,
2. Two-stage WR schemes,
3. Three-stage WR schemes,
4. Multisplitting WR schemes: Jacobian-, Gauss-Seidel-Types,
where, we simplifiy the inversion of the matrices between the one-stage to the
three-stage method, means the simplification of the inversion is done via addi-
tional inner iteratvie stages. Further the multisplitting approach allows to be
more flexible in the parallelization of a multi-stage method.
2.1 One-stage WR method
In the following, we discuss the one-stage WR method.
1. We have the following WR method (in parallel):
1
(MA+hM1)ynk++11 =h(N1+ hNA)ynk+1+MAynk+1 (11)
−NAynk +hfn+1,
y0k+1 =y0, k =0,1,...,K, n=0,1,2,...,J,
2. We have the following WR method (in serial):
1
(MA+hM1)ynk++11 =h(N1+ hNA)ynk+1+MAyn (12)
−NAyn+hfn+1,
y0 =y , k =0,1,...,K, n=0,1,2,...,J,
n+1 n
4
where we apply the algorithm with p = 50,q = 6,J = 20,h = 0.1, here we
have ∆x=1.0 and D =1.0.
∆x2
We have 2 possible stopping criteria:
(a) Error bound:
We have an stopping error norm with ||yk+1−yk ||≤10−3.
n+1 n+1
(b) Fix number of outer-iterative steps: K =20.
2.2 Two-Stage WR method
The two-stage WR method is given as:
1. We have the following Two-Stage WR method (in parallel):
1
(MA+hM2)znν++11 =hN2znν+1+h(N1+ hNA)ynk+1+MAznν+1 (13)
−NAynk+hfn+1,
z0ν+1(t0)=yk(t0)=y0, k =0,1,...,K,
ν =0,1,...,ν , n=0,1,2,...,J.
k
2. We have the following Two-Stage WR method (in serial):
1
(MA+hM2)znν++11 =hN2znν+1+h(N1+ hNA)ynk+1+MAyn (14)
−NAyn+hfn+1,
z0 (t )=yk , ν =0,1,...,ν ,inner iteration
n+1 n n+1 k
yk+1 =zνk , k =0,1,...,K,outer iteration
n+1 n+1
z0 (t )=y , initialization
n+1 n n
n=0,1,2,...,J,.
Two-Stage WR algorithm (serial) 1 is the given as:
Given the initial vector y0 =y(0) ,
zn0+1(0)=y0,
for n=0,1,...,J do
for k =0,1,...,K do
for ν =0,1,...,ν do
k
(MA+hM2)znν++11 =
hN2znν+1+h(N1+ h1NA)ynk+1+MAyn−NAyn+hfn+1,
end
yk+1 =zνk ,
n+1 n+1
z0 =yk+1,
n+1 n+1
end
yn+1 =znνK+1,
zn0+1(tn+1)=yn+1,
end
Algorithm 1: Two-Stage WR algorithm (serial)
5
where we apply the algorithm with p = 50,q = 6,J = 20,h = 0.1. Here we
have ∆x=1.0 and D =1.0.
∆x2
We have 2 possible stopping criteria:
(a) Error bound:
We have an stopping error norm :
for the outer iteration with ||yk+1−yk ||≤10−3,
n+1 n+1
and for the inner iteration ||zν+1−zν ||≤10−3
n+1 n+1
(b) Fix number of outer-iterative steps: K = 5 and inner iterative steps
ν =4.
k
2.3 Three-Stage WR method
The three-stage WR method is given as:
1. We have the following Three-Stage WR method (in parallel):
1
(MA+hM3)z˜nµ++11 =hN3z˜nµ+1+hN2znν+1+h(N1+ hNA)ynk+1
+MAz˜nµ+1−NAynk+hfn+1,
z˜0µ+1(t0)=z0ν+1(t0)=yk(t0)=y0, k =0,1,...,K,
ν =0,1,...,ν ,µ=0,1,...,µ , n=0,1,2,...,J.
k νk
2. We have the following Two-Stage WR method (in serial):
1
(MA+hM3)z˜nµ++11 =hN3z˜nµ+1+hN2znν+1+h(N1+ hNA)ynk+1
+MAyn−NAyn+hfn+1,
z˜0 (t )=zνk , µ=0,1,...,µ ,first inner iteration
n+1 n n+1 νk
z0 (t )=yk , ν =0,1,...,ν ,second inner iteration
n+1 n n+1 k
yk+1 =z˜µνk, k =0,1,...,K,outer iteration
n+1 n+1
z˜0 (t )=y , initialization
n+1 n n
n=0,1,2,...,J,.
6
Three-Stage WR algorithm (serial) 2 is the given as:
Given the initial vector y0 =y(0) ,
zn0+1(0)=y0,
for n=0,1,...,J do
for k =0,1,...,K do
for ν =0,1,...,ν do
k
for µ=0,1,...,µ do
νk
(MA+hM3)z˜nµ++11 =hN3z˜nµ+1+hN2znν+1+h(N1+
h1NA)ynk+1+MAyn−NAyn+hfn+1,
end
zν+1 =z˜µν ,
n+1 n+1
z˜0 =zν+1,
n+1 n+1
end
yk+1 =zνk ,
n+1 n+1
z0 =yk+1,
n+1 n+1
end
yn+1 =znνK+1,
zn0+1(tn+1)=yn+1,
end
Algorithm 2: Three-Stage WR algorithm (serial)
where we apply the algorithm with p = 50,q = 6,J = 20,h = 0.1. Further
we have ∆x=1.0 and D =1.0.
∆x2
We have 2 possible stopping criteria:
(a) Error bound:
We have an stopping error norm :
for the outer iteration with ||yk+1−yk ||≤10−3,
n+1 n+1
and for the inner iteration ||zν+1−zν ||≤10−3
n+1 n+1
and for the second inner iteration ||z˜µ+1−z˜µ ||≤10−3
n+1 n+1
(b) Fix number of outer-iterative steps: K = 5 and inner iterative steps
ν =2, µ =2.
k νk
2.4 Multisplitting WR method
We have the following Multi-splitting WR method (in serial/parallel):
L
1
(MAl +hM1,l)ynl,+k+11 =h(N1,l+ hNAl) El,mynl,+k1 +MAyn (15)
!
m=1
X
−NAyn+hfn+1,
yl,0 =y , l =1,...,L, k =0,1,...,K, n=0,1,2,...,J,
n+1 n
where we apply the algorithm with p =50,q = 6,J =20,h=0.1 and the error
norm ||yk+1−yk ||≤10−3, here we have ∆x=1.0 and D =1.0. Further L
n+1 n+1 ∆x2
are the number of the processors.
7
Withoutloosingthegeneralityofthemethod,weconcentrateonthefollowing
to L=2.
Jacobian-Method The first processor is computing:
1
(MA1 +hM1,1)yn1,+k1+1 =h(N1,1+ hNA1) E1,1yn1,+k1+E1,2yn2,+k1 +MAyn
−NAyn+hfn+1, (cid:16) (cid:17) (16)
y1,0 =y , l=1,...,L, k =0,1,...,K, n=0,1,2,...,J,
n+1 n
The second processor is computing:
1
(MA2 +hM1,2)yn2,+k1+1 =h(N1,2+ hNA2) E2,1yn1,+k1+E2,2yn2,+k1 +MAyn
−NAyn+hfn+1, (cid:16) (cid:17) (17)
y2,0 =y , l=1,...,L, k =0,1,...,K, n=0,1,2,...,J.
n+1 n
where we decide if we have to switch of the mixing means: means if we have
fulfilled:
|| E1,1yn1,+k1+E1,2yn2,+k1 −yn1,+k−11||≤||yn1,+k1−yn1,+k1−1|| (18)
||(cid:16)E2,1yn1,+k1+E2,2yn2,+k1(cid:17)−yn2,+k−11||≤||yn2,+k1−yn2,+k1−1|| (19)
(cid:16) (cid:17)
we do not switch of the mixing, but if the mixing has a larger error we have:
yk =y1,k , (20)
n+1 n+1
Remark 1. The multisplitting is switched off, if one partial solution is much
more accurate, than the other partial solution. Then we only apply the best
approximation.
Gauss-Seidel-Method (decoupled version, serial with 2 processors) In
this version, we apply the wel-known standard Gauss-Seidel method, which has
the drawback of the serial treatment with the results.
The first processor is computing:
1
(MA1 +hM1,1)yn1,+k1+1 =h(N1,1+ hNA1) E1,1yn1,+k1+E1,2yn2,+k1 +MAyn
−NAyn+hfn+1, (cid:16) (cid:17) (21)
y1,0 =y , l=1,...,L, k =0,1,...,K, n=0,1,2,...,J,
n+1 n
Remark 2. Here, we can apply the result of the first processor, if we assume,
that part is much more faster to the second processor.
8
The second processor is computing:
1
(MA2 +hM1,2)yn2,+k1+1 =h(N1,2+ hNA2) E2,1yn1,+k1+1+E2,2yn2,+k1 +MAyn
−NAyn+hfn+1, (cid:16) (cid:17) (22)
y2,0 =y , l=1,...,L, k =0,1,...,K, n=0,1,2,...,J.
n+1 n
where we decide if we have to switch of the mixing means: means if we have
fulfilled:
|| E1,1yn1,+k1+E1,2yn2,+k1 −yn1,+k1−1||≤||yn1,+k1−yn1,+k−11|| (23)
||(cid:16)E2,1yn1,+k+11+E2,2yn2,+k(cid:17)1 −yn2,+k−11||≤||yn2,+k1−yn2,+k1−1|| (24)
(cid:16) (cid:17)
we do not switch of the mixing, but if the mixing has a larger error we have:
yk =y1,k , (25)
n+1 n+1
Remark 3. Themultisplittingisswitchedoff,ifonepartialsolutionismuchmore
accurate,thantheotherpartialsolution.Otherwise,weapplythemixtureofthe
results based on the multisplitting method.
Gauss-Seidel-Method (decoupled version) The first processors compute:
1
(MA1 +hM1,1)−h(N1,1+ hNA1)E1,1 yn1,+k1+1
(cid:18) (cid:19)
1
=h(N1,1+ hNA1) E1,2yn2,+k1 +MAyn−NAyn+hfn+1, (26)
y1,0 =y , l=1,..(cid:16).,L, k =0(cid:17),1,...,K, n=0,1,2,...,J,
n+1 n
The second processors compute:
1
(MA2 +hM1,2)−h(N1,2+ hNA2)E2,2 yn2,+k1+1
(cid:18) (cid:19)
1
=h(N1,2+ hNA2) E2,1yn1,+k1 +MAyn−NAyn+hfn+1, (27)
y2,0 =y , l=1,..(cid:16).,L, k =0(cid:17),1,...,K, n=0,1,2,...,J.
n+1 n
Gauss-Seidel-Method(coupledversion) Here,weapplythecoupledversion
of the GS method, which means one processor is faster with the computation
and the other processor can profit from the improved computations.
9
The processors compute:
1
(MA1 +hM1,1)−h(N1,1+ hNA1)E1,1 yn1,+k1+1
(cid:18) (cid:19)
1
=h(N1,1+ hNA1) E1,2y˜n2,+k1+1 +MAyn−NAyn+hfn+1, (28)
(cid:16) 1(cid:17)
(MA2 +hM1,2)−h(N1,2+ hNA2)E2,2 yn2,+k1+1
(cid:18) (cid:19)
1
=h(N1,2+ hNA2) E2,1y˜n1,+k1+1 +MAyn−NAyn+hfn+1, (29)
y1,0 =y , l =1,..(cid:16).,L, k =0,(cid:17)1,...,K, n=0,1,2,...,J,
n+1 n
y2,0 =y , l =1,...,L, k =0,1,...,K, n=0,1,2,...,J.
n+1 n
where we have two cases:
– If Processor 1 is faster than Processor 2:
y˜1,k+1 =y1,k+1, (30)
n+1 n+1
y˜2,k+1 =y2,k , (31)
n+1 n+1
– If Processor 2 is faster than Processor 1:
y˜1,k+1 =y1,k , (32)
n+1 n+1
y˜2,k+1 =y2,k+1. (33)
n+1 n+1
Remark 4. For all multisplitting methods, we can also extend the one-stage
waveform-relaxationmethod to a multi-stage waveform-relaxationmethod.
3 Numerical Experiments
Inafirstexperiment,weapplyapartialdifferentialalgebraicequation(PDAE),
which combines partial differential and algebraic equations.
We choose an experiment, which is based on the following two equations:
∂tc1+∇·Fc1 =f1(t), inΩ×[0,t], (34)
∇·Fc2 =f2(t), inΩ×[0,t], (35)
F=−D∇, (36)
and we have the following DAE problem:
A∂ c+Bc=f(t), in[0,t], (37)
t
The analytical solution is given as:
y =[cos(t),sin(t),t,cos(t),sin(t),t,...,cos(t),sin(t),t]∈IRm, (38)
10
and we have to calculate f(t) as:
f(t)=A∂ y(t)+By(t), in[0,t], (39)
t
where y is the analytical solution.
We apply the error in L2 or Lmax-norm means:
I
1
err (t)= ( (y (x ,t)−y (x ,t))2)1/2, (40)
L2 ∆x ana i num i
i=1
X
I
err (t)=max||y (x ,t)−y (x ,t)||. (41)
max ana i num i
i=1
In the following we deal with the semidiscretized equation given with the
matrices:
I
 ... 
A= ∈IRm×m, (42)
I
 
 0 
 
 0
 
 
where I,0∈IRp×p
We have the following two operators for the splitting method:
4−1
−1 4−1
 
B1 = ... ... ... ∈IRp×p (43)
 
 −1 4−1
 
 −1 4
 
 
B1 −I
−I B1 −I
D  
B = · ... ... ... ∈IRm×m (44)
∆x2  
 −I B1 −I
 
 −I B1
 
 
with pq =m, where we assume D =1.
∆x2
Means A,B are m×m block-matrices.
We have the following splitting:
0
..
 . 
N = , (45)
A 0
 
 1 I 
 100 
 0
 
 