151 82 18MB
English-French-Italian Pages Ed. Cremonese, Roma 1968. [446] Year 2011
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