UN METODO POLINOMIALE PER IL CALCOLO DI FUNZIONI DI ...alvise/TESI_STUDENTI/CACCO.pdf · loro...

UN METODO POLINOMIALE PER IL CALCOLO DIFUNZIONI DI MATRICI NON SIMMETRICHE BASATO

SUI PUNTI DI LEJA

Jacopo Cacco

Indice

1 Introduzione 5

2 Funzioni di matrici 72.1 Tecniche per il calcolo approssimato di funzioni di matrici . . . . . . . . . . 8

2.1.1 Serie di Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1.2 Diagonalizzazione di matrici . . . . . . . . . . . . . . . . . . . . . . . 92.1.3 Decomposizione di Schur . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2 Alcuni esempi di funzioni di matrice . . . . . . . . . . . . . . . . . . . . . . 102.2.1 Radice quadrata . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.2 Esponenziale di matrice . . . . . . . . . . . . . . . . . . . . . . . . . 112.2.3 Logaritmo di matrice . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Approssimazione con i punti di Leja 133.1 Weakly Admissible Meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . 143.2 Algoritmo DLP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4 Descrizione del metodo e risultati 214.1 Il metodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

4.1.1 Costruzione dell’ellisse . . . . . . . . . . . . . . . . . . . . . . . . . . 214.1.2 Costruzione della mesh . . . . . . . . . . . . . . . . . . . . . . . . . 224.1.3 Algoritmo DLP: estrazione dei punti di Leja . . . . . . . . . . . . . . 234.1.4 Interpolazione delle funzioni di matrici con il metodo di Newton . . 25

4.2 Risultati: un problema di convezione-di!usione . . . . . . . . . . . . . . . . 26

A Codici MATLAB 29A.1 Codice 1: Costruzione dell’ellisse . . . . . . . . . . . . . . . . . . . . . . . . 29A.2 Codice 2: Costruzione della mesh . . . . . . . . . . . . . . . . . . . . . . . . 30A.3 Codice 3: Estrazione dei Punti di Leja . . . . . . . . . . . . . . . . . . . . . 31A.4 Codice 4: Interpolazione con il metodo di Newton . . . . . . . . . . . . . . 33

3

4 INDICE

Capitolo 1

Introduzione

Data una matrice A ! Rn!n, e una funzione di matrice f analitica in una certa regionedel piano complesso, ci poniamo il problema di calcolare f(A) ottenendo l’interpolazionepolinomiale di f con il metodo di Newton. I nodi di interpolazione saranno i cosiddettipunti di Leja, ottenuti tramite l’algoritmo DLP (Discrete Leja Points). Tali punti vengonodefiniti asintoticamente ottimi in quanto asintoticamente garantiscono che la convergenzauniforme sia la stessa dei polinomi di miglior approssimazione. Inoltre hanno una costantedi Lebesgue che cresce al piu linearmente con il grado di interpolazione. Utilizzeremo ipunti estratti dall’algoritmo DLP come nodi per il polinomio interpolante, e ricaveremouna ricorrenza a tre termini che permettera di calcolare in maniera semplice il valore delpolinomio interpolante calcolato in A. Produrremo infine alcuni risultati numerici con dellematrici e delle funzioni test. Le funzioni sui quali e!ettueremo i test saranno note funzionianalitiche (esponenziale, radice quadrata, logaritmo), sulle quali esistono gia degli altrimetodi per il loro calcolo numerico, e sulle quali sara quindi possibile confrontare i nostririsultati.Questo metodo, che ovviamente produce un calcolo approssimato di f(A), ha il vantaggioche, dal punto di vista computazionale, ha un costo in operazioni pari a quello di un pro-dotto di matrici ad ogni passo dell’iterazione (mn3, dove n e la dimensione della matrice em e il numero delle iterazioni del metodo).

Nel Capitolo 2 introduciamo la teoria necessaria sulle funzioni di matrici, ricordando illoro legame con l’Integrale di Cauchy per funzioni di una variabile complessa, e successiva-mente descrivendo le principali tecniche gia esistenti per il calcolo approssimato di radicequadrata, esponenziale e logaritmo.

Nel Capitolo 3 vediamo in dettaglio l’approssimazione tramite i punti di Leja discreti. Inparticolare mostriamo come sono ottenuti numericamente a partire da una discretizzazio-ne di un insieme compatto, detto Mesh debolmente ammissibile, e perche sono importanti

5

6 CAPITOLO 1. INTRODUZIONE

nell’ambito della teoria dell’approssimazione.

Nel Capitolo 4 mostriamo come si possono valutare funzioni di matrici con interpolantinei punti di Leja, evidenziando la formula ricorsiva ottenuta dal metodo di Newton perl’approssimazione polinomiale.Successivamente esibiremo i risultati dei test numerici, sulla base di un problema fisico diconvezione-di!usione.Osservermo che la convergenza del metodo e garantita per le funzioni che utilizzeremo cometest, e che dal punto di vista computazionale il metodo e indubbiamente e"ciente (sia perla costruzione della mesh, sia per la formula ricorsiva); noteremo tuttavia che a gradi diinterpolazione alti, e con una scelta della base polinomiale poco accurata, si possono veri-ficare problemi di malcondizionamento che possono rallentare notevolmente la convergenza.

Capitolo 2

Funzioni di matrici

Data una matrice A ! Rn!n e una funzione scalare f , definiamo l’immagine di f tramite Ain modo tale da generalizzare la definizione di funzione in una variabile complessa, f(z), z! C. Cio ci fornisce numerosi vantaggi, come la possibilita di sostituire l’operazione di di-visione in ambito scalare con la moltiplicazione per la matrice inversa in ambito matriciale.Si consideri ad esempio la funzione scalare

f(z) =1 + z2

1" z2, t #= 1 (2.1)

Risulta naturale definire la corrispondente funzione di matrice

f(A) = (I "A)"1(I +A2), 1 /! !(A) (2.2)

dove !(A) corrisponde allo spettro di A, ossia all’insieme degli autovalori di A.Diamo quindi la definizione di una funzione di matrice come la generalizzazione del Teoremadell’intergrale di Cauchy:

Definizione 1. Data una matrice A ! Cn!n di spettro !(A) e una funzione f che siaanalitica in un insieme compatto K, e tale che K $ !(A) , definiamo

f(A) :=1

2"i

!

Kf(z)(zI "A)"1dz (2.3)

Ci poniamo quindi il problema di calcolare numericamente l’immagine tramite f di unagenerica matrice A. Dal punto di vista computazionale e un problema molto complesso,sia per la dimensione di A, che puo essere molto grande, sia perche non sempre e possibilecalcolare esplicitamente l’intergrale di Cauchy.

E’ bene ricordare che esistono altri modi per definire la funzione di matrice. I piu noti,oltre a quello appena descritto, sono quelli che sfruttano la forma canonica di Jordan el’interpolazione di Hermite (cf. [1]). Queste ulteriori definizioni sfruttano delle proprieta

7

8 CAPITOLO 2. FUNZIONI DI MATRICI

algebriche delle matrici che permettono di ricondurre il calcolo di f(A) al calcolo di f(z),cioe al calcolo di una funzione scalare. Tali definizioni sono comunque equivalenti alla (2.3),pertanto i risultati che otterremo in questo lavoro sono indipendenti dalla definizione dif(A).Dal punto di vista applicativo, il calcolo di f(A) trova applicazione in numerosi campi, nonsolo prettamente matematici.Si pensi ad esempio alla funzione esponenziale, e al ruolo fondamentale che essa giocanell’ambito della risoluzione delle equazioni di!erenziali lineari. Oppure, nel campo dellaprobabilita e statistica, di notevole interesse e il calcolo di esponenziale e logaritmo dellematrici stocastiche (particolari matrici con somma di righe e colonne pari a 1). O anco-ra, in algebra lineare, capita spesso il problema di dover calcolare la radice quadrata diuna matrice Hermitiana definita positiva (decomposizione di Cholewsky), o di trovare lasoluzione di equazioni non lineari del tipo XAX = B.

2.1 Tecniche per il calcolo approssimato di funzioni di ma-trici

Vediamo brevemente alcune delle tecniche piu comunemente utilizzate per approssimare lefunzioni di matrici.

2.1.1 Serie di Taylor

Quando la funzione f e analitica, ossia derivabile un numero infinito di volte nel campocomplesso, essa e sempre sviluppabile come serie di potenze:

f(z) =#"

k=0

ak(z " #)k, ak =f (k)(#)

k!(2.4)

Vale inoltre il seguente teorema di convergenza, che riportiamo per il caso matriciale:

Teorema 1. [6, p. 76] Sia f una funzione analitica con espansione in serie di potenze (odi Taylor) come in (2.1) e con raggio di convergenza r. Data allora una matrice A ! Cn!n,e definita la f(A) come serie di potenze, cioe da

f(A) =#"

k=0

ak(A" #I)k, (2.5)

se e solo se ognuno degli autovalori $i di A soddisfa una delle seguenti condizioni:

• |$i " #| < r

• |$i " #| = r e la serie fni"1($) (dove ni e l’indice di $i) e convergente nel punto$ = $i

2.1. TECNICHE PER IL CALCOLO APPROSSIMATO DI FUNZIONI DI MATRICI 9

Il problema dell’approssimazione tramite la serie di Taylor e dovuta al fatto che siha a che fare con delle somme infinite che sono, ovviamente, incalcolabili. Pertanto siutilizzano le somme parziali, utilizzando un numero finito ma su"cientemente grande ditermini. A tal proposito, il seguente teorema fornisce un limite superiore all’errore dovutoal troncamento della serie di Taylor:

Teorema 2. [6, p. 77] Sia f sviluppabile in serie di potenze come in (2.1), con raggio diconvergenza r. Se A ! Cn!n, con %(A " #(I)) < r, allora definita una qualunque normamatriciale vale la seguente disuguaglianza:

%f(A)"s"1"

k=0

ak(A" #I)k% & 1

s!max0$t$1

%(A" #I)sdsf

dts(#I + t(A" #I))%. (2.6)

2.1.2 Diagonalizzazione di matrici

Se la matrice A e diagonalizzabile, cioe se esiste una matrice diagonale D = diag($i),tale che A = ZDZ"1, dove i $i sono gli autovalori di A e la matrice Z ha per colonnegli autovettori di A, allora il calcolo di f(A) si riduce al calcolo di f($i), poiche valef(A) = Z diag(f($i))Z"1. Tuttavia questo metodo funziona molto bene solo con una certaclasse di matrici, ossia le matrici normali (matrici che commutano con la propria traspostaconiugata, tra cui le matrici simmetriche) altrimenti ci possono essere dei problemi dicondizionamento.

2.1.3 Decomposizione di Schur

Data una matrice A, la sua decomposizione di Schur e definita come QTAQ = T , con Qunitaria e T triangolare superiore, con in diagonale gli autovalori di A. Tale decomposizionepuo essere prodotta con l’algoritmo QR. In questo modo, si ha che f(A) = Qf(T )QT , equindi il problema del calcolo di f(A) si riduce a quello di calcolare f(T ). Ci viene allorain aiuto il seguente risultato:

Teorema 3. [6, p. 84] Sia T ! Cn!n una matrice triangolare superiore, e sia la funzionef definita sullo spettro di T. Allora f(T ) = F e triangolare superiore con fii = f(tii) e

fij ="

so,...,sk%Sij

ts0,s1ts1,s2 . . . tsk!1skf [$s0 , . . . ,$sk ], (2.7)

dove $i = tii, Sij e l’insieme di tutte le successioni crescenti di interi che iniziano in i efiniscono in j, e f [$s0 , . . . ,$sk ] e la di!erenza divisa di ordine k di f in $s0 , . . . ,$sk .

Il teorema appena enunciato produce, tra l’altro, delle formule esplicite per il calcolodelle fij nelle forme canoniche di Jordan e dei blocchi triangolari 2x2 (Teorema di Kenny-Laub, [6, p. 76]).


Il problema di questo metodo consiste nel fatto che ha un costo computazione elevatissimo(addirittura O(sn)); questo problema e stato parzialmente risolto dall’algoritmo di Parlett(cf. [6, p. 86], che riduce il numero di operazioni a O(2n3/3). Tuttavia, per come e costruitol’algoritmo, esso non funziona nel caso in cui la matrice T abbia autovalori ripetuti. Perovviare questo problema si usa utilizzare l’algoritmo suddividendo T in blocchi diagonali2x2.

2.2 Alcuni esempi di funzioni di matrice

In questa sezione vedremo brevemente come vengono utilizzati alcuni di questi metodi perle funzioni di matrice quadrata, esponenziale e logaritmo.

2.2.1 Radice quadrata

La radice quadrata di una matrice non e, in genere, unica. Generalmente tuttavia, datauna matrice A si e interessati a calcolare l’unica radice quadrata X definita positiva. Unatale matrice non e sempre definita, infatti notiamo che dalla definizione di funzione dimatrice si ha '

A =2

"A

! #

0(t2I +A)"1dt. (2.8)

Vediamo subito quindi che una delle condizioni per poter calcolare questa funzione dimatrice e che lo spettro di A non contenga autovalori con parte reale negativa.Vale inoltre il seguente Teorema riguardante l’esistenza di radici quadarate reali:

Teorema 4. [6, p. 138] Sia A ! Rn!n non singolare. Se A possiede autovalori negativi,allora non possiede radici quadrate reali definite positive. Se A non ha autovalori realinegativi, allora esistono esattamente 2r+c radici quadrate di A, dove r e il numero diautovalori reali distinti di A, mentre c e il numero di coppie distinte di autovalori complessiconiugati.

I principali metodi utilizzati per calcolare la radice quadrata di A sono basati sulladecomposizione di Schur, gia vista nel precedente paragrafo, e il metodo di Newton.Il metodo di Schur, in particolare, prevede la suddivisione della matrice in blocchi 2x2.

Nel metodo di Newton, invece, si suppone che Y sia una soluzione approssimata dell’e-quazione X2 = A, e che esista una qualche matrice E tale per cui Y + E = X. AlloraA = (Y + E)2 = Y 2 + Y E + EY + E2. Trascurando il termine di secondo ordine E2

abbiamo la seguente ricorrenza:

X0, XkEk + EkYk = A"X2k , Xk + 1 = Xk + Ek, k = 0, 1, 2, . . . (2.9)

2.2. ALCUNI ESEMPI DI FUNZIONI DI MATRICE 11

In questa forma, il metodo di Newton e decisamente dispendioso, dovendo ad ogni iterazionerisolvere un equazione nell’incognita Ek. Con il seguente lemma tuttavia siamo in gradodi ridurre notevolmente il numero di conti:

Lemma 1. [6, p. 139] Supponiamo che nell’iterazione di Newton X0 commuti con la ma-trice A, e che le iterate Xk siano tutte ben definite. Allora, per ogni k, Xk commuta conA e Xk+1 =

12(Xk +X"1

k A).

Questo ci permette di riscrivere in maniera piu semplice il metodo sopra descritto,ponendo X0 = A, che commuta ovviamente con se stessa.

Osservazione 1. Per poter usare questo metodo, che sfrutta l’inversa della matrice Adobbiamo richiedere che lo spettro non contenga l’autovalore 0; questa sara l’ipotesi cherichiederemo anche nel metodo polinomiale che spiegheremo nei capitoli successivi. Esisto-no comunque degli altri metodi iterativi che garantiscono la convergenza anche per matricisingolari.

2.2.2 Esponenziale di matrice

Ricordiamo che la matrice esponenziale e definita nella seguente maniera:

eA =#"

k=0

Ak

k!(2.10)

e che la serie di potenze relativa a tale funzione ha raggio di convergenza infinito. Quello del-lo sviluppo di Taylor e il metodo usato dallo stesso MATLAB per il calcolo (approssimato)della matrice esponenziale tramite il comando expm, che utilizzeremo successivamente perla verifica dei nostri risultati numerici. Un altra possibilita e la rappresentazione tramiteil limite

eA = lims&#

(I +A

s)s. (2.11)

Anche per questa funzione, inoltre, esistono degli algoritmi basati sulla decomposizione diSchur. Uno di questi in particolare utilizza le di!erenze divise di Newton per il calcolo dif(T ) (si veda il paragrafo 2.1.3). Il problema che si presenta nel calcolo delle di!erenzedivise si ha quando, ad esempio, i nodi di interpolazione $i e $i+1 sono molto vicini tra diloro. Dato quindi che la di!erenza divisa f [$i,$i+1] si calcola come f [$i,$i+1] = (f($i+1)"f($i))/($i+1"$i), possono verificarsi problemi di accuratezza nel calcolo del denominatore,che quindi si amplificano nel valore finale della frazione. Tuttavia, per ottenere un calcolopiu accurato delle di!erenze divise, si puo utilizzare un risultato dovuto a Opitz [1,pag.250].


2.2.3 Logaritmo di matrice

Il logaritmo di A ! Rn!n e qualunque matrice X tale che eX = A. Se A e non singolare, essapossiede infiniti logaritmi, questo perche nel campo complesso il logaritmo e una funzioneperiodica di periodo 2"i. Quello che calcoleremo sara quindi il logaritmo principale di A,ossia la matrice X il cui spettro e contenuto nella striscia {z : "" < ((z) < "}. Perpoter definire il logaritmo di A, inoltre, e necessario che i suoi autovalori abbiano partereale non negativa. Un primo metodo per il calcolo del logaritmo sfrutta le serie di Taylor.Ricordiamo che valgono le seguenti relazioni:

log(I +A) = A" A2

2+

A3

3" . . . , %(A) < 1 (2.12)

e la cosiddetta Serie di Gregory (1668), che si costruisce a partire dal caso scalare sottraendogli sviluppi di Taylor di log(1 + x) e di log(1" x):

log(A) = "2#"

k=0

1

2k + 1((I "A)(I +A)"1)2k+1, min

i)($i(A)) > 0 (2.13)

Notiamo subito che la Serie di Gregory ha un raggio di convergenza molto piu granderispetto alla serie di Taylor, e in particolare essa converge per qualunque matrice hermitianadefinita positiva. Tuttavia la convergenza puo essere molto lenta se A non e abbastanzavicina all’identita.Altri metodi per il calcolo del logaritmo sfruttano la decomposizione di Schur. Come nelcaso dell’esponenziale, anche qui si ripresenta il problema del calcolo della di!erenza divisaf [$1,$2], quando $1 e $2 sono vicini tra loro. Per evitare questo problema si opera unasostituzione di questo tipo:

log $2 " log $1 = log$2

$1+ 2"i Avv(log $2 " log $1) (2.14)

dove Avv e l’indice di avvolgimento. Poniamo Avv(log $2 " log $1) = U e z = ($2 "$1)/($2 + $1) e otteniamo che la (2.11) equivale a

log

#1 + z

1" z

$+ 2U (2.15)

Infine, ricordando che arctanh(z) = 12 log(

1+z1"z ), concludiamo che

f [$1,$2] = t122 arctanh(z) + 2U

$2 " $1(2.16)

Capitolo 3

Approssimazione con i punti diLeja

In questo capitolo descriveremo un algoritmo e"ciente per l’interpolazione polinomiale, daapplicare poi alle funzioni di matrici. Come visto nel Capitolo 2, possiamo ricondurci alcaso delle funzioni scalari complesse in una variabile.L’algoritmo che utilizzeremo ci fornira i Discrete Leja Points (di seguito abbreviato conDLP), i cui vantaggi sono i seguenti:

• la costante di Lebesgue #n cresce moderatamente con n (vedi paragrafo successivo);

• se {zj}j=1,...,n e un set di punti di Leja di grado n"1 per il dominio $ * C, allora le co-stanti di Lebesgue del sottoinsieme {zj}j=1,...,m crescono anch’esse “moderatamente”per m & n.

Vediamo subito come la crescita moderata di #n renda i DLP dei buoni punti di ap-prossimazione. Ricordiamo che definito p'n il polinomio di miglior approssimazione per lafunzione f , per un noto teorema (cf. [9]) vale la relazione:

%f " pn% & (1 + #n)%f " p'n% (3.1)

per ogni polinomio pn ! Pn(K). Sia quindi pn il polinomio interpolante nei punti di Leja;estraendo la radice n-esima e passando ai limiti otteniamo

limn&#

%f " pn%1/n & limn&#

(1 + #n)1/n%f " p'n%1/n (3.2)

che, vista la crescita moderata della costante di Lebesgue con n, equivale a

limn&#

%f " pn%1/n & limn&#

(1 + n)1/n%f " p'n%1/n = limn&#

%f " p'n%1/n (3.3)

Dato che p'n e il polinomio di miglior approssimazione, non puo verificarsi la disuguaglianzastretta, quindi il polinomio pn, che interpola sui DLP, e asintotico al polinomio di migliorapprossimazione.

13

14 CAPITOLO 3. APPROSSIMAZIONE CON I PUNTI DI LEJA

3.1 Weakly Admissible Meshes

Dato un insieme compatto K * C, i polinomi di grado massimo n ristretti a K formanouno spazio vettoriale Pn(K) di dimensione n + 1. Preso quindi un sottoinsieme di puntidistinti di K, A= {a1, . . . , an+1}, e una funzione f : K + C, vogliamo trovare un polinomiopn ! Pn(K) tale che

pn(ai) = f(ai), i = 1, . . . , n+ 1 (3.4)

Cio equivale a risolvere il seguente sistema di n+ 1 equazioni lineari:

pn(ai) =n+1"

j=1

cjPj(ai) = f(ai), i = 1, . . . , n+ 1 (3.5)

dove le costanti cj ! C sono le incognite del sistema, e i Pj sono i polinomi che costituisconouna base di Pn(K). In forma matriciale questo sistema si scrive nella forma

V (a;P )c = F, (3.6)

dove V (a;P ) = [Pj(ai)]1$i,j$n+1 e una matrice di Vandermonde, c e il vettore dei coe"-cienti e F e il vettore dei valori f(ai). Se i punti di A sono a due a due distinti, la matricedi Vandermonde e non singolare, e quindi possiamo definire il polinomio interpolatore diLagrange come

li(z) =detV (a1, . . . , aj"1, z, aj+1, . . . , an+1;P )

detV (a;P ), i = 1, . . . , n+ 1, (3.7)

con li(aj) = &ij , e da cui otteniamo il polinomio interpolante i punti dell’insieme A

pn(z) =n+1"

i=1

f(ai)li(z). (3.8)

Ora, se l’insieme dei punti A* K, di cardinalita n + 1, e tale da massimizzare il valoreassoluto del determinante della matrice V (a;P ) — tali punti sono detti punti di Fekete, (cf[2]) — dato che |li(z)| & 1 e li(ai) = 1, la costante di Lebesgue (cf. [11]), per definizionesara

#A := maxz%K

n+1"

i=1

|li(z)| & maxz%K

n+1"

i=1

1 = n+ 1 (3.9)

I punti di Fekete, quindi, data la crescita moderata della costante di Lebesgue, costituisco-no un buon insieme di punti di interpolazione, e si nota che gli stessi punti e la costantedi Lebesgue sono indipendenti dalla base polinomiale scelta. Il problema e che, dal pun-to di vista computazionale, questi punti nel continuo sono di"cili da calcolare, perche lamassimizzazione del determinante della matrice di Vandermonde e un problema di otti-mizzazione non lineare da applicare ad una matrice di dimensione arbitrariamente grande

3.1. WEAKLY ADMISSIBLE MESHES 15

e che coinvolge n + 1 variabili complesse. Tale problema e classificato come NP-hard, eper questo motivo si usa ricavare i punti con degli algoritmi greedy (letteralmente“golosi”)che lavorano su delle discretizzazioni del compatto K, e che ricavano, quindi, dei puntiapprossimati. E’ chiaro che queste discretizzazioni non possono essere scelte arbitraria-mente, ma dovranno essere su"cientemente grandi e dense tali da approssimare i puntidi Fekete, o quantomeno - e nello specifico sara questo il nostro caso - da avere la lorostessa distribuzione asintotica al crescere del grado di interpolazione (si veda [1]) . D’altraparte, per ragioni di tipo computazionale, richiederemo che questi insiemi siano di cardi-nalita minima indispensabile per raggiungere un buon livello di approssimazione. Un tipodi discretizzazioni che raggiungono lo scopo sono le mesh debolmente ammissibili (WeaklyAdmissibile Meshes), di cui ricordiamo la definizione di seguito:

Definizione 2. Dato un compatto K * C, una mesh debolmente ammissibile (di seguitoabbreviata con WAM), e una sequenza di sottoinsiemi discreti Mn * K, tale che valga laseguente disuguaglianza polinomiale

%p%K & Cn%p%Mn , ,p ! Pn(K) (3.10)

dove Cn > 0 e card(Mn) cresce al piu in maniera polinomiale in n.Una mesh si dice ammissibile (di seguito abbreviata con AM) se e una WAM in cui lasuccessione Cn e limitata, ossia Cn & C, ,n ! N.

Si possono dimostrare una serie di proprieta interessanti per le (W)AM, di cui ricor-diamo le seguenti quattro:

1. per mappe a"ni di del dominio, limmagine di una (W)AM e ancora una (W)AM,con la stessa costante C;

2. ogni sequenza di insiemi contenenti Mn, con cardinalita a crescita polinomiale, e una(W)AM con la stessa costante Cn;

3. ogni sequenza di insiemi unisolventi di interpolazione la cui costante di Lebesguecresce al piu in modo polinomiale con il grado n e una WAM, con Cn la costante diLebesgue stessa;

4. un’unione finita di(W)AM e una(W)AM per la corrispondente unione di compatti,con costante Cn il massimo delle corrispondenti costanti.

Prima di descrivere l’algoritmo di estrazione dei punti di Leja, vediamo come costruiredelle WAM di cardinalita minima da utilizzare come input. Il punto di partenza, dato uncompatto K ! C e la disuguaglianza di Markov, che e soddisfatta se esistono delle costanti' e µ , entrambe positive, tali che

%p(%K & µ(deg(p))!%p%K , ,p ! Pn(K) (3.11)

L’esistenza di tali WAM e conseguenza del seguente risultato di Pommerenke:


Teorema 5. (Pommerenke, [8]) Ogni K ! C compatto connesso, che non sia un singolopunto, soddisfa la disuguaglianza di Markov per µ = e

2cap(K) , e ' = 2, dove cap(K) e la

capacita logaritmica1 dell’insieme K.

Ricordando che il diametro di un insieme e la massima distanza tra tutte le coppie dipunti appartenenti a quell’insieme, e dato che per un insieme compatto K vale la relazionecap(K) - diam(K)

4 , (cf. [10]), e su"ciente utilizzare un limite inferiore per il diametro diK, e applicare la seguente

Proposizione 1. (Piazzon, Vianello, [7]) Sia C > 1 e sia K ! C compatto connesso chesoddisfa la disuguaglianza di Markov come nel Teorema 1. Sia la frontiera di K data dauna curva C1 descritta da ( : [a, b] + )K, con la condizione che maxt%[a,b] |(((t)| & # .Allora gli insiemi

Mn =

%(

#a+

k(b" a)

M " 1

$: k = 0, . . . ,M " 1, n ! N'

&, (3.12)

dove

M =

'C(b" a)#en2

(C " 1)diam(K)

(, (3.13)

formano una mesh ammissibile per K.

Da questo risultato discende un corollario che utilizzeremo per i nostri esperimentinumerici. In buona sostanza non facciamo altro che applicare la Proposizione 1 al casodell’ellisse, che come noto puo essere scritta come curva parametrica tramite la relazione

((t) = a sin(t) + ib cos(t), a, b ! R, t ! [0, 2"] (3.14)

dove a e b sono le lunghezze in modulo dei suoi semiassi.

Corollario 1. Data un’ellisse parametrizzata come in (3.11) essa soddisfa alle ipotesi dellaProposizione 1 e pertanto la successione di insiemi Mn, definita come sopra, forma unamesh ammissibile di costante C, e con diam(K) = max{a, b}, e # = 2 ·max{a, b}.

Dimostrazione. La curva ( definita in (3.11) che descrive l’ellisse, non solo e di classe C1

ma addirittura di classe C# in quanto somma di funzioni di classe C#. Inoltre

|(((t)| = |a sin(t) + ib cos(t)| & |a sin(t)|+ |i||b cos(t)| & |a|+ |b| & 2max{a, b} =: # (3.15)

Possiamo quindi applicare la Proposizione 1, ricordando che diam(K) altri non e che l’assemaggiore, e quindi risulta

M =

'4"Cen2

C " 1

((3.16)

1Per la definizione della capacita logaritmica di un insieme, si veda [5].

3.2. ALGORITMO DLP 17

Osservazione 2. Per ridurre ulteriormente la cardinalita della mesh senza comprometter-ne la funzionalita, nei codici MATLAB useremo la seguente parametrizzazione equivalentedell’ellisse: ((t) = a

2 sin(t) + i b2 cos(t). Questo ci permettera di ridurre # alla quantitamax{a, b} e quindi di ridurre ulteriormente la cardinalita M della mesh.

Osservazione 3. Grazie alla Proposizione 1 abbiamo costruito una mesh che cresce inmaniera quadratica rispetto al grado massimo del polinomio. Sperimentalmente si verificache e possibile creare delle mesh di cardinalita inferiore, ma che discretizzano comunquebene il compatto K. Nei nostri programmi faremo un confronto fra una mesh di cardinalitache cresce con n2, e una che invece cresce con n log n. Vi e inoltre evidenza sperimentaledel fatto che e possibile costruire delle buone mesh di cardinalita 2n e 3n.Per il caso n2, invece, un’ulteriore riduzione dei punti della mesh che utilizzeremo nei codicisara inoltre motivata dal fatto che, nel caso specifico dell’ellisse, vale 2cap(K) > diam(K),e quindi si porra

M =

'4"Cnmax{a, b}en2

(Cn " 1)(a+ b)

((3.17)

3.2 Algoritmo DLP

Siamo finalmente pronti per definire l’algoritmo di estrazione dei punti discreti di Leja, chedi seguito chiameremo algoritmo DLP (Discrete Leja Points). Come detto in precedenza,tale algoritmo prende in input la nostra WAM, crea la matrice rettangolare di VandermondeA = V (a;P ) nella base scelta P , e di seguito estrae i punti di Leja in un vettore *.Ricordiamo che il nostro scopo e massimizzare il determinante di tale matrice; per farlol’algoritmo DLP e!ettua la decomposizione LU della matrice, ossia decompone la matrice Anel prodotto di una matrice triangolare inferiore L con diagonale unitaria e di una matricetriangolare superiore U . Questa decomposizione avviene tramite l’eliminazione di Gauss,che determina i pivot nel seguente modo:

• Il primo pivot e quell’elemento Ar11 tale che |Ar11| sia massimo con r1 = 1, . . . , n

• Ogni riga j #= r1 viene rimpiazzata da se stessa meno Aj1/Ar11 volte la r1-esima riga.In questo modo otteniamo una matrice equivalente alla prima con il pivot in posizione{r1, 1}, e le entrate aj1 = 0 per j #= r1. Infine con una semplice permutazione, siscambia la prima riga con quella contenente il pivot.

• Si applica quindi lo stesso procedimento alla sottomatrice A1 = A2$i,2$j , e allesottmatrici successive fino ad esaurire la dimensione di A.

Si e dunque soliti scrivere la decomposizione tramite la formula PA = LU , dove P ela matrice di permutazione utilizzata per e!ettuare gli scambi, mentre U contiene nelladiagonale i pivot precedentemente trovati. Passando ai determinanti, e ricordando che


le matrici di permutazione hanno determinante ±1, mentre le matrici triangolari hannodeterminante pari al prodotto degli elementi in diagonale, otteniamo che

| det(PA)| = | det(A)| = | det(LU)| = | det(U)|, (3.18)

cioe il prodotto dei pivot. Cerchiamo di capire cosa abbiamo ottenuto: ad ogni passodell’eliminazione di Gauss determiniamo il pivot scegliendo l’elemento di modulo massimonella prima colonna della (sotto)matrice considerata. Questo pivot viene inserito nelladiagonale della matrice U , ed essendo questa triangolare contribuisce a renderne massimoil determinante. Al termine dell’iterazione avremo dunque massimizzato il determinantedella matrice A.Scriviamo quindi ordinatamente l’algoritmo DLP, usando uno pseudo-codice in MATLAB:

Algoritmo DLP

• a = An = {a1, . . . , aM};

• A = V (a;P );

• [L,U,!] = LU(A,’vector’); % ! vettore di permutazione

• ind = !(1 : N); * = a(ind) % Estrazione dei punti

Osservazione 4. Si presti attenzione alla possibilita che ci siano dei problemi di malcondi-zionamento della matrice di Vandermonde. Cio e dovuto all’utilizzo della base monomiale,che rischia di far terminare l’algoritmo per singolarita della matrice. Per risolvere il pro-blema, normalmente, e su"ciente eseguire un paio di fattorizzazioni QR sulla matrice A,che sostituiscono la base monomiale con una base ortogonale. Questo sistema, tuttavia,funziona bene solo fino a gradi di interpolazione non troppo elevati, e in ogni caso l’algorit-mo perde parte della sua e"cienza. Pertanto si e soliti (lo vedremo anche nei nostri codici)cambiare la base monomiale standard (la indichiamo ad esempio con {1, z, . . . , zn}) con unabase monomiale shiftata nel baricentro b della mesh e scalata del valore S = max |zi " b|(gli zi sono i nodi della mesh), ossia

%1,

z " b

S, . . . ,

#z " b

S

$n&(3.19)

Lo shift della base, in particolare, ha una valenza fondamentale nel calcolo dei punti e nelgarantire la loro bonta. Infatti, se il dominio non e centrato nell’origine, mentre inveceparte della sua frontiera le si avvicina, si rischia di compromettere la simmetria nel risul-tato dell’estrazione e quindi la bonta dei punti di interpolazione. Questo avviene perchel’algoritmo estrae un maggior numero di punti lontano dall’origine.

3.2. ALGORITMO DLP 19

Osservazione 5. Un’altra importante osservazione da fare riguarda la stretta dipendenzadei punti di Leja estratti dall’algoritmo con la base P da noi scelta. Questo o!re il vantaggioche, se alla base P aggiungiamo un ulteriore elemento, i punti gia ottenuti con la base dipartenza rimangono comunque validi.

Capitolo 4

Descrizione del metodo e risultati

4.1 Il metodo

Nei capitolo precedente abbiamo visto come si calcolano i DLP. Vediamo ora in che modopossono essere utilizzati nell’ambito dell’approssimazione delle funzioni di matrici. Datauna funzione di matrice f e una matrice A di spettro !(A), tale per cui la definizione dif(A) abbia senso, il nostro metodo di approssimazione polinomiale si articola nel modoseguente:

1. Determinazione del compatto K contenente lo spettro di A. Tale passaggio e ne-cessario dalla definizione di funzione di matrice. Quello che faremo, quindi, datauna generica matrice quadrata A, sara approssimare gli autovalori e costruire la piupiccola ellisse che li contiene.

2. A partire dal bordo di questa ellisse costruiremo una WAM da utilizzare per ricavarcii nodi di interpolazione.

3. Estraiamo dalla WAM i punti di Leja con l’algoritmo DLP.

4. Con il metodo di Newton determiniamo il polinomio interpolante, e il valore dellafunzione con una ricorrenza a tre termini.

Per ognuno dei passi del metodo abbiamo prodotto delle routine in MATLAB, il cuicodice e riportato in appendice, con le quali abbiamo e!ettuato i test e ricavato i risultatinumerici, che descriveremo di seguito.

4.1.1 Costruzione dell’ellisse

In questo passo, l’unica cosa cui bisogna stare attenti e la scelta della matrice A, datoche le funzioni radice quadrata, esponenziale e logaritmo richiedono delle caratteristicheparticolari sullo spettro. In particolare:

21

22 CAPITOLO 4. DESCRIZIONE DEL METODO E RISULTATI

• per la radice quadrata e il logaritmo di matrice sceglieremo la matrice A con lospettro !(A) contenuto tutto nel semipiano complesso destro, ossia tutti con tutti gliautovalori con parte reale positiva.

• per l’esponenziale di matrice sceglieremo la matrice A con %(A) < 1. Se fosse infatti%(A) - 1 l’esponenziale tenderebbe, in modulo, ad esplodere, e quindi con il metodopolinomiale non otterremmo una buona approssimazione.

Una volta ricavato lo spettro, si considerano gli autovalori con massima e minimaparte reale, e quello con massima parte immaginaria, e si determinano cosı assi e centrodell’ellisse.Costruiamo poi l’ellisse con una semplice function (vedi Codice 1) che prende in input isemiassi e le coordinate del centro, e utilizza la parametrizzazione definita in (3.11).

4.1.2 Costruzione della mesh

Il programma afek_ellipse.m (vedi Codice 2) richiede quali input i parametri dell’ellisse(semiassi, centro) e il grado di interpolazione e costruisce la mesh da cui poi verrannoestratti i punti di Leja. Come anticipato nel precedente capitolo, sono presenti due metodiper il calcolo delle mesh: il primo (selezionabile ponendo la variabile extraction_type=1),utilizza la formula ottenuta nel Corollario 1, e quindi genera una mesh che ha crescitaquadratica con n; il secondo invece crea una mesh di cardinalia che cresce con n log n. Lecostanti C sono state ricavate sperimentalmente.

Il seguente grafico mostra come cresce la costante di Lebesgue in funzione del grado diinterpolazione rispetto alla crescita di tipo lineare, sia nel caso del metodo n2, sia nel casodel metodo n log n.

(a) Metodo n2 (b) Metodo n log n

4.1. IL METODO 23

Notiamo che nel caso (a) la precisione e sicuramente maggiore, ma nei nostri test cisiamo dovuti arrestare al grado di interpolazione 10, perche i tempi di esecuzione eranotroppo elevati. Nel caso (b) invece l’andamento della costante di Lebesgue e oscillante,ma comunque in media e sub-lineare. Inoltre i tempi di esecuzione sono notevolmente piulenti. Quindi il metodo n log n assicura non solo una crescita moderata della costante diLebesgue, ma anche un notevole risparmio in termini di operazioni macchina.

4.1.3 Algoritmo DLP: estrazione dei punti di Leja

L’estrazione dei punti viene eseguita tramite la routine approx_fekete.m. (vedi Codice 3)Si noti che con tale programma e possibile estrarre anche i punti approssimati di Fekete,che non sono oggetto di questo lavoro, ma che ci tornano utili per confrontare i risultatinumerici. I punti approssimati di Fekete sono dei punti che approssimano i punti di Feketenel continuno, e, come i punti di Leja, sono considerati dei buoni punti di interpolazione.Come si puo vedere, il programma prende la mesh dell’ellisse che abbiamo creato con laroutine precedente, costruisce la matrice di Vandermonde gia nella base riscalata e shif-tata (si veda il Capitolo 3, Oss. 4), e!ettua due di fattorizzazioni QR, e poi in base allaextraction_type scelta applica l’algoritmo DLP (oppure l’algoritmo AFP per gli Approxi-mate Fekete Points) e ci restituisce i punti che utilizzeremo per l’interpolazione.

Vediamo ora i risultati. Abbiamo e!ettuato il test su una mesh di grado 30 per un’el-lisse di assi a = 0.1 e b = 1, centrata nell’origine. Una volta ottenuti i punti sia conl’algoritmo AFP sia con l’algoritmo DLP, testiamo i nodi di interpolazione su una funzio-ne. Ci limiteremo per ora al caso scalare; per i risultati nel caso matriciale rimandiamo latrattazione alla sezione successiva. La funzione da interpolare che utilizzeremo come teste una funzione esponenziale definita da y(z) = e5z.

(a) Costanti di Lebesgue, rispetto alla crescita lineare:gli asterischi rappresentano gli AFP, i triangoli i DLP

(b) Errore di interpolazione in scala logaritmica: gliasterischi rappresentano gli AFP, i triangoli i DLP


Di nuovo, notiamo come la crescita della costante di Lebesgue sia molto piu lenta negliAFP che non nei DLP, pur mantenendosi questa, in media, al di sotto della crescita ditipo lineare (si confronti con la linea rossa). Osserviamo pero che, in ogni caso, l’errore diinterpolazione dei DLP a gradi alti e per la nostra funzione test quasi indentico.

Facendo riferimento a quanto visto nel paragrafo 3.2, in cui spiegavamo il funzionamentodell’algoritmo di estrazione, mostriamo anche che se fissiamo un grado (sia questo m) perl’interpolazione, i punti DLP estratti per il grado m sono tali che un suo sottoinsieme deiprimi k elementi e buono per l’interpolazione. Vediamo di seguito un esempio con m = 30,in cui calcoliamo le costanti di Lebesgue dei set di punti {x1, . . . , xk} per k = 1, . . . m. (Siveda a tal proposito il Codice 4).

N log(N) METHOD.>> DLP.DEG: 0 PTS: 1 LEB: 5.00000e+04 WAM CARD.: 511DEG: 1 PTS: 2 LEB: 2.99765e+00 WAM CARD.: 511DEG: 2 PTS: 3 LEB: 1.28561e+00 WAM CARD.: 511DEG: 3 PTS: 4 LEB: 2.98982e+00 WAM CARD.: 511DEG: 4 PTS: 5 LEB: 2.35168e+00 WAM CARD.: 511DEG: 5 PTS: 6 LEB: 3.96907e+00 WAM CARD.: 511DEG: 6 PTS: 7 LEB: 2.75498e+00 WAM CARD.: 511DEG: 7 PTS: 8 LEB: 4.82396e+00 WAM CARD.: 511DEG: 8 PTS: 9 LEB: 4.81389e+00 WAM CARD.: 511DEG: 9 PTS: 10 LEB: 4.66208e+00 WAM CARD.: 511DEG: 10 PTS: 11 LEB: 4.71705e+00 WAM CARD.: 511DEG: 11 PTS: 12 LEB: 9.52893e+00 WAM CARD.: 511DEG: 12 PTS: 13 LEB: 7.10738e+00 WAM CARD.: 511DEG: 13 PTS: 14 LEB: 6.81177e+00 WAM CARD.: 511DEG: 14 PTS: 15 LEB: 6.79881e+00 WAM CARD.: 511DEG: 15 PTS: 16 LEB: 6.01624e+00 WAM CARD.: 511DEG: 16 PTS: 17 LEB: 9.85031e+00 WAM CARD.: 511DEG: 17 PTS: 18 LEB: 8.10743e+00 WAM CARD.: 511DEG: 18 PTS: 19 LEB: 7.69111e+00 WAM CARD.: 511DEG: 19 PTS: 20 LEB: 1.29537e+01 WAM CARD.: 511DEG: 20 PTS: 21 LEB: 5.84122e+00 WAM CARD.: 511DEG: 21 PTS: 22 LEB: 1.03854e+01 WAM CARD.: 511DEG: 22 PTS: 23 LEB: 9.12034e+00 WAM CARD.: 511DEG: 23 PTS: 24 LEB: 1.06420e+01 WAM CARD.: 511DEG: 24 PTS: 25 LEB: 9.07345e+00 WAM CARD.: 511DEG: 25 PTS: 26 LEB: 1.20676e+01 WAM CARD.: 511DEG: 26 PTS: 27 LEB: 1.01271e+01 WAM CARD.: 511DEG: 27 PTS: 28 LEB: 7.88695e+00 WAM CARD.: 511DEG: 28 PTS: 29 LEB: 3.20247e+01 WAM CARD.: 511DEG: 29 PTS: 30 LEB: 1.62282e+01 WAM CARD.: 511DEG: 30 PTS: 31 LEB: 1.81662e+01 WAM CARD.: 511

4.1. IL METODO 25

L’andamento delle costanti ci dice che, anche in questo caso, i punti di interpolazionesono buoni.

Osservazione 6. Si presti attenzione a non confondere l’esito di questo test con quellomostrato nel Paragrafo 4.1.2, in cui per m = 1, . . . , 30 viene calcolata la costante di Lebe-sgue dei DLP di grado m. In quest’ultimo test, invece, la cardinalia della WAM e fissata,mentre varia il grado della costante di Lebesgue.

4.1.4 Interpolazione delle funzioni di matrici con il metodo di Newton

Ricordiamo che il nostro scopo e l’ approssimazione delle funzioni test con dei polinomi.Calcoleremo cioe

Ym = pm"1(A) .= f(A) (4.1)

dove pm"1 ! Pm"1(K). L’interpolazione avverra tramite il metodo di Newton, che ricor-diamo essere definito nella seguente maniera:

pm"1(z) = f0 + f1(z " *1) + · · ·+ fm"1(z " *1) . . . (z " *m"1) (4.2)

dove gli {*i}i=1,...,m sono punti distinti del piano complesso (ossia i nostri nodi di in-terpolazione: i punti di Leja), e dove le fk sono le note di!erenze divise del metodo diNewton:

f0 = f [*1], fk = f [*1, . . . , *k+1], k = 1, . . . ,m" 1 (4.3)

Con queste notazioni, si mostra che la (4.2) soddisfa la seguente ricorrenza a tre termini:

)***+

***,

Y0 = O(m+1)!(m+1)

Y1 = f0I

Ym = Ym"1 +fm"1

fm"2(A" *m"1I)(Ym"1 " Ym"2), m - 2

(4.4)


che sara quella che useremo nei nostri codici. La routine che implementa la ricorrenzae riportata nel Codice 4. Ai risultati numerici a"ancheremo una stima dell’errore, os-sia confronteremo il valore di Ym con il valore ottenuto usando le funzioni gia presenti inMATLAB per il calcolo di esponenziale, logaritmo e matrice quadrata di matrice, ossia lefunzioni expm, logm, sqrtm. Questi programmi sfruttano le tecniche di calcolo spiegate nelcapitolo 2.

4.2 Risultati: un problema di convezione-di!usione

Produrremo le matrici test a partire da un classico problema fisico di convezione-di!usione.Tale matrice si ricava dalla discretizzazione dell’operatore

L(U) = "U (( + c U ( (4.5)

dove U ! C([0, 1]), e c e una costante. In forma matriciale otteniamo una matrice tridia-gonale A definita positiva che utilizzeremo per il calcolo di radice quadrata e logaritmo dimatrice. Per l’esponenziale di matrice, invece, utilizzeremo "A, che avra tutti gli autova-lori negativi (questo per evitare che la matrice diventi troppo grande e dia quindi probleminell’approssimazione polinomiale).I test sono stati e!ettuati su matrici di dimensione 50 / 50, mentre il grado m dell’inter-polazione si ferma a 45÷ 47 in quanto a gradi superiori si verificano problemi di malcondi-zionamento che, probabilmente, potrebbero essere risolti con una scelta migliore della basepolinomiale della matrice di Vandermonde. Questo ci permetterebbe di raggiungere unamigliore convergenza del metodo che, soprattutto per la radice quadrata e il logaritmo,osserviamo essere un po’ lenta.Riportiamo di seguito i risultati con c = 0 e c = 4.

4.2. RISULTATI: UN PROBLEMA DI CONVEZIONE-DIFFUSIONE 27

(a) Esponenziale di matrice, c=0 (b) Esponenziale di matrice, c=4

(c) Radice quadrata di matrice, c=0 (d) Radice quadrata di matrice, c=4

(e) Logaritmo di matrice, c=0 (f) Logaritmo di matrice, c=4

Appendice A

Codici MATLAB

Alleghiamo di seguito i codici scritti in linguaggio MATLAB utilizzati dal Metodo polino-miale descritto nel Capitolo 4.

A.1 Codice 1: Costruzione dell’ellisse

%-------------------------------------------------------------------------% Input: matrice quadrata A% Output: autovalori di A, semiassi e centro dell ’ellisse% che contiene lo spettro di A%-------------------------------------------------------------------------

spect=eig(A);

re_min=min(real(spect ));re_max=max(real(spect ));im_max=max(imag(spect ));im_min=min(imag(spect ));

a=abs((re_max -re_min )/2)*2;b=abs((im_max -im_min )/2)*2;x0=( re_max+re_min )/2y0=0;

%--------------------------------------------------------------------------% Ellipse%--------------------------------------------------------------------------

function z=ellipse(theta ,a,b,x0 ,y0)z=(a/2)* cos(theta )+(b/2)*i*sin(theta )+x0+i*y0;

29

30 APPENDICE A. CODICI MATLAB

A.2 Codice 2: Costruzione della mesh

function [fek ,mesh_ellipse ,C,ind]= afek_ellipse(a,b,x0 ,y0 ,deg ,extraction_type)

wam_type =2;

switch wam_typecase 1

fprintf(’\n \t N^2 METHOD .’);mesh_ellipse=WAM_Ellipse(a,b,x0,y0,deg);

case 2fprintf(’\n \t N log(N) METHOD .’);C=5;mesh_ellipse=WAM_Ellipse_nlogn(a,b,x0,y0,deg ,C);

end

fprintf(’\n \t WAM DIMENSION: %5.0f’,length(mesh_ellipse ));[fek ,C,ind]= approx_fekete(mesh_ellipse ,deg ,extraction_type );

%--------------------------------------------------------------------------% ADDITIONAL ROUTINES.%--------------------------------------------------------------------------

%--------------------------------------------------------------------------% WAM_Ellipse%--------------------------------------------------------------------------

function Mn=WAM_Ellipse(a,b,x0 ,y0 ,n,C)

% a,b: ELLIPSE AXIS.% x0 ,y0: ELLIPSE CENTER.% n: DEGREE.% C: WAM FACTOR.

% Mn: WAM ON THE ELLIPSE DEFINED BY a,b,x0,y0,n.

if nargin < 6C=200;

end

ell_cap =(a+b)/2;mu=exp (1)/ ell_cap;

A.3. CODICE 3: ESTRAZIONE DEI PUNTI DI LEJA 31

bv=2*pi; av=0;alpha=max(a,b);

NN=(n^2)*(bv -av)* alpha*mu/2;

M=ceil(NN*C/(C -1))+1;

k=0:M-1;v=k*2*pi/(M-1);Mn=ellipse(v,a,b,x0 ,y0);

%--------------------------------------------------------------------------% WAM_Ellipse%--------------------------------------------------------------------------function Mn=WAM_Ellipse_nlogn(a,b,x0 ,y0 ,n,C)

% a,b: ELLIPSE AXIS.% x0 ,y0: ELLIPSE CENTER.% n: DEGREE.% C: WAM FACTOR.% Mn: WAM ON THE ELLIPSE DEFINED BY a,b,x0 ,y0,n.

if nargin < 6C=200;

end

M=ceil(C*n*log(n));

k=0:M-1;v=k*2*pi/(M-1);Mn=ellipse(v,a,b,x0 ,y0);

A.3 Codice 3: Estrazione dei Punti di Leja

function [fek ,C,ind] = approx_fekete(mesh ,deg ,extraction_type)

%numero di iterazioni QR per ridefinire la baseniter =2;

%shift e riscalamento della base monomiale

b=sum(mesh)/ length(mesh); %baricentro della mesh S=max(abs(mesh -b));S=max(abs(mesh -b)); %coefficiente della riscalaturaz=(mesh -b)/S; %mesh shiftata e scalata


%costruzione matrice di Vandermonde

for j=1:deg+1V(j,:)=z.^(j-1); %riga j-esima della matrice di Vandermonde

end;

% fattorizzazioni QR

B=V’;P=eye(length(V(: ,1)));

% matrice unitaria , memorizza i cambiamenti di base

for it=1:niter ,[Q,R]=qr(B,0);U=inv(R); % matrice di cambiamento di baseB=B*U; % matrice di Vandermonde nella nuova baseP=P*U; % memorizzazione cambiamento di base

end;

%soluzione sistema lineare

C=B’; m=ones(length(C(: ,1)) ,1);

switch extraction_type

case 1

fprintf(’\n \t >> AFP.’)

%vettore dei termini notiw=C\m;

%indice dei pesi non nulliind=find(abs(w)>0);

case 2

% algoritmo LUfprintf(’\n \t >> DLP.’)

[L,U,Plu] = lu(B,’vector ’);ind=Plu(1: size(U,1));

A.4. CODICE 4: INTERPOLAZIONE CON IL METODO DI NEWTON 33

end

%punti di Fekete approssimati corrispondenti

fek=mesh(ind);

A.4 Codice 4: Interpolazione con il metodo di Newton

A titolo di esempio, vediamo la routine m_leja.m nel caso dell’esponenziale.

function Y = m_leja(A)

deg =30;

% Estrazione punti di fekete / leja

[fek ,mesh_ellipse ,C,ind]= afek_ellipse(a,b,x0,y0,deg ,extraction_type );figure (3)plot(real(fek),imag(fek),’r*’)

% Interpolazione: calcolo delle differenze diviseXX=fek;YY=feval(’exp ’,XX);

[d, Df , c, f] = newton_interpolation_vett(XX,YY,XX);

% ricorrenza a tre termini

m=length(XX);Y(:,:,1)= zeros(size(A)); Y(:,:,2)=f(1)* eye(size(A));for i=3:m

Y(:,:,i)=Y(:,:,i-1)+(d(i-1)/d(i -2))*(A-XX(i-2)*eye(size(A)))*(Y(:,:,i-1)-Y(:,:,i-2));

end

eA = expm(A);for i=1:mE(:,:,i)=Y(:,:,i)-eA;end

for i=1:mnorma(i)=norm(E(:,:,i));

endplot(log10(norma))

Bibliografia

[1] L. Bos, J.-P. Calvi, N. Levenberg, A. Sommariva, M. Vianello, Geometric WeaklyAdmissible Meshes, Discrete Least Squares Approximations and Approximate FeketePoints, preprint, 2009.

[2] L. Bos, S. De Marchi, A.Sommariva, M.Vianello Computing multivariate Fekete andLeja points by numerical linear algebra, 2010, preprint.

[3] L. Bos, S. De Marchi, A. Sommariva, M. Vianello, On Multivariate NewtonInterpolation at Discrete Leja Points, 2011, to appear.

[4] L. Bos, S. De Marchi, A. Sommariva, M. Vianello Weakly Admissible Meshes andDiscrete Extremal Sets, 2011, preprint.

[5] W. Dijkstra, M.E. Hochstenbach, Numerical approximation of the logarithmic capacity,preprint

[6] Nicolas J. Higham, Function of Matrices: Theory and Computation, Siam, 2008.

[7] F. Piazzon, M. Vianello Analytic transformations of admissible meshes

[8] Ch. Pommerenke, On the Derivate of a Polynomial, Michigan Math. J. 6 (373-375),1959.

[9] A. Sommariva, Apprssimazione di funzioni, Dispense del corso di Analisi Numerica(last rev. 23/10/2010)

[10] M. Tsuji, Potential Theory in Modern Function Theory, Chelsea Publishing Co , 1975.

[11] Wikipedia, (Lebesgue Constant interpolation) http://en.wikipedia.org/wiki/

Lebesgue_constant_(interpolation)

35

UN METODO POLINOMIALE PER IL CALCOLO DI FUNZIONI DI ...alvise/TESI_STUDENTI/CACCO.pdf · loro...

Documents

Transcript of UN METODO POLINOMIALE PER IL CALCOLO DI FUNZIONI DI ...alvise/TESI_STUDENTI/CACCO.pdf · loro...