Multiple grid methods for equations of the second kind with applications in fluid mechanics

(1)

M U L T I P L E

G R I D M E T H O D S

F O R E Q U A T I O N S

O F T H E S E C O N D K I N D

W I T H A P P L I C A T I O N S

I N F L U I D M E C H A N I C S

H . S C H I P P E R S

(2)

(3)

o Ni M Oí O CS O -4

l i l i

lili!

i t un •o o * Ul 0s BIBLIOTHEEK TU Delft P 1696 4035 C 512586

(4)

(5)

M U L T I P L E

G R I D M E T H O D S

F O R E Q U A T I O N S

O F T H E S E C O N D K I N D

W I T H A P P L I C A T I O N S

I N F L U I D M E C H A N I C S

(6)

(7)

M U L T I P L E

G R I D M E T H O D S

F O R E Q U A T I O N S

O F T H E S E C O N D K I N D

W I T H A P P L I C A T I O N S

I N F L U I D M E C H A N I C S

TER VERKRIJGING VAN DE GRAAD VAN DOCTOR IN DE TECHNISCHE WETENSCHAPPEN AAN DE TECHNISCHE HOGESCHOOL DELFT, OP GEZAG VAN DE RECTOR MAGNIFICUS PROF.IR. B.P.TH. VELTMAN,

VOOR EEN COMMISSIE AANGEWEZEN DOOR HET COLLEGE VAN DEKANEN TE VERDEDIGEN OP

DONDERDAG 22 APRIL 1982 DES MIDDAGS PROEFSCHRIFT OM 14.00 UÜR DOOR H A R M E N S C H I P P E R S Wiskundig ingenieur ST. NICOLAASGA GEBOREN TE

MATHEMATISCH CENTRUM, AMSTERDAM 1982

(8)

door de promotor

(9)

(10)

(11)

Tab Ze of contents ^ Acknowledgements t t i 1. INTRODUCTION 1 1.1. Compact operators 1 1.2. P r o j e c t i o n methods 3 1.3. M u l t i p l e g r i d methods 4 1.4. Scope of the study 7 2. MULTIPLE GRID METHODS FOR THE SOLUTION OF FREDHOLM INTEGRAL

EQUATIONS OF THE SECOND KIND 11

2.1. Introduction '1 2.2. Basic assumptions and r e s u l t s 12

2.3. I t e r a t i o n schemes with Nyström i n t e r p o l a t i o n 15 2.4. I t e r a t i o n schemes with p r o j e c t i o n i n t o f i n i t e dimensional

subspaces 21 2.5. Numerical r e s u l t s 26

3. AUTOMATIC NUMERICAL SOLUTION OF FREDHOLM INTEGRAL EQUATIONS OF THE

SECOND KIND 37 3.1. Introduction 37 3.2. T h e o r e t i c a l foundations 38 3.3. Automatic program 43 3.4. Numerical r e s u l t s 49 Appendix 54 4. MULTIPLE GRID METHODS FOR INTEGRAL EQUATIONS IN POTENTIAL THEORY. 59

4.1. R e g u l a r i t y r e s u l t 60 4.2. E r r o r a n a l y s i s 70 4.3. M u l t i p l e g r i d methods 80

4.4. Numerical r e s u l t s f o r smooth contours 87 4.5. A p p l i c a t i o n of m u l t i p l e g r i d methods to the c a l c u l a t i o n of

c i r c u l a t o r y flow around an a e r o f o i l 93

5. OSCILLATING DISK FLOW I l l

5.1. Introduction 112 5.2. A n a l y t i c a l r e s u l t s 114 5.3. Numerical approach 117 5.4. Numerical r e s u l t s 122 Summary 127 Samenvatting 129 Curriculum v i t a e 131

(12)

(13)

The author wishes to express h i s g r a t i t u d e to the governing board of the " S t i c h t i n g Mathematisch Centrum" f o r g i v i n g him the opportunity to carry out the i n v e s t i g a t i o n s presented i n t h i s t h e s i s , and f o r p u b l i s h i n g t h i s book.

During h i s stay a t the Mathematical Centre the author has g r e a t l y appreciated the pleasant cooperation w i t h the colleagues of the department of Numerical Mathematics. S p e c i a l thanks go to Dr. P.W. Hemker f o r the pleasure of several i l l u m i n a t i n g discussions and to Mr. P.M. de Zeeuw f o r assistance w i t h the programming.

For the t e c h n i c a l r e a l i z a t i o n of t h i s book the author wishes to thank Mrs. C.J. K l e i n Velderman-Los f o r her c a r e f u l typing of the manuscript, Mr. D. Zwarst f o r o r g a n i z i n g the production process, Mr. T. Baanders f o r the design of the front-cover, and Mr. J . Schipper and h i s colleagues who took care of the p r i n t i n g and the b i n d i n g .

(14)

(15)

CHAPTER 1

INTRODUCTION

In t h i s t h e s i s m u l t i p l e g r i d methods are studied f o r s o l v i n g the a l g e b r a i c systems that occur i n numerical methods f o r Fredholm equations of the second k i n d :

(1.1) f = Kf + g.

Let g belong to a Banach space X. We assume that the operator K i s compact from X i n t o X, and that (1.1) has a unique s o l u t i o n f e X. For the d e s c r i p -t i o n of -the m u l -t i p l e g r i d -techniques we use p r o j e c -t i o n me-thods. Equa-tions of the type (1.1) o f t e n a r i s e i n a p p l i c a t i o n s : p o t e n t i a l problems, s o l i d mechanics, d i f f r a c t i o n problems, s c a t t e r i n g i n quantum mechanics, water waves, e t c . ( c f . [63, [ 7 ] , [ 1 0 ] , [11]). I n the f o l l o w i n g sections we b r i e f -l y review the concepts of compactness, p r o j e c t i o n methods and m u -l t i p -l e g r i d methods.

1.1. COMPACT OPERATORS

H e u r i s t i c a l l y speaking, compact operators are operators that i n general possess some kind of "smoothing property". The f o l l o w i n g d e f i n i t i o n i s bor-rowed from ATKINSON [ 3 ] :

DEFINITION 1.1.1. Let X and Y be Banach spaces. A l i n e a r operator K: X •+• Y i s c a l l e d compact i f f o r every bounded sequence { fn} i n X, there i s a sub-sequence { fRj } such that {Kfn^} converges i n Y.

E q u i v a l e n t l y , K maps bounded subsets of X i n t o subsets w i t h compact closure i n Y.

A simple example of a compact operator i s given by the f o l l o w i n g i n t e g r a l operator on X = Y = C[0,1]:

(16)

(1.2) K f ( x ) = ( fe(x,y)f(y)dy, 0 where sup 0<x<l |fc(x,y)|dy < and l i m sup 6-H) 0<x<l |fe(x+6,y) - fe(x,y)|dy = 0. 0

By the A r z e l a - A s c o l i theorem i t follows that K i s compact on C[0,1] ( c f . [ 3 ] ) . A more abstract example i s given by Schauder i n v e r s i o n . Let Au = g

(g belongs to a Banach space Y) be a q u a s i l i n e a r e l l i p t i c operator equation w i t h :

Au = Au + a(u)-Vu + b ( u ) .

Define the l i n e a r operator

Avu = Au + a(v)»Vv + b(v) f o r f i x e d v e Z,

where the Banach space Z i s a subspace of the Banach space X with Z compact-l y imbedded i n X. I f A^u = g has one and oncompact-ly one s o compact-l u t i o n u = Kv s a t i s f y i n g the a p r i o r i estimate

Bullz £ C(R)llgllv, Vg e Y,

i f IIvllx < R, then the non-linear operator K: X •* X i s compact ( c f . [ 4 ] ) . The f i x e d points of

(1.3) u = Ku

coincide w i t h the s o l u t i o n s of Au = g. In t h i s t h e s i s examples w i l l be studied that are s i m i l a r to (1.2) and (1.3).

(17)

1.2. PROJECTION METHODS

A d e s c r i p t i o n of our m u l t i p l e g r i d techniques can be given by i n t e r -polatory p r o j e c t i o n operators. For the approximate s o l u t i o n of l i n e a r or nonlinear equations of the type (1.1) o r (1.3) p r o j e c t i o n methods are w e l l -known.

Let {X^} be a sequence of f i n i t e - d i m e n s i o n a l subspaces of X with dimen-s i o n N and adimen-sdimen-sume that f o r every N we have a l i n e a r p r o j e c t i o n operator T„ which p r o j e c t s X i n t o X^. The p r o j e c t i o n method f o r s o l v i n g (1.1) i s :

solve the approximate equation

(1.4)

fN " TNK fN + TN « '

f.T £ X . N N

One can solve (1.4) by reducing i t to a G a l e r k i n system ( c f . [ 9 ] ) , which i s a matrix equation of order N of the form

(1.5)

( W 4 = % •

On the assumption that (1.1) has a unique s o l u t i o n , i t can be proven that both ( I - T j j K )-' and (J -£)""' e x i s t f o r N l a r g e enough [ 9 ] .

EXAMPLE (Piecewise l i n e a r i n t e r p o l a t o r y p r o j e c t i o n s i n C[0,1]).

Consider (1.1) - (1.2) i n X = C[0,1]. Let I IN: 0 = xQ < Xj < ... < x^^ = 1 be a r e g u l a r g r i d on the i n t e r v a l [0,1]. Choose to be the space of piece-wise l i n e a r functions on spanned by (UQ,UJ,...,u^) w i t h :

u0( x ) ( X - X j ) / ( X Q - X J ) XQ < X < X J . elsewhere, ^ N - I ^ V W * V l " x - XN ' 0 , elsewhere, u.(x) ' ( x - xi_]) / ( xi- xi l) , ( x - xi + 1) / ( xi- xi + 1) , x. , < x < x., 1-1 1' X. < X < X . j . , 1 1+1 elsewhere,

f o r i 1,. . ,N-1. For the p r o j e c t i o n operator TN we choose the i n t e r p o l -atory operator defined by:

(18)

T f ( x ) = I f ( x , ) u . ( x ) . i=0

The G a l e r k i n system reads

(1.6)

g ( xQ) '

8 ( ^ )

where (^N)£j = !q fe(x^,y)Uj(y)dy. The s o l u t i o n of (1.4) i s given by the piecewise l i n e a r f u n c t i o n t a k i n g the value /\ at x^, i = 0,...,N, where /N i s the s o l u t i o n of (1.6). For other examples of p r o j e c t i o n methods Fredholm i n t e g r a l equations of the second k i n d we r e f e r to IKEBE [ 9 ] .

1.3. MULTIPLE GRID METHODS

The matrix i s i n general a f u l l m a t r i x , i . e . a l l elements of are non-zero. I f N i s small enough (1.5) can be solved by Gaussian e l i m i n a t i o n . For large values of N i t e r a t i v e methods are needed. I n t h i s t h e s i s m u l t i p l e g r i d methods are a p p l i e d i n order to solve the large non-sparse system (1.5) e f f i c i e n t l y . These methods can be seen as an extension of the i t e r a t i v e schemes given by BRAKHAGE [5] and ATKINSON [ 3 ] , who only use two g r i d s (one coarse and one f i n e g r i d ) . M u l t i p l e g r i d methods work w i t h a sequence of grids of i n c r e a s i n g refinement, which are used simultaneously to o b t a i n an approximation to (1.4), i . e . to compute e f f i c i e n t l y the s o l u t i o n on the f i n e s t g r i d .

Let Np, p = 0,1,2,...,1, be an i n c r e a s i n g sequence of integers (NQ < NJ < < N^) and l e t X be a short n o t a t i o n f o r XN ( fp, Tp and K are s i m i l a r l y d e f i n e d ) . In the context of m u l t i p l e g r i d i t e r a t i o n the

subscript p i s c a l l e d level. We need the f o l l o w i n g assumption f o r {Kp):

XQ c X, c ... c Xt c X.

Suppose we want to solve (1.4) w i t h N = N , i . e . :

(1.7) ( I - T ^ K ) ^ = T^g.

An approximation to f i s obtained by the f o l l o w i n g i t e r a t i v e process. Let (0)

the i n i t i a l guess be zero, i . e . =0. Perform the f o l l o w i n g steps f o r i = l ( l ) o :

(19)

(1.8a) f f+ i ) = T ^ K f ^ + TAg ,

( 1

.

8 b ) f

(i

+

D .

f

^i>

+ T i ( I

-

V l K )

- '

V l { T

^ l )

+ T

^ i )

}

.

The f i r s t step (1.8a) of t h i s scheme i s a Jacobi o r P i c a r d i t e r a t i o n ; the second step i s c a l l e d a coarse grid correction. As i s shown i n Chapter 2, the scheme (1.8) converges i f N j i s l a r g e enough. In order to understand the s t r u c t u r e of (1.8b) we f i r s t consider the second part of the right-hand s i d e . A f t e r (1.8a) the term between braces i s equal to the residue of the Jacobi i t e r a t i o n :

( i+ 1)

This residue r ^ 2 e i s projected to the subspace X^_j and we proceed with s o l v i n g

We note that t h i s equation i s of the same type as (1.7), but a l l s u b s c r i p t s are decreased by 1. In a m u l t i p l e g r i d method the s o l u t i o n of (1.9) i s approximated by y steps of the i t e r a t i v e process (1.8), except when i t has to be solved on the lowest l e v e l . In that case the l i n e a r system correspond-ing to (I-TQK)VQ = TQI'J*' i s solved by some d i r e c t method. In Chapter 2 i t i s shown that i t i s s u f f i c i e n t to take y = 2 f o r the number of i t e r a t i o n s .

The proposed m u l t i p l e g r i d method to solve (1.7) i s of a r e c u r s i v e type. Recursion takes place w i t h respect to the l e v e l number p. For the p r e c i s e d e s c r i p t i o n of a m u l t i p l e g r i d algorithm the programming language ALGOL 68 [12] i s convenient, because t h i s language can e a s i l y handle both the r e c u r s i o n mentioned and the data s t r u c t u r e s that appear i n a m u l t i p l e g r i d algorithm. In order to describe our m u l t i p l e g r i d method i n a concise, modular and readable form we f i r s t introduce the f o l l o w i n g ALGOL-68 proce-dures i n which MODE VEC = REF [] REAL.

(20)

PROC solve directly = (VEC f,g) VOID: determines t t

by means of Gaussian-elimination # # determines the s o l u t i o n of (I —K )f

P P

PROC zero = (INT p) VEC:

" d e l i v e r s the vector corresponding to the zero-element of X 9.

P PROC restrict = (VEC y ) VEC:

°P

" d e l i v e r s the vector corresponding to T , y , p-1 'p' i . e . restrict i s a representation of the operator T . : X •+ X , #.

p-1 p p-1 PROC prolongate = (VEC y ) VEC:

# d e l i v e r s the vector corresponding to T , y ,

r ° P+l P

i . e . prolongate i s a r e p r e s e n t a t i o n of the operator T , : X •*• X , #.

r p+1 p p+l

The f o l l o w i n g ALGOL 68 program describes our m u l t i p l e g r i d process:

PROC mulgrid - (INT p,a,VEC f,g) VOID: IF p = 0

THEN solve directly (f,g) ELSE TO a

DO f := g + Kp*fj

VEC residue = g-f+K^*f; VEC v := zero(p-l);

mulgrid(p-l, y, v, restrict (residue) ); f := f + prolongate (v)

OD FI;

The a c t u a l implementation of the procedures zero, prolongate and restrict depends on the choice of ^X^} and the p r o j e c t i o n operators {T }. Using uniform g r i d s and piecewise l i n e a r i n t e r p o l a t o r y p r o j e c t i o n s i n C[0,1], we give an implementation i n Chapter 3. I n [8] HACKBUSCH a l s o studied the above m u l t i p l e g r i d method f o r Fredholm i n t e g r a l equations of the second k i n d . In Chapter 2 we introduce an a l t e r n a t i v e m u l t i p l e g r i d method, which can deal w i t h a l a r g e r c l a s s of problems than the above method. The

(21)

implementation of t h i s new method i s given i n TEXT 3.3.2 (Chapter 3).

1.4. SCOPE OF THE STUDY

In t h i s t h e s i s we are concerned mainly with m u l t i p l e g r i d methods f o r the f a s t s o l u t i o n of equations (1.4). In Chapter 2 various m u l t i p l e g r i d methods are studied f o r these equations. For these i t e r a t i v e methods the reduction f a c t o r s , which determine the r a t e of convergence, are derived using the c o l l e c t i v e l y compact operator theory by ANSELONE [1] and ATKINSON [ 3 ] . T h e o r e t i c a l and numerical i n v e s t i g a t i o n s show that m u l t i p l e g r i d

2

methods give the s o l u t i o n of (1.4) i n 0(N ) operations as N •* •», whereas 2

other i t e r a t i v e schemes take at l e a s t 0(N l o g N) operations. In p r a c t i c e t h i s r e s u l t s i n algorithms f o r the s o l u t i o n of (1.4) that are s i g n i f i c a n t -l y more e f f i c i e n t than the other schemes.

For the automatic s o l u t i o n of Fredholm i n t e g r a l equations of the second kind a new code, c a l l e d solve int eq i s presented i n Chapter 3. The l i n e a r system obtained from the d i s c r e t i z a t i o n of the i n t e g r a l equation i s

i t e r a t i v e l y solved by a m u l t i p l e g r i d method. For a v a r i e t y of problems the performance of solve -int eq i s compared with Atkinson's program iesimp [ 2 ] . Using the number of k e r n e l evaluations as a basis f o r comparison, the cost of the new algorithm i s about 2/3 the cost of the algorithm iesimp; and i t appears to be equally w e l l or even more r e l i a b l e .

In Chapter 4 we discuss the numerical s o l u t i o n of a two-dimensional D i r i c h l e t problem f o r Laplace's equation. We use the c l a s s i c a l approach of representation of the s o l u t i o n by means of a doublet d i s t r i b u t i o n on the boundary of the domain. From the boundary c o n d i t i o n we obtain a Fredholm i n t e g r a l equation of the second k i n d f o r the doublet d i s t r i b u t i o n . We introduce a m u l t i p l e g r i d method which makes use of a sequence of g r i d s , that are generated by d i v i d i n g the boundary i n t o an i n c r e a s i n g number of smaller and smaller panels. On these g r i d s the doublet d i s t r i b u t i o n i s assumed to be constant over each panel. Assuming the boundary to s a t i s f y a c e r t a i n smoothness c o n d i t i o n we prove that the reduction f a c t o r of the

mul-1+a

t i p l e g r i d process i s l e s s than Ch , where h and a are a measure of the mesh-size and of the smoothness of the boundary, r e s p e c t i v e l y . We i l l u s t r a t e

t h i s t h e o r e t i c a l convergence r e s u l t w i t h the c a l c u l a t i o n of p o t e n t i a l flow around a Karman-Trefftz p r o f i l e .

In Chapter 5 we deal w i t h the nonlinear problem concerning the r o t a t i n g flow due to an i n f i n i t e d i s k performing t o r s i o n a l o s c i l l a t i o n s at an angular

(22)

v e l o c i t y fi s i n mt. This problem i s described by the Navier-Stokes and c o n t i n u i t y equations. By means of the Von Karman s i m i l a r i t y transformations the equations are reduced to a nonlinear system of p a r a b o l i c equations with p e r i o d i c conditions i n time. Applying these conditions we can reformulate the problem as a nonlinear operator equation of the type (1.3). This equa-t i o n i s solved by a m u l equa-t i p l e g r i d meequa-thod.

This t h e s i s i s based on the f o l l o w i n g p u b l i c a t i o n s of the author:

[A] Multiple grid methods for the solution of Fredholm integral equations of the second kind.

This paper has been w r i t t e n j o i n t l y with P.W. Hemker and has been published i n Mathematics of Computation 36 (1981), pp.215-232.

In f a c t , Chapter 2 of t h i s t h e s i s i s an adapted v e r s i o n of t h i s paper.

[B] The automatic solution of Fredholm equations of the second kind. This paper w i l l appear i n the SIAM Journal on S c i e n t i f i c and S t a t i s -t i c a l Compu-ting and i s based on Chap-ter 3.

[C] On the regularity of the principal value of the double layer potential. This paper contains the t h e o r e t i c a l r e s u l t s of Section 4.1, which are applied to the aerodynamic problem of c a l c u l a t i o n of p o t e n t i a l flow around 2-D bodies. I t w i l l appear i n the Journal of Engineering Mathematics.

[D] Analytical and numerical results for the non-stationary rotating disk flow.

In t h i s paper the a n a l y t i c a l r e s u l t s of Section 5.2 and the f i n i t e d i f f e r e n c e schemes given i n Section 5.3 have been published. I t has appeared i n the Journal of Engineering Mathematics JJ3 (1979), pp.173-191.

[E] Application of multigrid methods for integral equations to two problems from fluid dynamics.

This paper i s based on Section 4.5 and Chapter 5. I t was presented at the NASA-symposium " M u l t i g r i d methods", October 21-22, 1981, Moffett F i e l d , C a l i f o r n i a . I t has been published i n the NASA Conference P u b l i c a t i o n 2202.

(23)

REFERENCES TO CHAPTER 1

[1] ANSELONE, P.M., Collectively compact operator approximation theory, P r e n t i c e - H a l l , Englewood C l i f f s , N.J., 1971.

[2] ATKINSON, K.E., An automatic program for linear Fredholm integral equations of the second kind, ACM Transactions on Mathematical Software 2:, 1976, pp.154-171.

[3] ATKINSON, K.E., A survey of numerical methods for the solution of Fredholm integral equations of the second kind, SIAM, P h i l a d e l -phia, Pa., 1976.

[4] BERGER, M.S., Nonlinearity and functional analysis, Academic Press, New York, San Francisco, London, 1977.

[5] BRAKHAGE, H., Über die numerische Behandlung von Integralgleichungen nach der Quadraturformelmethode, Numerische Mathematik 2, 1960, pp.183-196.

[6] BREBBIA, CA., Ed., Recent advances in boundary element methods, Pentech Press, London, Plymouth, 1978.

[7] BREBBIA, CA., Ed., New developments in boundary element methods, CML P u b l i c a t i o n s , Southampton, 1980.

[8] HACKBUSCH, W., Die schnelle Auflösung der Fredholmsahen Integralglei-chung zweiter Art, Beiträge zur Numerischen Mathematik 9^ 1981, pp.47-62.

C9] IKEBE, Y., The Galerkin method for the numerical solution of Fredholm integral equations of the second kind, SIAM Review Ut, 1972, pp.465-491.

[10] JASWON, M.A. & G.T. SYMM, Integral equation methods in potential theory and elastostatics, Academic Press, London, 1977.

[11] TE RIELE, H.J.J., Ed., Colloquium Numerical treatment of integral equations, MC-Syllabus 41, Mathematisch Centrum, Amsterdam, 1979. [12] VAN WIJNGAAEDEN, A. et a l , Eds, Revised report on the algorithmic

language ALGOL 68, Springer V e r l a g , New York, Heidelberg, B e r l i n , 1976.

(24)

(25)

CHAPTER 2

MULTIPLE GRID METHODS FOR THE SOLUTION OF

FREDHOLM INTEGRAL EQUATIONS OF THE SECOND KIND

2.1. INTRODUCTION

M u l t i p l e g r i d methods have been advocated by BRANDT [5,6] f o r s o l v i n g sparse systems of equations that a r i s e from d i s c r e t i z a t i o n of p a r t i a l d i f f e r e n t i a l equations. Convergence and computational complexity of such m u l t i -p l e g r i d techniques have been studied since by HACKBUSCH [7,8] and WESSELING [12,13]. We intend to show that m u l t i p l e g r i d methods can also be used ad-vantageously f o r the non-sparse systems that occur i n numerical methods f o r i n t e g r a l equations.

In [10] we have a p p l i e d the m u l t i p l e g r i d technique to the s o l u t i o n of Fredholm i n t e g r a l equations of the second k i n d

1

(2.1.1) f ( x ) - fe(x,y)f(y)dy = g ( x ) , X 6 [ 0 , 1 ] , 0

where g belongs to a Banach space X. At the same time, HACKBUSCH [7] also used a m u l t i p l e g r i d technique f o r these problems. Moreover, he gave a proof of convergence. In t h i s proof he assumed the operator K, associated w i t h the kernel fc(x,y) to map from X i n t o a "smooth" subspace X c X, which has a stronger topology. I n t h i s chapter, f o r Hackbusch's method we give another proof, which f i t s i n t o the t h e o r e t i c a l framework developed by ANSELONE [1] and ATKINSON [2,3] f o r Fredholm equations. We assume that K i s compact from X i n t o X. In contrast to Hackbusch's a n a l y s i s , t h i s approach enables us to consider a l s o Nystrom i n t e r p o l a t i o n as a p e r m i s s i b l e i n t e r p o l a t i o n method. In a d d i t i o n , we introduce a new m u l t i p l e g r i d method f o r Fredholm i n t e g r a l equations, which can deal w i t h a l a r g e r c l a s s of problems than the method proposed by Hackbusch.

In 1978 STETTER [11] introduced the Defect C o r r e c t i o n P r i n c i p l e f o r the formulation of various i t e r a t i v e methods. We s h a l l apply t h i s p r i n c i p l e

(26)

because i t also appears to be an expedient t o o l to formulate m u l t i p l e g r i d techniques.

In Section 2.2 we c o l l e c t some r e s u l t s from papers by ATKINSON [2,3] and PRENTER [ 9 ] . In Section 2.3 we cast the i t e r a t i v e schemes of BRAKHAGE [4] and ATKINSON [2,3] into the context of the Defect C o r r e c t i o n P r i n c i p l e and m u l t i p l e g r i d i t e r a t i o n . Furthermore we give a proof of convergence of the m u l t i p l e g r i d schemes with Nystrom i n t e r p o l a t i o n . In Section 2.4 we t r e a t other i n t e r p o l a t i o n methods and we extend the i t e r a t i v e schemes of Section 2.3 f o r subspaces X of X of f i n i t e dimension N . These schemes are

* P P used as a basis f o r the c o n s t r u c t i o n of a general algorithm f o r the s o l u t i o n of Fredholm equations of the second k i n d . This algorithm i s more e f f i c i e n t than the algorithms by BRAKHAGE [4] and ATKINSON [2,3] because these schemes

3 2

take 0(Np) and 0Q$^ log N^) operations, r e s p e c t i v e l y , whereas the m u l t i p l e g r i d schemes r e s u l t i n OCN^) operations. In Section 2.5 we i l l u s t r a t e the t h e o r e t i c a l r e s u l t s of the previous sections by some numerical examples and we comment on the computational complexity.

2.2. BASIC ASSUMPTIONS AND RESULTS

Equation (2.1.1) can be w r i t t e n s y m b o l i c a l l y as

(2.2.1) Af = g, g e X,

where X i s a Banach space and A = I - K, w i t h I the i d e n t i t y operator on X and K the l i n e a r operator associated w i t h the kernel fe(x,y). A i s assumed to have a bounded inverse on X. We s h a l l discuss the convergence of a sequence of approximations to the unique s o l u t i o n of (2.2.1).

Let Xp, p = 0,1,2,..., be f i n i t e - d i m e n s i o n a l subspaces of X and l e t T , p = 0,1,2,..., be a bounded p r o j e c t i o n operator from X onto X^, i . e . T f = f f o r a l l f e X . We need the f o l l o w i n g assumptions f o r {X } and {T }:

P P P P A l . X_ c X. c ... c X c . . . c X, 0 1 p A2. l i m llf-T fH = 0 f o r a l l f e X. p p-X» c LEMMA 2.2.1. •C. = sup IIT II 1 p>6 P

(27)

PROOF. The lemma follows from the p r i n c i p l e of uniform boundedness, see ATKINSON [ 3 ] , p.18. •

The sequence {X^} i s thought to be associated w i t h i n t e r p o l a t i n g s p l i n e functions on a sequence of p a r t i t i o n s {n } of the i n t e r v a l [0,1] w i t h mesh-s i z e mesh-s {hp}. We amesh-smesh-sume that + 0 amesh-s p + ». Corremesh-sponding with the mesh-sequence {II } we approximate K by a sequence of operators {K }, K : X •* X. Analogous

P P P

to A = I - K, we also w r i t e A = I - K . I n the context of m u l t i p l e g r i d

P P r e

i t e r a t i o n , the subscript p i s c a l l e d level.

We use the f o l l o w i n g assumptions on K^, p = 0,1,2,..., A3. K i s a l i n e a r operator: X •+ X.

P

A4. {Kp} i s a c o l l e c t i v e l y compact family of operators, i . e . , the set S = {Kjjf I P s 0 an d II f II S 1} has compact closure i n X.

A5. l i m UK f-Kfll = 0 f o r a l l f e X. p-X» P

A6. K = K T . P P P

LEMMA 2.2.2. From the Assumptions A3 - A5 follow: ( i ) K is compactj

( i i ) the sequence {K } is uniformly bounded, i.e. C„ = sup UK II < °°;

P 1 p>0 P

( i i i ) l i m II (K-K )Mil = 0 for any compact operator M: X->-X; p-x» P

( i v ) let a = sup sup II (K-K )K„H. then l i m a = 0 . P q>p l>0 q S- p-«o P PROOF. See ATKINSON [ 3 ] , p.96 and p.138. •

LEMMA 2.2.3. Let the finite dimensional subspace X^ c x be sufficiently large (i.e. the mesh-width of the coarsest discretization is sufficiently small). From the existence of a bounded inverse of A = I - K and the Assump-tions A3 - A5 follow:

( i ) (I-K )_1 exists on X for p > 0 and C, = sup II (I-K )~'ll < <»; P p>0 P

( i i ) l l f - f II < C,llKf-K f i t w h e r e f is the solution of (2.2.1) and f of

p i p P

(2.2.2) tt-Kp)fp =

8-PROOF. See ATKINSON [ 2 ] , p.18. •

(28)

LEMMA 2.2.4. From the Assumptions A2 - A6 follow:

( i ) For any compact operator M on X i n t o X; l i m II (I-T )Mll = 0: p-x» p

( i i ) if XQ is sufficiently large, then (I-TMC^) 1 exists on X for p > 0 and C. = sup II (I-T K )_ 1l l < °°.

4 p>0 P P Let f be a s o l u t i o n of

P

(2.2.3) (I-T K ) f = I g . P P P P

According to Lemma 2 . 2 . 4 ( i i ) , f e x i s t s and i s unique; i t follows from (2.2.3) that f e X .

P P LEMMA 2.2.5. Let

then

b = sup sup II (I-T )K It 3 P q>p l>0 q

l i m b = 0. p

PROOF. Let >f = { Kpf | p > 0 and II fII < 1}. By Assumption A4, V has compact closure i n the Banach space X. R e f e r r i n g to Lemma 1 of ATKINSON [ 3 ] , p.53, the convergence of T^f i s uniform on compact subsets of X. Then

sup II (I-T )zll -»-0 f o r q -+• <=° ze"f q

and therefore b^ 0 as p •* <». •

LEMMA 2.2.6. Let the subspaae X_ c x be sufficiently large; then

I f -T fll S C,C.llKf-K f l l . p p 1 4 p

where f is the solution of (2.2.1).

PROOF. See PRENTER [ 9 ] , Theorem 6.3. •

As a consequence of Assumption AÍ, the f o l l o w i n g lemma i s t r i v i a l .

LEMMA 2.2.7. Let q < p, i.e. dim(Xq) < d i m ( Xp) ; then

T T = T . p q q

(29)

2.3. ITERATION SCHEMES WITH NYSTROM INTERPOLATION

In t h i s s e c t i o n we use the Defect Correction P r i n c i p l e (cf. STETTER [11]) to formulate a c l a s s of i t e r a t i v e methods f o r the s o l u t i o n of (2.2.2), This equation i s w r i t t e n as

(2.3.1) Apfp = g, g e X,

with Ap = I - Kp. The Defect Correction P r i n c i p l e defines the f o l l o w i n g i t e r a t i v e process:

(2.3.2)

f n =0, _p,0

f . . = B g + (I-B A ) f .. p,i+l p6 p p p , i

Here Bp denotes some approximate inverse of A , which i s b i j e c t i v e and continuous i n X. The s o l u t i o n f of (2.3.1) i s a f i x e d point of (2.3.2) and

(2.3.2) w i l l converge to f i f the reduction f a c t o r

II I-B A II < 1. P P

Several well-known i t e r a t i v e schemes f o r s o l v i n g Fredholm i n t e g r a l equations of the second kind can be formulated w i t h i n t h i s framework. The i t e r a t i v e scheme of BRAKHAGE [4] i s obtained by taking the f o l l o w i n g approximate inverse

(2.3.3) B( I ) = I + (I-K J "1 K . P P-1 P

Here we n o t i c e that the operator (I-K^) ', q > 0, as a mapping on X i n t o X describes the process of d i s c r e t i z a t i o n , s o l u t i o n of the d i s c r e t e problem

( i . e . the s o l u t i o n of a square l i n e a r system) and subsequent Nystrom i n t e r -p o l a t i o n (see e.g. [10]). Other kinds of i n t e r -p o l a t i o n are treated i n the next s e c t i o n .

The second i t e r a t i v e scheme of ATKINSON [ 2 ] , p.19, a r i s e s from

(2.3.4) B<2 ) = I + ( I - KQ) 1 Kp.

The reduction f a c t o r s of the corresponding i t e r a t i v e processes are estimat-ed i n the f o l l o w i n g theorem.

(30)

THEOREM 2.3.1.

( i ) l l l - B ^ A II + 0 as p •*• °°; P P

( i i ) l l l - B( 2 )A II < C(Xn) a s p + » ,

p p 0

C(XQ) < 1 for XQ sufficiently large.

PROOF. ( i ) S u b s t i t u t i o n of the e x p l i c i t expressions f o r A^ and Bp'^ y i e l d s : I - B ^ ° An = I - {1+ (I-K . ) " ' K }(I-K ) = P P P-1 P P = K - (I-K K (I-K ) = P P-1 P P = (I-K ,) 1 (K -K ,)K . P-1 P P-1 P

From Lemmas 2.2.2 and 2.2.3 we get the f o l l o w i n g bound f o r the norm

(I-K ) ' (K -K )K I < C,(a+a„ .) _{p-1 p p-1 p 3 p p-1}

(2)

( i i ) Analogously, we get f o r B w i t h XQ s u f f i c i e n t l y l a r g e ,

II I-B( 2 )A II < C,(a +an) = C(X„). p p 3 p 0 0

From Lemma 2.2.2 i t f o l l o w s that C(XQ) < 1 f o r a l l s u f f i c i e n t l y large XQ. •

We remark that the approximate inverses and use only two

(1) P (2) P

l e v e l s : B uses the l e v e l s p-1 and p, whereas B uses the l e v e l s 0 and

P. (3) P (4)

p. We now introduce approximate inverses B and B , which use p+1 l e v e l s . P P

They are defined r e c u r s i v e l y as f o l l o w s :

(2.3.5) ,(3) 0 (3) Bo " ( I"Ko) • I + Q(3>K , P-1 P P = 1.2,..., and (2.3.6) ( I - K0) ',

Q ^ ( i - V i

+

V '

p = 1

>

2

- - - >

w i t h Q^3), j = 3,4, p = 0,1,2,..., given by

(31)

Q

(j)

=

V

( I

_

B

( J )

A

)

m B ( j ) P m=0 P P P

f o r some p o s i t i v e integer y. From the f a c t that s a t i s f i e s the e q u a l i t y ( I - Q( j )A ) = ( I - B( j )A )Y

P P P P

we see that Q_*^ i s an approximate inverse of Ap and i t s a p p l i c a t i o n i s equivalent to the a p p l i c a t i o n of y i t e r a t i o n steps of (2.3.2) w i t h the use of the approximate inverse B ^ . I n f a c t , t h i s i s the motivation f o r t h i s

( i ) P

d e f i n i t i o n of Qp and i t i s the b a s i s f o r the a c t u a l (recursive) implemen-t a implemen-t i o n of implemen-the meimplemen-thod.

In the f o l l o w i n g d e f i n i t i o n we give a short n o t a t i o n f o r the reduction f a c t o r f o r the various i t e r a t i v e processes.

DEFINITION. = »I-B( j )A II , j = 1,2,3,4. P P P THEOREM 2.3.2. p p p - i p p p p p - i p PROOF, ( i ) By d e f i n i t i o n and Hence I - B(1) A = I - { I + A_ 1. ( I - A ) } A p p p - P py p I - B(3) A = I - {1+ [ I - ( I - B(3> A , )YJ A ' , ( I - A ) U = p P p - l p - r p - i p p = I - B( 1 )A + ( I - B( 3} A . j V ' . d - A l A = P p p - i p - r p - r p' p = I - B(1)A + ( I - B <3>A , )Y( B0)- I ) A . p P p-1 p - r p p I I I - B( 3 )A II < « I - B( 1 )A I + I I I - B( 3} A JY( J I - B(,)A I + I I I - A I } , p p P P p-1 p-1 P P P c( 3 ) , (1) + C( 3 ) \C( 0 + „K | | ) p P p-1 p P

(32)

( i i ) S i m i l a r to the case ( i ) . Now we have I - B( 4 )A = I - [ I - ( I - B ^ An , )yl A ~ \ ( A . - À + I U p p p-1 p-1 p-1 p-1 p p Hence o I - B( 1> A + ( I - B( 4> A , )YB( , )A . p P p-1 p-1 p p I I I - B(4)A I < | | I - B(1)A II + I I I - B(4} A JYI I I - B(1) A - I I I , p P P P p-1 p-1 P P ç( 4 ) (1) + Ç( 4 )Y (1) + P P P"' P

By Theorem 2.3.1 we know that + 0 as p + •; conditions f o r (1)P (4) P to vanish depend on y, UK « and ç , whereas the conditions f o r ç

(1) P P . P depend on y and only. In order to study t h i s dependence f u r t h e r we

prove the f o l l o w i n g lemma.

LEMMA 2.3.3. Let k e K. and y e U be given, let { vp | vp> 0 , p=l,2,...} be a non-increasing sequence with d = i n f - P and let {w^} be defined by

h = v i Iw = v + wY (v +k). If either ( i ) ( i i ) P P P-1 P

I

y > \, 0 < k < d < l and vi < *c» (y i 2 and l v , < i ( / i + (|)2 - | h

tfcen a C > 0 exists such that v < w £ Cv . P P P

PROOF, ( i ) We define c = 4^r; then c > 1 and the conditions on {v } are _{d-k P} w r i t t e n as

(33)

and

(2.3.8) (l+c)v. < 1.

We show that the lemma i s true f o r C = 1+c. From (2.3.8) we see Wj = Vj < < (l+c)Vj < 1. Now we show by induction that wp < ( l + c ) vp < 1 assuming that w , < (l+c)v , < 1. From (2.3.7) follows v k c + r U + i v Ky c + r p N p ' ( H O ^1 vY- ; f v , • ^ z i k ) P-1\ P-1 v / V c p wY (v +k) < ( l + c )YvY (v +k) < cv p-1 p p-1 p p w = v + wY . (v +k) < (l+c)v . p p p-1 p p / k 2 c k

( i i ) We assume v. < / (-TTT) + o _ "^r f o r some 0 < c 2 1 and we show 1 2 d (1+c)2 2d

wp < (l+c)Vp; then the lemma i s proven by t a k i n g c = 1. For any v e [0,VjD we have

v2 + k _ _ c < Q > d ( H e )2 k c

Hence (l+c)v(v+-j) < ~y+^' ^ assumption we know

W] =V] < ( l + c j v j < /{Mg£l}2 + C _ M g £ l < /c < l .

Now we show by i n d u c t i o n that Wp < (l+c)Vp < 1, assuming that w . < (l+c)v , < 1: p-1 p-1 ( H O V p _1( vp_]+ | ) P ( l+c ) V p _]( vp_1+ ^ i k ) ( 1 + c )

V > (

1 +

f )

< c

>

* p' ( l + c )YvY ƒ 1 + A.) < c, y - 2,3 P 1 \ v_/

(34)

w = v + wr .(v +k) < (l+c)v . • P P P-1 P P

THEOREM 2.3.4. Let y > 2 and let satisfy

P P 1

for some 0 < d < 1. Then

( i ) if Vj < -^-[/d2+C2 - w£t& C2 as defined in Lemma 2.2.2, it follows that c<3) < 2v j and

P 1 P C4-)

( i i ) •£ƒ v, <; + dz- 1} it follows that c ^ 2v . " 1 2d P P

PROOF, ( i ) Let {w }be defined as i n Lemma 2.3.3 w i t h k = C„ = supIlK II , then

P 2 p^O P

i t f o l l o w s from the proof of Lemma 2.3.3 that w < 2v ,

P P

Therefore we show £ w by i n d u c t i o n : from the d e f i n i t i o n of £ ^ we

P P P derive - Vj Wj and by Theorem 2.3.2 Y (3) (1) + (3) ( (1) + ||K ||) < v + WY (v + c„) = w . p P P-1 P P P P-1 P 2' p ( i i ) S i m i l a r l y , with {w } defined as i n Lemma 2.3.3 w i t h k = 1, we prove

(4) P C < w and hence P P C( 4 ) < w < 2v . P P P (3)

REMARK. I f B i s defined with y = 1, then a s i m i l a r proof y i e l d s t h a t , P

f o r any decreasing sequence fvp } w i t h v

sup UK II = C, < d = i n f — 2 - < 1, p>0 p 2 p vp - l

(35)

C0 ) * v , P P we have

By f ^ we denote the r e s u l t of a a p p l i c a t i o n s of the Defect Correction

P,a (.}

Process on l e v e l p with approximate inverse B , j = 3,4, when we take zero as the i n i t i a l approximant.

With the a i d of the previous theorem and Lemma 2.2.3 the f o l l o w i n g theorem i s immediate.

THEOREM 2.3.5 (Approximation theorem). Under the hypotheses of Theorem 2.3.4 the multiple grid process yields approximate solutions for which the follow-ing error estimate holds:

Hf_f( j ) | l < c ||K f-Kf| + ( 2 dp" ' v1)al l f II, j = 3,4,

p,a 3 p V p J ' '

where f and f are the solutions of (2.2.1) and (2.2.2), respectively. PROOF. The proof follows immediately from

l l f - f ( j )l l < l l f - f II + I f - f ( j )l l . • p,o p p p,o

2.4. ITERATION SCHEMES WITH PROJECTION INTO FINITE DIMENSIONAL SUBSPACES

In t h i s s e c t i o n we expand the technique used i n Section 2.3 to f i n d the s o l u t i o n i n Xp of the equation (2.2.3):

(2.4.1) A f = g , g e X , P P BP BP P

where A = I - T K i s a mapping on X i n t o X. We assume that Xn i s s u f f i

-P -P -P v v 5 _ j 0

c i e n t l y large such that (I - T Kp) e x i s t s f o r a l l p > 0.

Analogous to the approximate inverse of Ap i n the previous s e c t i o n , we now introduce:

(36)

B B (1) P (2) I + T A ',T ,T K , P p-1 P-1 P P = I + T A_ TnT K , p 0 0 p p' B0 S( 3 ) 1 P A0 T0 ' I + T Q ^ T ,T K , psp - l P-1 P P B0 S( 4 ) P A0 0' = 1 - 1 1 ,T t T q ^ ' l , ( I - K , + T K ) T , p p - i p

P

V-I

p - i p - i P P P with ;(j) Y - l _{I j = 3,4,} m=0 f o r some p o s i t i v e integer -y.

The operators B^^, j = 1,2,3,4, are a l l mappings on X i n t o X. The s o l u t i o n fp e X^ of (2.4.1) i s approximated by a defect c o r r e c t i o n process of the form

(2.4.2)

'p,0 = 0,

f - ^ , = B g + ( I - B A ) f .. p,i+l p6p p p p , i

We n o t i c e that B ^ and B ^ y i e l d i t e r a t i v e processes that are equivalent

P P

r e s p e c t i v e l y with the "One Step Method" and the " M u l t i Grid Method" discus-~(2)

sed i n HACKBUSCH [ 7 ] . B ' y i e l d s an i t e r a t i v e process analogous to p ~(4)

Atkinson's method, whereas B y i e l d s a new m u l t i p l e g r i d method with bet-_{t e r convergence p r o p e r t i e s than B}(3) Analogously to Section 2.3, but r e s t r i c t i n g the domain of the operators to X , we see that here Q ^ : X •* X i s an approximate inverse of A : X •* X

p' P P P P P P

and the a m p l i f i c a t i o n operator on Xp i n t o X^ of a defect c o r r e c t i o n step w i t h i s

P

I - Q( j )A = ( I - B( j )A )Y.

~ ( i ) • • • Thus, one a p p l i c a t i o n of i s equivalent w i t h the y times a p p l i c a t i o n of

(37)

Q( j ) = CI - ( I - B( j )A )Y] A_ 1. P P P P

The convergence of the process (2.4.2) depends on the L i p s c h i t z constant o the operator I - B A as a mapping X •+ X . Therefore i t s reduction f a c t o r

p p P P i s given by

»Tp(I-SpAp)||.

This reduction f a c t o r i s studied i n the remainder of t h i s s e c t i o n .

THEOREM 2.4.1.

( i ) l l Tp( I - B( 1 )Ap) l l -y 0 as p -»• »;

( i i ) l l Tp( I - B( 2 )Ap) i l < C(XQ) as p C ( XQ) < 1 for XQ sufficiently large.

PROOF, ( i ) S u b s t i t u t i o n of the e x p l i c i t expressions f o r B^'^ and A y i e l d s P P w i t h % = p-1:

T ( I - B( 1 )A ) = T X7'{A -T T A }K .

p p p p I I i p p p

The expression between braces i s r e w r i t t e n as

(I- YTI(I- VTP VKP} -Therefore we have IIT ( I - B( 1 )A )ll < P P P S I T H A T ' i r i ( I - T J K I + p I L % p

+ IIT£1 {II ( I - T ) K l + II (I-TP)K^H«Kpll +

+ I T I I (K -K)K II + I T IIII(K-K„)K 1 } ] .

P P P P V p J

Using the lemmas 2.2.1 to 2.2.5 we obtain the proof of ( i ) by the same arguments as used f o r the proof of lemma 2.3.1.

(38)

D + oo II d-T )K II , II (I-T )KJI and II (K -K)K II vanish, whereas the other terms

* ' v p p p 0 p p

tend to a constant value depending on X~. •

DEFINITION.n( j ) = IIT ( I - B( j )A )H, j = 1,2,3,4. p p p p

THEOREM 2.4.2.

p p p-1 p p

PROOF. We use the n o t a t i o n = I - B ^ A , j = 1,3,4. From Assumption A6 „ P P P

and the d e f i n i t i o n s of A and B , J = 1,3,4, i t i s c l e a r t h a t : P P and Hence and a l s o T A = A T , P P P P I = B(^ T , j = 1,3,4. p p p p' I M'^ = M ^ ' l , j = 1,3,4, P P P P J p p p ( i ) From the d e f i n i t i o n of Bp ' ^we ge t> w i t h I = p-1: M( 1 ) = I - A -T AT'T„T K A = T K -T AT'T.T K A . p p p i I p p p p p P I I P P P

Use lemma 2.2.7 and the r e l a t i o n A ^ T ^ = T^A^1 to prove that (1) — ] ~

MV

'

= T K - A. T T K A . p p p H i p p p

These r e l a t i o n s are used to prove that

M( 3 ) = I - U + T (I-MP^ JXI'T.T K }A = p p' 1 1 I p p J • T M„(3) p P * = M(1> + T MJ3> P P «• Y

(39)

Hence

I t i s e a s i l y v e r i f i e d that M = T M and KT ^ = T M . By means of

P P P P P P ' Lemma 2.2.7 we get M( 3 ) = M( 1 ) +MP )Y( T K - M( 1 )) . P P % P P P n(3), n(1) +

,

( 3

r(n

( 1 ) + H T l l l l K l l ) . P P P-1 P P p ( i i ) M^4 ) = I - + T ( I - M .( A ) Y) A T ' T . ( A . + T K )T }A = p p i p p a i nv & p py p p = M( 1 ) + T M.(4)YT. (1 + A T V T K )A T . p p J l I I I P P P p

From t h i s expression we conclude that M = T M ^ = M ^ T . Using Lemma

(4) (4) p p p p p 2.2.7 we obtain X M„ ' = M„v ' and p I I M( 4 ) = M( 1 ) + M <4)V M(1> )T . P P * P P Hence n(4> , n ( 0+n( 4 )Y (1) + | | T | | P P p-1 P P

THEOREM 2.4.3. Let y > 2 and let satisfy np'5 < vp = dp - Iv] ; > for some 0 < d < 1; then

( i ) if v, < •^•{/d2+C2c2 - C^C^, it follows that np3) £ 2 vp J ( i i ) if Vj < -^-{/dz+c2 - C,}, i t follows that < 2vpl

where Cj and C£ are defined as in Section 2.2.

PROOF• ( i ) Use Lemma 2.3.3 w i t h k = and Theorem 2.4.2. ( i i ) Analogously w i t h k = C,. •

By f ^ we denote the r e s u l t of a a p p l i c a t i o n s of the Defect C o r r e c t i o n

P'a . ~ ( i )

Process on l e v e l p w i t h approximate inverse Bp , j = 1,3,4, when we take zero as the i n i t i a l approximant.

THEOREM 2.4.4 (Approximation theorem). Under the hypotheses of Theorem 2.4.3 the multiple grid process yields f " , for which the following error

esti-p, c mates hold:

l l f - f( j )l l < llf-T fll + C.C.IlKf-K fll + ( 2 dP _ 1v f l l f II, p,a p 1 4 p \> p '

(40)

where f and f are t/ze solutions of (2.2.1) and (2.2.3)., respectively.

PROOF. For j = 1,3,4 we have

l l f - f ^ l l < llf-T fil + Il T f - f « + II f - f ^ l l p,a p P P P P»o

and the proof follows from Lemma 2.2.6 and Theorem 2.4.3. •

We notice that the usual d i s c r e t i z a t i o n methods e a s i l y s a t i s f y the f i r s t condition of Theorem 2.4.3 as i s i l l u s t r a t e d i n Section 2.5. The other condition of Theorem 2.4.3, which requires an upperbound on v., i s essen-t i a l l y a requiremenessen-t on essen-the coarsesessen-t d i s c r e essen-t i z a essen-t i o n used i n essen-the mulessen-tiple g r i d algorithm. This condition i s also discussed i n the next section.

2.5. NUMERICAL RESULTS

In t h i s section we i l l u s t r a t e the t h e o r e t i c a l convergence r e s u l t s from the previous sections and we make some remarks about the computational com-p l e x i t y of the various methods. We s h a l l only show numerical r e s u l t s obtain-ed with the methods that appear to be the most e f f i c i e n t . These methods are

(2) . ~(3) defined by the approximate inverses B (Atkinson's method), B

~(4) P P

(Hackbusch's method) and B (a new method with better convergence proper-P

t i e s ) .

As an example, the i n t e g r a l equation 1

(2.5.1) f(x) - A cos(Trxy)f(y)dy = g(x) 0

i s solved for various values of the parameter X (cf. HACKBUSCH [7] who gives

r e s u l t s for the same equation); g(x) i s chosen such that

f(x) = eXcos(7x).

The operators Kp are defined by means of the repeated trapezoidal r u l e :

X.

VI

f(x) = f w fe(x,x.)f(x.),

P j=0 J

where the nodal points {x.} are uniformly d i s t r i b u t e d ( xo = 0, xN = 1 and

J P - i the weights {w^ } are given by { £ hp, hp, hp hp, J hp} , with hp = (Np) . The

(41)

p r o j e c t i o n operators are defined by piecewise l i n e a r i n t e r p o l a t i o n at the nodal points {x.}. The various g r i d - l e v e l s are r e l a t e d by N = 2N

J P p-1 For the operators {K } and {T } we know ( c f . ATKINSON [3] and PRENTER

P P [9]) t h a t , i f f i s the s o l u t i o n of (2.5.1): (2.5.2) IlK f-Kfll = 0 ( h2) , P P (2.5.3) IfT f - f II = 0 ( h2) , P P (2.5.4) a = 0 ( h2) P P and (2.5.5) b = 0 ( h2) f o r p -* », P P

with ap and bp defined as i n Lemmas 2.2.2 and 2.2.5, r e s p e c t i v e l y . Using these estimates, we e a s i l y derive (see the proof of Theorem 2.4.1)

, v = Ch2 ,. p P p-1

Because the successive mesh-sizes are r e l a t e d by h = h„2 p we have P 0 (2.5.6) n 0) < v = 4Ch24"P.

p p 0

Comparing t h i s expression with the assumption on i n Theorem 2.4.3 we P ( i )

see that d = 1/4. In the same theorem conditions on n j are formulated for the m u l t i p l e g r i d methods to converge. Comparing these conditions we conclude that the c o n d i t i o n on n, i n the process defined by B ^ i s

Pru / *3 \ independent of C0 = sup IlK II , whereas i n the process defined by B the

. m p-o p p

c o n d i t i o n on n, ' becomes stronger as sup IlK II increases. I n Figure 2.1

. p>0 P „,(3) ~(4)

we sketch the regions of convergence induced by Bp and Bp as derived from Theorem 2.4.3 w i t h d = 1/4 and Y = 2.

Hence, from Theorem 2.4.3 one may expect that both m u l t i p l e g r i d methods y i e l d s i m i l a r r e s u l t s as IlKll & I whereas they d i f f e r f o r IlKll » 1. For the i n t e g r a l equation (2.5.1) IlKll » 1 holds f o r X» 1.

In Tables 2.5.1 -2.5.3 we give the observed reduction f a c t o r n(N ;Nn) = [llf . ^ , - f . l / l f , - f nI U1 / l,

p" 0' p,i+l p , i p , l p,0

for the i t e r a t i v e methods defined by B ( 2 ) , ]J(3)an(j B , r e s p e c t i v e l y , p P P

(42)

Figure 2.1. The m u l t i p l e g r i d convergence regions. The coarsest g r i d reduction f a c t o r rij versus

C~ = sup IK II ; C. = sup IIT II . 2 p>0 P 1 psb P

w i t h y = 2. The dependence of liCN^;!^) on Np, the number of mesh i n t e r v a l s i n the f i n e s t g r i d , and on NQ, the number i n the coarsest g r i d , i s shown. The value of i i s s u i t a b l y chosen and II-II denotes the maximum norm. From Table 2.5.1 we see that the reduction f a c t o r s of Atkinson's method tend to a constant value as Np °°. As was expected, i t decreases as NQ increases. In the case of convergence, the Tables 2.5.2 and 2.5.3 a s y m p t o t i c a l l y show s i m i l a r r e s u l t s . However f o r l a r g e r values of X the new m u l t i p l e g r i d

method needs fewer s u b i n t e r v a l s i n the coarsest g r i d . The quotients n(Np;NQ)/ri(Np_j ;NQ) approximate the value d = 1/4, which i s i n agreement

(43)

> o X N \ p \ 2 4 8 16 32 64 4 8 .23 1 0_ 1 .28 1 0_ 1 .58 10~2 1 16 .30 10~ .72 10~2 -2 .15 10 32 .30 1 0_ 1 .76 10"2 .18 10"2 .38 10~3 64 .30 10~ .78 10_ .19 10"2 .47 10~3 -4 .83 10 128 .30 1 0_ 1 .78 10"2 .19 10"2 .52 10~3 .10 i o "3 -4 .24 10 4 .11 i o+1 8 .16 1 0+ 1 .18 10° 10 16 .17 1 0+ 1 .22 10° .36 1 0_ 1 32 .17 1 0+ 1 .23 10° .45 1 0_ 1 .86 10~2 64 .17 1 0+ 1 .24 10° .48 1 0_ 1 .11 i o_1 .21 10~2 128 .17 1 0+ 1 .24 10° .48 1 0_ 1 .11 i o_1 .27 10~2 .38 10"3 4 .64 1 0+ 1 8 .11 i o+ 2 .14 1 0+ 1 .40 10° 100 16 .14 1 0+ 2 .16 1 0+ 1 .40 10° 32 . 1 5 ! 0+ 2 .16 1 0+ 1 .42 10° .99 1 0_ 1 64 .15 1 0+ 2 .16 1 0+ 1 .45 10° .15 10° .33 1 0- 1 128 .15 1 0+ 2 .16 1 0+ 1 .49 10° .16 10° .41 1 0_ 1 .68 10~2

Table 2.5.1. Reduction f a c t o r s f o r the two-grid method (2)

defined by B (Atkinson's method). P

(44)

V o X N \ P \ 2 4 8 16 32 64 4 8 .31 1 0_ 1 .98 10~ .94 i o "2 1 16 .24 10"2 .24 10 2 .23 i o "2 32 .62 10"3 .62 10"3 .62 i o '3 .62 10"3 64 . 14 10~3 . 14 i o "3 . 14 i o "3 .14 IO"3 .14 10" 3 128 -4 .35 10 .35 i o "4 .35 i o "4 -4 .35 10 .34 10" 4 .35 10" 4 4 .32 10° 8 .12 1 0+ 1 .10 10° 10 16 .42 1 0+ 1 .12 10° .25 i o '1 32 .20 1 0+ 3 .18 10 .13 i o "1 .62 10~2 64 .23 1 0+ 6 .91 i o "2 .24 10~2 .23 IO"2 .19 10" 2 128 .40 1 0+ 1 2 .46 i o "3 .57 10"3 .52 10"3 .53 10" 3 .51 10 -3 4 8 .43 1 0+ 1 .11 1 0+ 4 .11 i o+ 1 100 16 +7 .66 10 .77 i o + 2 .29 10° 32 .17 1 0+ 1 7 .79 i o+ 4 .51 i o+ 1 .10 10° 64 .82 1 0+ 3 4 .46 i o+n .15 i o+ 3 .33 10° .29 10" 1 128 .80 1 0+ 7° .86 1 0+ 2 3 .96 i o+ 7 .43 10° .24 10" 1 .85 10" 2

Table 2.5.2. Reduction f a c t o r s f o r the multiple g r i d method ~(3)

(45)

V o

X N \

p \

2 A 8 16 32 64 4 .31 1 0- 1 8 .95 10"2 .94 IO"2 1 16 .23 10~2 .23

io"

2 .23 10"2

32 .62 10~3 .62

io"

3 .62 IO"3 .62

io"

-

.3899 32 - - .0694 64 - - .0194 128 - - .0050

Table 2.5.4. The r a t i o : i t e r a t i o n e r r o r a f t e r 2 sweeps/ approximation e r r o r s , I f , - f II i . e . P»2 P»°° . II f - f I _p,» Number of s u b i n t e r v a l s : NQ = 2 (a divergent i t e r a -t i o n process i s deno-ted by - ) .

(47)

X N P B< 2 ) P S( 3 ) P ~(4) P 1 16 .0003 .0001 .0001 32 .0017 .0000 .0000 64 .0075 .0000 .0000 128 .0306 .0000 .0000 10 16 .0936 .3202 .3202 32 .6088 .2000 .0692 64 ' 2.7111 .0341 .0194 128 11.1310 .0056 .0050 100 16 91.0760 34.0160 34.0160 32 563.4392 - 34.7089 64 2480.5082

-

24.3138 128 > i o4 - 0.1220

Table 2.5.5. The r a t i o : i t e r a t i o n e r r o r a f t e r 2 sweeps/ approximation e r r o r s .

As Table 2.5.4, but with N = 8.

We conclude t h i s s e c t i o n w i t h some remarks about the asymptotic compu-t a compu-t i o n a l complexicompu-ty. I compu-t i s our purpose compu-to e s compu-t a b l i s h compu-the f a s compu-t convergence of the m u l t i p l e g r i d methods rather than to construct e f f i c i e n t computer codes. Therefore, i n our implementation k e r n e l - f u n c t i o n s are re-evaluated whenever they are used. Then, the number of kernel evaluations i s equal to the number of m u l t i p l i c a t i o n s involved i n the matrix * vector computations d e f i n i n g the operation counts. The overhead costs (a.o. a r i t h m e t i c operations used f o r the i n t e r a c t i o n between the g r i d s ) are neglected. Asymptotically f o r N ->-<*>, the operation counts per i t e r a t i o n f o r the various approximate inverses are:

(1) 2 2 Bv ;: 3.25 N , B : 2.5 N , P P P P B<2>: 2 N2, B<2>: 2 N2, P P P P B( 3 ): 2 N2 2l o g N , 5( 3 ): 3 N2, P P P P P B( 4 ): 2.5 N2 2l o g N , B( 4 ): 3.5 N2. P P P P P

(48)

Here we ignored the d i r e c t s o l u t i o n on the coarsest g r i d and we applied the m u l t i p l e g r i d methods w i t h y = 2 on a l l l e v e l s .

2 .

NOTE. The number of k e r n e l - f u n c t i o n evaluations i s Np i n the l i n e a r case when they are computed once and stored. However, f o r large values of Np one may need e x t e r n a l devices and the e f f i c i e n c y of the m u l t i p l e g r i d algorithms may depend on the 1 / 0 - f a c i l i t i e s of the computer system. I n the nonlinear case or i n the case when kernel-functions are re-evaluated whenever they are used, the number of k e r n e l - f u n c t i o n evaluations per i t e r a t i o n i s equal to the number of m u l t i p l i c a t i o n s given i n the table above.

A s y m p t o t i c a l l y a l l methods need only two i t e r a t i o n s to obtain a r e s u l t (2) ~(2) of the order of the t r u n c a t i o n e r r o r , except the methods w i t h B and B

(1) ~ ( DP P

which need 0 ( l o g N ) sweeps. For the methods w i t h B and B the coarsest P P P

g r i d s t i l l has N /2 mesh-intervals; on t h i s g r i d the problem i s solved by a

P . . . 1 3 d i r e c t method (e.g. Gauss-elimination) and therefore we have to add T T N to

12 p the t o t a l computational complexity. Thus, f o r the total amount of asymptotic computational work we get the f o l l o w i n g t a b l e :

B _(1): ±K3 • 6.5 N2, B( 1> p 12 p p' p 8P2 ) ; iNo + (?(Np 108V ' 5P 2 ) B <3) : | N3+4 N22l o g N , B( 3 ) P 3 0 p 6 p' p ( 4 ): | Nn 3 + 5 N2 2l o g N , B<4> P 3 0 P 6 P' P B 1 3 2 12 Np + 5 Np> H + fl^logNp), 2 3 2

I

N3 + 6N2, 3 0 p

From these tables we see that the m u l t i p l e g r i d methods become cheaper than Atkinson's method whenever the l a t t e r needs more than three i t e r a t i o n s .

In order to get an impression of the q u a l i t i e s of the various methods we suggest to measure by experiments the f o l l o w i n g r a t i o (which shows the amount of computational work per d i g i t accuracy obtained):

Number of m u l t i p l i c a t i o n s to obtain f p,a ° N2 * 1 0l o g l l f - f II

P P»a

For the m u l t i p l e g r i d methods we choose a = 2 because t h e i r reduction f a c -tors tend to zero as N •*• ». For Atkinson's method we determine a such that

P

K i s minimal. Better methods are now characterized by a smaller K .

a a Table 2.5.6 shows f o r the m u l t i p l e g r i d methods that small values of

(49)

reasonable range of small NQ, i t seems not worthwhile to determine an optimal Nn. B( 2 ) P 5(3J P ~<4) P 7.55 (11) 2.00 (2) 16 3.46 (5) -8.42 (2) 1.86 (2) 32 64 2.73 (3) 6.20 (2) 2.97 (2) 5.95 (2) 2.20 (2) 6.10 (2)

Table 2.5.6. For problem (2.5.1) w i t h X = 100 and N = 128

P the experimental r a t i o s K^, where a i s given between parentheses; f o r t h i s problem

1 0l o g l l f - f II = -3.5.

6 p,oo

The asymptotic work estimates and the convergence property discussed i n Section 2.4 lead us to p r e f e r f o r IlKll « 1 and f o r IlKl! » 1.

P P F i n a l l y , we remark that the same m u l t i p l e g r i d techniques can be applied to nonlinear problems and the s t r u c t u r e of m u l t i p l e g r i d algorithms y i e l d s estimates f o r the approximation and t r u n c a t i o n errors i n a n a t u r a l way. A l l these features together can be used to construct an automatic program f o r s o l v i n g Fredholm i n t e g r a l equations of the second k i n d . In Chapter 3 such a program i s constructed.

REFERENCES TO CHAPTER 2

[1] ANSEL0NE, P.M., Collectively compact operator approximation theory, Englewood C l i f f s , New Jersey, P r e n t i c e - H a l l , 1971.

[2] ATKINSON, K.E., Iterative variants of the Nyström method for the numerical solution of integral equations, Numerische Mathematik 22 (1973), pp.17-31.

[3] ATKINSON, K.E., A survey of numerical methods for the solution of Fredholm integral equations of the second kind, SIAM, 1976. [4] BRAKHAGE, H., Über die numerische Behandlung von Integralgleichungen

nach der Quadratur-formelmethode, Numerische Mathematik 2 (1960), pp.183-196.

(50)

[5] BRANDT, A., Multi-level adaptive solutions to boundary-value problems, Mathematics of C o m p u t a t i o n a l ('977), pp.333-390.

[6] BRANDT, A., Multi-level adaptive techniques for singular perturbation problems, i n : Numerical A n a l y s i s of Singular P e r t u r b a t i o n

Prob-lems (P.W. Hemker & J.J.H. M i l l e r , eds), Academic Press, London, 1979.

[7] HACKBUSCH, W., Die schnelle Auflösung der Fredholmschen Integralglei-chung zweiter Art, Beiträge zur Numerischen Mathematik 9^ (1981), pp.47-62.

[8] HACKBUSCH, W., An error analysis of the nonlinear multi-grid method of second kind, Aplikace Matematiky ^6 (1981), pp.18-29.

[9] PRENTER, P.M., A collocation method for the numerical solution of integral equations, SIAM J . Numer. Anal. J_0 (1973), pp.570-581. [10] SCHIPPERS, H., Multi-grid techniques for the solution of Fredholm

integral equations of the second kind, Colloquium Numerical Treatment of I n t e g r a l Equations, MC-Syllabus 41, Mathematisch Centrum, Amsterdam (1979).

[11] STETTER, H.J., The defect correction principle and discretization methods, Numerische Mathematik 29 (1978), pp.425-443.

[12] WESSELING, P. & P. SONNEVELD, Numerical experiments with a multiple grid and preconditioned Lanczos type method, Procs of the IUTAM-Symposium on Approximation Methods f o r Navier-Stokes Problems, 1980, Lecture Notes i n Mathematics 771, pp.543-562, Springer.

[13] WESSELING, P., The rate of convergence of a multiple grid method, Procs of the Dundee B i e n n i a l Conference on Numerical A n a l y s i s ,

(51)

CHAPTER 3

AUTOMATIC NUMERICAL SOLUTION OF

FREDHOLM INTEGRAL EQUATIONS OF THE SECOND KIND

3.1. INTRODUCTION

In t h i s chapter we describe an algorithm f o r the automatic numerical s o l u t i o n of Fredholm i n t e g r a l equations of the second k i n d :

b

(3.1.1) f ( x ) - | fc(x,y)f(y)dy = g ( x ) , x e [ a , b ] . a

The a l g o r i t h m i s an improvement of Atkinson's automatic program iesimp [3] i n the sense that a new i t e r a t i v e method i s used f o r the s o l u t i o n of the system of equations that a r i s e s from the approximation of (3.1.1). Our i t e r a t i v e methods are m u l t i p l e g r i d methods that work w i t h a sequence of g r i d s of decreasing mesh-size. These g r i d s are used simultaneously to ob-t a i n an approximaob-tion ob-to ob-the conob-tinuous problem (3.1.1). The m u l ob-t i p l e g r i d methods used can be seen as extensions of Atkinson's i t e r a t i v e scheme, that uses only two g r i d s : a coarse and a f i n e g r i d . Convergence and computational complexity of the m u l t i p l e g r i d methods have been studied i n Chapter 2. The program has been w r i t t e n i n the a l g o r i t h m i c language ALGOL 68, because i n t h i s language we can e a s i l y and e f f i c i e n t l y handle the data s t r u c t u r e s and the r e c u r s i v e procedures that appear i n m u l t i p l e g r i d methods.

In Chapter 2 a d e s c r i p t i o n of our m u l t i p l e g r i d methods has been given by means of c o l l e c t i v e l y compact operators and i n t e r p o l a t o r y p r o j e c t i o n s onto subspaces of piecewise continuous f u n c t i o n s . I n Section 3.2 some relevant r e s u l t s are r e c o l l e c t e d . Based on the t h e o r e t i c a l foundation of Section 3.2, the program f o r the automatic s o l u t i o n of Fredholm equations, solve int eq, i s described i n Section 3.3. Numerical examples i l l u s t r a t i n g the method are given i n Section 3.4, where comparisons are made with Atkinson's automatic program iesimp. Further a p p l i c a t i o n s of solve int eq are described i n Chapters 4 and 5.

(52)

3.2. THEORETICAL FOUNDATIONS

In t h i s s e c t i o n we w r i t e equation (3.1.1) i n operator notation as f o l l o w s :

(3.2.1) (I-K)f = g, g e X,

where X i s a Banach space and K: X •* X the l i n e a r operator associated with the kernel fe(x,y). I t i s assumed that (I-K) has a bounded inverse on X. We approximate the s o l u t i o n of (3.2.1) by a sequence of i n t e r p o l a t i n g s p l i n e functions fp w i t h knots a t the points GP = { t ^ | a = tg< t j ... < tjj = b}. The g r i d s { G ^ } , P = 0, 1,2,.. . ,5., are constructed such that G Q C G J C... <=G^. Let h be a measure of the mesh-size defined by:

P '

h = max | t . - t . ,1. P i 1 i 1-1'

In our algorithm we take the sequence of g r i d s { GP} uniform with NP = 2PNQ, so that

B l . h = 2~Ph„. P 0

Let Xp, p = 0,1,2,...j& be the f i n i t e - d i m e n s i o n a l subspaces of i n t e r p o l a t i n g s p l i n e functions on GP and l e t Tp, p = 0,1,2,...,£ be the corresponding i n t e r p o l a t i n g operators. With t h i s choice of {X } and {Tp} the Assumptions A1-A2 of Section 2.2 are s a t i s f i e d . As i n Section 2.2 we approximate K by a sequence of approximating operators {Kp}, Kp: X -*- X, s a t i s f y i n g the Assumptions A3 - A5 of Section 2.2. Moreover, we make some assumptions about

the smoothing p r o p e r t i e s of {K } and the order of approximation of the operators {Kp} and {Tp}. Let S be the f o l l o w i n g subset of X:

S = (K f I p > 0 and I f I < 1}. P '

By Assumption A4 t h i s set has compact closure i n X. In t h i s chapter we take the stronger assumption that the functions i n S are s u f f i c i e n t l y smooth. Furthermore, f o r p -*• « we assume:

B2. sup II (K-K )fll 5 C, h", a > 0,

fe§ P 1 P

B3. sup || (I-T )f|| < C0 hB, g > 0.

(53)

The constants a and g are associated with the order of the i n t e g r a t i o n r u l e (e.g. t r a p e z o i d a l r u l e ) and of the i n t e r p o l a t i n g s p l i n e f u n c t i o n . These assumptions are i l l u s t r a t e d by the f o l l o w i n g examples. Let X be the Banach space C[a,b] of continuous f u n c t i o n s , provided with the supremum norm. Assume that the kernel f u n c t i o n fe(x,y) i s twice continuous d i f f e r e n t i a b l e with respect to x. For t h i s c l a s s of k e r n e l functions the Assumptions B2 and B3 can be v e r i f i e d f o r the f o l l o w i n g example.

EXAMPLE 1. The operator i s defined by the repeated t r a p e z o i d a l r u l e and the operator Tp by continuous, piecewise l i n e a r i n t e r p o l a t i o n . In t h i s case a = g = 2.

Analogously we obtain the f o l l o w i n g example i f fe(x,y) i s four times continuously d i f f e r e n t i a b l e w i t h respect to x.

EXAMPLE 2. Kp i s defined by the repeated Simpson r u l e and Tp by continuous, piecewise cubic i n t e r p o l a t i o n . In t h i s case a = g = 4.

EXAMPLE 3. Finite element methods for integral equations from potential theory. Let D be a simply connected f i n i t e plane region bounded by a smooth contour S with continuous curvature. S i s given by the parametric equations x = X ( s ) , y = Y ( s ) , s e [0,1]. The k e r n e l f u n c t i o n i s given by

For the i n t e r p o l a t i n g s p l i n e functions we take the piecewise constant func-tions ( i . e . s t e p - f u n c t i o n s ) . The space X must be chosen such that f o r each N i t contains t h i s c l a s s of f u n c t i o n s . Furthermore, pointwise operations should be defined on X. Therefore we choose X to be the Banach space of regulated functions (see Section 4.2), provided with the e s s e n t i a l supremum norm. The operator Tp i s defined by piecewise constant i n t e r p o l a t i o n at the midpoints. For t h i s p a r t i c u l a r example i t i s not necessary to approximate the i n t e g r a l operator because T Kf (f e X ) can be computed a n a l y t i c a l l y . However, the t h e o r e t i c a l r e s u l t s of Chapter 2 apply i f Kp i s defined by KTp. From the r e s u l t s given i n Chapter 4 i t follows that a = 1 + p, where p i s a measure f o r the smoothness of S (0 < p < 1) and (5=1.

For these examples we obtain the f o l l o w i n g estimates f o r ap and bp, which were defined i n Lemma 2.2.2 and 2.2.5, r e s p e c t i v e l y .