REND.SEM. MAT. UNIVERS. POLITECN. TORINO

Vol. 42*. 1 (1984)

G. Di Lena - D. Trigiante

IL METODO DELLE LINEE IN ANALISI NUMERICA

Introduzione.

La vastissima rivoluzione culturale e tecnologica indotta daH'awento dei calcolatori elettronici ha permesso tra l'altro il trattamento di problemi ma-tematici molto complicati e nello stesso tempo ha esteso I'applicazione della matematica ad altre discipline tradizionalmente lontane da essa.

L'Analisi Numerica, anch'essa enormemente sviluppatasi in questi anni, rappresenta il fulcro intorno al quale ruotano tutte le nuove applicazioni della Matematica. Essa ha pero bisogno di metodi matematici di tipo costruttivo per poterne fare degli algoritmi da implementare sui calcolatori.

Non sempre pero i metodi matematici che dimostrano ad esempio l'esi-stenza della soluzione di un problema sono di tipo costruttivo ed in tal caso, pur non disconoscendo l'importanza teorica del risultato a cui si perviene mediante tali metodi, la loro utilita pratica rimane scarsa.

Vi e certamente una differente preferenza (anche perche differenti sono gli obiettivi) sui metodi di indagine tra Analisti e Analisti Numerici.

Il metodo delle linee, che qui presenteremo dal punto di vista dell'Ana-lisi Numerica, e pero un metodo che sta avendo un successo sempre crescente presso entrambe le categorie di ricercatori e cio perche pur essendo un potente

Questo articolo contiene gli argomenti esposti dal secondo autore in una conferenza tenuta a Torino il 16 giugno 1983.

26

mezzo di indagine teorica, ha caratteristiche spiccatamente "costruttive". Esso nacque circa cinquanta anni fa [1], come un metodo per dimostrare

l'esistenza di soluzioni di equazioni a derivate parziali di tipo parabolico e fu continuato in seguito per gli stessi scopi, specialmente ad opera di autori in lingua russa (ved. [2] e la estesa bibliografia ivi citata). Recentemente molti autori l'hanno usato con successo (ved. [3], [4], [5]). In Analisi Numerica e stato introdotto in epoca molto piu recente, ottenendo molto successo (ved. [6].[7],[8],[9]ecc).

II metodo delle linee.

Questo metodo consiste nell'approssimare una equazione a derivate par­ziali mediante un sistema di equazioni differenziaii'ordinarie, ottenuto sosti-tuendo le derivate parziali rispetto ad una o piu variabili mediante rapporti incrementali.

Nel caso di equazioni paraboliche od iperboliche, a seconda che si di-scretizzi la variabile tempo o la variabile spaziale si ottengono due differenti tipi di metodi ognuno/ dei quali porta ad un problema matematico diverso.

Consideriamo ad esempio il seguente prbblema:

ut = uxx + f(x,t)

(1) w(0,0 = "(1 ,0 = 0 > 0<t<T

w(#,0) = <t>(x)

Possiamo discretizzare la variabile t, ponendo ti = ih, i = 0,1, . . . ,7V, ed h = T/N. Approssimiamo la derivata ut con

U\X,tn + i) U\X,tn)

con 0 < n < N- 1. Si ottiene cosi un sistema di equazioni differenziali ordi-narie

u{x, tn+i) — u{x, tn)

h

con u" = d2 u/dx2. Introducendo i vettori

u \x, tn+i) + f{x, tn+i)

27

u(x) =

u(x,tx)

u{xit2) ' / ( * , * ! )

; F(*) 0(*)

«(X, t^) / ( * » * * ) ,

e la matrice

^AT ~

'+1 0

- 1 +1 0

o '•.. ' •.

0 0

si ottiene il sistema:

(2)

con le condizioni:

u -" = ~r Au~ F(x) — <f>(x)

«(0) = »(1) = 0

<l>(x)

0

Si ottiene cosi un problema ai valori ai limiti. Ogni componente del vet-tore u{x) rappresenta il valore approssimato di u(x,tj) lungo la linea £, parallela all'asse x. Cio spiega 'For igine del nome.

Invece di discretizzare la variabile f, si potrebbe discretizzare la varia­bile x, ponendo Ax = UN e xi -...i'Ax (i = 0,1, . . . ,N}•. + 1), ed approssi-mando la derivata seconda mediante il rapporto

u(xi + 1, t)-2u(xj, t) + u{xi.l,t)

Ax2

Si ottiene cosi il sistema

dii(Xjf t) u(Xj+l, t) — 2U(XJ, t) + U(XJ-1 , t)

dt

u{Xi,0) = <f>(Xi)

Ax'4 + M.0

28

Introducendo i vettori

'«(*i,0 u(x2t t) f(oc2,t)

u(t) ; F(*)

0(^2)

0 =

u(xN, t) f(*N-> *)l 4>(xN)

e la matrice

il sistema puo scriversi nella forma:

(3)

dU 1

s

U(0) =-*

Si ottiene quindi un problema ai valori iniziali. Le componenti del vettore U(t) rappresentano le approssimazioni delia

soluzione u(x, t) lungo le linee 'Xj. Poiche, come si e visto, a seconda cjlel tipo di discretizzazione usato, si perviene a due problemi diversi, si e soliti distiri-guere il metodo delle linee in trasverso (e il primo caso da noi considerato) e longitudinale. II vantaggio nell'usare questo ultimo e chiaro sia dal punto di vista analitico che numerico. Dal punto di vista numerico infatti gli studi dei metodi per approssimare le soluzioni di equazioni differenziali ordinarie svi-luppati negli ultimi venti anni hanno chiarito sufficientemente il problema a valori iniziali e cio spiega quindi la preferenza degli analisti numerici per il metodo longitudinale.

Come abbiamo visto, tale metodo approssima il problema iniziale me-diante un sistema di equazioni differenziali il cui numero di equazioni tende a diventare infinito quando Ax tende a zero. Cio introduce una difficolta in

29

piu rispetto al problema di stabilita dei metodi numerici studiati per un siste-ma di equazioni diff erenziali con un numero fissato di equazioni.

Come e noto per ogni metodo numerico si definisce una regione, detta di Assoluta stabilita, tale che se, il metodo viene applicato al sistema test

y - Ay (4)

y(0)=y0 , .

con A matrice n x n con autovalori a parte reale negativa, il metodo pro­duce una successione di vettori y,- (i = 1,2,...) tale che ||y,-1| < ||y,--ill (ri­spetto ad una norma opportuna) se, indicati con X, (i = 1,2,..., N) gli auto­valori di A, si ha che h-\ siano interni alia regione di Assoluta stabilita.

Nel caso in considerazione in cui le matrici ^ possono essere di di-mensione crescente, la naturale generalizzazione di questa condizione sarebbe quella di richiedere che tutte le quantita b\i(n) (i = 1,2,..., n\ n = 1,2,...) appartengano alia regione di Assoluta stabilita. Cio pero e esatto solo nel caso parabolico, mentre nel caso iperbolico si ottengono condizioni errate (ved. [8]).

Per spiegare tale diversita analizziamo in qualche dettaglio i vari casi.

Caso parabolico.

Per semplicita limiteremo la nostra analisi al caso lineare. Sia dunque

(5)

ut = Au

u(x,0) = <t>(x)

w(0,t) = u(ir,t) = 0

0(0) = .0Or) = O

ove A e l'operatore b/bx2. Come e noto le autofunzioni di questo operatore nella classe delle funzioni definite in [0,7r] ed assumenti valori nulli agli estre-misono sen^^c: (^ = 1,2,...) con autovalori Xfc=-&2.

La soluzione del problema e

(6) u(x,t)= 2 Cb sen kxti. k=i

con oo

(7) <t>{x)= 2 Cksenkx

-kh

30

Discretizziamo ora la variabile spaziale. Sia

7T (8) Ax = Af + 1

e Xj = i Ax, i = 0 ,1 , . . . , N + 1. Usando le stesse notazioni del paragrafo precedente si perviene al proble­

ma semidiscreto:

dU 1 _ (9) dt ~ Ax2 N

U(0) = <f> .

Per poter confron tare la soluzione del problema (5) con quella del pro­blema approssimato (9), introduciamo Toperatore' PN cosi definito: sia v(x) unafunzione definitain [0,7r].

Poniamo

I v(x2) (10) PNv= :

\ v(xN)

Usando opportune norme nei diversi spazi si ha lim H?^ v|| = IMI •

La condizione di consistenza risulta essere

(11) \\m\\BNPNu-PNAu\\ = 0.

Gli autovalori di BN sono

N

Vediamo ora come Toperatore BN approssima Toperatore A sono

'*k = -k2pk„ * =-1,2,... ,AT

kn S € n 2(AT + 1)

Pk" ~ \ kit

2(N + 1)

avendo posto

Poiche

lim Pk = 1

31

si ha che

lim nk = \k . N-+ °°

Cioe lo spettro di BN tende alio spettro dell'operatore d/dx2. Una cosa simile awiene per gli autovettori. Infatti gli autovettori della

matrice BN sono

sen

sen Z , =

sen W + l

e immediatamente si verifica che

Zk = PN sen&#;.

Si puo concludere quindi che la matrice BN approssima bene Toperato-re A nel senso che lo spettro di BN tende alio spettro di A e cosi pure gli autovettori.

II risultato precedente puo essere generalizzato. Si ha infatti il seguente teorema:

Teorema. Sia A un operatore definito in uno spazio Hilbertiano H e sia inoltre:

1) {AN} una famiglia di matrici

2) {PN} una famiglia di operatori di proiezioni tra lo spazio H e lo spazio euclideo iV-dimensionale EN tali che VxGH lim ||/\,a:|| = llffll

N-* °° A/'

3) XN Ed una successione convergente a X tale che esiste una successione di vettori yN G EN in corrispondenza dei quali si abbia:

(12) lim \\AMyN-\NyN\\ = 0

4) VxeH, lim \\PNAx-ANPNx\\ = 0

5) Hx EH tale che P^x = yN

32

Allora si ha che X e autovalore di A corrispondente all'autofunzione x.

Dim. Infatti si ha

\\PN(Ax~- \x)\\<\\PNAx-ANPNx\\ + \\ANPNX- XNyN\\ + \\(XN - \)%\\

e per le ipotesi 3) e 4) si ha che

lim \\PN(Ax-Xx)\\ = 0

da cui per l'ipotesi 2) discende la tesi. In questo caso si puo quindi dire che Toperatore discreto AN rappre-

senta bene l'operatore continuo e non meraviglia che l'uso dei soli autovalori delle matrici ^4N basti a garantire, come vedremo, la stabilita dei metodi numerici.

Caso iperbolico.

Distinguiamo due sottocasi:

a) Problema iperbolico del primo ordine. b) Problema iperbolico del secoridb ordine.

Nel caso a) il problema tipico che considereremo sara:

(13) u(Q,x)=g(x) . x>0 u(t,0) = 0

mentre il problema tipico del caso b) sara

0<x<ir (14) u(0,x) = g(x) u{t,0)-- = U{t, 7T) = 0

Iprobleihi discreti sono rispettivamente • t

dU 1

(15) > dt Ax

U(0) = 4>

ANU

33

d2U = J__ (16) dt2 ~ Ax2

N

U(0) = $ .

Lo spettro dell'operatore - "3— , anche nei casi in cui esiste, non carat-

terizza la soluzione del problema (13). II problema approssimato (15) pero e

earatterizzato fortemente dagli autovalori della matrice -~r-AN. Non dob-

biamo quindi aspettarci che l'uso dei soli autovalori di questa matrice possa caratterizzare la stabilita e la convergenza dei metodi discreti usati per appros-simare la (15) e difatti una condizione di stabilita ottenuta per tale via e meno restrittiva di quella ottenuta da C.F.L. (vedi [8]).

82

Nel caso b) invece gli autovalori dell'operatore -r—j" caratterizzano la

soluzione ed un ragionamento analogo a quello gia fatto nel caso parabolico puo far si anche in questo caso.

Famiglie di matrici.

Abbiamo visto come in alcuni casi la condizione ottenuta usando la no-zione di Assoluta stabilita e gli autovalori delle matrici non e quella corretta. L'uso della nozione di Assoluta stabilita e pero molto comodo perche rende facilmente intuibili le proprieta di stabilita dei metodi e quindi in definitiva diventa un aiuto prezioso nello scegliere a priori un metodo adatto al parti-colare problema da risolvere. E possibile, come ora vedremo, salvare questo concetto, introducendo lo spettro ampliato di una famiglia di matrici.

Sia dunque {A„}nEN+ una famiglia di matrici.

Definizione. Dicesi spettro della famiglia {An} l'insieme dei numeri com-plessi X tali che, per ogni e > 0, esistono nEN+ e x G P , x i= 0 in cor-rispondenza dei quali si ha:

(17) H ^ x - X x I K e | |*| | .

Lo spettro di {An} sara indicato con S({An}) (oppure 5).

Owiamente l'insieme degli autovalori di tutte le matrici An, che indicheremo con 2 , e contenuto in S. Gli elementi S\X saranno detti quasiautovalori.

Elenchiamp alcune proprieta, rimandando a [8] per le dimostrazioni:

34

Pl. S c chiuso.

P2 : se le matrici An sono normali, allora S = 2 .

P3 : se {An} e {Bn} sono due famiglie di matrici tale che _ .

lim \\An-BH\\ = 0

allora le due famiglie hanno gli stessi quasi-autovalori.

P4 : se {An}t {Bn} e {Tn} sono tre famiglie di matrici tale che VnGN+

si abbia

ese F«€'N*, 3r>0 tale che

. I ir ;11| || Tn\\<r

allora

S({An}) = Sa'Bn})>

P5 : se P e un compatto nel piano complesso non contenente alcun elemento di S({An}), allora

sup HLz-AJ-'WK + oo .

zEP

P6 : sia S({An}) compatto e contenuto in un insieme aperto e semphcemen-te connesso £2. Iholtre sia / una funzione analitica in £2, allora

S({f(An)})=f(S({An})).

P7 : sia q — sup j X |. Se q < 1 si ha Vn, m EN+

xes

I U O < *

essendo k indipendente da m ed n.

P8 : se q>l, Vn>0, 3m,nEN+ tale che \\A™\\>k.

Come esempio di spettro di una famiglia di matrici consideriamo il seguente caso, il quale comprende entrambe le famiglie gia considerate. Sia

35

(18) Ln

a 7 0

P . a 7 0

0

Sussiste il seguente teorema:

0

0 'P nxn

Teorema. Lo spettro della famiglia di matrici Cn e formato dai numeri com-plessi' X tali che ie radici dell'equazione

(19) yr2 + ( a - X ) r + 0 = O

siano in modulo minori od uguali a 1.

Dimostrazione. Non e restrittivo supporre |0/7l < 1. Sia x il vettore forma­to dalle prime n componenti della successione definita recursivamente dalla equazione aile differenze del secondo ordine:

(20)

con xl = 1, x2 =

Siha:

yxi+i -(a-X)Xj 4-/3^-1 = 0

X ~ a

w-l .n-\ \\Cnx-\x\\ = \$xn-l + (*-\)xn\ = !7*w+il = MciTi • + c2r2 l)\

ove rx ed r2 sono le radici della (19). Consideriamo il caso in cui rxi=r2

ed entrambe minori di 1. Si ha

1 f\-a-yr2\

rx - r2 \ 7 1 f\-a.-yrx

c2 u-r e quindi

36

(21) i 7 * » + i l =

1

ry-r2

.n-\ (\-Qi-yr2)r" * 4- (k-a-yr^r^ n~\

< .n—\ « - i ( X - a - 7 r 3 ) r p + ( X - a - 7 r i ) r r |H*II .

essendo ||#|| > 1. Nel caso in cui r2=rx e | r t | < l si ha:

(22) l7*»+il< |[7^i + « (X-a -7) ] r w ~ 1 | l l ^ l l .

E ovvio che se una delle due radici ha modulo maggiore di 1, allora il corrispondente X non fa parte dello spettro.

Diamo ora una espressione esplicita dello spettro. Dalla equazione (19) risulta che lo spettro S({Cn}) e dato dai numeri

complessi /

(23) \ = a + (yr + —)

con | r | < l e 18 yr

< 1. Le ultime condizioni discendono dai fatto che le

radici della (19) devono essere di modulo minore od uguale ad 1. Quindi ~~~

0 (24)

Poniamo

Sostituendo si ha

< l r | < l .

r = q(cos6 4- isend) .

(25) X = a + [yq 4- — Jcos-0 + [yq - ) sen.0

con 6<6<n

<q<l

Ponendo X = x 4- iy si ha che la frontiera delrinsieme S e data da

(26) x = a 4- (P + 7) cos 0

y = (0 ~ 7) sen 0

37

Alio stesso risultato si perviene se -r- < 1.

Tale frontiera e una ellisse di centro (a, 0) e semiassi (j3 4- 7) e (0 - 7). Distinguiamo tre casi importanti:

1) a = - 2, 0 = 7 = 1 . In questo caso abbiamo proprio la matrice Bn consi-

derata nel paragrafo precedente ottenuta discretizzando l'operatore -r-~T • L'ellisse precedente degenera nel segmento reale-.

(2,7) * + 2 = 2cos0

0<6<n

se si tiene conto del fatto che gli autovalori della matrice Bn sono dati da

2.kn (28) Xk = - 2 + 2/3 cos , k = l,2,...,n

si ottiene una conferma della proprieta P2.

2) a = + 1, j3 = — 1. In questo caso abbiamo proprio la matrice An ottenuta

discretizzando l'operatore -r— . L'ellisse degenera in un cerchio di centro

1 e raggio 1. Da notare che in questo caso tutte le matrici della famiglia hanno autovalore 1, mentre lo spettro e un insieme di valori di misura non nulla. Si pub usare la P6 per ricavare lo spettro di {~An} oppure di {zl + A„} con 2 qualunque.

3 ) a = 0,.j3 = - l , 7 = 1 . Queste famiglie di matrici si ottengono quando si

discretizza l'operatore -z— con il metodo leap-frog. In questo caso l'ellisse

degenera nel segmento dell'asse immaginario compreso tra —2i e 2i.

Stabilita mediante lo spettro di una famiglia.

Lo spettro di una famiglia di matrici risulta essere la nozione che sosti-tuisce gli autovalori per ottenere la esatta condizione di stabilita. Richiedendo infatti che tutto lo spettro sia contenuto nella regione di Assoluta stabilita del metodo si ottiene la stabilita uniforme rispetto alio spazio e al tempo (ved. [4]).

Dalla proprieta P2 si ha che se le matrici sono normali le condizioni ot-

38

tenute in questo modo sono le stesse di quelle che si otterrebbero usando i soli autovalori ed infatti alcuni autori (ved. [9]) avevano gia ottenuto tali condizioni che si accordavano con le classiche condizioni di Courant-Fried-richs-Lewy. _,

Nel caso di matrici non normali (ad es. caso iperbolico del primo ordi-ne) lo spettro e notevolmente piu ampio deirinsieme degli autovalori. La condizione ottenuta e quindi piu restrittiva e coincide con quella di Courant, ottenuta direttamente.

Analizziamo in dettaglio i risultati per ciascun esempio presentato nei paragrafi precedenti.

a) Caso parabolico. La famiglia di matrici {BN} ha, spettro dato da:

(29) S= {x = 2(-l + cos0), O < 0 < T T }

Se si usa ii metodo di Eulero esplicito per integrare l'equazione differenziale, poiche questo metodo ha come regione di Assoluta stabilita il cerchio

(30) D = { z e C | | z + l | . < l }

si deve richiedere che

(31)

da cui risulta

(32)

che e la condizione di Courant-Friedrichs^Lewy. Owiamente se si usa il metodo di Eulero implicito, essendo quest'ulti-

• At • . ' mo /4-stabile non vi e da imporre nessuna condizione su ~r~r . Si ha cosi la stabilita incondizionata per questo metodo.

b) Caso iperbolico del primo ordine. La famiglia di matrici - ^ na» come abbiamo visto, spettro dato da

(33) 5 = {2GCH2 + I K I } .

Se si usa il metodo di Eulero esplicito, bisogna imporre che

At (34) ^ S CD

At Ax2 SCD

At Ax2 .2

39

e quindi A*

<35) A ^ 1

ottenendo ancora una volta la esatta condizione di stabilita.

c) Caso iperbolico del secondo ordine. Questo caso richiede un ulteriore ap-profondimento, perche la condizione di stabilita non e direttamente applica­ble essendo il sistema di equazioni differenziali ordinarie (16) del secondo ordine.

Introducendo i vettori w e \frEE2N cosi definiti

dU (V \ / * \

(36) v - - , W«(J, •..(J e la matrice a blocchi D2N

(37) Dw = \BN

il sistema (16) puo porsi nella forma:

dw _ 1 (38) ~d7"Ax2^2N D-KTW

w(0) = \l/

La famiglia di matrici {D2N} non e norrnale, ma e normale la famiglia {D\N} essendo

/AX2BN 0

(39) D\N=[ V 0 Ax2Bl

Poiche la famiglia {DlN} ha spettro St dato da Ax2S ove 5 e dato da (33), applicando la P6 si ha che la famiglia {D2N} ha spettro

(40) s\/2=iAx\S\in

essendo

\S\1'2 = {x = ± V 2 ( l - c o s 0 ) , 0 < 9 < IT]

Essendo lo spettro un segmento dell'asse immaginario, e immediato stabilire quali metodi scegliere per discretizzare la variabile tempo.

40

E necessario infatti che il metodo scelto abbia l'asse immaginario sulla frontiera della regione di Assoluta stabilita. Tra i metodi con questa proprie­ty, i piu noti sono il metodo leap-frog e il metodo dei trapezi. Quest'ultimo,

essendo A -stabile, non pone condizioni sul rapporto ~r—. Il metodo leap-

-frog (ved. caso 3 del paragrafo precedente) invece richiede che

At 2 < 2 Ax

e cioe

At ( 4 1 ) A^< 1-

B I B L I O G R A F I A

[1] E. Rothe: Zweidimensionale parabolische Randwertaufgalen als Grenzfall eindimen-sionaler Randwertaufgaben. Math. Ann. 102 (1930) 650-670.

[2] O.A. Liskovets: The Method of Lines (Review). Differential Equations 1 (1965) 1308-1323.

[3] C. Gorduneanu: Some Applications of Rothe's Method to Parabolic and Related Equations. Applied Non Linear Analysis (1979) 111-121.

[4] W.Walter: The Line Method for Parabolic Differential Equations. Lecture Notes in Mathematics, N. 430, Springer-Verlag (1974).

[5] Kannan, D. Trigiante: Esistenza di soluzioni periodiche per problemi parabolici ed iperbolici con il metodo delle linee. In preparazione.

[6] R.F. Albrecht: Approximation to the solution of partial differential-equations by the solutions of Ordinary Differential Equations. Num. Math. 2 (1960) 245-262.

[7] R.F. Sincovec, N.K. Madsen: Software for Non linear Differential Equations. ACM transactions on Math. Soft. 1 (1975) 232-260.

[8] G. Di Lena, D. Trigiante.- On the stability and convergence of lines method. Rendi-conti di Matematica, 3 (1983) 113-126.

41

[9] G. Dalquistt Problems Related to the Numerical Treatment of Stiff Differential Equations. In A. Giinter editor - International Computing Symposium. (1973) 307-314, North Holland.

G.DI LENA - D. TRIG I ANTE - Dipartimento di Matematica, Universita di Bari, Via Nicolai, 2 - Bari.

Lavqro pervenuto in redazione il 16/V1/1983