implicit variable time-step Reuse

WORDT NIET UITGELEENrJ

Reuse of preconditioners in implicit variable time-step methods

Jawad Al-Temimi

Master's Thesis

Reuse of preconditioners in implicit variable time-step methods

Jawad Al-Temimi

University of Groningen Informatica

Postbus 800

Rijksunjve,sjtejt Groningen BibI,othk Wiskunde &Infornjç

Postbus 800 9700 AV Groningen 'ef. 050 - 3634001

Contents

1

Introduction

²

2 Problem Description

2

2.1 Model Problem 2

2.2 Space Discretization of the Problem ^{. .} 3

2.3 Time Discretization 5

2.4 The variable time-step algorithm ^{. . .} 6

3 Various rewritings of the linear system and results

⁷

3.1 Introduction 7

3.2 Method 1 8

3.3 Method 2 ¹¹

3.4 Method 3 ¹³

3.5 Method 4 ¹⁴

3.6 Condition numbers and iteration counts ¹⁶

4

Discussion and Conclusion

18

1

Introduction

In this report we study the reuse of a LU-decomposition in an implicit method for time-integration with a variable time step. The model problem is a parabolic equation and the 9-method is used for the time-integration. The discretization of this problem leads in every time step to a system of the form Ax

= f.

In general, the work required to compute the LU-decomposition is 0(n3) and the work to backsolve with computed LU factors is 0(n2). Though these amounts are for full matrices, also for sparse matrices the factorization and solve is expen- sive. In any case, it is efficient to reuse the LU-decomposition as much as possible.

If L\t is constant, the system will not change and we can solve Ax =

f

^{easily by}

using the preconditioner or LU factorization of A. If it will take various values, the diagonal of the matrix will change. We will study several rewritings of the system in order to be able to reuse the preconditioner of previous time steps.

2 Problem Description

In this section we will describe the continuous model problem and its discretiza- tion. For the discretization we follow the method of lines approach, hence we first descretize in space and next apply a suitable method to the resulting system of ODE's. In one case we apply the 9-method

2.1 Model Problem

We want to solve the following parabolic equation in two dimensions

ut=u+u+g(x,y,t), O<x,y<1, t>O

subject to the initial and boundary conditions

u(x,y,O)=f(x,y),

for

O<x,y<1

u(O,y,t) = fL(y,t)

=

^0,

u(1,y,t) =

fR(Y,t)

=

⁰

u(x, 0, t) = fB(x,

t) =

sin((fiexp(—0.lt) + f2)t) sin(irx), for

x =

0, ¹

u(x,l,t) =

fT(x,t)

=

^0,

and g(x,y,t) =

^0,

for t>0

2

2.2 Space Discretization of the Problem

For the space discretization we use a cartesian grid with meshsize h=—1---_N+1'

x=ih,y3=jh, 1i,jN

On this grid we use the following approximations

021t(x2,y3,^t) U÷1,,(t)—2U1,(t)+U1_i,(t)

Ox2 — h2

y3,t) U,,i(t)—2U, (t)-i-U,,_1(t)

-

where U2,3(t) = U(x,yj, t)

This results use a system of ODE's

— F •(U,t) = ^HU1,3

+g,3(t),

dt ^{— 2,3}

where

g3(t) =

g(xi,y3, t), and the difference operator H is implicitly defined by

— h ⁺ U1,3÷1(t)—2U, (t)+U1,_1(t)

HU1 2(t) — U1i,,(t)—2U, (t)+U1_1,3(t

___________________________

h

Hence we have the following system

— U1,3(t)—2t4, (t)+U_1,, (t) Ui,j+i (t)—2U;, (t)+U,_i (t)

dt — ^h1

Along the boundaries we have

=

U2,,(t)—2U1,3(t)

+ U1,÷i(t)—2Ui,(t)+Ui,_i(t)

+ ^(Ubot) + gi,3(t))

dt h2 h2

— —2UN,3(t)+UN_l,j(t)

+ UNi+1(t)—2UNg(t)+UNi_1(t)

+ (Uitt)

_{+ 9N,j(t))}

dt — h2

= (Ji,j(t)—2U,i(t)+U1_1,1(t) U1,(t)—2U1,1 (t) (Ui(t)

—F gi,1 (t))

dt h2 h2

•j=N,i=1,2,...,N:

= U+l,N(t)2U,N(t)+Ui_1,N(t) —2U,N(t)+U1,N_1(t)

+

^(Ut.N1(t) + g;,p,r(t))

dt + _h2

Where

U0,,(t) = U(x0,y,t) = fL(Yj,t), ^{UN+1J(t) =} U(XN+l,yj,t)

=

fR(yj,t),

U2,0(t) = U(x,yo,t) =

fB(x,t),

^{and U,N+1(t) =} U(x2,yN+1,t) = fT(xI,t),

where i = ^1,

2,..

^., N,

j

^{= 1,}

2,...

^, N,

and t> 0

For the implementation it is convenient to define

G,3=g, for

G1, g1,3(t) + ^Uo(t) = g1,,(t) + for

j

^{= 2,3,. . .,N}

GN,3 = gN,3(t) + ^UN+,(t)

=

^gN,3(t)

+ ftt),

for

j

^{= 2,3,. .. ,N — 1}

C2,i = g,1(t)

+

^Ui,o2(t)

= g,i(t)

+

1B,t),

_{for i = 2,3,. .. ,N — 1}

G,N = g,N(t) + = g,N(t)

+ 1T,t),

_for

_{i =}

_{2,3,. .. ,N — 1}

and

C1,1 = g1,(t) + [fB(x1, t) + fL(Y1, t)]

G1,N = g1,N(t)

+ [fT(xl,t) + fL(yN,t)1

GN,1 = gN,1(t) + [fB(xN, t) + fR(Y1, t)]

CN,N = gN,N(t) + [fT(xN, t) + fR(YN, t)]

By using this C instead of g we can solve the same problem using homogeneous Dirichiet boundary conditions. In the following we write the system in matrix vector notation.

First define the following matrices

—41 0...

⁰

LI 0...

⁰

L=(0 T

:: .

where I is the identity matrix of order N. Furthermore, we define the vectors

U1 U1,j

U2 U2,j

U3 U3,j

UN_i UN_1,

UN UN,3

and likewise

C1(t) G1,,(t)

G2(t) G2,3(t)

G3(t) C3 .(t)

G(t) = ^. ,G,(t)

'

GN_1(t) GN_1,2(t)

GN(t) GN,,(t)

Herewith the system of ordinary differential equation can be written as

=fiu+c

^(2.1)

U(O) = F1 (2.2)

2.3 Time Discretization

We can apply any suitable time-step method to (2.1). We choose the 9-method which yields

= O(uiU' + G(t')) + (1 —

^O)(11U

+ G(t')), where Lt

denotes the current time step.

This is in fact a reuse of the form

A1U' ⁼

^{+ V (2.3)}

where

= —I

^{— OH,}

P= -i+(i-O)k

V =

9G(t')

^{+ (1 —O)C(t')}

For 9 = ^we have a second order account discretization. For other values of 9 E [0, 1] it is first order only. For the analysis it is convenient to rewrite the system to be solved.

Define

2 2

a

= F= ^—h2H,

f+i = -(PU ^{+ V), A1 =} aI + F

^(2.4)

then (2.3) is equal to

A TTfl+l —

— Jn+1

5

F

2.4 The variable time-step algorithm

In this section we will sketch an algorithm to control the step size Lt in our time integration problem.

Algorithm 1.

Choose Lt =

(a) We will take two steps (algorithm 2.) to solve (2.3) from t till t + 2it, First we begin from t to t + t, Second from t + tXt to t + 2t with step size At,

we have the initial state then we will compute U,

(b)

We compute fni, f, f÷i and D, where D =

^where

f = iIu

+ g(t), [Dt3 approximates local truncation error].

(c) If not <DLt3 <E/.t, where e is the global error, then we will compute

new it,

^by

'

^The new value of At has to be not more than two time the old value if we want to avoid rapid changes in tt.

(d) If Lt not decreased, then increase t by two times the old value of At.

(e) go to (a).

Algorithm 2.

Step: (Perform one time step).

(a) Compute

A1 from (2.4).

(b) Solve A1U'' =

fn+i from (2.5).

We assume that we have a good preconditioner PA for A, so the equation

AU ₌ f

can be solved easily. The goal is to find a way to solve the equation

A+iU' =

^fn+1 without constructing a new preconditioner for the matrix AM4, but rewrite the equation such that we can reuse PA. A+iU44 = fn+1 is solved in four ways to be described in the next section.

The result of this algorithm Figure 1 is the solution of parabolic equation (two dimension) using Crank-Nicolson method (0 =

)

with variable time-stepping, for two time values (top) and of a specific point in space (bottom left), time-step during integration (bottom right).

0.4 0.2 0

-02

—04 -0.6

Solii T- 6.3

0

^{50 100} _T ^{150 200 250}

Figure 1: solution of the model problem variable time-stepping, for two time values (bottom left), time-step during integration

0 ^{50 100 150} I

20

200250

using Crank-Nicolson (9 =

)

^with

(top) and of a specific point in space (bottom right)

3 Various rewritings of the linear system and results

For convenience we define the systems By = g and Ax =

f

which are equivalent

to A+iU1 = f+i ^{and AU =}

^In, respectively. We will discuss four ways to solve the linear system By = g without constructing a preconditioner for B, but by rewriting B to an expression such that we can reuse the preconditioner of A. In this report we will take as preconditioner of A the inverse of A. The application of this preconditioner means that we have to solve a system with A.

For that we use the LLT factorization of A.

3.1 Introduction

We need to know how we can find the eigenvalues for our matrices and also the condition numbers because the convergence depends on the condition num- ber. First we state an important theorem (Gerschgorin theorem) to locate the eigenvalues of the matrix A:

Theorem 1. (Gerschgorin's circle theorem) Let A be a square complex matrix.

Around every element a12 on the diagonal of the matrix, we can draw a circle with radius the sum of the absolute value of the other elements on the same

Sdution at T- 136.3

row E3

^a3I. The interior of such circles are called Gerschgorin discs. Every eigenvalue of A lies in one of these Geschgorin di.scs.

Theorem 2. The eigenvalues of a real symmetric matrix are real.

Corollary 1. The eigenvalues of the matrix F defined in (.4) are real and on the interval 0 A <8.

Definition 1. If V is a vector space over C, the spectrum of a linear mapping

T: V —* V is the set

o(T) =

{A E C: I T — Al I is not invertible},

where I denotes the identity mapping. If V is finite dimensional, the spectrum of T is precisely the set of its eigenvalues. For infinite dimensional space this is not generally true, although it is true that each eigenvalue of T belong to ci(T).

The spectral radius of T is

p(T) = sup{A : A E a(T)}.

In the following we rewrite the linear system By = gsuch that the preconditioner of A can be employed and analyse it. The first is simply preconditioning by A.

The second is a rewriting of the system. The third is the second with A as

preconditioner. In the fourth method we try to minimize the condition number using a free parameter.

Lemma 1. Let h(A) be a rational function. Suppose A is an eigenvalue of F then h(A) is an eigenvalue of h(F).

Proof:

Suppose U*FU = A, where U is the matrix of all eigenvectors of F, and A is a diagonal matrix with the eigenvalues of F on the diagonal. Then U*(h(F))U = h(U*FU) = h(A). Hence h(diag(A1, A2,..., As)) = diag(h(Ai), h(A2),..^{. ,}h(A)), giving the eigenvalues of h(F) on the diagonal.

3.2 Method 1

We will use the preconditioner of A to solve the equation By =

, so

^{we have to}

solve

(aI + F)'[bI + F]y = , where =

(al

+ F)'g

We can rewrite this system in the following way:

hi(F)y =

where

hi(F)=(aI+F)'[bI+F]

Using Lemma 1 the eigenvalues of h1 (F) are:

h1(A) =(a+A)'[b+AJ

(3.1)

where A is an eigenvalue of F.

In Figure 3 we have drawn h1(A) on the interval [0,8] where the eigenvalues according to Corollary 1 are located. It seems that h1(A) is decreasing function for a < b, and increasing function for a> b.

a=O.l.b—O.2 a..O.I.b.O.06

0.6

04 ¹

) 2 4 6 8 0 2 4 6 8

a=1,b=2 a.lb.0.5

040;8

a=10b—20 a—10b—5

0768 168

^{(e) (e)}

Figure 2: Plots of h1(A) for various values of a and b

This property is convenient to get an impression of the condition number is(h1 (F)) which is bounded by:

,c(h1(F)) ki(a,b), where Ici(a,b) = (3.2)

mlno<A<8 Ihi(A)I

In order to show that indeed h1 (A) is monotonous, we consider the first derivative of hi(A):

h' 1A' — (a+A)—(b+A) — (a—b) 1k 1 — (a+A)2 — (a-I-)2

Since we take positive time-steps, a and b are positive. Then we have two cases:

• if a <

b, then h (A) 0, we conclude that h1 (A) is a decreasing function, and condition number is bounded by:

< h(O) — b(a+8) h1 F

— h1(8) — a(b+8)

• if a> b, then h (A) 0, and the condition number will be bounded by:

, ^, < ^{hi(8) —} a(b+8)

icihiiFjj

— h1(0) — b(a+8)

The upperbound ki(a, b) is also a good approximation since, for the standard discretization employed, we know that the extremal eigenvalues of F tend to 0 and 8 when the mesh is refined.

In Figure 4, k1 the estimate of the condition number of h1 (F), is depicted for some values of a and b < 2a. We bound b in a since usually in adaptive time-step algorithms one allows only a change by a factor 2.

1.0.1 a.0.1

0.05 0.06 0.07 0.08 0.09 0.1 0,1 012 014 016 0.18 02

____

1

a—b a=10

(b) (b)

Figure 3: The estimated condition number k1 (a, b) for some values of a,

b <

2a

The estimated condition numbers for b = ^2a and b = (the end points in the graphs of Figure 4) are:

a

a+16 a+8

ki(a,

&

= a+ 8 ^{E [1,2]} and K1(a, 2a)

=

a+ 4 ^{[1,2] (3.5)}

3.3 Method 2

In this section we will try to rewrite the equation By =^g such that precondition- ing with A-' has more effect. We will consider the rewriting in this section and we will study the preconditioned version in the next section. We will rewrite B in an expression of A as follows:

B=()[A—F]+F=()A+F—()F= ^(A+F-F)

Define C = ( — 1)F, then By =

()(A

_{+ C)y =} ^g

Multiplying the left- and right-hand side by (I —

CA')

yields the second system we will study:

()(A—CA1C)y=2

where

(I —

CA')g

Again the system matrix can be expressed in F:

()(A

^—

CA'C) = ()[(aI + F)

_— ( — 1)F(aI +

F)'(

— 1)F]

So we can rewrite the above system as follows:

h2(F)y

=

where we can rewrite h2(F) to

h2(F) =

()[(aI

+ F) — ( — 1)F(aI + F)'( — 1)F]

=

()(aI + F)'[(aI + F)2 — F2(

^{— 1)2]}

= ()(aI + F)'[a2 + 2aF + F2 — F2 _{+ F2}

- F2]

=

()(aI

⁺

(al +

So using Lemma 1, the eigenvalues of the matrix h2(F) are given by:

h2(A) = (a+ ^{))_l[(2b —a)A2} ± 2b2A + ab21

(3.6) where A is an eigenvalue of F.

Figure 5 shows h2(A) for some values of a = ^0.1, 1, 10 and b = 2a and b = , and we see that for these cases h2(A) is an increasing function.

We can also show that h2(A) is in general a monotonic increasing function by considering the derivative:

..0.1bO2 a.O.1b.0.05

___

a=lb=2 a—lb—O.5 8

68

^{a1O.b=20 0'}

²⁴⁶³

^a.l0,b.5

a (e)

02468

(e)

Figure 4: Plots of h2(A) for various values of a and b

h'

f\\

— b(a+)[2A(2b—a)+2b2]—b[(2b—a)A2+2b2A+ab21 —

____________________

2k 1 — b2(a+A)2 (a+A)2

We consider in our experiments only cases where (2 _—

)

0, which makes the counter positive for al ). 0. Hence, since the denominator is positive, we have Hence h2(A) is increasing function. So h2(8) is the maximum of h2(F), and h2(0) is the minimum:

h "8' — ⁽ 8\_iI(2a)(8)2+2(8)th2l ^— 64(2—)+16b+ab — 64(2—)+b(a+16)

2k 1 — a + _{b J — (a+8) — (o+8)}

h2(0)

= (a)'[ç]

^{= b}

An estimate of the condition number is now given by:

'

_b'— ^{h2(8) —}

64(2—)+b(a+16)

^— 64(2—f)

a+16

ic2a, j—

— — + 3.7

h2(0)

b(a+8)

b(a+8)

a+8

In Figure 6, k2 the estimated condition number depicted for some values of a and b with

<b

^{< 2a.}

For b = ,

^{2a we have:}

a

a+16

⁴⁸

a+16

ic2(a,

& = a + 8 ^{E [1,2] and K2(a, 2a)}

= a(a +8) +

a +8

^{e [1,oo) (3.8)}

Note that k2(a, )

=

ki(a, )

^and that k2(a, 2a)ki(a, 2a). So no advantage with this approach is expected, it is only a starting part for method 3.

a..O.1 a.1 a.1O

005 01 _(b)0Th

02 51l52 '5W15

_{(b) (b)}

Figure 5: The estimated condition number k2(a, b) for some values of a,

b

2a

3.4 Method 3

We will study in this section the system of the previous section with precondi- tioner A:

A-'()(A

^—

CA'C)y = A'2

Hence, up to a factor, the system matrix is of the form (I — ^{(A—'C)2). So for} (A-'C) small we will have a condition number close to 1. Again we write the system matrix as a function of F.

Clearly:

h3(F) = (al

+ F)-'h2(F)

h3(F) = (al

+

So using Lemma 1, the eigenvalues of h3(F) are given by:

h3(A) = (a + A)_2[2b — a)A2 +2b2A + ab2]

(39) where A is an eigenvalue of F.

In Figure 7, h3(A) is drawn for some values of a, b =^2a, and b = , and we see that this function is decreasing.

We can show that h3(A) is in general a monotonic decreasing function by consid- ering the derivative:

h' — b(a+A)2 [2A(2b—a)+2b2] —[2b(a+A)I [(2b—a)A2+2b2A+ab2j 3k 1 —

— —A(2a2—4ab+2b2) — —2A(a2—2ab+b2) — —2A(a—b)2

— b(a-f.A)3 — b(a+A)3 — b(a+A)3

Hence, for the considered values of a, b and A (all being positive) h'3(A) <0, and therefore h3(A) is a decreasing function. So h3(0) is the maximum of h3(F) and h3(8) is the minimum:

h "8' ^— 8\_2[(2b—t1)(8)2+2b2(8)+(th2l ^— 64(2—)+16b--ab

— a + I L b — (a+8)2

h3(0) =

(a)_2[] = b

a.0.1.b.02 a—O.l.b—0.06

0[

2 0.6

18 04

6 8

.—1.b—2 a—1.b.0.5

____

02468

a—l0.b—20 a—1Ob—5

___

(e) (•)

Figure 6: Plots of h3(7t) for various valuesof a and b An estimate of the condition number is now given by:

b — ^{h3(0) —}

_________

— b(a+8)2

k3(a, 3

— h3(8) ^— 64(2-)+16b+ab — 64a(2_—

)

+ l6ab + a2b ^{( .10)}

(a+8)

We see in Figure 8 the estimated condition number depicted for some values of a and b with

b 2a.

For b = 2a we have:

a

(a+8)2 (a+8)2

⁴

E[1,00) and

(3.11)

Note that k3(a, 2a) <ki(a, 2a).

3.5 Method 4

In this section we will study (using Mathematica) the following preconditioned system:

a— a—l0

45 ^{• 01245}

4 0124

01235 3

0123 2.5

2 01225

1.5 0122

1

___ ______

05 ¹ 15 2 5 10 15 20

(b) (b)

b+A

(a—b)(b+A)

h4(A)= —

________

a+A

a(8+a)

a=O.l 40

35

30 25

r20

15 10

005 01 015 02

(5)

+

Figure 7: The estimated condition number 1c3(a, b) for some values of a,

b

2a

(A—' — yI)By =

where

=

^(A—' —'yI)g

Again we can express the system matrix in F:

(A-'

^— 'yI)B = h4(F) where

h4(F)=((aI+F)' —'yI)(bI+F)

The freedom 'y is chosen such that h4(A) is as constant as possible on the interval

[0, 8]. This leads to the requirement that h4(O) = h4(8) which results in 'y =

a(8±a) So h4(A) assumes the form:

(3.12)

In Figure 8, h4 (A) is shown and we see that there is only one internal extreme which is a maximum for a> b and a minimum for a < ^b. The extreme is found for A = —a

+

+ a2 with value:

(—a+b+s/8a+a2)2

8a + a2 313

For a > ^b we can see that h4(0) is the minimum eigenvalue of h4(F), then the condition number for a> b is bounded by:

ic(h4(F)),c4(a,b)= (-a+b+/8a+a2)2

(3.14)

b—a/2 b—2a

oe

__

03C0246

(0) (e)

Figure 8: h4(A) for various values of a, b 2a and b =

For

a <

^b the estimated condition number is just the reciprocal of the one for

a>b.

In figure 9, k4 the estimated condition number is depicted for some values of a and b in

b 2a. Note that k4(a, )

< k(a,

)

^and that k4(a, b) < ^{k1(a, b)}

for a> 0 and <b

< 2a.

3.6 Condition numbers and iteration counts

In this section, we compare condition numbers and iteration counts. We see in Table 1 the estimated condition number for all methods for some values of a and b. In Table 2 we find the iteration numbers of all methods using MATLAB's CG to solve the system By = g. The initial guess is zero, g is a random vector, and the iteration is stopped if the residual is less than 1.10—6.

We see the estimated condition number of Method 2 is always larger than that of Method 1, as already derived in general in Section (3.3). This is expressed in the number of iterations shown in Table 2. Method 3 gives condition numbers smaller than those of Method 1 if a < b as noticed in Section (3.4). However, for a > b

it is much worse than that of Method 1 for a, b small, so it is not appropriate to use it in that case.

The estimated condition number of Method 4 is smaller than that of Method 1,

a0,I ^a.0.1

0.05 0.06 007 0.08 0.09 0.1 0.1 012 0.14 0,16 018 02

a—I a.1

IC ^. 1.4

08 09

16l82

a—b a—10

Figure 9: The estimated condition number k4(a, b) for some values of a,

b

2a

which is not surprising since Method 1 is a special case in the class of considered methods to derive Method 4. Method 4 will be the one with the smallest condition number in this class. We see that smaller condition numbers do not always result in a smaller number of iterations. This might be due to the initial solution and right-hand side.

a

^b

^{k(hi(F)) [}

^{ic(h2(F)) ic(h3(F)) ic(h4(F)) El}

0.1 0.2 1.9756 61.2469 1.3225 1.64 0.1 0.05 1.9877 1.9877 40.7516 1.7950

1 2 1.8000 7.2222 1.2462 1.25

1 0.5 1.8889 1.8889 4.7647 1.4705

10 20 1.2857 1.7111 1.0519 1.0212

10 5 1.4444 1.4444 1.2462 1.0897

Table 1: The condition number ic of

h1(F), h2(F), h3(F) and h4(F) for a =

0.1,1,10 and

<b2a

fta

0.1

b yi(iter)

[ y2(iter) ] ^{y3(iter) y4(iter)}

0.2 5 25 5 7

0.1 0.05 4 4 16 7

1 2 6 15 5 5

1 0.5 6 6 12 6

10 20 5 7 4 3

10 5 6 6 5 4

Table 2: The iteration number for n = 10

4 Discussion and Conclusion

Method 4 shows that if we want to reduce the condition number then that can be obtained by using a free parameter. It would be of interest to generalize this idea. For example, the preconditioner of Method 3 is simply the product

A'(I

^{— CA—')} arid, since C is just a factor times F, it belongs to the class of

preconditioners of the form A—' —

'yA'

FA'. Method 1 is in this class for 'y = 0,

the method with the smallest condition number in this class will be better than both Methods 1 and 3.

It is doubtful whether the extra effort that is put in a preconditioner is payed back in a sufficient decrease of the number of iterations. We see in Method 3 an extra solve, which often dominates the cost, is necessary per iteration. We hope that the number of iterations is reduced by a factor 2. This is in the cases where it functions well, a < ^b clearly not the case.

Maybe an ideal method in the above mentioned class will do better. Method 4 gives a reduction of the condition numbers with hardly no extra costs, so this will do better in general than Method 1. In practice we deal with an approximation of the inverse of A due to an incomplete factorization. Then the optimization process based on the inverse of A only leads to reasonable results if the incom- plete factorization is accurate. In that case it does not make sense to consider sophisticated approximation with more free parameters. So only simple forms

are relevant.

Our conclusion in this study is that computation time can be brought down by doing something just slightly less trivial than just preconditioning with the previous matrix. It remains however to be investigated whether this holds also if we consider incomplete factorization.

References

[1] John C. Strikwerda, Finite Difference Schemes & Partial Differential Equa- tions, University of Wisconsin-Madison, 1989.

[2] Anne Greenbaum, Iterative Methods for Solving Linear Systems, University of Washington, Seattle, Washington, 1997.

[3] Kendall E. Atkinson, An Introduction to Numerical Analysis, University of Iowa, 1988.

[4] Joel N. Franklin, Matrix Theory, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1968.