WaveForm Relaxation Methods Using Low-order Composite Quadrature Schemes

Numerical ODEs-Final Project
Team : Ahmad + Mohammad

Introduction

Consider the following system of ODEs:  
	y'=f(y),     IC: y(0)=(Subscript[y, 0])
Picard-Lindelof iteration is given by:
	(y^(n+1))(t)=y(0)+(∫_0^t)f((y^n)(τ))dτ
From ODE Theory: under reasonable assumptions on f, Picard iteration is guaranteed to converge;
	(y^n)(t)-> y(t)  as  n-> ∞       (however convergence can be slow).
	
Example 1.  Consider the 1D IVP:  
	y'=λ y,             y(0)=1.
Analytic solution is
	y(t)=(e^(λt))=(∑_(k=0)^∞)((((λ t))^k)/(k!)).
Picard iteration with (y^0)=(Subscript[y, 0])=1 gives
	(y^1)(t)=1+(∫_0^t)λdτ=1+λ t
	(y^2)(t)=1+(∫_0^t)λ(1+λ t)dτ=1+λ t+((((λ t))^2)/(2!))
	         ⋮
	(y^n)(t)=(∑_(k=0)^n)((((λ t))^k)/(k!))
Obviously
	(y^n)(t)-> y(t)  as  n-> ∞,
but very slowly if |λ| is large:

λ=10;
LogPlot[Abs[(E^(λ t))-(∑_(k=0)^200)((((λ t))^k)/(k!))],{t,0,5},PlotRange->All,AxesLabel-> {"t","Error"}]

Waveform Relaxation

Waveform relaxation (or dynamic iteration method) is a highly parallel iterative method for solving very large systems of ordinary differential equations. Over the years this method has been applied almost exclusively for the systems that model VLSI (Very Large Scale Integrated) electronic circuits. In this project, we will examine a waveform relaxation method based on Picard-Lindelӧf iteration. The idea is to approximate the integral above rather than computing it exactly. We will focus on (low-order) quadrature rules to approximate the integral.

Numerical Integration

Suppose that the interval [a,b] is is subdivided into r subintervals of equal width h, and set (Subscript[t, i])=a+i h, for i=0,1,...,r, and (Subscript[t, i+1/2])=(1/2)((Subscript[t, i]) + (Subscript[t, i+1])) for i=0,1,...,r-1.

1.  Composite Trapezoidal Rule takes the form

(∫_a^b)f(t)dt=(h/2)[f((Subscript[t, 0]))+2f((Subscript[t, 1]))+···+2f((Subscript[t, r-1]))+f((Subscript[t, r]))]+O((h^2))

2.  Composite MidPoint Rule takes the form

(∫_a^b)f(t)dt=h[f((Subscript[t, 0+1/2]))+f((Subscript[t, 1+1/2]))+···+f((Subscript[t, r-1+1/2]))]+O((h^2))

3.  Composite Simpson Rule takes the form

(∫_a^b)f(t)dt=(h/6)[f((Subscript[t, 0]))+4f((Subscript[t, 1]))+2f((Subscript[t, 2]))+4f((Subscript[t, 3]))+···+2f((Subscript[t, r-2]))+4f((Subscript[t, r-1]))+f((Subscript[t, r]))]+O((h^4))

Waveform with Low-Order Rules

Let (Subscript[y, i]^n)≈(y^n)((Subscript[t, i])),  (Subscript[y, i+1/2]^n)≈(y^n)((Subscript[t, i+1/2])), and set
	(Piecewise[{{(Subscript[Y, 0])=(([(Subscript[y, 0]),(Subscript[y, 0]),...,(Subscript[y, 0])])^T),  }, {(Y^n)=(([(Subscript[y, 0]^n),(Subscript[y, 1]^n),...,(Subscript[y, r]^n)])^T),  }, {(Y^(n+1/2))=(([(Subscript[y, 1/2]^n),(Subscript[y, 3/2]^n),...,(Subscript[y, r-1/2]^n)])^T),  }, {(F^n)=(([f((Y^n))])^T),  }, {(F^(n+1/2))=(([f((Y^(n+1/2))),0])^T),  }}])
then resulted WR iteration reads as

	(Y^((n + 1))) = (Subscript[Y, 0]) + α ( A (F^n) + B (F^((n + 1/2))) )

where:
1. If the trapezoidal rule is used: then α=h/2,B=0, and
	A=(Subscript[(({{0, 0 , … ,  ,  , 0}, {1, 1, ⋱,  ,  , ⋮}, {1, 2, 1,  ,  ,  }, {1, 2, 2, 1,  ⋱,  }, {⋮,  ,  , ⋱, ⋱, 0}, {1, 2, …,  , 2, 1}})), (r+1)×(r+1)])
2. If the midpoint rule is used: then α=h, A=0, and
	B=(Subscript[(({{0, …,  ,  ,  , 0}, {1, ⋱,  ,  ,  , ⋮}, {1, 1,  ,  ,  ,  }, {1, 1, 1,  ,  ,  }, {⋮,  ,  , ⋱,  ,  }, {1, …,  ,  , 1 , 0}})), (r+1)×(r+1)])
3.  If the simpson rule is used: then α=h/6, and
	A=(Subscript[(({{0, 0 ,  …,  ,  , 0}, {1, 1, ⋱,  ,  , ⋮}, {1, 2, 1,  ,  ,  }, {1, 2, 2, 1,  ⋱,  }, {⋮,  ,  , ⋱, ⋱, 0}, {1, 2, …,  , 2, 1}})), (r+1)×(r+1)]),          B=4 (Subscript[(({{0, …,  ,  ,  , 0}, {1, ⋱,  ,  ,  , ⋮}, {1, 1,  ,  ,  ,  }, {1, 1, 1,  ,  ,  }, {⋮,  ,  , ⋱,  ,  }, {1, …,  ,  , 1 , 0}})), (r+1)×(r+1)])

Remark 1.  Example 1 above suggests that (Subscript[y, i]^n)->(Subscript[y, i]^∞) faster as i decreases.
Remark 2.  We have used quadratic interpolation to approximate (y^n)((Subscript[t, i+1/2])) in Midpoint and Simpson Rules as follows: let (Subscript[y, I]) be the quadratic interpolant of y over the interval [(Subscript[t, i]),(Subscript[t, i+1])] determined by the conditions (Subscript[y, I])((Subscript[t, i]))=y((Subscript[t, i])), (Subscript[y, I])((Subscript[t, i+1]))=y((Subscript[t, i+1])), (Subscript[y, I]^('))((Subscript[t, i]))=y'((Subscript[t, i])), and (Subscript[y, I]^('))((Subscript[t, i+1]))=y'((Subscript[t, i+1])). Thus we can write
 	(Subscript[y, I])(t)=a+b (t-(Subscript[t, i]))+c (t-(Subscript[t, i]))(t-(Subscript[t, i+1])), 
 where a, b, and c are constants to be determined. Imposinge above conditions yield
 	a=y((Subscript[t, i])) ,
 	b=(y((Subscript[t, i+1]))-y((Subscript[t, i])))/h ,
 	c=((f(y((Subscript[t, i+1])))-f(y((Subscript[t, i]))))/(2h)),
where h=((Subscript[t, i+1])-(Subscript[t, i])). Consequently, we have
	(Subscript[y, I])(t)=y((Subscript[t, i]))+((y((Subscript[t, i+1]))-y((Subscript[t, i])))/h) (t-(Subscript[t, i]))+((f(y((Subscript[t, i+1])))-f(y((Subscript[t, i]))))/(2h)) (t-(Subscript[t, i]))(t-(Subscript[t, i+1])),
and so at the midpoint (Subscript[t, i+1/2])=(1/2)((Subscript[t, i])+(Subscript[t, i+1])) we have the approximation (after simplification)
	 (y^n)((Subscript[t, i+1/2]))≈(Subscript[y, I])((Subscript[t, i+1/2]))=(1/2)(y((Subscript[t, i]))+y((Subscript[t, i+1])))+(h/8)(f(y((Subscript[t, i+1])))-f(y((Subscript[t, i])))).

Fixed Window Waveform Scheme

In the fixed window waveform approach, we iterate over the entire time window. 
• We adopt the following criteria to measure convegence:
	iterate until  ‖(Y^(n+1))-(Y^n)‖<IterTol    ...(1)
for a given iteration tolerance, IterTol.

Dynamic Windowing (Experimental)

As we noted above, waveform approach suffers from slow convergece especially if
• we iterate over a large time window,
• Eigenvalues are “large”.

To overcome these difficulties, we used a new approach we called a Dynamic Windowing. Here is the idea:
for simplicity let us consider the scalar problem above y'=λ y,   y(0)=1. The exact solution and Picard iterations (see above) are
	y(t) = (e^(λt)) = (∑_(k=0)^∞)((((λ t))^k)/(k!))
	(y^n)(t)=(∑_(k=0)^n)((((λ t))^k)/(k!)).
Notice that (y^n)(t) is just the (n^th) degree Taylor polynomial of y(t) around t=0. From calculus, we know that
	|y(t)-(y^n)(t)| <((Max|(y^((n+1)))|)/((n+1)!))|t((|)^(n+1)).
Now given n and a tolerance Tol, find the window width H so that 
	((Max|(y^((n+1)))|)/((n+1)!))(H^(n+1))≤ Tol.
Estimating Max|(y^((n+1)))| by ρ:=|max Evalue| of the linearized system at the end of the previous window, and solving for H, yields
	H≤ sf*(((((n+1)! Tol))^(1/(n+1)))/ρ)   ... (2)
where sf is a safety factor.
• Numerical examples suggest a very good results using this approach in terms of timing, number of function evaluations.

Possible Parallelization/ Acceleration (not implemented)

Gauss-Jacobi Iteration

Gauss-Seidel Iteration

Numerical Tests / Conclusions

We have tested both the Fixed and Dynamic windowing approaches on Gravitational and Kepler problems from the book. The codes are attached.

• In the dynamic windowing WR, the stepper and driver have the following syntax
	 DWR_StepperName[n_,r_,y_,H0_,f_]
	DWRDriver[t0_,tmax_,y0_,r_,n_,Tol_,stepper_,f_]
where
n: desired number of iterations on the window H0
r: number of nodes in the window
y: the initial value at the current window
H0: Width of the time window as suggested by the Driver, computed using (2)
f: name of the function on the right hand-side of the ODE
stepper: is the stepper name, such as DWRSIMP
t0: inital time
tmax: ending time
y0: y@t0
Tol: the user define tolerance as in (2)

• In the fixed window approach, the stepper and driver are compined in a single code and it takes the format
	FWR_Name[n_,r_,Tol_,y_,H0_,f_]
where
n: desired number of iterations (over all time interavl)
H0: the width of the time interval
Tol: the iteration tolerance as defined in (1)
y: the initial value (not nessecary y0)
Name: Name of the quadrature method to use, such as SIMP

We compared both approaches in tems of number functions evaluations and timings. The results are shown in the following tables:

({{"Gravitational", , }, {"Method", "# Func Evals", "Timings"}, {"DWR-SIMP", 12600, 0.561601`7.200972865933416}, {"DWR-TRAP", 18200, 0.8112015`7.360673738543819}, {"DWR-MID", 50760, 2.3868042`7.829361786885906}, {"FWR-SIMP", 159400, 7.8000137`8.343640358985065}, {"FWR-TRAP", 130000, 6.9420122`8.293030366067898}, {"FWR-MID", 359400, 16.8480296`8.678094110342208}})  ({{"Kepler", , }, {"Method", "# Func Evals", "Timings"}, {"DWR-SIMP", 2720, 0.0936002`6.42282177021275}, {"DWR-TRAP", 3200, 0.0780001`6.3436401529738955}, {"DWR-MID", 10260, 0.3120005`6.945700283498601}, {"FWR-SIMP", 17174, 0.5304009`7.1761492458170855}, {"FWR-TRAP", 80000, 10.0152176`8.452205383116842}, {"FWR-MID", 134274, 5.4444096`8.187495784591581}})
D = Dynamic Windowing,  F = Fixed Window

Comments/concluding remarks:
1. We can notice the superiority of the simpson method over the other inegration rules, TRAP and MID.
2. Even the Dynamic windowing approach is built on a crude analysis, it is proved to be efficient and extremely fast.
3. WaveForm relaxiation methods are most commonly used in special applications such as electrical circuits simulations, the main advantages are its simplicity and the ease of parallelization.
4. Drawbacks: large memory requirement, slow (as it is) unless they combined with some sorts of acceleration techniques.
5. We have tried hard to study the stability but we failed to come up with usefull and interesting results.

References

[1] Yanan Yu and Ashok Srinivasan, Dynamic Iterations for the Solution of Ordinary Differential Equations on Multicore Processors, IPDPS, 2009.
URL: http : // www.cs.fsu.edu/~asriniva/papers/ipdps09.pdf
[2] J.White, F.Odeh, Alberto L.Sangiovanni - Vincentelli and A.Ruehli, Waveform Relaxation : Theory and Practice, EECS Department, University of California, Berkeley, Technical Report No.UCB/ERL M85/65, 1985.
URL: http : // www.eecs.berkeley.edu/Pubs/TechRpts/1985/ERL - 85 - 65. pdf
[3] Jacob K.White, Multirate Integration Properties of Waveform Relaxation with Application to Circuit Simulation and Parallel Computation, EECS Department, University of California, Berkeley, Technical Report No.UCB/ERL M85/90, 1985.
URL: http : // www.eecs.berkeley.edu/Pubs/TechRpts/1985/ERL - 85 - 90. pdf
[4] A.Bellen, Z.Jackiewicz, M.Zennaro, Contractivity of Waveform Relaxation Runge - Kutta Iterations and Related Limit Methods for
Dissipative Systems in the Maximum Norm, SIAM Journal on Numerical Analysis, Vol.31, No.2 (Apr., 1994), pp.499 - 523
URL: http : // www.jstor.org/stable/pdfplus/2158013. pdf