Available online at www.sciencedirect.com
Mathematics and Computers in Simulation 77 (2008) 179–188
Superlinear convergence of asynchronous multi-splitting waveform relaxation methods applied to a system of nonlinear ordinary differential equations M. El-Kyal a,b,∗ , A. Machmoum b,a b
a LIP2E, ENSA, Universit´ e Ibn Zohr, B.P. 1136 Agadir, Morocco LAMA, Facult´e des sciences, Universit´e Ibn Zohr, B.P. 8106 Agadir, Morocco
Available online 31 August 2007
Abstract We prove the superlinear convergence of asynchronous multi-splitting waveform relaxation (MSWR) methods applied to a system of nonlinear ordinary differential equations. This study is based on the technique of nested sets. It allows to specify the class of the convergence. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Numerical analysis; Waveform relaxation methods; Multi-splitting; Asynchronous algorithms
1. Introduction Multi-splitting (MS) algorithms were introduced by O’Leary and White [8] to solve in parallel linear systems of the form Ax = b where A is a nonsingular matrix. Several splitting A = Ml − Nl , l = 1, . . . , L, where Ml is nonsingular, are used together with a set of diagonal nonnegative (weighting) matrices El such that they add up to the identity. The number L is called number of splitting. (p+1) (p) The resulting algorithm of multi-splitting consists in calculating L x l according to x by solving L systems l (p) (p+1) Ml y = Nl x + b, l = 1, . . . , L, and while imposing x = l=1 El y . Waveform relaxation (WR) is a technique to solve large system of initial-value problems (IVP’s) consisting of ODE’s in parallel. The right hand side of the system is split into subsystems. One then solves iteratively all the subsystems in parallel and exchanges the information. Multi-splitting waveform relaxation (MSWR) method is a generalization of the method of waveform relaxation by using techniques of multi-splitting in the sense that the multi-splitting allows overlap between subsystems. This technique was introduced by Frommer and Pohl [5] in the case of linear ordinary differential equations. Our goal is to generalize this formulation to the nonlinear case and to show its superlinear convergence in the asynchronous mode. The main characteristic of an asynchronous mode is that the local algorithm not have to wait at predetermined messages to become available. We allow some processors to communicate more frequently than others, and we allow the communication delays to be substantial and unpredictable. Note that synchronous algorithms in the computer science sense are particular cases of our formulation of asynchronous one. Our goal is to prove the superlinear convergence of ∗
Corresponding author. E-mail address:
[email protected] (M. El-Kyal).
0378-4754/$32.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2007.08.012
180
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
AMSWR methods when one applies them to the nonlinear differential systems, this study is based on the technique of nested sets. It permits to specify kind of the convergence using the theory introduced by Bertsekas and Tsitsiklis [2]. The paper is organized as follows: Section 2 is devoted to succinct background about asynchronous algorithms formulation and asynchronous convergence theorem of Bertsekas and Tsitsiklis with a characterization of superlinear convergence. In Section 3, we first give the definition of multi-splitting waveform relaxation methods applied to fixed point problems. Then, we prove the superlinear convergence of these algorithms for ODE’s. 2. Preliminaries Let Bi , i ∈ {1, . . . , α} be Banach spaces equipped with the norms .i and B = αi=1 Bi equipped with the norm . = max1≤i≤α (.i /γi ), γi > 0 ∀i ∈ {1, . . . , α}. We consider a mapping T from D(T ) (definition domain of T) in B such that T (x) = (T1 (x), . . . , Tα (x)). One calls fixed point problem the problem of the form: find x∗ such that x∗ = T (x∗ ). 2.1. Asynchronous algorithms formulation Define J = {J(p)}p ∈ N a sequence of nonempty subsets of {1, . . . , α} such that: ∀i ∈ {1, . . . , α},
{p ∈ N, i ∈ J(p)} is infinite.
(1)
and S = {s1 (p), . . . , sα (p)}p ∈ N a sequence of Nα such that: ∀i ∈ {1, . . . , α},
∀p ∈ N, si (p) ≤ p.
(2)
and ∀i ∈ {1, . . . , α},
lim si (p) = +∞.
p→+∞
The asynchronous algorithm denoted (T, x0 , J, S) is defined by, ⎧ Given x0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ for l = 1, . . . , α, for p ≥ 1, ⎪ ⎪ ⎪ Tl (x1s1 (p) , . . . , xαsα (p) ), if l ∈ J(p), ⎪ p+1 ⎪ ⎪ = ⎩ xl p xl , if l ∈ / J(p).
(3)
(4)
- J(p) is the set of components to be updated at step p. - si (p) is the delay due to ith processor when it computes the ith block at the pth iteration. 2.2. Superlinear convergence of asynchronous algorithms Theorem 1. [2] Let {Xp }p ∈ N be a sequence of nonempty closed sets satisfying: (i) Nested sets condition: . . . ⊂ Xp+1 ⊂ Xp ⊂ . . . X0 ⊂ B, (ii) Synchronous convergence condition: T (x) ∈ Xp+1 ∀x ∈ Xp and ∀p ∈ N, Furthermore, if {xp }p ∈ N is a sequence satisfying xp ∈ Xp for any p, then every limit point of {xp }p ∈ N is fixed point of T, p p p (iii) Box condition: for any p ∈ N, there exists α sets Xi such that Xp is their Cartesian product: Xp = X1 × · · · × Xα . If x0 = (x10 , . . . , xα0 ) ∈ X0 , then every asynchronous algorithm (T, x0 , J, S) converges to fixed point of T. To define the superlinear convergence of asynchronous algorithms, denote by p the diameter of Xp : p = supx,y ∈ Xp x − y.
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
181
Definition 2. Let us write that the synchronous convergence property of the successive approximation iterations associated to T and initialised by x0 ∈ X0 is superlinear if limsupp→+∞ (p )1/p = 0. Let us consider a sequence of delays S to which corresponds {s(p)}p ∈ N ⊂ N defined by: s(p) = min sl (p).
(5)
1≤l≤α
Under hypothesis (2), (3) and (5), we have lim s(p) = +∞
p→+∞
and ∀p ∈ N
s(p) ≤ p.
Definition 3. Let us consider J and S satisfying (1)–(3) and the sequence {s(p)}p ∈ N associated to S defined by (5). Let us consider the sequence of natural integers {pl }l ∈ N strictly increasing defined in the following way: (1) p0 = 0, (2) p1 is the smallest natural integer such that: 0≤s(p)≤p
3. Asynchronous multi-splitting waveform relaxation method for nonlinear ODE’s Let us consider the system of ordinary differential equations: x (t) = F (x(t), t) x(0) = x0
(6)
with t ∈ [0, tmax ] and F = (F1 , . . . , Fα )t : Rα × [0, tmax ] → Rα Let us assume that: x∈C =
i=α
C1 ([0, tmax ], R) i=1
with for i ∈ {1, . . . , α} and t ∈ [0, tmax ], |Fi (x(t), t) − Fi (¯x(t), t)| ≤ L1 |x(t) − x¯ (t)|∞ where |.|∞ is the norm defined by |x(t)|∞ = max1≤i≤α |xi (t)|.
(7)
182
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
For any i ∈ {1, . . . , α}, the space C1 ([0, tmax ], R) is endowed with the norm: |xi |K = maxt ∈ [0,tmax ] e−Kt |xi (t)| where K satisfies 1 − e−Ktmax 1 . < K 2L1
(8)
The space C is endowed with the norm xK = max1≤i≤α |xi |K . We write: |xi |K,t = sup0≤τ≤t e−Kτ |xi (τ)| and xK,t = max1≤i≤α |xi |K,t . To problem (6), we associate the fixed point mapping T as follows: T : D(T ) → C u→x
(9)
such that for i = 1, . . . , α, xi is the solution of the system xi (t) = Fi (u1 (t), . . . , ui−1 (t), xi (t), ui+1 (t), . . . , un (t), t), xi (0) = x0i
∀t ≥ 0
Proposition 7. [3] Under hypothesis (7) and (8), T is .K contractive of constant β and its fixed point x∗ coincides with the solution of the problem (6). Remark 8. The resolution of our problem is reduced then to the resolution of a fixed point problem. In all the sequel and for any p ∈ N, ρp is the function defined on [0, tmax ] with values in R by: p p ⎧ KL1 t ⎨ , if t < tp K − L1 p! ρp (t) = ⎩ 1, if t ≥ tp where tp =
K−L1 1/p . KL1 (p!)
Lemma 9. (Technical lemma [3]) (1) (2) (3) (4)
KL1 t For t < tp , K−L < 1, 1 p For p = 0, 1, . . . tp+1 ≥ tp , limp→+∞ tp = +∞, There exists a positive integer p0 such that:
∀p ≥ p0 ,
ρp (t) =
KL1 K − L1
p
tp p!
∀t ∈ [0, tmax ].
3.1. Superlinear convergence property of the asynchronous multi-splitting waveform relaxation method Let L be a fixed integer which will be called the number of splitting. A multi-splitting of the problem is a collection of L fixed point mappings T l which fixed point coincide with the solution of the problem (6) and L diagonal nonnegative (weighting) matrix El such that they add up to the identity. The resulting algorithm of multi-splitting is the extended fixed point mapping defined by: T¯ : X0 → CL U = (u1 , . . . , uL ) → X = (x1 , . . . , xL )
L m . where for all l ∈ {1, . . . , L}, xl = T l m=1 Em u
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
183
To build these mappings, let Il , l ∈ {1, . . . , L}, some nondisjoined subsets from {1, . . . , α} such that L l=1 Il = {1, . . . , n} and Ilc their complementary. For v and u in C and l ∈ {1, . . . , L}, we define the vector σil by σil (v, u) = w and σil (v, u)(t) = w(t) such that wj = vj , if i, j ∈ Il or i, j ∈ Ilc ∀j ∈ {1, . . . , α} wj = uj , otherwise. Definition 10. For all l ∈ {1, . . . , L}, we define an implicit mapping T l associated to T by: ∀u ∈ C,
Til (u) = Ti (σil (T l (u), u)).
∀i ∈ {1, . . . , α},
(10)
Proposition 11. [3] For all l in {1, . . . , L}, the fixed point mapping T l is .K contractive, its constant βl is less than or equal to the constant β of T and its fixed point x∗ coincides with that of T. Definition 12. For all p ∈ N, let us define the truncated fixed point mapping T¯ p by T¯ p : D(T¯ p ) ⊂ CL → CL
(11)
(u1 , . . . , uL ) → (x1 , . . . , xL ) with for all l ∈ {1, . . . , L} ⎧ ⎪ xl (t) = ul (t) ,if t ≥ tp+1 ⎪ ⎨ ⎪ xl (t) ⎪ ⎩
=
vl (t)
,if t < tp+1
with vl (t)
=
Tl
L
Em u
m
m=1
where for all m ∈ {1, . . . , L}, Em is a matrix which satisfies: ⎧ Em is diagonal ⎪ ⎪ ⎪ ⎪ ⎨ Em ≥ 0
(12)
L ⎪ ⎪ ⎪ Em = I (identity matrix in Rα×α ). ⎪ ⎩ m=1
Define a sequence of nonempty sets (Xp )p ∈ N in CL where Xp = p Thus for all l ∈ {1, . . . , L}, the set Xl is defined as follows: p
χl = {xl ∈ C/||xl − x ||K,t ≤ ρp (t)||xl − x ||K
l=L
p l=1 Xl .
∀t ∈ [0, tmax ] and xl (0) = x0 }.
Lemma 13. ∀p ∈ N, Xp+1 ⊂ Xp . Proof. It will be sufficient to prove that for all t ∈ [0, tmax ], ρp+1 (t) ≤ ρp (t). Let t ≥ 0, (1) if t < tp , then
ρp+1 (t) = =
KL1 K − L1
p+1
p+1 ˜ L ˜ ptp t p+1 Lt ˜ p+1 t =L = (p + 1)! (p + 1)! p + 1 p!
˜ ˜ Lt Lt ρp (t) ≤ ρp (t) ≤ ρp (t)(the first part of Lemma 9), p+1 p
184
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
(2) if tp ≤ t < tp+1 , then p+1 ˜ p+1 t ρp+1 (t) = L (p + 1)!
and ρp (t) = 1.
However t < tp+1 , we have (Lemma 9) p+1 ˜ p+1 t L <1 (p + 1)!
then ρp+1 (t) ≤ ρp (t), (3) if t ≥ tp+1 , then ρp+1 (t) = 1 = ρp (t). then for all t ≥ 0, ρp+1 (t) ≤ ρp (t).
/ Xp+1 , this is equivalent Let U ∈ X0 and let p be the greatest integer such that U ∈ Xp . With this definition of p, U ∈ to p
U ∈X X
p+1
p
where X X
p+1
=
U ∈ Xp U∈ / Xp+1
.
(13)
Taking into account the mappings defined in (10) and (11) and the condition (13), let us define the fixed point mapping T¯ associated to the problem (6) by: T¯ : X0 → CL U = (u1 , . . . , uL ) → X = (x1 , . . . , xL ) where X = T¯ p (U) if U ∈ Xp Xp+1 .
Lemma 14. The synchronous iterations associated to T¯ and initialised by U 0 ∈ X0 converge superlinearly. Proof. (I) We prove that for all p = 0, 1, . . ., T¯ (Xp ) ⊂ Xp+1 . Let U = (u1 , . . . , uL ) ∈ Xp and X = (x1 , . . . , xL ) = T¯ (U). Two cases arise:
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
(a) t < tp+1 ∀l ∈ {1, . . . , L}, then
l
x =T
l
L
⎧ ⎪ xl (t) = v(t) ⎪ ⎨ ⎪ ⎪ ⎩ v(t) =
Tl
L
185
Em u
m
m=1
Em u
m
m=1
then
xil = Ti
∀i ∈ {1, . . . , α}, then we have xil (t) = x0,i + then |xil (t) − xi (t)|
0
t
Fi
σil
σil
xl ,
Tl
L
0
t
E m um
,
m=1
Elm um
L
E m um
m=1
= Ti
σil
(s), s
Em u m
m=1
ds
L t l l m (s), s − Fi (x (s), s) ds = F i σi x , Em u 0 m=1 t L l l ≤ Em um (s), s − Fi (x (s), s) ds F i σi x , 0 m=1 t L L l l m l m (s), s − Fi σi x , (s), s ds Em u Em u ≤ F i σi x , 0 m=1 m=1 t L l m + (s), s − Fi (x (s), s) ds Em u F i σi x , 0 m=1 t t L m L1 Em u (s) − x (s) ds + L1 |xl (s) − x (s)|∞ ds. ≤ 0 0 ∞
∞
L L t m L1 Em (u (s) − x (s)) ds Em = Iα×α = 0 m=1 m=1 ∞ ⎛ ⎞ L t j=α m ⎝ ⎠ L1 max1≤i≤α (Em )ij (uj (s) − xj (s)) ds = 0 m=1 j=1 t L m = L1 max1≤i≤α ((Em )ii (ui (s) − xi (s))) ds (Em is diagonal) 0
xl ,
L
L L1 Em um (s) − x (s) ds m=1
m=1
m=1
But
L
m=1
186
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
t
≤
L1
0
m=1 t
≤
L
0
(Em )ii max1≤i≤α (max1≤m≤L |(um i (s) − xi (s))|) ds
L1 max1≤i≤α (max1≤m≤L |(um i (s) − xi (s))|) ds
then −Kt
e
|xil (t) − xi (t)|
t
≤ 0
L
Em = Iα×α
,
m=1
L1 e−K(t−s) max1≤i≤α (max1≤m≤L | e−Ks (um i (s) − xi (s))|) ds
+
t
0
L1 e−K(t−s) max1≤i≤α (|e−Ks (xil (s) − xi (s))|) ds
k k KL1 1 − e−Kt s ds · x0 − x K +L1 xl − x K,t K − L k! K 1 0 p t p KL1 s L1 l ds · x0 − x K + x − x K,t , ≤ L1 K − L1 K 0 p!
≤ L1
t
e−K(t−s)
then max1≤i≤α e Thus
K − L1 K
−Kt
|xil (t) − xi (t)|
≤ L1
KL1 K − L1
p
t
0
sp L1 l ds · x0 − x K + x − x K,t . p! K
p p+1 KL1 t .x0 − x K K − L1 (p + 1)! p+1 p+1 t KL1 ≤ .x0 − x K . K − L1 (p + 1)!
xl − x K,t ≤ L1
xl − x K,t
i.e.
(b) t ≥ tp+1 . for t ≥ tp+1 and i = 1, . . . , α, we have xil (t) = uli (t) ⇒ e−Kt |xil (t) − xi (t)| = e−Kt |uli (t) − xi (t)| then |xil − xi |K,t = sup e−Kr |xil (r) − xi (r)| = max( sup 0≤r≤t
0≤r
e−Kr |xil (r) − xi (r)|,
−xi (r)|) ≤ max(xil − xi ,
sup
sup
tp+1 ≤r≤t
tp+1 ≤r≤t
e−Kr |uli (r)
e−Kr |uli (r) − xi (r)|).
But sup
tp+1 ≤r≤t
e−Kr |uli (r) − xi (r)| ≤ |uli − xi |K,t ≤ ρp (t)x0 − x K = ρp+1 (t)x0 − x K
moreover sup 0≤r
e−Kr |xil (r) − xi (r)| ≤ sup ( sup e−Kr |xil (r) − xi (r)|) ≤ sup |xil − xi |K,δ δ
δ
≤ sup ρp+1 (δ)x0 − x K ≤ sup |xil − xi |K,δ ≤ x0 − x K δ
δ
(because ρp+1 (δ) ≤ 1).
Then |xil (r) − xi (r)|K,t ≤ max(x0 − x K , ρp (t)x0 − x K ) ≤ x0 − x K = ρp+1 (t)x0 − x K . p+1
In the two cases, we proved that for i = 1, . . . , α, xi ∈ Xi
.
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
p+1 Finally, x ∈ i=α . i=1 Xi Consequently, ∀l ∈ {1, . . . , L}, (II) Let p the diameter of Xp . p ≤ 2
sup
t ∈ [0,tmax ]
p+1
x l ∈ Xl
187
.
ρp (t)x0 − x K .
According to Lemma 9, there exists p0 such that p ≥ p0 ⇒ ρp (t) =
KL1 K − L1
p
tp . p!
Then p ≤ 2
KL1 K − L1
p
p
tmax 0 x − x K p!
therefore limsupl→+∞ (p )2
1/p
KL1 K − L1
tmax (p!)
1 p
x0 − x K = 0 1/p
then limsupl→+∞ (p )1/p = 0.
Theorem 15. Any asynchronous algorithm (T¯ , X0 , J, S) a ssociated to T¯ converges superlinearly towards X fixed point of T¯ where X = (x , . . . , x ). L
times
Proof. (1) By Lemma 13, we have proved that ∀p ∈ N, Xp+1 ⊂ Xp . p p (2) By construction, ∀p ∈ N, Xp = X1 × · · · × XL . (3) By Lemma 14, the synchronous successive approximation iterations associated to T¯ and initialised by X0 ∈ X0 are superlinearly convergent. Then by Definition 2 and Lemma 5, we are sure that the asynchronous algorithm (T¯ , X0 , J, S) associated with T converges superlinearly towards the solution of the problem (6). References [2] D.P. Bertsekas, J.N. Tsitsiklis, Parallel and Distributed Computation, Prentice Hall, Englewood Cliffs, New Jersey, 1989. [3] M. El-Kyal, Les m´ethodes parall`eles asynchrones de multi-d´ecomposition appliqu´ees aux syst`emes diff´erentiels alg´ebriques d’indice inf´erieur ou e´ gal a` 1. Th`ese de l’universit´e de Franche-Comt´e de Besanc¸on, 1999.
188
M. El-Kyal, A. Machmoum / Mathematics and Computers in Simulation 77 (2008) 179–188
[5] A. Frommer, B. Pohl, A comparison result for multisplittings and waveform relaxation methods, Numer. Linear Algebra Appl. 2 (4) (1995) 335–346. [7] J.-C. Miellou, P. Cortey-Dumont, M. Boulbrachˆene, Perturbation of Fixed Point Iteration Methods, Adv. Para. Comput. 1 (1990) 81–122. [8] D.P. O’Leary, R.E. White, Multi-splittings of matrices and parallel solution of linear systems, SIAM J. Alg. Disc. Math. 6 (1985) 630–640.