Numerical Analysis of Partial Differential Equations 3642110568, 9783642110566 [PDF]

Zitiervorschau

Jacques Louis Lions ( E d.)

Numerical Analysis of Partial Differential Equations Lectures given at a Summer School of the Centro Internazionale Matematico Estivo (C.I.M.E.), held in Ispra (Varese), Italy, July 3-11, 1967

C.I.M.E. Foundation c/o Dipartimento di Matematica “U. Dini” Viale Morgagni n. 67/a 50134 Firenze Italy [email protected]

ISBN 978-3-642-11056-6 e-ISBN: 978-3-642-11057-3 DOI:10.1007/978-3-642-11057-3 Springer Heidelberg Dordrecht London New York

©Springer-Verlag Berlin Heidelberg 2010 st Reprint of the 1 Ed. C.I.M.E., Ed. Cremonese, Roma 1968 With kind permission of C.I.M.E.

Printed on acid-free paper

Springer.com

CENTRO IfJTER!~JAZIONALE PIATEMATICO CSTIVO

-

I'

S.

ALBERTONI

-

ALCUNI METODI D I CALCOLO NELLA TEORIA DELLA D I F F U S I O H C D E I 1JEUTRONI"

Corso tenuto ad Ispra

d a l 3-11 Luglio 1967

ALCUIII METDDI DI CALCOLO IJELLA TCORIA A MULTIGRIJPPI DELLA DIFFUSIONF: D C I I.JCUTROEI1

PARTE la

5 1

-

-

S , Albertoni

Teoria s t a z i o n a r i a a multif!ruppi,

(I,?)

- F,lotazioni e Problemi

.

Sia ncRn un aperto l i m i t a t o e xr ( x l , x 2 , . .xn) t R seguito considereremo g l i s p a z i ~ ' ( 0 1 , H~ (Dl, H A ( R ) , i n t e r 0 >o,

( L ~ ) , ( ) ~ , ( saranno ) i g

- prodotti

.

FIel

e sa g & diretti

d i L 2 , H I , Hi. Volendo considerare principalmente problemi d i "trasmissione" r = U e s i a n o t = 30, r i = a n 1 . rispettivanente supporremo 1

",

l e frontiere d i R ,

ni:

yrs

- a n r n ans

cio& l e p a r t i d i frontiers conuni a

sono l e " i n t e r f a c c e " , e

?r,zs,

I n f i n e vr s i a l a normale, d i r e t t a verso l ' e s t e r n o , d i aQr. Pre messo c i b supponiamo assegnate l e funzioni r e a l i : D~ ( x l , tli(x), ?q ci(x) c o ~~(x),fl(xl.(~,q=

.

= 1 , 2 , . , n ; i = 1 , 2 , , , ,g) I problemi da r i s o l v e r e sono due,

e cio8: Problema A

-

Trovare l a soluzione ~ ( x d) e l sistema:

soddisfacente a u[, = 0, e a l l e condizioni d i "trasrnissioneu

valide per x c y P s , essendo ur, us l e r e s t r i z i o n i d i u a

d

l a d e r i v a t a conormale a as2

nr e

QS

r

e

rispettivanente.

-

Osservazione I Dal punto d i v i s t a f i s i c o l e (1) forniscono l a d i s t r i b u z i o n e s t a z i o n a r i a d e l f l u s s o neutronico e n t r o un r e a t t z r e R (composto da r e g i o n i d i m a t e r i a l i d i v e r s i R i ) ripartito i n g gruppi d i energia decrescente ( u i 5 il f l u s s o d e i neutron i d i iE e n e r g i a ) ove f i sono Le s o r g e n t i neutroniche e s t e r n e e: Di

Ai

un "coef f i c i e n t e " d i d i f f usione un "coef f i c i e n t e " d i assorbirnento

~i un " c o e f f i c i e n t e " d i rirnozione d a l gruppo (Ccco perch6 cl:O) B~ un " c o e f f i c i e n t e " d i f i s s i o n e .

-

i-1

a 1 gruppo i

Osservazione I1 E ' u t i l e per il s e g u i t o dare pi2 e s p l i c i t a m e ~ t e l a s t r u t t u r a formale d e l sistema (1):

Ipotesi sul problema A: 1) Positivita dei coefficienti D ~ A, ~ ,gi, Ci (supposti misurabili a limitati in R) e cioB: ~~(x))a>o, ~ ~ ( x ) s b > o~~(x))c>o , i (x)2d>o; xi, v sono coefficienti numerici >o. (i=2,3...),D Pq Inoltre si suppone C 1 (x) : 0. 2) Ellitticit; uniforme e cio8: 4 {A ; p=1,2.. .n, A 6 Rl,x c T ) P P esistono m, H>o tali che:

per ogni i=1,2...g. Inoltre a sari abbastanza regolare da garantire che tutte le operazioni di traccia abbiano un senSO.

-

Problema B Trovare il massimo autovalore >o, lo, e la corrispondente autosoluzione positiva l+,(x)>o, per il Problema A omogeneo, e cioB per fl=o e ulr=o; tal caso fisicamente corrisponde ad una ripartizione di neutroni autosostenentesi (reattore critic01

.

§

2

- Soluzione del problema A

Cercheremo u in IIo1 (n) (u 1 r=o): le condizioni di trasmissione appariranno come "condizioni naturali", automaticamente soddisfatte nella formulazione variazionale, come B ben noto nel caso d'una equazione sola (g-1). Per g>l il Pb non 6 autoaggiunto. Si k allora prcferito seguire un ragionamento elementare basato sul fatto che c 1(x)~o. Introdotte le forne lineari:

a(u,v)=

I n

(cpq 'D

Ei 1

Q

1

- - + Aiu iv i)dn

P9 . axp axq

siano A , B, C i corrispondenti operatori (matriciali) penerati (nel senso delle distribuzionil, Per le ipotesi fatte sui coefficienti essi risultano def initi su tutto ( H : ) ~ con codominio c ( ~ Gli ) elenenti d i natrice corrispondenti sono:

ove :

Introdotti A, 3 , C il Problema A si riconduce alla.risoluzione dell'equazione in u:

I1 procedimento esistenziale k allora il seguente: Se u=o allora la: Aou=(A-C)u=f 5 costituita (si osservi la struttura delle ( 2 1) da g equazioni disaccoppiate (cl~ 0 del ) tipo: (6)

Assumendo fe(f1-I)g

(o pin senplicemente

( L Z ) ~ )per

le ipotesi

1) e 2 ) ) esiste allora un operatore G: di Green(') che ci fornisce u1 = G: f 1 ( 6 : b un isomorfisrno di H-'(O) su H:(Q)). Trovato u l , u2 b fornita da: u2=G~(~'u1+f2) = G:(C~G:~~+~), e in i-2fi-1 nerale si ha: ui= C io ( C iGo

ge-

Ne consegue llesistenza (u=o) di un operatore matriciale di Green, Go, per il nostro problema, che realizza un isomorfismo di ( H " ) ~

su (i~:)~,

e u=Gof r ( H : ) ~ soddisfa "naturalnente" alle

condizioni di trasmissione. I1 caso p # o si tratta riducendo la (5) nella forma:

Essendo Go un operatore limitato da 2

catori in L (01, e llimmersione di

( L ~ ) ~ +(H:)~, Bij moltipli(I~:)~+(L~)~compatta, allo-

ra T b compatto da (II:)~+(H:)~, e pertanto il Problema A 2 ricondotto ad un classic0 Problema di Riesz-Fredholm. Lo spettro puntuale { p : } amnette llunico punto di accumulazione all'infini to, e i l i i L- si accumulano solo verso zero. LJ :

5 3

-

a)

In tal caso si ha:

Soluzione del Problema I3

risolvere la (10) si 6 trovata una rappresentazione a nucleo (funzione di Green) dell'operatore L-l. Se g.1 la funzione di 1 Green G (x,y) 6 stata trovata, nel caso dei coefficienti discontinui, e per la la volta, da ~tam~acchia(') come soluzione di: Pep

(11)

A'G'(X,~)

=

6

x IY

(6x,y misara di Dirac in (xiy))

Allora u =

1

~'(x,y) f(y)dy ? soluzione di A 1 u=f e si ha

n

-

Osservazione Questo fatto di Q'CQ non ci permette (bench? plausibile) di concLudere che G' t L' (QxQ) (e i suoi iterati). Allbra sfruttando il fatto che per l'operatore L si possono (vedi (6)) fare dei ragionamenti per singole equazioni (ez sendo queste disaccopiabili) si pu6 costruire subito una matrice "fornale" d i Green di elenenti ti;:

Si osservi che Habetler e Martino (8) gii nel 1958 avevano considerato il Problema B assumendo perb formalmente l'esisteq za della funzione di Green. (O)

Ad esempio:

se i Gi fossero t L'(QXR), Si verifica poi subito che ~ - ~ z 1{ t2 ~tale ~ ) che

-

L-~LUEU, e ~( H : ) ~ ,

e pertanto L - ~ I { ~ 5 ~la~matrice } di Green del

problema Lu = 4 .

-

Osservazione Per le ipotesi di positivits fatte sulle ci l'opee tore L-' lascia invariato il cono, in (L*)~,dei vettori >o, come pure, per le ipotesi sulle B ~ l'operatore , L-~B. (compatto in (L2)E). Ne consegue subito, per noti risultati di Krein-Rutman, (l), che esiste un autovalore massimo (dominante e semplice) X o del problena, . . = L - ~ B ~positivo , e maggiore del valore assoluto di ogni aitro autovalore, a1 quale corrisponde un'autosoluzione uo pure >o in Q.

roue

b)

-

Determinazione iterativa di A,.

Richiamiamo il Teorema di I. Marek, (6) Se II,K sono due operatori lineari con DH, D K C K , spazio di Hilbert, e se K B limitato e B 1 esiste limitato da X-DH, e H - ~ K 2 limitato con autovalore dominante X, (e autosoluzione x,) all2 ra il process0 iterative di Rayleigh-Kellog descritto dalle: u (0) = x(~)

converge:

(approssimazione zero!

Per avere una determinazione i t e r a t i v a d i X o basta a p p l i c a r e t a l e Teorema n e l 7ostr0 caso K=B, H=L, X ; ( L ~ ) ~ . S i o s s e r v i ora che il calcolo d e l l ' i t e r a z i o n e m + l s i deve r i s o l v e r e un'equazione d e l t i p o :

e l e (131, come a 1 s o l i t o disaccoppiandosi, permettono d i t r o v a r e l e u(m+1)9i una per v o l t a ( i = 1 , 2 , . .g) risolvendo problemi d e l

.

tipo:

A ciascuna i t e r a z i o n e s i devono p e r c i s r i s o l v . r e problemi d i t i p o ben noto. Qcesto consente p e r c i s d i applicare tecniche svai r i a t e s t u d i a t e per problemi d e l t i p o A u = f .

Ad esempio l a u ("' il funzionale:

) 9i

puh e s s e r e t r o v a t a

'L

,

(9)

, rnininizzando

essendo f i nota d a l l a precedente i t e r a z i o n e . Ne consegue che adattando questo metodo s i ha un c i c l o doppiamente iterative-variazionale per trovare uo, X o . Evidentemente l a u ) 9 i pus e s s e r e t r o v a t a anche c o l metodo d e l l e d i f f e r e n z e f i n i t e , e appunto s i sono f a t t e d e l l e esper i e n z e numeriche comparative a 1 riguardo.

5 4

- Esperienze numeriche

a ) Risoluzione d e l Problema ( 1 5 ) a t t r a v e r s o il metodo d i Riesz. Se s i assune una'lbase" f i n i t a F v ( x ) ; v = 1 , 2 , . . . sentano l e

(

i cone

M

=

i a F 1

,

M, e s e s i rappre-

l e condizioni d i ninino per

il funzionale (15) diventano:

i = 1 , 2 , ...g

- a= oI

;

aa.

v = 1,2,.,.M

1,v

Queste c i danno un sistema algebrico l i n e a r e d e l t i p o ( p e r ogni i ) : R ia (m+l),i-B(ra),i -

...* aN(m+ 1 1 ,i

;

(i=1,2,...~)

essendo R,S,Q m a t r i c i NXN d i elementi n o t i . Ad esempio:

L'inversione d e l l e m a t r i c i R 2 s t a t a f a t t a con l'algoritmo

di

Gauss. b) Confronti t r a il metodo iterative variazionale d e s c r i t t o e quello d e l l e d i f f e r e n z e f i n i t e (assunto cone elemento d i confronto). Caso monodimensionale: t r a t t i ; 1.1=20.

g=2, r = 5 , D,A,3,C funzioni c o s t a n t i

X o (con d i f f e r e n z e f i n i t e )

F 1

v

0

=

Xdf

sono l e autosoluzioni d e l l ' o p e r a t o r e d i Laplace nonodin?ensionale (con il n o s t r o metodo) = X N

a

Ecco una t a b e l l a r i f e r e n t e s i a i v a r i c a s i : L'errore % 2 a 1 pih a t t o r n o a110 0,1%. Circa l'andamento d e l l e s o l u z i o n i u l , u2 n e i c a s i s p e r i rnentati s i ha un accord0 d e l nos t r o metodo con q u e l l o a l l e d i f ferenze f i n i t e s i n o a '3 c i f r e s i gnif i c a t i v e n e l l e zone c e n t r a l i , e uno meno buono (2 c i f r e ) n e l l e a l t r e zone.

-

Caso bidimensionale R i s u l t a t i analoghi a i precedenti, perch2 l ' e y r o r e 5 (!I eguale s i a FET l ' a s s e xl che per l ' a s s e x 2 ) n e i n o s t r i e s p e r i n e n t i non ha n a i superato l o 0,38. (Fv sono p r o d o t t i d i autosoluzione d e l l f o p e r a t o r e precedente). Caso tridimensionale: E - 2 , r=3; D , A , B , C

c o s t a n t i a p e z z i , simmetria

r i s p e t t o a i p i a n i x , = o , x2=0. (Fv p r o d o t t i d i autosoluzione come prima). X):

L

Risultati a ) NZ = M

Y

=

I1

1,

X

= 3 e cio2 N t o t a l e abbastanza piccolo

()I=FJ t1 N 1: con 6 i t e r a z i o n i (30'' IB?? 7090) A H approssima X X Y Z

e n t r o il 3,3%. b ) aurnentando II

na),

Y

df

da 3 a 10 l ' e b s i r i d u c e a 1 3 4 (1'47" d i rnacchi

,

In generale X = XI,[ 6 d i t i p o monotono (crescente i n N) e i n 15 it: r a z i o n i a 1 p i h , n e i c a s i c o n s i d e r a t i , s i ha l ' a u t o v a l o r e can l a approssimazione c e r c a t a (1%) mentre 6 ben noto che con il netodo dell e d i f f e r e n z e f i n i t e il nunero d e l l e i t e r a z i o n i s a l e , i n genere, a l meno a c i r c a 50460. u l ,u2 ( s u l l e r e t t e y=6, Z=8 i n f i g u r a l sono i n buon accord0 con i valori ottenuti a differenze f i n i t e , tranne n e l l e i n t e r f a c c e ove l o scart o 4 ?. 2 , 3 % ,

-

Osservazione I1 netodo i t e r a t i v o variazionale f o r n i s c e l'autovalore massimo i n modo a s s a i soddisfacente, s i a p e r precisione che pep ternp i d i c a l c o l a t o r e ( r i s p e t t o a 1 metodo d e l l e d i f f e r e n z e f i n i t e ) , ma invece 4 i n f e r i o r e a quest'ultimo metodo p e r l a precisione d e l l a t i bulazione d e l l a soluzione s p e c i e per quel che riguarda l'andamento d e l l a u2 che pu6 presentare d e i "picchi" n e l l e zone non c e n t r a l i ( e vicino a l l e i n t e r f a c c e ) , picco a v o l t e a s s a i ma1 d e s c r i v i b i l e c o l metodo variazionale.

PARTE 2a

- Teoria a multigruppi dipendente dal tenpo

-

1 E' ben noto che nella teoria della diffusione dei neutroni nell'approssimazione a pih gruppi g di velocith, che supporremo due per sempliciti, l'evoluzione del tempo dei flussi veloce e lento, rappresentati da u l , u2 e della concentrazione dei cosi detti "neutroni ritardati" rappresentata da C, 8 retta dal seeue; te sistema: 5

essendo assegnate: 1) i coefficienti (funzioni nisurabili e limitate essenzialmente

>O) e le funzioni di "sorgente"

f1,2,3

'

2) le condizioni (Dirichlet) per le ul, uz a1 contorno O o una matrice diagonale NXN a c o e f f i c i e n t i >o una matrice diagonale NXN a c o e f f i c i e n t i >o una matrice tIIXN a c o e f f i c i e n t i >o

+ + +

H

+

a22 h

El

+

X

A22+

,

,

una matrice NX?I1 a c o e f f i c i e n t i

>o

, ,

,

,

.

-

Osservazione I Lo s t u d i o d e l problerna d i s c r e t o ( 2 ) $ egualmente n o l t o importante i n a n a l i s i numerica, perch$, a d i f f e r e n z a d i quanto accade n e i c a s i p a r a b o l i c i , c ~ n c e r n e n t ii n generale p i c o meno l a d i f f u s i o n e d e l c a l o r e , l a matrice pub avere a u t o v a l o r i > O , il che cornporta a v o l t e una crescenza n o l t o r a p i d a d i $ cosa sempre d e l i c a t a da con-:rollare d a l punto d i v i s t a numerico, Inoltre l e v ( i n generale c o s t a n t i ) sono 5 l o 6 , mentre 192 g l i a l t r i c o e f f i c i e n t i sono 2 1 , e p e r t a n t o d a l punto d i v i s t a numerico s i incontrano d i f f i c o l t a s i m i l i a q u e l l e che s i hanno n e i prohlemi d i "boundary layer" connessi con equazioni d i f f e r e n z i a l i contenenti p i c c o l i parametri n e l l e d e r i v a t e p i h a l t e , Usando schemi i m p l i c i t i ( p e r r a g i o n i d i s t a b i l i t s ) s i generano da ( 2 ) "grossi" s i s t e m i l i n e a r i per i q u a l i occorrono netodi it: r a t i v i l a cui convergenza, che e r a da indagare, B s t a t a v e r i f i c z t a i n , (3)

- Proprieth del sistema (2).

5 2

In ( 3 ) sono s t a t i o t t e n u t i i seguenti r i s u l t a t i :

5 irriducibile; 6 essenzialmente >o; 3) Q possiede un autovalore wo>-A cui corrisponde un autovettore v>o, t a l e che se ai 5 un qualsiasi a l t r o autovalore 8: R a.w 0'

5 3

- Metodi d i risoluzione d i ( 2 1 ,

a ) Metodo Esplicito: + ( t )= (I+At Q ) ((t-At). Tale metodo 8 s t a t o s c a r t a t o nei n o s t r i c a s i perch8 ha una soglia d i s t a b i l i t a t r o ~ po bassa (At troppo piccolo). b) Metodo Implicito: 4 ( t ) = ( I - ~ t~ ) " + ( t - A t ) ,

Ad ogni passo temperale c ' 6 da risolvere un sistema del tipo:

Se A t w o c l a l l o r a , i n base a l l a p r o p r i e d 5 del 5 2 , i metodi d i Jacobi e Gauss-Seidel r e l a t i v i sono convergenti.

-

Osservazione Quando s i ha un "transiente" molto rapido s i 6 trovat0 che anche il metodo implicit0 (e pure q u e l l i d i Crank-Flicolson e s i Saulyev ( 4 ) r i s u l t a n o molto imprecisi, Pertanto s i pone il problema d i trovare qualche metodo meno imprecise. La valutazione ( 4 ) ha f o r n i t o l ' i d e a base per il seguente metodo che chiamereno metodo U.

Nel netodo w s t ; pensato d i esprimere l a soluzione n e l l a forma:

Allora l a ( 2 ) s i trasforrna i n :

In ogni caso per6 c t & il problema d i determinare w o che non B determinabile con procedimento t i p o "metodo d e l l e potenze Rayleigh-Kellog" non essendo l t a u t o v a l o r e q u e l l o d i modulo massimo. Questa questione 6 abbastanza d i f f i c i l e d a l punto d i v i s t a numer i c o , perch; & vero che s i pu6 t e n t a r e d i " t r a n s l a r e " l o s p e t t r o a1 f i n e d i condursi ad un problema d i autovalore d i massimo mod2 l o , ma c o s i facendo, dovendo poi s o t t r a r r e il passo d i t r a s l a z i c ne, s e questo k molto grande s i pu6 perdere ogni s i g n i f i c a t o . Allora posto B=M-N :

c i s i r i d u c e , come equazione a g l i a u t o v a l o r i per w (essendo w o 0 autovalore d i -VB) a l l a ( 5 ) : (M+uo v-' x = Mx. I n t r o d o t t o un parametro f i t t i z i o v s i dimostra che l t e q u a z i o n e 1 (M+ w ~ - ' ) x Nx possiede un autovalore d i massimo modulo c u i v corrisponde un a u t o v e t t o r e > o o t t e n i b i l e c o l metodo iterative d e l l e potenze. Q u e l l o che s i dimostra k che p = p ( w ) & monotona decrescente, e che v = l individua w,. I1 sistema ( 4 ) & poi r i s o l t o con il metodo i m p l i c i t o CrankNicolson ed i r e l a t i v i metodi i t e r a t i v i r i s u l t a n o convergenti. ( 3 )

-

§

4

- Metodo dei passi frazionari

Recentemente, (5) , abbiamo pure speri~entatoil metodo di Narchuck ed altri per la (1) (g= 1) assurr.endo C =o e come deconposizione dellfoperatore A una del tipo (vedi § 5 ) :

a a) -

Al. ax (D ax + ib a A2= a (D -1 + ;b aY aY

(7 1

Lo schema alternato per il passaggio da tn+ tntlltr[tn,tntl[ 2 il seguente: n+i

-v - tdt AIU

- 0

fornente u(tn+,) = untl (tntl). La discretizzazione della ( 8 ) ci da poi sistemi del tipo: "-1 (9)

& = dt

v-' & dt

-

a1U

;

V-'

a2U

;

al,a2 matrici tridiagonali

matrice diagonale;

Osservazione Un primo vantaggio del metodo k che abbiamo ora a che fare con matrici tridiagonali (invece di pentadiagonali) invertibili anche con metodi diretti.

L1w0 c a l c o l a t o come i n d i c a t o a l l a f i n e d e l § 3 r i s u l t a oo = 63,69, mentre il valore e s a t t o B 63,21, (Le d i f f e r e n z e f i n i t e sovrastimano l ' a u t o v a l o r e e , a n o s t r a conoscenza, c i sono r i s u l t a t i t e o r i c i per questa stima s o l o n e l caso d i c o e f f i c i e n t i c o n t i n u i i n T I , I r i s u l t a t i numerici sono r i p o r t a t i i n t a b e l l a dove

;; &

il valore nedio su 0 d e l l a soluzione e s a t t a , h l a soluzione approssimat a o t t e n u t a con l a trasformazione w ed il metodo d i Marchuck, 3, l'analoga soluzione senza trasformazione w , u quello o t t e n u t o I Mu con il metodo i m p l i c i t o che f a seguito a l l a trasformazione w , me5 tre la 6 q u e l l a che s i r i f e r i s c e a 1 metodo i m p l i c i t o d i r e t t o . D

cIM

Come s i vede i n ogni caso l a trasforrnazione w , conunque s i a assoc i a t a ad a l t r e tecniche, d2 i r i s u l t a t i m i g l i o r i , ed il metodo d e i p a s s i f r a z i o n a r i s i B r i v e l a t o superiore, n e i c a s i f a t t i , a 1 metodo i m p l i c i t o .

t

S o l u z i o n e esatta

-

-u

Mu

-

-

U~

U

~

~

u

U~~

Bibliograf i a (1)

S.Albertoni: "Metodi v a r i a z i o n a l i per c e r t i s i s t e m i d i equazioni a derivate parziali"

- 1st.

Lombardo Scienze e L e t t e r e

Vol. 100-1966.

(2)

-

-

S .Albertoni M.Lunel1i G.Mangioni: "Metodi i t e r a t i v i var i a z i o n a l i per problemi e l l i t t i c i n e l l a t e o r i a d e i r e a t t o r i nucleari" A t t i Seminario Mat. e Fis. Univ,Modena Vol.XIV ,1965.

(3)

A.Daneri

- A.Daneri - 1 , G a l l i ~ a n i : "A

numerical approach t o

t h e time dependent neutron d i f f u s i o n equations" EUR 3742e

-

1968.

(4)

1,Galligani: "Numerical s o l u t i o n s of t h e time dependent d i f fusion equations using t h e a l t e r n a t i v e method of Saul'yev" Calcolo Vo1.2, Suppl. 1, 111

(5)

S.Albertoni

-

- 1965,

- A.Daneri - G.Geymonat:

"Existence and approxi-

mation theory f o r general d i f f e r e n t i a l equations of t h e multigroup d i f f u s i o n r e a c t o r theory" ( i n corso d i pubblicazione), (6)

1,Marek: " I t e r a t i o n s of l i n e a r unbounded o p e r a t o r i n non s e l f - a d j o i n t eigenvalores problems and Kellog i t e r a t i o n processes" (Cech.Meth.Journa1 1 2 (1962)).

(7)

G.Birkoff

-

R,S.Varna:

"Reactor C r i t i c a l i t y and non negative

matrices" J.Soc.Indust.App1.Math. (8)

G. I. Habetler

- M.A.

6 , 354-377 (19581.

Martino: "The multigroup d i f f u s i o n

equations of r e a c t o r physics, KAPL, 1886, J u l y 28, 1958 (9)

G.Stampacchia: "Su un problema r e l a t i v o a l l e equazioni d i tipo e l l i t t i c o d e l 2 O ordine" Ricerche d i Mat., Vol. V (19561,

CENTRO INTERNAZIONALE MATEMATICO ESTIVO (C. I. M. E. )

"PROBLEMS OF OPTIMIZATION AND NUMERICAL STABILITY' IN COMPUTATIONSfI

Corso tenuto ad

Ispra dal 3-11 Luglio 1967

PROBLEMS OF OPTIMIZATION AND NUMERICAL STABILITY IN COMPUTATIONS

)

by I. ~ a b u g k a(Praga)

Computer Science is a new scientific discipline. An of this discipline

is the numerical

important part

mathematics. The "Art of Computationn

i s becoming science ; new questions and problems become important.

A typical problem i s the problem of the creation of numerical methods, the determination

of

their

"worthn and, in general, the choice of

the most suitable method for the given purpose. For example, the program-library in a computing centre centains mostly many algorithms for solving single mathematical problerfis. Opinions on the expedience of these algorithms a r e usually quite different and subjective. This statement i s still more apparent when a method of applied mathematics

is to

be

appreciated, especially in the field of scientific-

technical computations. These scientific part

of the computer science My paper will

in

I think that

these

to

the

computations

which

are

more o r

computations may be characterized a s way

required it

I have some experience.

deal with questions

this kind of computations.

t.on

technical computations a r e that

which

ciated with

and constructive

-

l e s s asso-

a mathematical

of processing (transformation) of the given informaone )

,

I

am

sure

that

in

scientific-techical

i s necessary to emphasize the knowledge of information

which we may collect and the appreciation

of its reliability. Further it is

In this paper some results obtained recently in Prague will be given defines numerical analysis a s the theory of constructive 2 ) ~ e n r i c i23 methods in mathematical analysis (with emphasis on the word llconstructiven) .

necessary to formulate clearly the required information on the given problem

. The

necessity of a mathematical and constructive way cf this pro-

cessing i s obvious here The "clarity part for a

fl

.

of the given and required information

successful solution

is

an important

of a technical problem. Numerical mathe-

matics a r e the rudiments of this constructive processing of information. Numerical method

generates (in a constructive manner) a mapping, from

the class (space) of the given information to the class of the required one. It i s important that this mapping i s defined on the entire class - of be the domain

information. This class will thod (mapping)

.

Numerical process creation

of definition of the given me-

of the given

is

an

exact constructive law (prescription) of

mapping.

.Computation is

a

given case. We

shall talk about

without

round-off

tion) in

a

concrete realisation

of the numerical process in the

- realisation

exact

e r r o r s and about a

real

when we compute

realization ( o r disturbed realisa-

computation.

It i s obvious that there a r e many different manners tive creation of

cess is It

a construc-

one given mapping, i. e. many processes exist which

transform the given information to the requested mathematical

of

problem. It

is

evident

one and solve the s a m e

that the question of choosing a pro-

very important. i s c l e a r that the choice and every optimization must necessarily be

relative to the given information. This does not me methods might

not be

advantageous

in

mean, however, that so-

a certain generality.

The manner in which we appreciate the method i s of great importance. My experience is that, from the practical point of view, important to respect an incredulity dulity can be of different kinds

of

the given information

. Some of

them

it

is

very

. This incre-

will be shown

in the next

v

I. Babuska

part of the paper clusion

-

. It i s

essential that the method -and in

general

be stable with respect to these incredulities. I think

stability i s one of

the most important points

when

all

con-

that this

choosing a method in

practice. In the next part I

shall point

out some aspects of these questions.

2 . T h e p r o b l e m of q u a d r a t u r e f o r m u l a s 1 )

In this section I shall show

some aspects of ideas, which I

men-

tioned previously, in a simple case of quadrature formulas. Let our task

be to

determine

numerically

Jo We shall function

suppose f(x)

1.

we

know

the following about 'the integrated

:

The function

the period 2.

that

2v We

can

f(x) i s

a continuous

periodic

function with

. evaluate only the

function

f(x) (i. e. compute the va-

lues of f(x) )

.

In this

the simpliest quadrature formula T (f) i s mostly used in n

case,

practice,

with

This formula i s the well known trapezoid formula. I will now analyse the question, if there a r e any reasons for selectins

') In this errors.

part we a r e

not dealing with the problems of the round-off

the trapezoid formula ; we may ask e. g. why the Simpson-formula isn't better than the formula previously mentioned. Some arguments for choosing the trapezoid formula (in this case of integration of a periodic function) a r e included in some papers,

e. g. Milne

(25) , Davis (18) and others.

The e r r o r bounds for the trapezoid formula a r e studied in many papers. See (4) (5), (21), (24)/ and others. We will now analyse the proF

blem of the choice of the quadrature formula according to the information we mentioned previously. In our considerations we shall confine the class of possible formulae to the linear one, The choice of the quadrature formula means, in our case, to determine of the sequence of linear functionals I in the form n

with the requirement

that

Jn(f)

+

J(f) (weak) for all functions f(x)

of the given class of functions. We shall measure the amount of work in using a formula by the number of evaluations of the integrated function. Let us

now assume that

B is a Banach space. Then we can de-

fine (2.4)

w (n, B) = inf

sup

(2

a j f (yj) - J(f)

1

and

w(n, B) is the minimal possible e r r o r under the assumption that we know only that I\ f

\\ B 4

1.

3

(11,B)

has analogous

meaning when we confine our-

selves to use equidistant points in the quadrature formula.

We shall further introduce

(2.6)

/\(n,B)

" ~ ~ 6 4

A (n, B) the space

sup

=

the error-bound

i s evidently

An objective measure of convenience of the given formula i s

B.

given here by the comparison of appreciation

of the trapezoid formula in

is

A (n, B) with

w(p, 8).s ( n , 8 ) reSP.

This

obviously relative to the space B.

The choice of the space

B

is

very problematical in practice

.

In majority of cases there is a large incredulity a s to whetherit is convenient

- B . If the

to take the integrated function a s an element of a certain space

conclusion on the suitability of a formula is strongly dependent on the choice of B , then to use that will

the conclusion

is

not "stablett and it is not advantegeous

formula in practive. Further we shall see that this ttunstabilitytt

appear in the Fase

formula whose

error

of the optimal formula, i.e. equals

9

(n, B) o r

will strongly depend on the space

B

when we use the

w (n, B) then

. Conversely

, a formula will be

3 (n,B)

advantageous in practice if its e r r o r is nearly equal to w(n, B)

but more o r l e s s

Later we shall an optimal

one,

see

has

that

this

only the trapezoid formula which

property. We now

of periodic

Definition 2.1

The Hilbert space

be

periodic

1) Every f E H

2) Let

3) If

11

f

is

introduce a class of

H (over complex numbers) will be said

if: a

2w

periodic,

continuous function.

signify the norm in the space

H will

i s not

functions.

f h H , then g(x) = f(x+c)e H

The space

or

independent of the space B.

Banach spaces

to

the results

be said to

C, then

for every real c and

be strongly periodic

if

!I f I!

=lip H H .

it

is periodic

and if :

1) e

ikx

eH, k =

5) If

jl

(

for

..., -

1, 0, 1,

& ( k l , then

...

~

(

ikx

(IH=

~

~H ~

and IIe e

~

-ikx Oe ~

\IH

.

~

e

l

~

~

~

O606d2

D

and

does

At

not

depend

t h e beginning

dic function. It

I think it

of

on

n.

t h i s section it w a s s a i d that f(x) i s a

i s obvious that t h i s information

It is evident

that

now

too

we

H

perio-

is insufficient. However,

is convenient t o a s s u m e that the function

of a periodic o r strongly periodic s p a c e

f(x) i s a n element

.

have

a l a r g e incredulity a s

r e g a r d s the concrete selection of the s p a c e H. The importance of this incredulity is well Theorem 2.1

seen

Let

i n t h e next t h e o r e m and example.

H

be a

strongly

periodic s p a c e with

the n o r m

20

\I~I\'-

(2.9)

JY-l2 +

A \fl\

2

dx,

A > 0.

Then the error-bound of t h e formula

where

i s equal

to Q

2

The t h e o r e m if

l

we

wing

(n, H) 2..1

.

affirms

that

the f o r m u l a

a r e using t h e equidistant

(2.10) is

net. Now we s h a l l

example

Example 2 . 1

Let

f(x) =

e

Wsin x

, = 3,lO.

an

optimal one

introduce

the follo-

l

Then

(f)=

1

2r f(x) dx = 4,88079258586502208..

1

.

. 2815,71662846625447.

resp

In Tab. 2.1 we show t h e r e s u l t obtained by the trapezoid formula R(A' f o r A = 1 . F r o m this table we see that a n optimal formula used in a n n inconvenient s p a c e may give bad r e s u l t s . We s e e that the conclusion of the convenience of the optimal formula i s v e r y llunstablell with r e s p e c t to the choic e of H)

. From

this table we a l s o s e e that t h e trapezoid formula ( C = 1 ) gives

v e r y good r e s u l t s ; however, the following t h e o r e m is t r u e : Theorem 2.2

F o r e v e r y periodic s p a c e

H

A (n. H) > g (n, H)

(2.11) This theorem periodic

shows

that the trapezoid formula cannot

space. Nevertheless t h i s

formula i s v e r y advantageous in

tice. The explanation of t h i s fact c a n be s e e n T h e o r e m 2.3

. Let

H

be a

b e optimal in a

periodic

in the following

for all

periodic s p a c e s

(except

statement:

space. Then

No o t h e r formula has the p r o p e r t y that the left-hand s i d e of ded

prac-

f o r a finite number

(2. 12) i s bounof

indices

of

n) This

t h e o r e m shows that the e f f i c ~ e n c yof the trapezoid

formula i s

I

Nunber of points n

O
0

.

w is indepen-

the optimal overrelaxation parameter.

introduce an example.

Example 4.2.

Let

us solve the one dimensional proble y n = 1, y(0) =y(l)=O

with the finite-difference

method and overrelaxation. Because of the

round-off e r r o r the iterations do not, in general

converge to the requi-

uquasi-converge" in a more o r l e s s well known 1. . sense. In Fig. 4.6 we see n = - ~n dependence on h. h We s e e a good agreement with theorem 4.3. It i s possible to

red solution. They will

formulate the theorem See

(12)

.

4.3 for

a

0

< w -< 1 in a more general form.

Further processes have also been

investigated. I shall mention

here the stability of the Kellog process for the determination of eigenvalue (see (27) ) and a the theory of reactors form mapping

numerical process for solving a problem of (44) and the process for computation of con-

. (see also

I have shown a

(12) )

.

few different aspects of incredulity with regard to

the given information which appear in computations. I think that this kind of investigations is very important when choosing an algorithm in general.

-

59

-

Fig. 4. 1. a.

I. b decadic Fig. 4.1. b

10

(0.5, 1 ) I . h decadic

n

-

Fig. 4 . 2 .

( 0.5

,

II. b

1

>

Fig. L. 3.

Fig. 4 . 4 .

(300,

300.5)

II. b

I. ~ a b u i k a References ----------

I11

A.A. A ~ ~ ~ M o6 BI I:e p e H O C e

PPaHMllAbIX Y C J I O B M ~AJlR

CMCTeM J ~ H H ~ ~ H H x

O ~ H K I ~ O B ~ H H ~UM @ @ e ? e H ~ Y I t v r b H bYl ~~ B B H ~ H M BblrI.MaT. ~ ~ 1$14a.

I4 MaT.

1 ~ 0 1 ,1; 542-545.

L2]

A.A.. A6pa~oe: B a p a a ~Mf2TOAa ~

131

A.A. Abramov : Transfer of boundary conditions for system of ordinary linear differential equations. Proc. of IFIP Congress 65, p. 420

[5]

I. Babuska : Uber die optirnale Berechnung d e r Fourierischen Koeffizienten. Apl. Mat. 11, 1966, 113-122.

[61

1. Babuska : U'ber universelloptimale Quadraturformeln. Apl. Mat. 1968

17J

I. Babuska , M. P r a g e r : Numerisch stabile Methoden zur ~ a s u n g von Randwertaufgaben. ZAMM 1961 , H. 4-6

181

i1.&36yd~a, C .J.C06one~:

IIpOrOHKEI.

, i t . ~ h l r l . ~ B T . YI MBT.

9143.

1961, 1, 349-351

Apl.mat.

~ ~ T H ; ~ I M ~ ~9MCJIeHHblX u M R MeTOAOB.

10, 1965, 96-129

[91

I. Babuska, M. ~ r a ' ~ e E. r , .VitBsek; Numerische ~ t a b i l i t 2 t von Rechenprozessen. Wiss. Z. Techn. Hochsch. Dresden 1963, 12, 101-110.

1101

I. Babuska, M. Prager, E. VitBsek; Numerick6 r*es'eni differenciilnich rovnic 1964 , SNTL

0 11 I. Babuska, M.PrBger, E. VitBsek: Stability of Numerical Processes, Proc. of IFIP 65 , 602-603

.

1121

I. ~ a b u g k a ,M. PrBger, E. VitBsek: Numerical Process in Differential Equations. Interscience Publishers 1966

(131

F. L. Bauer: Numerische ~ b s c h z t z u nund ~ Berechnung von Eigenwerten nichtsymmetrischen Matrizen. Apl. Mat. 10, 1965, 178-189.

1141

F. L. Bauer et al. Moderne Rechenanlagen , Stuttgart 1965, p. 64.

[151

F.L. Bauer: Genauigkeitsfragen bei der Losung linearer Gleichungssy-

b63

F.L. Bauer, H. Rutishauser, E. Stiefel : New Aspect in Numerical Quadrature. Proc. of Symp. in Appl. Mat. 1963 , XV, 199-218.

steme. ZAMM 46, 1966, 409-421.

I. Babuska

1171 R. Bauman : Algol Manual d e r Alcor-Gruppe, Sonderdruck a u s Elektronischen Rechenanlagen H 516 (1961) H 2 (1962) R. 01denburg ,Munchen. 1181 P . J.Davis: On the Numerical Integration of Periodic Analytic Functions. Proceedings of Symposium Madison 1959. [191 G. G. Dahlquist : On Rigorous E r r o r Bounds in the Numerical Solution of Ordinary Differential Equations. Numerical Solution of Nonlinear Differential Equations. Wiley 1966, 89-96 . b0] G. G. Dahlquist: P r i v a t e communication. [21] H. Ehlich : Untersuchungen z u r numerischen F o u r i e r analyse. Math. Zeitschr. 91 (1966) , 380-420

.

1221 G. Hammerlin: Uber ableitungsfreie Schranken f:r Quadraturfehler. Numerische Mathematik 5, 1963 , 226-233; 7, 1965 , 232-237. 1231 P. Henrici: Elements of numerical Analysis. J . Wiley New York-London-Sydney, 1964.

Sons, Inc

1241 D. J a g e r m a n : Investigation of Modified Mid-Point Quadrature F o r m u ? l a , Math. of Comp. 20 1966 , 78-89. I,

1251 G. Kowallewski : Interpolation und genahreQuadratur. Leipzig 1 9 3 0 , p. 130

1271 I. Marek: Numerische ~ t a b i l i t z td e r P r o z e s s e vom Keloggschen Typus. Liblice 1967 Apl. Mat. 13, 1968. [28] J. Milota ; Universal Almost Optimal F o r m u l a e Solutions of Boundary Value P r o b l e m s f o r Ordinary Differential Equations. Liblice 1967 , Apl. Mat. 13, 1968.

-

1291 R. E. Moore: The automatic Analysis and Control of E r r o r in Digital Com putation. Vol. 1, 61-130 . Proceedings of a s e m i n a r University of Wisconsin, Madison Octobre 5-7, 1964. 1301 R. E. Moore; Interval Analysis. P r e n t i c e Hall 1966 f31]

.

R. E. Moore: P r a c t i c a l Aspect of Interval Computation. Liblice 1967. Apl. Mat. 13, 1968 .

[321 M. P r a g e r : Numerical Stability of t h e Method of Overrelaxation. Liblice 1967. Apl. Mat. 13, 1968. \331

P . P r i k r y l : On Computation of F o u r i e r Coefficients in Strongly Periodic Spaces. Liblice 1967, Apl. Mat. 13, 1968.

I. Babuska

A. Sard : L i n e a r Approximation. Providence 1963.

K. Segeth : On Universally Optimal Quadrature F o r m u l a e Involving Values of Derivatives of Integrand. Liblice 1967 , Apl. Mat. 13, 1968.

H. J. Stetter : Numerical Approximation of F o u r i e r - T r a n s f o r m . Num. Math. 8 , 1966 , 235 - 249. J. Taufer: On Factorization Method. Apl. Mat. 11, 1966,427-452

.

I

J. Taufer : Faktorisierungsmethode f u r ein Randwertproblem e i n e s l i -

nearen Systems von Differentialgleichungen, Liblice 1967 Mat. 13, 1968.

. Apl.

I

J. Taufer: Faktorisierungsmethode fur ein Eigenwertproblem eines l i n e a r e n Systems von Differentialgleichungep , Liblice 1967 E. VitBsek: Numerical. Stability in Solution of ordinary Differential Equations of Higher O r d e r , Liblice 1967, Apl. Mat. 13, 1968.

J. H. Wilkinson : Rounding e r r o r s in algebraic H. M. S. 0. 1963

p r o c e s s e s . London

J. H. Wilkinson : A. Survey of E r r o r s Analysis of Matrix Algorithms. Liblice 1967, Apl. Mat. 13, 1968. 1

R. Zezula : Numerische ~ t a b i l i t L t eines Algorithmus z u r Berechnung d e s E i g e n p a r a m e t e r s e i n e s Matrizenoperator mit Hilfe d e r Reduktionsmethode und d e r Banachschen Iterationene. Liblice 1967, Apl. Mat. 13, 1968. [45]

R. Zurmuhl : Matrizen und i h r e technische Anwendungen. Berlin 1964.

CENTRO INTERNAZIONALE MATEMATICO ESTIVO (C. I. M. E. )

J. H. BRAMBLE

ERROR ESTIMATES IN ELLIPTIC BOUNDARY VALUE PROBLEMS

C o r s o t e n u t o a d I s p r a d a l 3 - 11 L u g l i o

1967

INTRODUCTION Error Estimates in Elliptic Boundary Value Problems J. H. Bramble (University of Maryland) In these lectures, I will discuss some methods of obtaining error estimates for finite difference approximations to solutions of elliptic differential boundary value problems. Because of the limited time, I shall restrict my attention to the Dirichlet problem, although some of the methods are easily carried over to other boundary conditions. The first part will be devoted to second order problems. In fact, in order to illustrate the methods, I will restrict my attention to the Dirichlet problem for ~oisson'sequation, with "zero boundary values" and the classical Dirichlet problem for Laplace's equation. The last part will be devoted to some results qn higher order elliptic equations. In most cases, I will not give more than a sketch of the proof, indicating the details. Instead of choosing to discuss a general class of operators and a corresponding general class of difference schemes, I shall tonsider a specific operator and certain specific difference schemes so as not to obscure the essential points of the method of analysis. I will not, during these talks, make extensive references to related work but shall include in the bibliography a number of closely related papers. I shall restrict the references to the specific results under discussion.

I choose to formulate the first problem in a weak form. Let R be a bounded open set in E with boundary n

aR and let A be the Laplace

operator

We need the following class of functions V defined on R :

where C~ is the class of infinitely differentiable functions with support 0

contained in R and 'H

0

is the Hilbert space obtained by completing cm 0

with respect to the norm

All functions for the present will be assumed to be real valued. We now state problem I. Problem I:

Given F

I

E

(L,)'

, find

u

E

uA$= ,\d$ E V

, 1

Lp

6 p < ~ _ 2 such that

.

R

Here (L,)'

is the space of continuous linear functionals defined on L,

The quantity

is the value of the functional F at the point 4

.

For example, if "F c L1" then = completion of the

cm

.

F

The space L is the P functions in R with respect to the norm

We will discuss this problem later but let me first remark that when F and aR are sufficiently smooth Problem I is just the classical problem

The second problem is the classical Dirichlet problem Problem 11:

Given f

E

c'(~R)

, find

u

E

c'(R)

such that

Au = 0 in R

We will refer to well known results (eg. on regularity) for this problem when needed. We can immediately state the following theorems. Theorem 1: There exists a unique solution to problem I. Theorem 2: If

aR is such that at each point of

aR there exists

a barrier then problem I1 has a unique solution. Theorem 2 is classical. As we will see Theorem 1 (existence part) can be proved by means of a difference method. Let us now formulate corresponding difference problems and investigate their properties. Let ENh be the set of mesh points in E

N

form (ilh,

', iNh) , for h

> 0 and

il,

, i.e., points of the

. ., iN

integers.

.

We make some definitions. a) Kh = R r) ENh

( ei

is the vector with 1 in the i s position and 0 in the others ) N

We can immeidately state the following. Lemmal: in E Nh

Suppose

.

- A V 3 0 in R h , V 2 0 in EN,, h

- R,, .

Then V 5 0

From this it immediately follows that the discrete problem

Ah V(x) = F(x)

, x E R,,

has always one and only one solution for any given F and' g can introduce the discrete Green's function Gh(x,y)

.

defined as

Thus we

We can now make some statements about G Lemma2: Lemma 3:

.

.

Gh(x,y)30

F or any V defined on

\

This follows from uniqueness in the discrete problem, We next introduce the following function. Let

L Then we can prove Lemma 4:

For suitably chosen yN

Gh(x,y) .' V(x-y)

,

Y

E

Rh

.

, a,

and do

(independent of h )

The proof of lemma 4 consists i n showing t h a t it i s possible t o choose yN

,a

and d

0

i n such a way that f o r y

A

h,x

A

n>x

E

-N ~(x-y)=h

,

x = y

V(x-y)aO

,

X # Y

V(x-y)

all x

0

Then i t follows from lemma 1 'applied t o V(x-y) true

- Gh(x,y)

t h a t lemma 4 is

From lemma 4 we can prove

Lemma 5: of

.

Let

h and x

N

1$ p < ~ _ 2

.

Then there i s a constant

C

P

independent

such that

After using lemma 4 the sum i s estimated by comparing the sum with corresponding analogou~i n t e g r a l s using the f a c t that X

# Yo

.

, for

For the case N 2 3

x

-

y

is subharmonic,

6 1 N-2,

example, we can obtain the estimate

where S i s a s u f f i c i e n t l y l a r g e sphere containing R and having center a t an a r b i t r a r y point y

o

E

R

.

The integral i s convergent i f

p


0

sup luR

where K(E)

6

sh1

L K(E)

is a constant which depends on

E,

The proof of this is based on the following two lemmas. Lemma 7:

Let aR

.

Let d(x)

in 'ii

E

cL

and u be harmonic in R and h

be the distance from the point x

E

R to the boundary

(which is well defined in a strip S6 of fixed width 6 ) every

E

> 0

there is a K(E)

for x

E

Rhfl

Ss

.

such that

- B6lder continuous

.

Then for

This lemma follows from the mean value theorem for harmonic functions. One obtains an estimate for the HBlder continuity of the second derivatives which depends on the distance to the boundary. The next lema is crucial. Lemma 8:

For every

E > 0

there exists K(E)

such that, if aR E 'C

The proof of this.lemma is tedeous, long and involved, but the motivation is the following. We consider formally

where G is the continuous Green's function and we have assumed that d has been suitably extended to R

.

Then, formally,

It is not difficult to see using the maximum principle that for aR E C2

$(XI

a c dE(x)

Hence the procedure for proving lemma 8 is based on constructing a suitable comparison function and using the maximum principle, lemma 1. We shall now turn our attention, for a time to problem I. First I will briefly sketch a proof of Theorem 1. For uniqueness it suffices to

IVA~= 0 , V ~ E V

show that if cm

$ E CQ

then there exists a 4

E

then v = 0

.

But one can show that if

V such that A4 = $

.

Uniqueness follows.

The existence follows once the inequality

is established.

(One used the Hahn-Banach theorem.)

Since we shall be interested in the convergence properties of related difference schemes, we shall sketch the proof of (8) by the difference method. For any integrable function f we define

where C (x) is the cube with center x h axes, defined previously. We let

4

h

be the solution of

, side h and sides parallel to

Then it follows from lemmas 3 and 5 and HBlder's inequality that

where

h

is the extension of

Oh as a constant in each cube Ch '

What we would like to show is Lemma 9:

Let

l N/2

.

Then ti, n

+

u

But from theorem 9 and the continuity of $ in R it follows that for N each fixed x E R ch(xy.) Gfx,.) weakly in Lp , 1 1 p < N-2 '

-

-+

Since F

E

L

P

, for some

q > N/2 the theorem follows,

In this same spirit we can prove Theorem 13: Let F

E

L

9

aR

for some q > N/2 and

E

c2

.

.

Then 6 + u uniformly as h -+ 0 h To show this we simply approximate F strongly in L by a sequence 9 with Fn E c~(R) for each n Then F

.

so that

where C does not depend on x and $n

E

V

.

For large n the second

term is small and by the previous remarks

n'

Onh Finally, if F is smooth and

uniformly as h

-+

O

.

aR is smooth we have the analog of

theorem 6. Theorem 14: Suppose F we have for every

where K(E)

E

coy'

and aR

E > 0

is a constant independent of h

.

E

c2

.

Then for problem I

This is similar to theorem 6 for 11. Take u such that 1 u 1

E

"'c

,

Then set u2

=

- U1

.

Aul

-

F

Theorem 14 then follows from

theorems 4 and 6. In order to obtain rate of convergence estimates which are of higher order it is clearly necessary to modify the difference scheme near the boundary. We will for the time being still be considering problem I but with various assumptions on aR and F We shall now redefine Rt, and 3%

a)

c)

Nh(x)

3%

.

.

is the set of "neighbors" of x with respect to

ah

is the set of points on aR which lie on "mesh lines".

For V defined on

aRh

we define

a

,

where uih is the distance between adjacent points of direction 0 < ai

2

-N

where m is an integer and

where k is a multi-index kl = (kl integers, lkl =

1 ki

- . , k.J

)

,

kl,

...,k.J non-negative

and

The fur u and uh the solutions of I and I2 respectively

In this case since we have quite specific knowledge of the behavior of the solution at the origin we obtain an estimate for the error which is point dependent. Thus even if the solution u has the form 2-Nt6 ~ ( x )= 1x1

t

uh

we get from the theorem that

+

regular function

,

6> 0

,

u uniformly on every compact subset

not containing the origin. Also note that under the assumption (13) if u

E

,6 > 0 (m=Z)

then the convergence is second order. This shows clearly that the usual sufficient condition that u

E

4 C (K) is far from necessary for 0(h2)

E

h C (m=0) we obtain a uniform rate of hh

convergence. Note also that when u or h

1-E

,

As might be imagined the theorem is proved by using the representation lemma 3 and estimating the resulting expressions. The details are long and technical and are found in

[7]

but I want to point out the crucial

points. First of all the essential ingredient is the majorant (lemma 4)

.

for Gh

This tells us that, as might be expected Gh behaves quite

like the fundamental solution for Laplace's equation in the neighborhood of the singularity. Having this we are led to proving a sum relation analogous to a well known integral expression Lemma 12:

If

-N < p ,q < 0

and x

,z E

Ix-yl .:ah , I Z - ~::ah, ~ a > O , V y ER,,

Rh

are such that

then

The proof is done by showing that the sum is majorized by C

I

R

Ix-yIP Iz-ylQdy For this formulation, however, we can get a sharper theorem than

either theorem 6 or 15. Theorem 17: Let aR

N

=

2).

E

C* (piecewise, with no reentrant cusps if

Suppose that u is the solution of I1 and uh is the solution

of 1I2 (the analog of I1 for the reformulation). 1 Then for

E

> 0

Let u

E

cPyA(K)

.

The proof is similar to that of theorem 6. Note that with u

E

c2"(F)

cl",

A

the theorem 4 would show only a rate of hA and for

< 1 we would conclude nothing. However, we get second order

-

convergence when the second derivatives are HBlder continuous in R , Clearly these methods are not restricted to these particular difference formulations for the Dirichlet problem for Poisson's equation, One can treat a)

More general operators (second order).

b)

Various boundary conditions.

c)

Eigenvalue problems.

d)

Various difference approximations.

I shall discuss an example of the last extension, since it brings out the fact that in the transition from the interior to a curved boundary one can (in the Dirichlet problem) take approximations which are of the order of accuracy (locally) worse by a factor of hL and still obtain as a global error that of the interior. For this example I choose N = 2 and consider the nine point approximation

so that Ah is now locally 4th order. Now we take A (14) in

R,,

at points of Rh

involved in (14). to

A

.

At

(say

$ = R,, - %

$

)

to be defined by

h

where only I$, points are

we take a second order approximation

aRh we take as in the second formulation. One can then show

by appropriate modifications of the previous Green's function method that if u,, is the solution of our new problem I1 (for Au = 0 ) then we have 3 Theorem 18: Let u be the solution of I1 and u uh

E

C6

(n .

Then if

is the solution of 113

But in fact we can lessen considerably the requirement that u

E

6 C (R)

and obtain by the methods of theorem 17 Theorem 19: Let aR r respectively. Then if u

E

c2

and u and u the solution of I1 and 113 h

CP" (T)

So far the approximations mentioned have all possessed a common property, i.e., that of being "of positive type." This means that if ~~j is the matrix of coefficients of the linear system then

t h e second condition possibly f a i l i n g near the boundary.

In f a c t ,

it is t h i s condition, together with

t h a t makes lemma 1 t r i v i a l . To show t h a t t h i s is j u s t a convenience I wish t o give another example,

Suppose instead of (15) we used the 9-point

It i s possible t o show t h a t the e r r o r is of the order h

and probably a theorem l i k e 19 i s also true.

4

O(h ) approximation

4

when u

E

6 C @)

Since the properties (15) and (16)

a r e not possessed by the resulting system the analysis is much more d i f f i c u l t . It i s i n t e r e s t i n g t o note, however, that a corresponding d i s c r e t e Green's

function w i l l s t i l l be positive, a f a c t which i s no longer completely t r i v i a l .

As i s evident the preceding discussion is i n many respects special f o r second order equations, since much use was made of the maximum principle o r , what is the same, the p o s i t i v i t y of the Green's function.

Thus i t appears

t h a t , i n attempting t o t r e a t higher order equations, we should work more with norms other than the maximum norm.

I would like first to sketch some results of Thornge[lO ] on higher order equations and difference approximations. Consider the differential operator

..

are multi-indices, i.e. B=(B1,. ,BN) , where the B 's j n The a are non-negative integers, ( 6 ( = Bj and similarly for y BY j=l are real constants and

where B and y

1

.

We assume that L is elliptic, i.e. for real 5 =

The Dirichlet problem (111)

has a smooth solution provided F and aR are sufficiently smooth.

Consider approximations of the form

where u = u(Sh) and the C 's are complex numbers defined for all a 5 a but zero except for a finite number of a's A point ([+a)h will be

.

called a neighbor of [h

if Ca f 0

.

This time Rh will be defined as those points of RflENh whose neighbors also lie in R

.

-

Define

\=

.

a\

The characteristic

polynomial of L is defined as the trigonometric polynomial h

where 0 = (O1,...,ON)

, (a,0)

=

1 aj 0j

.

Because of periodicity 0

can

be taken in the set

We say that Lh is consistent with L if at an a&trary

point (which we

take as the origin)

Lh uo

=

+

L u(0) (

Now it can be shown

o(1)

when h

+ O(hk) -

+

0

.

consistant of order k )

Lemma.13:

Lh is consistent with L if and only if

p(B) = L(0)

+

o

( [elzm)

when 0

+

0

.

Now we shall denote the set of complex valued mesh functions on

R,

by Dh and

The sum will always be finite since all functions considered will vanish outside some bounded set. Define

and

Now we call the difference operator Lh elliptic if

p(9) > 0 for 0 # 0

In particular if Lh is elliptic p(9)

E

S

p(9)

satisfies

.

is real so that C-a =

ra .

With this definition, Thom6e then gives two a priori inequalities which we state as the next two theorems. Theorem 20: Let Lh be consistent with L

.

Then Lh is elliptic

if and only if there is a constant C independent of u and h such that

The main tool in the proof is the "Fourier transform." For this reason only constant coefficients are treated.


o) there are solutions of ( 2 2 ) , ( 3 2 ) exponentially increasing with time.

(ii)when

T

G. Capriz

9. The non-linear problem.

We consider here t h e non-linear problem, whose s o l u t i o n represents t h e Taylor v o r t i c e s within our approximation :

with tile boundary conditions ( 2 2 ) . We give p r e c i s e sense t o t h i s problem by s t a t i n g t h e ~ c t where we seek a n o n - t r i v i a l solution: i t i s a s u b s e t a o f a Sanach space E of v e c t o r s (ly ,V ) obtained thus. Consider t h e s e t of functions 'f ( I 5 ) cleflined i n a s t r i p S1 l a r g e r than t h e s t r i p S :06 $ 1, of c l a s s C" i n S1 , periodic i n f with period 2q; introduce i n tile norm

5

3,

5

Y m

and l e t in

-

-

5

3

be t h e closure of with reference t o t h i s norm. Then E-is t h e Banach space of v e c t o r s ( ) V I 'V) with ?y , v i n y l and t h e norm

We consider a l s o t h e s e t of a l l functions f ( 5 , y ) , C m i n S , p e r i o d i c i n S w i t h p e r i d 2q, which vanish i n a s t r i p along t h e boundary of G. The space obtained by closure of t h e s e t with

G. Capriz

yIIE,

.

1:

11 will bc indicated with reference t o t h e norm El' Then i s t h e s e t of v e c t o r s (Y),V ) of C :iitii y 6 H2 e

and V t iil ; i n f a c t vectors such t h a t y 6 V L H ~s a t i s f y , i n a generalized scnsc, t h e boundary conditions at = 0, = 1. \re w i l l look thcn f o r s o l u t i o n s of our 2roblcrn (33) i n

G2,

7

3

I t i s p o s s i b l c t o silov f i r s t of a l l t h a t t l ~ e r ca r e no

e::cept

t h e t r i v i a l one, f o r