Universita degli Studi di Venezia

Corso di Laurea in Informatica

Appunti per le lezioni diCalcolo Numerico

Prof. Flavio Sartoretto

2 giugno 2009

Indice

1 Introduzione 61.1 Notazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.1.1 Esempi . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2 Costo Computazionale . . . . . . . . . . . . . . . . . . . . . . 71.3 Documentazione dei risultati . . . . . . . . . . . . . . . . . . . 81.4 Concetti di base. . . . . . . . . . . . . . . . . . . . . . . . . . 111.5 Definizioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2 Rappresentazione dei numeri. 132.1 Pentium FDIV flaw . . . . . . . . . . . . . . . . . . . . . . . . 132.2 Esempi di errori di arrotondamento. . . . . . . . . . . . . . . . 13

2.2.1 Errori di accumulazione . . . . . . . . . . . . . . . . . 132.3 Modi di arrotondamento . . . . . . . . . . . . . . . . . . . . . 15

3 Schemi ricorrenti. 16

4 Risoluzione di equazioni non lineari. 184.1 Esercizio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

5 Risoluzione di sistemi lineari. 225.1 Matrici e vettori . . . . . . . . . . . . . . . . . . . . . . . . . . 22

5.1.1 Autovettori e autovalori . . . . . . . . . . . . . . . . . 235.1.2 Norme . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

5.2 Sistemi lineari . . . . . . . . . . . . . . . . . . . . . . . . . . . 275.3 Decomposizione a valori singolari . . . . . . . . . . . . . . . . 27

5.3.1 Un esempio . . . . . . . . . . . . . . . . . . . . . . . . 285.3.2 Spazi vettoriali . . . . . . . . . . . . . . . . . . . . . . 295.3.3 Rango di matrici . . . . . . . . . . . . . . . . . . . . . 295.3.4 Sistemi sottodeterminati . . . . . . . . . . . . . . . . . 305.3.5 Sistemi sovradeterminati . . . . . . . . . . . . . . . . . 315.3.6 Filtraggio del rumore . . . . . . . . . . . . . . . . . . . 32

1

5.3.7 Compressione di immagini . . . . . . . . . . . . . . . . 325.4 Esercizio sull’ eliminazione di Gauss. . . . . . . . . . . . . . . 34

6 Costo Computazionale 376.1 Operazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 376.2 Accessi alla memoria . . . . . . . . . . . . . . . . . . . . . . . 396.3 Esempi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 416.4 Costo del metodo di Jacobi . . . . . . . . . . . . . . . . . . . 50

7 Interpolazione. 527.1 L’ algoritmo di Neville. . . . . . . . . . . . . . . . . . . . . . . 527.2 Operatori lineari. . . . . . . . . . . . . . . . . . . . . . . . . . 537.3 Interpolazione trigonometrica e FFT . . . . . . . . . . . . . . 547.4 Calcolo di splines . . . . . . . . . . . . . . . . . . . . . . . . . 56

8 Approssimazione. 598.1 Approssimazione di funzioni. . . . . . . . . . . . . . . . . . . . 59

9 Integrazione numerica. 649.1 Estrapolazione di Romberg . . . . . . . . . . . . . . . . . . . . 659.2 Integrazione di Gauss . . . . . . . . . . . . . . . . . . . . . . . 66

10 Equazioni differenziali ordinarie. 7010.1 Cenni sulla risolubilita delle ODE. . . . . . . . . . . . . . . . . 7010.2 Errori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7610.3 Sistemi di ODE. . . . . . . . . . . . . . . . . . . . . . . . . . . 7610.4 Runge-Kutta per sistemi di ODE. . . . . . . . . . . . . . . . . 7710.5 Risoluzione di un BVP con shooting. . . . . . . . . . . . . . . 7810.6 BVP e sistemi alle DF . . . . . . . . . . . . . . . . . . . . . . 8010.7 Strato limite . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

10.7.1 Formulazione del problema . . . . . . . . . . . . . . . . 8210.7.2 Risoluzione numerica . . . . . . . . . . . . . . . . . . . 8310.7.3 Conclusioni . . . . . . . . . . . . . . . . . . . . . . . . 86

11 Metodi Multigrid 8711.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8711.2 Nested iterations . . . . . . . . . . . . . . . . . . . . . . . . . 9111.3 Griglie triangolari . . . . . . . . . . . . . . . . . . . . . . . . . 9311.4 Multigrid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9511.5 Convergenza e Costo . . . . . . . . . . . . . . . . . . . . . . . 9711.6 Osservazioni interessanti . . . . . . . . . . . . . . . . . . . . . 10111.7 Multigrid adattivo . . . . . . . . . . . . . . . . . . . . . . . . 103

2

11.7.1 Interfacce fra sottogriglie . . . . . . . . . . . . . . . . . 10411.7.2 Adattamento della griglia . . . . . . . . . . . . . . . . 117

A Progetto di un contenitore 119A.1 Il problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119A.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

B Modelli di crescita 121B.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121B.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

C Flusso di calore 124C.1 Il problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124C.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

D Calcolo di angoli 127D.1 La Cefalometria . . . . . . . . . . . . . . . . . . . . . . . . . . 127D.2 Descrizione del problema . . . . . . . . . . . . . . . . . . . . . 129D.3 Approccio risolutivo . . . . . . . . . . . . . . . . . . . . . . . 129D.4 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

E Schemi di rilassamento 131E.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131E.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

F Valutazione di sistemi 133F.1 Descrizione del problema . . . . . . . . . . . . . . . . . . . . . 133F.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

G Information Retrieval 137G.1 Formulazione del problema . . . . . . . . . . . . . . . . . . . . 137

G.1.1 Latent Semantic Indexing . . . . . . . . . . . . . . . . 138G.1.2 Valutazione delle prestazioni . . . . . . . . . . . . . . . 139

G.2 Obiettivi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139G.3 Parte facoltativa . . . . . . . . . . . . . . . . . . . . . . . . . 139

H Costo computazionale 142H.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142H.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143H.3 Approfondimenti . . . . . . . . . . . . . . . . . . . . . . . . . 143H.4 Memorizzazione in forma sparsa . . . . . . . . . . . . . . . . . 143H.5 Facoltativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

3

I Deformazione di progetto 144I.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144I.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

J Grafica 2D 147J.1 Formulazione del problema . . . . . . . . . . . . . . . . . . . . 147J.2 Obiettivi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148J.3 Rappresentazioni . . . . . . . . . . . . . . . . . . . . . . . . . 149J.4 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149J.5 Parte facoltativa . . . . . . . . . . . . . . . . . . . . . . . . . 150

K Analisi di carico 151K.1 Il problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151K.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151

L Moto di un giocattolo 153L.1 Formulazione del problema . . . . . . . . . . . . . . . . . . . . 153L.2 Procedimento risolutivo. . . . . . . . . . . . . . . . . . . . . . 153L.3 Obiettivi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155L.4 Suggerimenti . . . . . . . . . . . . . . . . . . . . . . . . . . . 155L.5 Esempio di risoluzione . . . . . . . . . . . . . . . . . . . . . . 155

M Onde di traffico 159M.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159M.2 Descrizione e compiti . . . . . . . . . . . . . . . . . . . . . . . 159

N Analisi della dinamica 167N.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167N.2 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167

O Sistemi Hamiltoniani 168O.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168

O.1.1 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . 170O.2 Sistemi Simplettici . . . . . . . . . . . . . . . . . . . . . . . . 171

O.2.1 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . 172

P Test dello Shock Tube 173P.1 Descrizione del problema . . . . . . . . . . . . . . . . . . . . . 173P.2 Sistema fluidodinamico ideale . . . . . . . . . . . . . . . . . . 175P.3 Metodi di integrazione numerica . . . . . . . . . . . . . . . . . 175P.4 Integrazione temporale . . . . . . . . . . . . . . . . . . . . . . 176

P.4.1 Metodo di Lax . . . . . . . . . . . . . . . . . . . . . . 177

4

P.4.2 Il metodo di Lax-Wendroff . . . . . . . . . . . . . . . . 177P.5 Compiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178

- Bibliografia 178

5

Capitolo 1

Introduzione

F. Sartoretto

Calcolo Numerico

Testi adottati:

• A. Quarteroni e F. Saleri. Introduzione al Calcolo Scientifico. SpringerVerlag Italia, III edizione, 2006.

• F. Sartoretto e M. Putti. Introduzione alla Programmazione per Ela-borazioni Numeriche. Edizioni Libreria Progetto, Padova, 2008.

• F. Sartoretto, Appunti di Calcolo Numerico, Dispense, WEB pagedel Docente, 2009.

Testi consigliati:

• G. Gambolati. Lezioni di Metodi Numerici per Ingegneria e ScienzeApplicate. Cortina, Padova, 1994.

Prerequisiti:

• Calcolo I, con cenni sui numeri complessi.

• Calcolo II: risoluzione di equazioni differenziali ordinarie lineari a coef-ficienti costanti, calcolo di derivate parziali.

• Algebra delle matrici.

In questi appunti, si fa riferimento all’area WEB dedicata al corso. Ladirectory principale viene chiamata $CN = www.dsi.unive.it/∼sartoret-/italian/didattica/CalcoloNumerico/.

6

1.1 Notazioni

• a ≃ b significa “a” e approssimativamente uguale a “b”.

• f(x) = O(g(x)) per x → a (si legge f(x) e di ordine “O-grande” dig(x) per x→ a), significa che la funzione

∣∣∣∣f(x)

g(x)

∣∣∣∣

e limitata per x→ a, ossia esiste un intorno Ia di a ed esiste M ∈ ℜ+

t.c. (∀x ∈ Ia) |f(x)/g(x)| ≤M .

1.1.1 Esempi

• sin(x) = O(1) per x→∞.

Infatti |sin(x)/1| ≤ 1.

Notare che non esiste il limx→∞ sin(x)/1. Notare anche che sin(x) =O(xk), k ≥ 0, quando x→∞!

• sin(1/x) = O(1) per x→ 0.

Per le stesse ragioni di prima.

• 2(n+ 1)n = O(n2) per n→ +∞. Infatti limn→+∞ |2(n+ 1)n/n2| = 2.

•

Una relazione piu forte e:

f(x) ≈ g(x)⇐⇒f(x) = O(g(x)) e g(x) = O(f(x)).

Si dice che f(x) e g(x) sono simili.

1.2 Costo Computazionale

Sia A una matrice n× n e v un vettore di n componenti.Il prodotto (algebrico)

w = Av

richiede n prodotti e n−1 somme per calcolare ogni elemento di w. In totalerichiede quindi n(2n− 1) operazioni floating point.

7

• Se consideriamo unitario il costo di una moltiplicazione e un’addizione,il suo costo e O(2n2 − n) = O(n2), infatti

limn→+∞

2n2 − nn2

= 2.

• Se consideriamo unitario il costo di una somma piu una moltiplicazione,il costo e O(n(n− 1)) = O(n2).

• Se consideriamo unitario il costo di un prodotto scalare tra vettori didimensione n, allora il costo e O(n).

1.3 Documentazione dei risultati

Non basta calcolare risultati numerici corretti, bisogna saperli anche presen-tare in maniera corretta.

Valgono le seguenti regole:Usare tabelle, ma con parsimonia: un tabulato con dieci pagine di risultati

numerici e illeggibile! Riportare in tabella solo valori importanti e affidare agrafici la presentazione di altri risultati.

No

x**2

Nei grafici

8

• segnare sempre con un simbolo opportuno i valori calcolati, per sepa-rarli da eventuali linee approssimatrici;

• tratteggiare le linee approssimatrici;

• usare sempre la scala semilogaritmica per disegnare l’andamento diquantita che presentano forti variazioni (errori, scarti, etc.);

• indicare sempre chiaramente le scale in ascissa e ordinata;

• indicare sempre chiaramente il significato delle ascisse e delle ordinatee le eventuali unita di misura usate.

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

No

x**2"esempi.dat"

Quando vengono forniti dei dati e viene chiesto di approssimarli, riportaresempre i primi assieme ai secondi, nello stesso grafico e/o tabulato.

9

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

Def

orm

azio

ne (

cm)

Pressione (Kg)

Approssimazione di una deformazione

Si

ApprossimazioneDeformazioni osservate

Etichette degli assi

Titolo

Legenda

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

Def

orm

azio

ne (

cm)

Pressione (Kg)

Approssimazione di una deformazione

ApprossimazioneDeformazioni osservate

... uso di scale semi-logaritmiche e logaritmiche ...

10

1 2 3 4 5 6 7 8 9 100

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Iterations

Err

or

error

1 2 3 4 5 6 7 8 9 1010

−5

10−4

10−3

10−2

10−1

100

Iterations

Err

or

error

1.4 Concetti di base.

Calcolo Numerico o Scientifico:

Disciplina che tenta di dare una risposta numerica ad un problemamatematico usando un elaboratore digitale.

11

Obiettivo del Calcolo Numerico o Scientifico e

Trovare gli algoritmi che risolvono un problema matematico

• nel minor tempo,

• con il minimo impiego di risorse di calcolo,

• con la massima accuratezza possibile.

1.5 Definizioni

Un procedimento numerico si dice instabile se gli errori associaticrescono fino a distruggere la soluzione.

Un problema si dice mal condizionato se una piccola variazione neidati (puo) provocare una grande variazione dei risultati.

12

Capitolo 2

Rappresentazione dei numeri.

2.1 Pentium FDIV flaw

Eseguendo sui primi Pentium l’ operazione

EXA := (5-1/1000000)/(15-1/1000000);

Si otteneva il valore 0.333329, con un errore relativo: (EXA-0.333329)/EXA=

1.287461e-5.Eseguendo la stessa operazione su un 486 si ottiene 0.33333329. L’ errore

relativo e’:(EXA-33333329/100000000)/EXA =-1666671/499999900000000,ossia circa -3.333343e-9.

2.2 Esempi di errori di arrotondamento.

2.2.1 Errori di accumulazione

La rappresentazione con un numero finito di cifre da luogo ad errori. Eccoi valori ottenuti calcolando le quantita xis = x

(s)i = xi−1 + h, x

(s)0 = x0,

xip = x(p)i = x0 + ih, i = 1, . . . , n. Il valore iniziale esatto e x0 = 1.3.

Risultati ottenuti sul VAX con VAX FORTRAN

i xis xip

0 1.29999995231628418 1.29999995231628418

1 1.34999990463256836 1.34999990463256836

2 1.39999985694885254 1.39999997615814209

13

−10

−5

0

5

10

−10 −5 0 5 10

x0

x1

=x

0+h

x2

=x

1+h

x3

=x

2+h

x4

=x

3+h

x1

=x

0+h

x2

=x

0+

2∗h

x3

=x

0+

3∗h

x4

=x

0+

4∗h

x0

Figura 2.1: Discretizzazione uniforme di un intervallo.

14

3 1.44999980926513672 1.44999992847442627

4 1.49999976158142090 1.50000000000000000

5 1.54999971389770508 1.54999995231628418

6 1.59999966621398926 1.60000002384185791

7 1.64999961853027344 1.64999997615814209

8 1.69999957084655762 1.69999992847442627

9 1.74999952316284180 1.75000000000000000

2.3 Modi di arrotondamento

Consideriamo un numero floating point normalizzato, in un’aritmetica di baseβ, con mantissa o significante (significand) a p cifre, esponenti s, emin ≤ s ≤emax,

x = ±a1.a2 . . . ap × βq = m× βq

Ricordiamo che troncare x a r ≤ p cifre da luogo al valore

trr(x) = ±a1.a2 . . . ar × βq.

Due modi tipici di calcolare l’ arrotondmento a r cifre di x sono i seguenti.

1. Arrotondamento (senza attributo) [Com91]

rdr(x) = sign(m) tr(|m|+ 1

2β−r)× βq.

2. Arrotondamento pari (round to even) [Gol91]

rd(e)r (x) =

rdr(x) se ar+1 6= β/2,trr(x) se ar pari,rdr(x) altrimenti.

altrimenti.

I due modi differiscono quando vi e equilibrio, ossia se β = 10 quando la(r + 1)-esima cifra e 5. In tal caso, rdp(x) e il piu vicino floating point che

in valore assoluto e piu grande di |x|, mentre rd(e)p (x) dipende dalla cifra ar.

Chiamiamo vicino minore, xn, di x il piu grande floating point minore di x;chiamiamo vicino maggiore, xg, di x il piu piccolo floating point maggiore dix. L’arrotondamento pari restituisce quello tra xn e xg che ha la cifra menosignificativa pari.

Esempio: sia β = 10, p = 3. Se x = 1.25, rd2(x) = 1.3, rd(e)2 = 1.2,

mentre se x = 1.35, rd2(x) = 1.4, rd(e)2 = 1.4. Se ar e diverso da 5, i due modi

coincidono: x = 1.24, rd2(x) = 1.2 = rd(e)2 , x = 1.26, rd2(x) = 1.3 = rd

(e)2 ,

ecc.

15

Capitolo 3

Schemi ricorrenti.

Esercizio 3.0.1 Per calcolare i valori dei polinomi ortogonali T ∗n(x), n =

2, . . . , 20, per x = 2 e x = 1/2, si vuole utilizzare lo schema ricorrente

T ∗n+1(x) = (−2 + 4x)T ∗

n(x)− T ∗n−1(x)

T ∗0 (x) = 1, T ∗

1 (x) = x

Supponiamo che T ∗0 e T ∗

1 siano noti con errori ǫ0 e ǫ1 che siano uguali inmodulo.

Nel caso x = 2 lo schema si puo ritenere stabile?Nel caso x = 1/2 lo schema si pio ritenere stabile?

Soluzione.Le seguenti tabelle riassumono i risultati che si ottengono calcolando

alcuni termini della formula ricorrente. Caso x = 2.

ǫ0 = ǫ1 ǫ0 = −ǫ1ǫ0 −ǫ05 ǫ0 −7 ǫ029 ǫ0 −41 ǫ0169 ǫ0 −239 ǫ0

Caso x = 1/2.

ǫ0 = ǫ1 ǫ0 = −ǫ1ǫ0 −ǫ0−ǫ0 −ǫ0−ǫ0 ǫ0ǫ0 ǫ0ǫ0 −ǫ0−ǫ0 −ǫ0−ǫ0 ǫ0ǫ0 ǫ0

16

Quindi nel caso x = 2 lo schema non e stabile, lo e nel caso x = 1/2.

17

Capitolo 4

Risoluzione di equazioni nonlineari.

4.1 Esercizio

Calcolare con uno scarto relativo minore o uguale a T = 10−3 il valore di ξt.c.

ln(ξ) = cos(ξ) (4.1)

Risolvere questo problema

1. determinando un intervallo in cui cade la soluzione.

2. Eseguendo una iterazione con il metodo della bisezione, in modo daottenere una stima x0 della soluzione,

3. eseguendo un numero opportuno di iterazioni con lo schema di Newton-Raphson, dopo aver posto x0 = x0.

Arrotondare i risultati intermedi a 5 cifre decimali.Soluzione: schizzando le funzioni y = ln(x) e z = cos(x) si vede facilmente

(cf figura) che esiste una sola soluzione dell’ equazione data, che cade adesempio nell’intervallo I = [1, π/2].

Trovare il valore ξ che soddisfa l’ equazione equivale a trovare lo zero dellafunzione

f(x) = ln(x)− cos(x).

Applicando il metodo di bisezione all’ intervallo I, posto ri = stima dellasoluzione all’ i-esima iterazione, ei = |(ti − si)/(ti + si)|, otteniamo:

i s_i t_i r_i e_i

0 1.0000e+00 1.5708e+00 1.2854e+00 2.2203e-01

1 1.2854e+00 1.5708e+00 1.4281e+00 9.9923e-02

18

0

0.2

0.4

0.6

0.8

1

1 1.1 1.2 1.3 1.4 1.5

y

x

cos(x)log(x)

Figura 4.1: Grafico delle due funzioni.

Posto ora x0 = 1.4281, ei = |(xi − xi−1)/xi| e usando lo schema diNewton-Raphson, otteniamo:

i x_i e_i

0 1.4281e+00 NaN

1 1.3014e+00 9.7357e-02

2 1.3030e+00 1.2279e-03

3 1.3030e+00 0

Poiche lo scarto relativo e decrescente e alla terza iterazione e minore dellatolleranza richiesta, alla terza iterazione si ha l’ approssimazione richiesta,che, arrotondata alla quinta cifra decimale e x = 1.3030.

Esercizio 4.1.1 Consideriamo lo schema del punto fisso

xn+1 = ln(1 + xn) + 1, n = 0, 1, . . . (4.2)

x0 = 2 (4.3)

1. Provare che questo schema converge ad un valore ξ.

2. Usando la stima a priori sull’ errore, valutare il numero di iterazioni knecessario perche xk sia affetto da un errore relativo ǫk ≤ 0.003.

3. Usando la stima a posteriori sull’ errore, eseguire il numero n di itera-zioni necessario perche l’ errore relativo sia ǫn ≤ 0.003.

19

-3

-2

-1

0

1

2

-1 -0.5 0 0.5 1 1.5 2 2.5

y

x

log(1.0+x)x-1.0

Figura 4.2: Grafico delle funzioni y = x− 1 e y = ln(1 + x).

Arrotondare i risultati intermedi a 5 cifre decimali.

Soluzione: Se lo schema converge, deve convergere alle soluzioni di

x = g(x) = ln(1 + x) + 1

che si vede (cf figura) sono una negativa ξ−, e una positiva ξ.Studiando la funzione, si puo vedere che e 2 < ξ < 2.5. Essendo g′(x) =

1/(1 + x), posto ad esempio I = [2, 3] risulta

• g(I) ⊂ I

• g(x) e contrattiva in I, dato che

m = maxx∈I|g′(x)| = 1

3< 1.

Dato che x0 ∈ I, lo schema converge a ξ.Si ha x0 = 2, x1 = g(x0) ≃ 2.0986. La stima a priori sull’ errore ci dice

che∣∣∣∣ξ − xnξ

∣∣∣∣ ≤∣∣∣∣ξ − xn

2

∣∣∣∣ ≤1

2

∣∣∣∣mn

1−m

∣∣∣∣|x1 − x0| ≃ (7.3950E − 2)(1/3)n

quindi deve essere(7.3950E − 2)(1/3)n ≤ 0.003

20

ossia, passando ai logaritmi,

n ≥ ln(0.003/7.3950E − 2)

ln(1/3)≃ 2.9171

ossia n ≥ 3.Utilizzando la stima a posteriori

ǫi =

∣∣∣∣ξ − xiξ

∣∣∣∣ ≤|ξ − xi|

2≤ (1/2)

m

1−m |xi − xi−1| = eps(i)

Otteniamo i risultati

i xi eps(i)

0 2.0000 ?

1 2.0986 2.4650E-2

2 2.1310 8.0855E-3

3 2.1413 2.5959E-3

ed essendo ǫ3 < 0.003, ci possiamo fermare alla terza iterazione e il risultatoe ξ = 2.1413.

21

Capitolo 5

Risoluzione di sistemi lineari.

5.1 Matrici e vettori

A =

a11 a12 . . . a1n

a21 a22 . . . a2n...

... . . ....

am1 am2 . . . amn

= A(m× n)

A = (aij) = (bij) = B,

C = (cij) = (aij + bij) = A + B,

P = (pij) = (

n∑

k=1

aikbkj) = AB = A ∗B.

Il prodotto algebrico AB si puo eseguire solo se A = A(m × n), B =B(n× p).

y : Rn → R

n, y : x 7→ f(x),

z : Rn → R

n, z : x 7→ g(x),

z = (g f)(x) = g(f(x))

se f e g sono lineari e si fissa una base di Rn...

y = Ax, z = By, z = B(Ax) = B ·Ax

Supponiamo che le matrici siano quadrate: m = n.

22

Matrice trasposta di A: AT = ATij = Aji. Se A = AT si dice che e

simmetrica.Matrice coniugata di A ∈ C

n: A∗ = A∗ij = Aij.

Matrice coniugata e trasposta (trasposta coniugata) di A ∈ Cn: A+ =A+ij = (Aji) = Aji. Se A = A+ si dice che e Hermitiana.Matrice inversa, X, di una matrice A = A(n× n):

AX = XA = I, X = A−1.

A−1 esiste sse la matrice e non–singolare, ossia

det(A) 6= 0

ovvero sse le colonne di A sono linearmente indipendenti, sse le righe di Asono linearmente indipendenti.

Il Determinante

det(A) =n∑

j=1

∆ijaij ,

∆ij = (−1)i+jdet(Aij),

Aij = A\riga i, colonna j.

Aij = minore di A relativo all’elemento aij .

5.1.1 Autovettori e autovalori

A ∈ Cn2

, λ ∈ C, 0 6= u ∈ Cn

Definizione 5.1.1

Au = λu

• λ e un autovalore (eigenvalue) di A, λ ∈ λ(A) = spettro di A,

• u e un autovettore (eigenvector) di A, u ∈ V (A) = autospazio di A,

23

Au = λu ⇐⇒ p(λ) = det(A− λI) = 0.

p(λ) e il polinomio caratteristico di A.

Definizione 5.1.2 A e simile a B sse

A = S−1BS.

λ(A) = λ(B).

Definizione 5.1.3 Raggio spettrale di una matrice,

ρ(A) = maxi|λi|, λi ∈ λ(A).

5.1.2 Norme

Norme di vettori

‖·‖ : v ∈ Cn → R

≥0

1. ‖v‖ ≥ 0, ‖v‖ = 0 sse v = 0;

2. ‖cv‖ = |c|‖v‖, c ∈ C;

3. ‖v + w‖ ≤ ‖v‖+ ‖w‖ (disuguaglianza triangolare).

Esempi:

• norma 2 o Euclidea:

‖v‖2 =

√√√√n∑

i=1

v2i =√

v · v =√

vTv.

• norma p:

‖v‖p = (

n∑

i=1

|vi|p)1/p

se p =∞,

‖v‖∞ =n

maxi=1|vi|

24

Norme di matrici

‖·‖ : M ∈ Cn2 → R

≥0

1., 2., 3. +

4. ‖AB‖ ≤ ‖A‖‖B‖ (sub–moltiplicativita).

Esempio. Norma di Frobenius:

‖M‖F =

√√√√n∑

i,j=1

M2ij =

√Tr(MM+) =

√Tr(M+M).

Traccia di una matrice: Tr(A) =∑n

i=1 aii.

Definizione 5.1.4 Una norma di matrice ‖·‖M e naturale, oppure compati-bile con, oppure subordinata a, una norma di vettore ‖·‖v sse

‖A‖M = supv 6=0

‖Av‖v‖v‖v

.

Norme subordinate:

• norma-1 (massima sulle colonne):

‖A‖1 = maxj

n∑

i=1

|aij|,

• norma-infinito (massima sulle righe):

‖A‖∞ = maxi

n∑

j=1

|aij |,

• norma Euclidea (o di Hilbert):

‖A‖E =√ρ(AA+).

Teorema 5.1.1 Se ‖A‖M e compatibile con ‖v‖v, allora ‖A‖M ≥ ρ(A).

Matrici Speciali: vedi tabella 5.1.

25

Tipo matrice Proprieta Elementi in Autovalori Autovettori

Hermitiana H+ = H C reali ortonorm.*

Hermitiana e H+ = H C reali ≥ 0 ortonorm.*SPD u+H u > 0

Simmetrica HT = H ℜ reali ortonorm.*

Simmetrica e HT = H ℜ reali ≥ 0 ortonorm.*SPD u+H u > 0

Antisimmetrica H+ = −H C nulli ortonorm.*o Re = 0

Antisimmetrica H+ = −H ℜ nulli ortonorm.*e reale o Re = 0

λi, λiUnitaria U+U = I = U U+ C |λ| = 1 ?

U+ = U−1

Ortogonale UT U = I = U UT ℜ |λ| = 1 ?UT = U−1 λ = e± iφ

Tabella 5.1: Alcune matrici speciali. * = Gli autovettori possono essere scelti in modo da essere ortonormali.

26

5.2 Sistemi lineari

Dato un sistema lineare

3x1 − 2x2 + x3 = 4x1 + x3 = −7

2x1 + 5x2 − 3x3 = 2

A =

3 −2 11 0 12 5 −3

, b =

4−72

,

x =

x1

x2

x3

.

La forma compatta e:

Ax = b.

5.3 Decomposizione a valori singolari

Mediante matrici ortogonali e possibile diagonalizzare una qualsiasi matriceA ∈ Rm×n.1 Vale il seguente

Teorema 5.3.1 (SVD) Sia A ∈ Rm×n. Esistono due matrici ortogonaliU ∈ Rm×m e V ∈ Rn×n tali che

UTAV = Σ = diag(σi). (5.1)

Σ ∈ Rm×n e una matrice diagonale con elementi σi tali che se r e il rango diA,

σ1 ≥ . . . ≥ σr > σr+1 = . . . = σn = 0. (5.2)

Per la dimostrazione, si veda [Gv89]. Le quantita σi vengono chiamate ivalori singolari di A. ATA e una matrice simmetrica semi-definita positiva,in quanto per ogni vettore x ∈ Rn

xT(ATA

)x =

(xTAT

)(Ax) = (xA)T (Ax) = ‖Ax‖22 ≥ 0. (5.3)

1Sezione sviluppata da M. Lizza

27

I suoi autovalori λi = λi(ATA) sono reali non negativi. Risultato analogo

vale per la matrice AAT . Le colonne ui di U sono dette vettori singola-ri sinistri, le colonne vi di V vettori singolari destri e rappresentano, ri-spettivamente, gli autovettori di AAT e di ATA. Risulta λi = σ2

i . Infatti,AATU = (UΣV T )(UΣV T )TU = UΣV TV ΣUTU = UΣ2 e ATAV =(UΣV T )T (UΣV T )V = V ΣUTUΣV TV = V Σ2.

Corollario 5.3.1 Ogni matrice A ∈ Rm×n ammette una fattorizzazione

A = UΣV T , Σ = diag(σi), (5.4)

dove U ∈ Rm×m e V ∈ R

n×n sono matrici ortogonali e Σ ∈ Rm×n e una

matrice diagonale tale che σ1 ≥ . . . ≥ σn ≥ 0.

Si puo dimostrare che [Gv89]

‖A‖F =(σ2

1 + . . .+ σ2r

)1/2, ‖A‖2 = σ1.

Infatti ‖A‖2F = traccia(ATA) = traccia(V Σ2V T ) =∑

i λi(VΣ2V T ) =∑

i λi(Σ2) =∑

i σ2i . Inoltre ‖A‖22 = ρ(ATA) = σ2

1 .La SVD di una matrice ha interessanti proprieta, alcune condivise da

altre fattorizzazioni, come la fattorizzazione QR con pivoting sulle colonne,altre peculiari della SVD.

5.3.1 Un esempio

Consideriamo la matrice quadrata

A =1 0 10 1 11 0 1

Il suo rango e 2. La decomposizione a valori singolari A = UΣV T e:

U =−0.6280 0.3251 −0.7071−0.4597 −0.8881 0−0.6280 0.3251 0.7071

Σ =2.1753 0 0

0 1.1260 00 0 0

V =−0.5774 0.5774 −0.5774−0.2113 −0.7887 −0.5774−0.7887 −0.2113 0.5774

28

Le matrici U e V sono ortonormali. Osserviamo che Σ ha due elementi nonnulli, pari al rango di A. Costruiamo la matrice

A1 = U1Σ1VT1 = |2.1753| ·

−0.6280−0.4597−0.6280

· −0.5774 0.5774 −0.5774

Risulta‖A− A1‖2 = 1.1260 = σ2.

Come vedremo piu avanti, questo risultato ha portata generale.

5.3.2 Spazi vettoriali

Dalle equazioni (5.1) e (5.4) si ottiene

AV = UΣATU = VΣT ,

ossia Avi = σiuiATui = σivi

, i = 1, . . . ,minm,n. (5.5)

Dalla (5.5) si ricava che

span vr+1, . . . , vn = ker(A), span ur+1, . . . , um = range(A)⊥,

span u1, . . . , ur = range(A), span v1, . . . , vr = ker(A)⊥,

dove ker(A) = x t.c. Ax = 0, range(A) = y t.c. y = Ax, S⊥ = y t.c. xTy =0 ∀x ∈ S.

5.3.3 Rango di matrici

La nozione di rango e di fondamentale importanza nell’algebra lineare. De-terminare numericamente il rango non e un problema banale. Piccole per-turbazioni lo modificano sensibilmente. Vale il seguente

Teorema 5.3.2 (Eckart-Young) Sia A = UΣV T . Se k < r = rank(A) e

Ak =

k∑

i=1

σi · ui · viT , (5.6)

alloramin

rank(B)=k‖A− B‖2 = ‖A− Ak‖2 = σk+1 . (5.7)

29

Per la dimostrazione, si veda [Gv89].In altri termini, Ak e la migliore approssimazione di rango k della matrice

A. Se σk > σk+1 essa e unica. Il teorema di Eckart-Young vale anche per lanorma di Frobenius ‖ · ‖F .

E possibile stimare l’errore relativo commesso quando la matrice A, dirango r > k, viene sostituita con una matrice di rango k:

ek =‖A− Ak‖2‖A‖2

=σk+1

σ1(5.8)

Definizione 5.3.1 Sia ǫ > 0. L’ ǫ-rango, rǫ, della matrice A e:

rǫ = min‖A−B‖2≤ǫ

rank(B).

Dal Teorema 5.3.2 segue che rǫ e pari al numero di valori singolari stret-tamente maggiori di ǫ, ossia

σ1 ≥ . . . ≥ σrǫ > ǫ > σrǫ+1 = . . . = σn = 0 .

Le relazioni (5.6) valgono anche quando si mette Ar al posto di A, r =rank(A).

Vediamo alcune applicazioni della decomposizione a valori singolari.

5.3.4 Sistemi sottodeterminati

Sia A ∈ Rm×n, m < n, di rango r, e sia b ∈ Rm. Il sistema lineare Ax = b esottodeterminato: ha piu incognite che equazioni, percio ha infinite soluzioni,oppure nessuna. Sostituendo ad A la sua SVD e ponendo y = V Tx, c = UT b,si ottiene il sistema equivalente

Σy = c,

la cui soluzione y ha componentiσj yj = cj, j ≤ m, σj 6= 0,0 yj = cj, j ≤ m, σj = 0.

Il secondo insieme di equazioni e vuoto se r = m. Quindi il sistema hasoluzione se e solo se cj = 0 quando σj = 0. Se esistono, le soluzioni y sono

yj = cj/σj , j ≤ m, σj 6= 0,yj = arbitrario, j ≤ m, σj = 0, opp. m ≤ j ≤ n.

Infine, le soluzioni del sistema lineare originario sono i vettori x = V y.

30

5.3.5 Sistemi sovradeterminati

Sia A ∈ Rm×n, m > n, di rango r, e sia b ∈ R

m. Il sistema lineare Ax = b esovradeterminato: ha piu equazioni che incognite. Procedendo come nel casoprecedente, abbiamo:

σj yj = cj , j ≤ n, σj 6= 0,0 yj = cj , j ≤ n, σj = 0,

0 = cj , n < j ≤ m.

Se cj = 0 quando σj = 0, oppure j > n, esistono soluzioni y e sonoyj = cj/σj, j ≤ n, σj 6= 0,yj = arbitrario, j ≤ n, σj = 0.

Data la norma euclidea di vettore, ‖ · ‖2, calcoliamo una soluzione “alter-nativa” minimizzando la quantita

‖r‖22 = ‖Ax− b‖22.

‖r‖22 e minima sse ∂‖r‖22/∂xi = 0, i = 1, . . . , n. Queste relazioni equivalgonoal sistema (vedi la sezione 8.1)

ATAx = AT b, (5.9)

che e detto sistema delle equazioni normali. L’accuratezza della soluzionenumerica risente molto del condizionamento di A, che peggiora in ATA. Uti-lizzando la SVD e possibile risolvere il problema di minimizzazione senzapassare per il sistema lineare (5.9). Sia A = UΣV T . Le matrici ortogonalipreservano la norma euclidea, ossia se H e ortogonale, ‖Hv‖22 = vTHTHv =vTv = ‖v‖22.

‖r‖22 = ‖UT(AV V Tx− b

)‖22

= ‖ΣV Tx− UT b‖22

da cui, ponendo y = V Tx e c = UT b, otteniamo

‖r‖22 = ‖Σy − c‖22.

Quindi ‖r‖22 = 0 se le componenti delle soluzioni y sonoyj = cj/σj, j ≤ n, σj 6= 0,yj = arbitrario, j ≤ n, σj = 0.

Le soluzioni di (5.9) si ottengono dalla relazione x = V y. Se A non ha rangomassimo, esistono infinite soluzioni. Di solito si pone yj = 0 per ogni jtale che σj = 0, ottenendo cosı la soluzione x di norma minima, dato che∑

i y2i = ‖y‖22 = ‖V y‖22 = ‖x‖22.

31

5.3.6 Filtraggio del rumore

Molti problemi nel campo dell’ elaborazione del segnale utilizzano un model-lo di tipo vettoriale, in cui le misurazioni vengono organizzate in una matriceA ∈ Rm×n [Car91]. Le colonne ai ∈ Rm sono ai = si + ri, dove si e rirappresentano rispettivamente la componente del segnale e quella del rumo-re. Il modello assume che l’errore possa essere separato dal dato, e che lacomponente rumore risieda in un sottospazio ortogonale a quello del segnale.

In assenza di rumore si puø calcolare il rango di A. Purtroppo esso vienedi solito sovrastimato a causa degli errori (non sistematici) di misurazione che“disturbano” i dati della matrice. Supponiamo che si risieda in uno spaziovettoriale di dimensione k. Calcoliamo la SVD di A. Dalle equazioni (5.6)risulta si ∈ span u1, . . . , uk , ri ∈ span uk+1, . . . , un . Possiamo costrui-re la matrice Ak e stimare le due componenti (segnale e rumore) sfruttandoil Teorema (5.3.2).

5.3.7 Compressione di immagini

Una immagine puo essere rappresentata mediante una matrice I (bitmap) didimensioni m×n [FvDFH90]. Supponiamo m > n. Il generico elemento, Ii,j,definisce il colore dell’elemento di immagine (pixel). Per le immagini in biancoe nero Ii,j e un valore di luminosita (tipicamente compreso nell’intervallo[ 0, 255 ]N cosicche puø essere memorizzato in un byte). Figure anche semplicirichiedono grandi matrici per essere rappresentate. Per memorizzare e/otrasmettere una matrice–immagine bisogna comprimerla. Gli algoritmi pi˘efficienti per la compressione di immagini sono quelli che comportano unacerta perdita di dettaglio (lossy algorithms).

Utilizzando la decomposizione a valori singolari (SVD) [Gv89] e possibiledefinire una rappresentazione di dimensione ridotta che preserva la maggiorparte della struttura del modello. Una volta calcolata la decomposizione avalori singolari I = UΣV T , si costruiscono le approssimazioni Ik = UkΣkVk

T .Sia U = [u1 . . . um], V = [v1 . . . vn], Σ = diag(σ1, . . . , σn), σ1 ≥ . . . ≥ σr >σr+1 = . . . = σn = 0. Poniamo Σk = diag(σ1, . . . , σk), Uk = [u1 . . . uk],Vk = [v1 . . . vk]. Calcoliamo la matrice ridotta di rango k

Ik = UkΣkVkT . (5.10)

Ci chiediamo se queste approssimazioni possono rimuovere i dettagli presentimantenendo l’“essenza” dell’immagine.

Esaminiamo due immagini. I dati relativi sono riportati nella Tabella 5.2.Nelle Tabelle 5.3 e 5.4 sono riportati i risultati delle prove effettuate. Il para-metro di maggior interesse e il numero di complessivo di valori (detti anche

32

Immagine Dimensione Pixel Rangoearth.eps 250× 257 64250 250gatlin.eps 480× 640 320720 480

Tabella 5.2: Caratteristiche delle immagini considerate

k Celle125 7900062 3527832 1724815 783010 5170

Tabella 5.3: Risultati quantitativi di compressione (earth.eps)

k Celle240 326400120 14880060 7080030 3450015 17025

Tabella 5.4: Risultati quantitativi di compressione (gatlin.eps)

celle o pixel), necessari per memorizzare le matrici Uk, Σk e VkT . Nella

Tabella 5.5 vengono mostrate alcune delle immagini ottenute con le approssi-mazioni. Ad esempio si vede che per k=30 si riduce di nove volte l’occupazio-ne di memoria in gatlin (cf. tabelle 5.4). e 5.5 ) e si ha un’approssimazione diqualita accettabile. Per dimensioni minori, la qualita dell’immagine peggioraconsiderevolmente ma e ancora possibile farsi un’idea dell’aspetto “globale”,soprattutto per quanto riguarda le regioni in cui compaiono linee orizzontalie/o verticali.

Il fattore di compressione non e paragonabile a quello offerto da algorit-mi di compressione lossy come JPEG [Mv96]. E’ confrontabile con quelloofferto da GIF [Mv96], che perø non e lossy. Abbiamo comunque un inte-ressante esempio di come la decomposizione a valori singolari possa essereteoricamente utilizzata per comprimere immagini.

33

gatlin.eps k = 60

k = 30 k = 15

Tabella 5.5: Risultati qualitativi di compressione (gatlin.eps)

5.4 Esercizio sull’ eliminazione di Gauss.

Esercizio.Risolvere il sistema lineare

10−2x+ y = 1

x+ y = 2

con il metodo di eliminazione di Gauss senza pivoting,

a) in aritmetica esatta, ottenendo cosıla soluzione s.

b) In aritmetica a tre cifre decimali con virgola mobile normalizzata.

c) Come nel caso precedente, ma scambiando la prima riga con la seconda.

Nei casi b) e c) calcolare inoltre

er =‖e‖2‖s‖2

(5.11)

Caso a) I risultati sono:

34

matrice A vett. b

0.01 1 1

1 1 2

Eliminazione di Gauss

0.01 1 1

0 -99 -98

La soluzione e quindi: x = 100/99, y = 98/99.

Caso b) Si ottiene

Matrice A Vettore b

1.00e-002 1.00e+000 1.00e+000

1.00e+000 1.00e+000 2.00e+000

Eliminazione di Gauss

1.00e-002 1.00e+000 1.00e+000

0 -9.90e+001 -9.80e+001

Calcolo di y

1.00e-002 1.00e+000 1.00e+000

0 1.00e+000 9.90e-001

Calcolo di x

1.00e+000 0 1.00e+000

0 1.00e+000 9.90e-001

La soluzione e quindi: x = 1.00E+0, y = 9.90E−1. er ≃ 7.1420e−3.

Caso c) Si ottiene

Matrice A Vettore b

1.00e+000 1.00e+000 2.00e+000

1.00e-002 1.00e+000 1.00e+000

Eliminazione di Gauss

1.00e+000 1.00e+000 2.00e+000

0 9.90e-001 9.80e-001

Calcolo di y

1.00e+000 1.00e+000 2.00e+000

0 1.00e+000 9.90e-001

Calcolo di x

1.00e+000 0 1.01e+000

0 1.00e+000 9.90e-001

35

La soluzione e quindi: x = 1.01E+0, y = 9.90E−1. er ≃ 1.4143e−4.

Come si vede, nel caso c) il risultato e piu accurato e corrisponde ad unaeliminazione con pivoting.

36

Capitolo 6

Costo Computazionale

Lo scopo della programmazione numerica e quello di fornire programmi ac-curati ed efficienti per risolvere problemi matematici.1 L’ efficienza di uncodice viene valutata tramite il suo costo computazionale.

Il costo computazionale e volto, in ultima analisi, a “prevedere” il tempodi esecuzione degli algoritmi, senza eseguirli materialmente. Per algoritminumerici, la stima viene data usualmente in senso asintotico [CLRS05]. Se ilcosto “esatto” di un algoritmo A e una complicata espressione C(A) = f(n),funzione del numero di dati, n, e f(n) = O(n), si dice che C(A) = O(n).2

Ad esempio, se

C(A) =a n3 + b n

g n + h+ d n1/2 + e log(n),

allora C(A) = O(n2) (verificarlo per esercizio).

6.1 Operazioni

La stima piu semplice del costo computazionale di un algoritmo numeri-co, consiste nel valutare il numero di operazioni floating point che verrannoeseguite, in funzione della quantita di dati da elaborare. In tal modo si tra-scurano i costi di acquisizione dei dati e produzione dei risultati, che in alcuniprogrammi numerici possono essere non trascurabili, ma l’ idea e che in unprogramma numerico, i costi maggiori sono quelli dovuti alle computazioni,

1Sezione sviluppata con la collaborazione di Andrea Albarelli, Stefano Fedato, LucaVisentin.

2Ricordate che f(x) = O(g(x)) per x → ∞ (si legge f(x) e di ordine “O-grande” di

g(x) per x→∞ [Ber94]), significa che la funzione |f(x)/g(x)| e limitata per x→∞, ossiaesistono due costanti M, N > 0 t.c. |f(x)| ≤M |g(x)|, quando x > N .

37

Tabella 6.1: Brano di programma prosca che esegue il prodotto scalare didue vettori.

1 begin2 s := 0;3 for i := 1, n do4 s := s+ v(i) ∗ w(i);5 endfor6 end

ossia in ultima analisi alle operazioni floating point (FLOP, floating pointoperation).

Consideriamo il brano di programma prosca in tabella 6.1, che calcola ilprodotto scalare dei vettori v e w. Per semplicita chiamiamo “B” questo bra-no. Se si considera unitario il costo di una qualsiasi operazione floating point,il costo computazionale puo essere stimato in base alla dimensione n dei duevettori, ed e C(B) ≃ 2 · n (l’ algoritmo esegue n somme e n moltiplicazioni).Se si considera unitario il costo di una somma (e una sottrazione) e maggioreil costo c > 1 di un prodotto, il costo computazionale e C(B) ≃ n + n · c.Se invece si considera unitario il costo di un’ operazione affine, ossia del ti-po [Nie03] z = a ∗ x + y, il costo e C(B) ≃ n. Se si considera unitarioil costo del prodotto scalare di due vettori di (valore indicativo) m = 128componenti, il costo e C(B) ≃ n/m. Quale di questi e il costo “vero”? Larisposta e che tutte queste stime sono significative, ma una e preferibile al-le altre a seconda dell’ architettura della macchina sulla quale si esegue l’algoritmo. Il fatto e che fino a pochi anni fa una somma veniva eseguitain un tempo molto inferiore a quello necessario per una moltiplicazione, chea sua volta era inferiore a quello speso per una divisione. Tuttavia era, ede, pratica comune semplificare le stime assumendo che tutte le operazioniabbiano lo stesso costo. Oggi i processori hanno unita apposite per eseguireoperazioni floating point, spesso di tipo pipeline [Sta04]. Inoltre, gli algoritmidi moltiplicazione e divisione sono stati molto migliorati. Per queste e altreragioni, il costo di tutte le operazioni floating point puo essere considerato,in pratica, uguale. Esistono inoltre processori che eseguono un’ operazioneaffine ad ogni passo e processori che hanno un’ architettura vettoriale, ossiache eseguono come “operazione atomica” il prodotto scalare tra due vettoridi una certa dimensione. Percio, se la macchina esegue operazioni floatingpoint una alla volta e il costo delle operazioni e piu o meno lo stesso, pervalutare il costo dell’ algoritmo “B” useremo la stima C(B) ≃ 2 ·n, se invecel’ architettura e vettoriale, useremo la stima C(B) ≃ n/m, e cosı via.

38

Chip

RAM

CPU

cache

Dischi

bus

bus

Figura 6.1: Rappresentazione schematica della gerarchia di memoria in unprocessore.

6.2 Accessi alla memoria

In tutte queste stime non abbiamo tenuto conto del tempo di accesso allamemoria, che e stato implicitamente pensato come trascurabile rispetto altempo di computazione. Questa assunzione non e oggi realistica, anzi leCPU elaborano dati in un tempo molto piu basso di quello necessario peracquisire le informazioni dalla memoria [PH05].

La figura 6.1 mostra la rappresentazione schematica della gerarchia dimemoria in un processore. La RAM e la memoria centrale dell’ elaboratore.L’ accesso ai dati su disco e molto piu lento che a quelli nella RAM, ma quan-do si spegne l’ elaboratore i dati nella RAM vanno persi, ecco perche tutte leinformazioni debbono sempre essere memorizzate su disco. Putroppo il costoper bit memorizzato nella RAM e molto piu alto che nei dischi, quindi laquantita di dati memorizzabili nel disco e molto piu grande che nella RAM;solo le informazioni che debbono essere elaborate in un dato istante vengonocopiate nella RAM. Una cache e una memoria tampone ad accesso rapido,che ha dimensione inferiore alla RAM. Le cache di livello 1 (L1) sono inseritenel chip in cui risiede anche la CPU. Lo schema in figura 6.1 e una sempli-ficazione: in realta la gerarchia e piu complessa e vi e piu di una memoriacache anche nei processori piu comuni (vedi ad esempio la figura 6.2). Perchiarirci la necessita delle memorie cache, facciamo un esempio. Consideria-mo un moderno processore a 64 bit, che lavora a 3 GHz, ossia esegue 3× 109

operazioni al secondo. Semplificando un po’, possiamo pensare che, se glioperandi sono disponibili, possa essere eseguita un’ operazione floating pointad ogni cliclo, quindi 3×109 al secondo. Tuttavia per eseguire un’ operazione

39

Figura 6.2: Livelli di cache in un processore AMD Athlon 64.

bisogna procurarsi i valori da elaborare, tramite i bus, prendendoli (in lineadi principio) dalla memoria centrale, o RAM. La RAM non risiede sul chip incui si trova la CPU. Per scambiare informazioni tra i due, bisogna usare uncanale di collegamento, detto bus. Un moderno bus puo lavorare a 333 MHz,ossia fornire 333 × 106 operandi al secondo, prelevandoli dalla memoria. Asua volta la RAM ha una velocita di accesso di 433 MHz, superiore a quelladel bus. Un’ operazione floating point richiede due operandi e fornisce unrisultato che va inviato alla memoria. Ora: 2 operandi + 1 risultato = biso-gna scambiare con la memoria 3 valori ad ogni operazione. A causa del bus,si otterrebbe quindi una performance di 333/3 = 111MHz ≃ 0.1 GFLOPs(GFLOPs = Giga–FLOP/secondo, 109 operazioni al secondo), decisamentescadente!

Le cache sono inserite nel chip assieme alla CPU, la quale puo accedere aidati in essa contenuti ad una frequenza molto maggiore di quella dei bus. Adesempio una cache di primo livello (L1) ha una frequenza di accesso tipica di3 GHz. Se i dati risiedono in cache, si puo quindi ottenere una performance di3 GHz / 3 valori = 1 GFLOPs, molto piu soddisfacente. La dimensione dellacache e piu piccola di quella della RAM, per problemi tecnologici [Sta04]. Adesempio, a fronte di una RAM di 1 GB (1 GB ≃ 109 Byte, 1 Byte = 8 bit),le cache L1 di un comune processore possono contenere 64 KB (1 KB ≃ 103

Byte) ciascuna. Per rendere piu efficiente lo scambio, vengono messe anchecache di livelli superiori, ad esempio cache L2 da 1 MB. La figura 6.2 mostrala gerarchia delle cache di un processore AMD Athlon 64 di tipo K8. Notate

40

ad esempio alla base del disegno, a sinistra la cache di livello 1 (L1) per leistruzioni; a destra la cache di livello 1 per i dati. I processori piu potentihanno cache L2 da 64 MB e L3 da 128 MB.

Le gerarchie di cache debbono essere gestite da algoritmi che, senza l’intervento dell’ utente, si incaricano di far in modo che ad ogni istante dell’elaborazione, siano disponibili nella cache L1 i dati che servono alla CPU.Questi algoritimi sono di vari tipi, ogni tipo ha un nome, ad esempio 2 wayassociative, oppure full associative, oppure 4–way, etc. (vedi figura 6.2). Perapprofondire l’ argomento, si consultino [PH05, Sta04]. Naturalmente accadeche le informazioni necessarie alla CPU non siano disponibili nella cache L1,in tal caso occorre trasferirli dalla cache L2; se mancano anche da questa,dalla cache L3. Ogni volta che un dato che serve non e presente in unacache, si parla di una situazione di cache miss. Le “cache miss” rallentanonotevolmente la velocita di elaborazione. Se i dati non si trovano in alcunacache, occorre recuperarli dalla RAM e se non si trovano nepppure lı, daldisco. Tutti i dati sono sempre memorizzati anche nel disco, perche allospegnimento dellla macchina i dati nelle cache e nella RAM vanno persi.

Il costo dell’ operazione di copiatura dei valori da RAM a Cache e registrie viceversa non e stato considerato nel brano B. Supponiamo che il costodi ogni FLOP sia unitario e il costo medio di una operazione di copiaturaRAM–cache–registri, χ, sia C(χ) = φ≫ 1. Allora C(B) ≃ 2n + 4φn, se adogni iterazione bisogna prendere i valori in v(i), w(i) ed s, e poi scrivere ilrisultato nella variabile s: 4 accessi alla memoria per ognuno degli n elementidei vettori.

Le stime attuali di costo computazionale non tengono conto del costo diaccesso ai vari livelli di memoria, perche dipende fortemente dall’ architetturae non e facile calcolarlo accuratamente. Purtroppo questo fa sı che le stimedi costo computazionale non rispecchino i tempi di esecuzione dei codici.

6.3 Esempi

Supponiamo di voler eseguire un algoritmo numerico su un laptop IBM X31,dotato di processore Intel(R) Pentium(R) M a 1400MHz, con cache di 1024KB. Usiamo il sistema operativo Fedora Core 3, con kernel Linux 2.6.9-1.667.Iniziamo scrivendo codici in ambiente GNU Octave, version 2.1.57. (vedihttp://www.octave.org). L’ ambiente Octave permette di programmarealgoritmi numerici ad alto livello, ma non e adatto a elaborazioni su grandiinsiemi di dati. Per brevita chiameremo “S” l’ esecutore costituito da questamacchina, con questo sistema operativo e l’ ambiente Octave.

Stimiamo il costo di un’ operazione floating point. Non e sufficiente ese-

41

function [csom,cdiv] = flopTime(Nmax);

% tempi di esecuzione di operazioni floating point

v = rand(1,Nmax); w = rand(1,Nmax); z = rand(1,Nmax);

for oper=1:2,

sprintf(’oper=%9d’,oper)

cstv = []; count = 0; nv = [];

n = 64

while (n < Nmax),

count = count + 1;

n = 2 * n;

if (n > Nmax), break, end

if (oper == 1),

[ti,ui,si] = cputime(); % octave

for i=1:n,

z(i) = v(i) + w(i);

end

[tf,uf,sf] = cputime(); % octave

else

[ti,ui,si] = cputime(); % octave

for i=1:n,

z(i) = v(i) / w(i);

end

[tf,uf,sf] = cputime(); % octave

end

tot = tf - ti; usr = uf - ui; sys = sf - si;

cst = usr / n; cstv = [cstv cst];

nv = [nv n];

end

if (oper == 1),

csomv = cstv; csom = cst;

else

cdivv = cstv; cdiv = cst;

end

end

loglog(nv, csomv, nv,cdivv); axis([1 Nmax 1.0e-7 1])

return

Tabella 6.2: Algoritmo Octave per la stima dei tempi di esecuzione dioperazioni floating point.

42

1e-05

1e-04

100 1000 10000 100000

tem

po m

edio

(se

cond

i)

numero operazioni eseguite

2.5E-5

sommedivisioni

Figura 6.3: Tempo medio di esecuzione di n somme o n divisioni sul sistema“S”, in Octave.

guire una sola operazione per stimare il tempo di esecuzione, perche i modernisistemi operativi sono di tipo time-sharing [PH05]. La misura del tempo diCPU impiegato da un processo puo essere influenzata dalla presenza di altriprocessi concorrenti, in particolar modo quando il tempo di esecuzione delprocesso che ci interessa e piccolo. La figura 6.3 mostra il tempo medio dell’esecuzione di una somma e una divisione in Octave, quando vengono eseguite100 ≤ n ≤ 100000 operazioni dello stesso tipo, usando il codice riportato intabella 6.2. Come si vede dalla figura 6.3, i tempi diventano praticamenteuguali quando il numero di operazioni e elevato. Assumiamo quindi che iltempo di esecuzione di un’ operazione floating point sia C⊙ ≃ 2.5e-5 secondi.

La figura 6.4 mostra i risultati se si usa il programma Fortran riportatonelle tabelle 6.3 e 6.4, compilato con g77, a sua volta basato su gcc versione3.4.4. Questa volta il tempo medio di esecuzione di una somma e minore diquello di una divisione. Il calcolo del costo computazionale dovrebbe tenerneconto. Nel seguito, ci limiteremo a considerare algoritmi in Octave, persemplicita, e assumeremo uguali i tempi di esecuzione di tutte le operazionifloating point.

Sia dato un sistema lineare Ax = b; la matrice dei coefficienti, A siastrettamente diagonalmente dominante. Consideriamo l’ algoritmo di elimi-nazione di Gauss senza pivoting, che calcola la matrice triangolare superioreU e il vettore c = L−1b, tali che il sistema LUx = b e equivalente a quellodato, ossia ha la stessa soluzione, e L e una matrice triangolare inferio-

43

program flopTime

c tempi di operazioni floating point

implicit none

c

integer Nmax

parameter (Nmax=10000000)

integer oper, i

real v(1:Nmax), w(1:Nmax), z(1:Nmax), drand

c inizializzazione casuale dei vettori v e w

do i=1,Nmax

v(i) = rand()

w(i) = rand()

end do

c calcolo dei tempi di esecuzione

do oper=1,2

write(6,1)’n’,’secondi’

if (oper .eq. 1) then

write(6,1)’# somme’

else

write(6,1)’# divisioni’

end if

call flopTsub(Nmax,v,w,z,oper)

end do

c

1 format(1x,a12,1x,a13)

2 format(1x,i12,1x,e13.6)

stop

end

Tabella 6.3: Programma Fortran per la stima dei tempi di esecuzione dioperazioni floating point.

44

subroutine flopTsub(Nmax,v,w,z,oper)

c subroutine per il calcolo di tempi

C parametri di scambio: Nmax,v,w,z,oper in input;

implicit none

c

integer Nmax,n, oper, count, i

real v(1),w(1),z(1)

real cst, ui, uf, usr

logical finito

c

count = 0

n = 2048

finito = .false.

do while ((n .lt. Nmax) .and. .not.(finito))

count = count + 1

n = 2 * n

if (n .gt. Nmax) then

finito = .true.

else

if (oper .eq. 1) then

call second(ui)

do i=1,n

z(i) = v(i) + w(i)

end do

call second(uf)

else

call second(ui)

do i=1,n

z(i) = v(i) / w(i)

end do

call second(uf)

end if

usr = uf - ui

cst = usr / n

write(6,2)n,cst

endif

end do

return

2 format(1x,i12,1x,e13.6)

end

Tabella 6.4: Sottorogramma Fortran per la stima dei tempi di esecuzione dioperazioni floating point.

45

1e-08

1e-07

1e-06

1000 10000 100000 1e+06 1e+07 1e+08

tem

po m

edio

(se

cond

i)

numero operazioni eseguite

somme (Fortran)divisioni (Fortran)

Figura 6.4: Tempo medio di esecuzione di somme o divisioni in Fortran.

1e-04

0.001

0.01

0.1

1

10

100

1000

1 10 100 1000

tem

po d

i ese

cuzi

one

(sec

ondi

)

numero incognite

rilevatoteorico

Figura 6.5: Tempo di esecuzione dell’ algoritmo di eliminazione di Gauss suuna matrice di ordine N .

46

function [tot,usr,sys] = gaussTime(n);

% metodo di eliminazione di Gauss, senza pivoting

% su matrice casuale

A = rand(n);

% modifico A per renderla strettamente diag. dom.

for i=1:n,

som = 0;

for j=1:n,

som = som + abs(A(i,j));

end

A(i,i) = 2 * som;

end

b = ((1:n)*0+1)’;

U = A;

c = b;

[ti,ui,si] = cputime(); % octave

for i=1:n,

fact = U(i,i);

for j=i:n

U(i,j) = U(i,j) / fact;

end

c(i) = c(i) / fact;

for j=(i+1):n,

fact = -U(j,i);

for k=i:n

U(j,k) = U(j,k) + fact * U(i,k);

end

c(j) = c(j) + fact * c(i);

end

end

[tf,uf,sf] = cputime(); % octave

tot = tf - ti;

usr = uf - ui;

sys = sf - si;

return

Tabella 6.5: Eliminazione di Gauss senza pivoting.

47

DL

SL

CL

PL

AW

1e−04

0.001

0.01

0.1

1

10

100

1 10 100 1000

tem

po d

i una

ese

cuzi

one

(sec

ondi

)

numero incognite

Figura 6.6: Tempo di una esecuzione dell’ algoritmo di eliminazione di Gausssu varie macchine.

re [Gam94]. Sia N l’ ordine della matrice; il costo asintotico dell’ algoritmodi Gauss [Dem97] e O(3N2/2). La figura 6.5 mostra i tempi di esecuzionerilevati in un solo run del programma riportato nella tabella 6.5. Viene anchemostrata la curva (3N2/2)C⊙, che rappresenta il tempo teorico di esecuzione(asintotico) dell’ algoritmo. Le due curve si distanziano all’ aumentare dell’ordine della matrice: il tempo di recupero dei dati dalla memoria altera itempi di esecuzione, allontanando i valori teorici da quelli rilevati.

La figura 6.6 mostra i tempi di una esecuzione su alcune macchine, lecui caratteristiche sono quelle fornite dai comandi “cat /proc/cpuinfo” e“uname -a” quando si usa Linux, recuperate sotto Windows con comandiopportuni. Per comodita nel seguito individueremo la macchine con unasigla, formata da due lettere:

• DL = model name: AMD Duron(tm); cpu MHz: 1000.214; cache size:64 KB; ram size: 1024 MB; operating system: Suse Linux 10;

• SL = model name: AMD Sempron(tm)2400; cpu MHz: 1665.32; cachesize: 256 KB; ram size: 1024 MB; operating system: Fedora LinuxCore 4;

• CL = model name: Intel(R)Celeron(R)CPU 1.70GHz; cpu MHz: 1693.94;cache size: 128 KB; ram size: 1024 MB; operating system: Suse Linux10;

48

AW

SL

CLPL

DL 1e−04

0.001

0.01

0.1

1

10

100

1 10 100 1000

tem

po m

edio

di e

secu

zion

e (s

econ

di)

numero incognite

Figura 6.7: Tempo medio di esecuzione dell’ algoritmo di eliminazione diGauss su varie macchine.

• PL = model name: Dual Intel(R)Pentium(R)III CPU family 1400MHz; cpu MHz: 1390.948; cache size: 512 KB; ram size: 1024 MB;operating system: Debian Linux;

• AW = model name: AMD Athlon XP 2000; cpu MHz: 1666; cachesize: n.a.; ram size: 480 MB; operating system: Windows XP SP2.

Le curve dei tempi riportate nella figura 6.6 hanno un andamento ab-bastanza variabile quando la dimensione del sistema da risolvere e piccola.Cio e sempre dovuto al fatto che nei moderni sistemi operativi time–sharing,i diversi carichi del sistema influenzano le misure dei tempi di esecuzione.Notate la maggiore variabilita delle curve nel Celeron, quando 10 < n < 100,probabilmente dovuta alla minor quantita di cache a disposizione. Quandola dimensione e “grande”, le variazioni puntuali sono smussate.

La figura 6.7 mostra i tempi medi su 8 esecuzioni. Le oscillazioni sonoleggermente smussate.

Esercizio 6.3.1 Verificare che il costo computazionale del brano di codiceche esegue l’ eliminazione di Gauss senza pivoting e O(3N2/2).

49

6.4 Costo del metodo di Jacobi

Consideriamo il metodo di Jacobi (J) per risolvere il sistema lineare di ordinen

Ax = b. (6.1)

Stimiamo il costo computazionale C[J ] di tale metodo, inteso come numerodi operazioni floating point da effettuare, ponendo C[⊕] = C[⊗] = C[⊘] = 1,dove ⊕, ⊗, ⊘ sono la somma, moltiplicazione e divisione floating point.

Lo schema di Jacobi puo essere messo nella forma [Gam94, QS06]

xk+1 = Exk +D−1b = Exk + c, (6.2)

La matrice di iterazione e E = −D−1(L+ U).Il costo di una iterazione consta di un prodotto matrice–vettore P e di

una somma di vettori S. Assumendo di eseguire il prodotto P con l’algoritmostandard, il costo e di n2 moltiplicazioni e n2 − n addizioni. Il costo di unasomma di vettori e di n − 1 addizioni. Quindi il costo di una iterazione eC[J(n)] = O(2n2)].

Supponiamo di eseguire m iterazioni per ottenere un’ approssimazionexm tale che ‖em‖/‖e0‖ < τ , essendo τ una tolleranza prefissata, ei l’errorealla i-esima iterazione. Allora C[J(n)] = m ·O(2n2).

Sia ‖·‖ la norma euclidea di vettori. Assumiamo che la matrice A siasimmetrica. Allora possiamo stimare m [Gam94],

m ≥ log τ

log ρ(E)= m0.

Esempio. Consideriamo la matrice tridiagonale A che ha gli elementitutti nulli, tranne

Aij =

8, se se i = j,−1, se |i− j| = 1

, i, j = 1, . . . , n.

Gli autovalori di A sono [Gam94]

λA

(k) = 8− sin((n− 1)kπ/(n+ 1))

sin(nkπ/(n+ 1))= 8− µ(k),

k = 1, . . . , n.

Poniamo A = A+10−6. Grazie alla disuguaglianza di Wielandt-Hoffman [Par80],si ha

maxj

∣∣∣λj(A)− λj(A)∣∣∣ < 10−6.

50

Percio ρ(E) = ρ(−D−1(L + U)) = ρ(−D−1(A − D)) = ρ((1/8)A − I) ≃maxk |(1/8)µ(k)| = (1/8) maxk |µ(k)| = (1/8)µ(n). Ora

limn→∞

(1/8)|µ(n)| = 0.25. (6.3)

Una stima asintotica e quindi

m0(τ) = log τ/ log 0.25.

51

Capitolo 7

Interpolazione.

7.1 L’ algoritmo di Neville.

Siano x0, x1, < . . . , xn ascisse date e sia I = [a, b], a = mini=0,...,nxi, b =maxi=0,...,nxi, il piu piccolo intervallo che contiene le ascisse date.

L’ algoritmo di Neville serve per calcolare in un punto dato x ∈ I i valoriP0(x), . . . , Pn(x) dei polinomi

Polinomio grado punti interpolatiP0(x) 0 (x0, y0)P1(x) 1 (x0, y0), (x1, y1). . . . . . . . .Pn(x) n (x0, y0), . . . , (xn, yn)

La tecnica consiste nel costruire la tabellax0 y0 = T00(x)x1 y1 = T10(x) T11(x)x2 y2 = T20(x) T21(x) T22

. . . . . . . . . . . .xn yn = Tn0(x) Tn1(x) . . . Tnn

Gli elementi della tabella si calcolano tramite le relazioni

Ti,0 := yi, 1 ≤ k ≤ i, i = 0, 1, . . . , n

Ti,k :=1

xi − xi−kdet

(Ti,k−1 x− xiTi−1,k−1 x− xi−k

).

Si puo dimostrare che Pi(x) = Tii(x).L’ algoritmo di Aitken (riportato nel testo di G. Gambolati) e una variante

dell’ algoritmo di Neville.

52

7.2 Operatori lineari.

Siano S, Z spazi di funzioni.Un operatore e un funzionale

P : S → Z. (7.1)

Sia D(P ) ⊆ S il dominio di P , ossia l’ insieme delle funzioni per le qualiP e definito e sia R(P ) ⊆ Z la sua immagine o codominio.

Esempio: D : C0 → C0 definito come

Df =df

dx= f ′ (7.2)

e un operatore. Risulta D(D) = C1.P, P1, P2, P3 siano operatori t.c. D(P ) ⊆ D(P1)∩D(P2)∩D(P3) = S 6= ∅,

R(P1) ⊆ D(P2) e sia f ∈ S, α uno scalare. Le operazioni tra operatori sidefiniscono come

(P1 + P2)f = P1f + P2f (7.3)

(P2P1)f = P2(P1f) (7.4)

(αP )f = α(Pf) (7.5)

Indichiamo l’ identita con 1: 1f = f , ∀f ∈ D(P ).Se esiste un operatore P−1

D t.c.

PP−1D = 1, (7.6)

P−1D viene chiamato inverso destro di P .

Se esiste un operatore P−1S t.c.

P−1S P = 1, (7.7)

P−1S viene chiamato inverso sinistro di P .

Si puo dimostrare che se esistono P−1D e P−1

S , allora

P−1S = P−1

D = P−1 (7.8)

Definizione 7.2.1 Un operatore L e detto lineare se

• D(L) e uno spazio lineare,

• per tutti gli f, g ∈ D(L) e per tutti gli scalari α, β risulta

L(αf + βg) = αLf + βLg.

53

7.3 Interpolazione trigonometrica e FFT

f : [0, 2π]→ R, f(0) = f(2π)

f(xk) = yk, xk =2π k

(N + 1), k = 0, . . . , N.

f(xk) =a0

2+

M∑

k=1

(ak cos(kx) + bk sin(kx)),

M = N/2 (N pari).

f(xk) = f(xk), k = 0, . . . , N.

f(x) =M∑

k=−M

ckeikx, i2 = −1.

ak = ck + c−k, bk = i(ck − c−k), k = 0, . . . ,M.

f e detta trasformata reale discreta di Fourier (DFT, Discrete FourierTransform).

Condizioni di interpolazione:

f(xj) =M∑

k=−M

ckeikxj =

M∑

k=−M

ckeikjh = f(xj),

xj = x0 + jh = jh, k = 0, . . . , N.

Tc = f ,

T =

1 1 ... 1ei(−M)h ei(−M+1)h ... ei(+M)h

ei(−M)2h ei(−M+1)2h ... ei(+M)2h

... ... ... ...ei(−M)Nh ei(−M+1)Nh ... ei(+M)Nh

c =

c−Mc−M+1

...c+M

f =

f(x0)f(x1)...

f(xN )

T matrice di Toeplitz: Tjk = φ(j − k).

54

cm =1

N + 1

N∑

j=1

f(xj)e−imjh.

f = Tc = O((N + 1)2)→ O(N log2N)

Trasformata Rapida di Fourier, Fast Fourier Transform, FFT (DFFT).Analogamente:

c = T−1f

Inverse Fast Fourier Transform (IFFT).Data una funzione con molte componenti in frequenza (molte armoni-

che). Supponiamo che la componente di frequenza massima abbia frequenzaµ, allora bisogna campionare con frequenza almeno 2µ (Teorema di Shan-non). Aliasing : campionando con punti “troppo lontani”, si ha una cattivaapprossimazione.

-1

-0.5

0

0.5

1

-pi/2 pi/2 pi 3 pi/2 2 pi0.0

y

x

sin(x)DFT, f=0, xi=k*pi

punti di campionamento

Figura 7.1: Aliasing provocato dal campionamento sui punti xi = kπ. DFT= Discrete Fourier Transform.

... vedi anchehttp://www.kostic.niu.edu/aliasing.htm Ultimo accesso: 23 marzo 2006.

Nyquist frequency may be defined in different ways, namely:1. maximum signal frequency that could be reliably measured to avoid

aliasing = sampling frequency / 2.

55

2. minimum required sampling frequency to avoid aliasing = 2 * signalfrequency.

Intervallo qualsiasi:

f : [a, b]→ R, f(a) = f(b)

z =2π(x− a)(b− a) , z : [a, b]→ [0, 2π]

x =(b− a)z

2π+ a, x : [0, 2π]→ [a, b].

Al posto di cos(k x), sin(k x), useremo

cos

(k

2π(x− a)b− a

), sin

(k

2π(x− a)b− a

)

7.4 Calcolo di splines

Sia ∆ una partizione dell’intervallo [a, b], ossia ∆ = a = x0 < x1 < . . . <xn = b. Sia Y un insieme di valori reali Y = y0, . . . , yn.

Definizione: Una funzione reale S∆ : [a, b] → R e una funzione spline(cubica) relativa a ∆ se

• S∆ e almeno due volte derivabile in [a, b], ossia S∆ ∈ C2[a, b].

• in ogni sottointervallo Vi = [xi, xi+1], i = 0, . . . , n− 1, S∆ coincide conun polinomio di terzo grado:

S∆|Vi∈ P3.

Supponiamo che ∆ e Y siano fissati. e scriviamo S(x) al posto di S∆(x).Consideriamo i seguenti tre tipi di splines:

naturale : S ′′(a) = S ′′(b) = 0;

periodica : S(i)(a) = S(i)(b), i = 0, 1, 2;

vincolata : S ′(a) = y′a, S ′(b) = y′b.

Vogliamo calcolare le spline dei tre tipi, che interpolano i valori Y , ossiatali che S(xi) = yi, i = 0, . . . , n.

Definiamo i momenti Mj = S ′′(xj), j = 0, . . . , n, poniamo poi hj =xj − xj−1, ∆yj = yj − yj−1.

56

Il polinomioSj(x) = S(x)|[xj ,xj+1]

si puø scrivere nel seguente modo [SB80]

Sj(x) = Mj(xj+1−x)3

6hj+1+

Mj+1(x− xj)3

6 hj+1

+ Aj(x− xj) + Bj.

Aj e Bj sono date da

Bj = yj −Mjh2

j+1

6,

Aj =∆yj+1

hj+1− hj+1

6(Mj+1 −Mj).

I momenti si possono calcolare risolvendo il sistema tridiagonale lineareTm = d, dove

T =

2 λ0

µ1 2 λ1

µ2 2 λ2

... ... ...µn−1 2 λn−1

µn 2

m = (M0, . . . ,Mn)T , d = (d0, . . . , dn)

T .

Le definizioni delle componenti sono

λj =hj+1

hj + hj+1, µj = 1− λj , (7.9)

dj =6

hj + hj+1

(∆yj+1

hj+1− ∆yj

hj

), (7.10)

j = 1, . . . , n− 1. (7.11)

Il sistema e tridiagonale: conviene risolverlo con l’algoritmo di Thomas.Servono delle condizioni aggiuntive, che per splinesnaturali, sono λ0 = d0 = µn = dn = 0;vincolate, sono λ0 = µn = 1,

d0 =6

h1

(∆y1

h1− y′0

), dn = − 6

hn

(∆ynhn− y′n

).

57

Le splines periodiche si ottengono risolvendo il sistema T ′m′ = d′, dove

T ′ =

2 λ1 µ1

µ2 2 λ2

µ3 2 λ3

... ... ...µn−1 2 λn−1

λn µn 2

m′ = (M1, . . . ,Mn)T , d′ = (d1, . . . , dn)

T .

Valgono ancora le relazioni (7.9), (7.10), (7.11), e

λn =h1

h1 + hn, µn = 1− λn,

dn =6

h1 + hn

(y1 − ynh1

− yn − yn−1

hn

).

58

Capitolo 8

Approssimazione.

8.1 Approssimazione di funzioni.

Il problema dell’ approssimazione di funzioni si puo schematizzare nel se-guente modo.

E’ data una funzione f(x) che vogliamo analizzare nell’intervallo I =[a, b]. Sia S = a0φ0(x) + . . . + anφn(x), ai sono scalari, lo spazio linearegenerato dalle funzioni φ0(x), . . . , φn(x). Si vogliono calcolare i parametric0, . . . , cn, della funzione f(x) = c0φ0(x) + . . .+ cnφn(x) tale che

d(f, f) = ming∈S

d(f, g) (8.1)

dove d(·, ·) e una opportuna misura di “distanza tra funzioni.La scelta delle funzioni φi(x) e importante. Spesso si pone φi(x) = xi, in

modo che S e lo spazio dei polinomi di grado minore o uguale a n e in talcaso si parla di approssimazioni polinomiali.

Tuttavia se e noto che la funzione da approssimare ha certe caratteri-stiche, sara conveniente scegliere le funzioni φi(x) che meglio rappresentanoqueste caratteristiche. Ad esempio in alcuni casi si prende φk(x) = ekx.

La scelta della “misura di distanza e anch’essa importante e varia a se-conda del tipo di modello da cui scaturisce f . Siano x0, . . . , xn valori dati esia

I = [a, b], a = minx0 . . . , xn, b = maxx0 . . . , xn,il piu piccolo intervallo che li contiene. Se f e nota nei punti x0, . . . , xn,ossia sono noti i valori yi = f(xi), i = 1, . . . , n, una scelta e

d1(f, g) =

n∑

i=0

|f(xi)− g(xi)|. (8.2)

59

In tal caso, se φi(x) = xi e S e lo spazio dei polinomi di grado minore ouguale a n, si parla di problema della interpolazione polinomiale. Sidimostra che il polinomio f(x) che rende minima la distanza (8.2) e unicoed e il polinomio per cui d1(f, f) = 0, ossia f(xi) = yi, i = 1, . . . , n.

Vengono definite anche distanze in cui entrano le derivate della funzione,ad esempio

d(1)1 (f, g) =

n∑

i=0

(|f(xi)− g(xi)|+ |f ′(xi)− g′(xi)|). (8.3)

e in tal caso si parla di interpolazione di Hermite.Se si usa la distanza d1 ma al posto dei polinomi si usano funzioni razionali

del tipo

Rn,m(x) =Pn(x)

Qm(x)=a0 + . . .+ anx

n

b0 + . . .+ bmxm(8.4)

abbiamo la cosiddetta approssimazione di Pade. In questo caso non edetto che vi sia sempre soluzione al problema, ossia che esista una funzioneRn,m(x) tale che

Rn,m(xi) = yi, i = 0, 1, . . . , n+m (8.5)

Se questo problema ha soluzione, la funzione che si ottiene e particolarmenteefficace nel rappresentare funzioni che hanno delle singolarita, che sono invecemal approssimate da polinomi.

Un altro caso importante e quello della interpolazione polinomiale atratti: mostreremo che non e conveniente interpolare punti con polinomi digrado elevato. Se i punti sono molti, conviene interpolarli (ossia calcolare lafunzione che rende uguale a zero la distanza (8.2)) con polinomi continui atratti, ossia funzioni che sono polinomi opportunamente “raccordati per darevita a una funzione definita in tutto l’ intervallo I.

Se d = d2, dove

d2(f, g) =

n∑

i=0

(f(xi)− g(xi))2, (8.6)

il problema di calcolare la funzione g che minimizza d2 in uno spazio di di-mensione minore di n, viene detto approssimazione discreta ai minimiquadrati. Se S e uno spazio di polinomi, viene detta approssimazione po-linomiale discreta ai minimi quadrati. Notare che se dim(S) = n, allora lasoluzione e quella per la quale d2(f, g) = 0, ossia il polinomio interpolatore!

Se d = d2, ma

S = 1, cos(x), sin(x), cos(2x), sin(2x), . . . , cos(kx), sin(kx),

60

Figura 8.1: I polinomi interpolatori di grado elevato presentano fortioscillazioni.

61

Figura 8.2: Esempio di approssimazione ai minimi quadrati.

62

si parla di approssimazione trigonometrica discreta ai minimi quadrati.Se d = d

(c)2 , dove

d(c)2 =

∫ b

a

(f(x)− g(x))2dx (8.7)

si parla di approssimazione continua ai minimi quadrati.Infine se d = d∞ dove

d∞(f, g) = supx∈[a,b]

|f(x)− g(x)|, (8.8)

si parla di approssimazione uniforme o approssimazione di Cheby-chev in [a, b]. Se S e uno spazio di polinomi, si parla di approssimazionepolinomiale uniforme e la soluzione e detta polinomio di Chebychev.

Il problema dell’ approssimazione di Chebychev e in generale piu difficileda risolvere di quello dell’ approsimazione ai minimi quadrati, sia nel casocontinuo che in quello discreto. In teoria l’ approssimazione di Chebychevsarebbe preferibile perche fornisce una stima buona su tutto l’ intervallo,mentre la stima con la distanza d

(c)2 puo fornire una soluzione f

(c)2 (x) che

differisce molto da f(x) in un ristretto sottointervallo di [a, b].

63

Capitolo 9

Integrazione numerica.

Esercizio 9.0.1 Calcolare le approssimazioni del valore di

I =

∫ 1

0

x2dx (9.1)

che si ottengono usando gli schemi di integrazione

1. dei trapezi semplice,

2. dei trapezi composta, con n = 2 suddivisioni.

3. di Gauss-Legendre a tre punti.

Calcolare le stime dei valori assoluti degli errori commessi, chiamiamolirispettivamente E1, E2 e E3, usando le stime teoriche e confrontarli con glierrori effettivamente generati.

Soluzione.Usando la regola dei trapezi semplice abbiamo:

I ≃ 10 + 1

2= 1/2.

Usando la regola dei trapezi composta, con n = 2 risulta:

I ≃ 1

2(0 + 1

2+

1

4) = 3/8.

Poiche l’ intervallo di integrazione non e [−1, 1], non possiamo applicaredirettamente la formula di Gauss-Legendre, tuttavia ponendo x = 1

2y + 1

2,

abbiamo

I =1

8

∫ 1

−1

(y + 1)2dx.

64

n pesi ascisse1 2 02 1 ±0.57735 026923 5/9 ±0.77459 66692

8/9 04 0.34875 48451 ±0.86113 63116

0.65214 51549 ±0.33998 104365 0.23692 68851 ±0.90617 98459

0.47862 86705 ±0.53846 931010.56888 88889 0

Tabella 9.1: Pesi ed ascisse per lo schema di integrazione di Gauss-Legendre.

Consultando una tabella con i pesi e le ascisse di Gauss-Legendre e usandolo schema a 3 punti otteniamo

I =1

8(0.5556(−0.7746+ 1)2 +0.8889(0+1) +0.5556(0.7746+1)2) ≃ 0.3333.

Per quanto riguarda gli errori, risulta

E1 ≤∣∣∣∣(b− a)3

12f ′′(ξ)

∣∣∣∣ = 1/6,

E2 ≤∣∣∣∣(b− a)3

1222f ′′(ξ)

∣∣∣∣ = 1/24,

e infineE3 = 0

perche f (5)(x) = 0.I valori esatti degli errori sono uguali a quelli stimati.

9.1 Estrapolazione di Romberg

T (h) =b− an

(f0 + fn

2+

n∑

i=1

fi

)(9.2)

f ∈ Cm+2[a, b]

65

T (h) =

∫ b

a

f(x)dx+ τ1h2 + τ2h

4 + . . .+

τmh2m + αm+1(h)h

2m+2

• |αm+1(h)| ≤Mm+1 per ogni h = (b− a)/n.

• τk non dipende da h.

T (0) =

∫ b

a

f(x)dx (9.3)

h0 = b− a, h1 = h0/n1, h2 = h0/n2, . . . (9.4)

n1 < n2 < n3 < . . .

Ti,0 = T (hi), i = 1, . . . , m

Tm,m(h) = a0 + a1h2 + a2h

4 + . . .+ amh2m

polinomio interpolatore

Tm,m(hi) = T (hi), i = 0, . . . , m

Estrapolare: Tm,m(0) ≃ T (0).Integrazione di Romberg: hi := b−a

2i .Neville:h0 T (h0) = T0,0

h1 T (h1) = T1,0 T1,1

h2 T (h2) = T2,0 T2,1 T2,2

. . .

Ti,0 = T (hi),

Ti,k = Ti,k−1 +Ti,k−1 − Ti−1,k−1[

hi−k

hi

]2− 1

Puo essere che estrapolazione = “vecchie formule”. Esempio:h0 = b− a, h1 = (b− a)/2.T0,0 = b−a

2(f(a) + f(b)) ← trapezi

T1,1 = b−a2

13(f(a) + 4f((a+ b)/2) + f(b)) ← Simpson

66

9.2 Integrazione di Gauss

I(f) :=

∫ b

a

f(x)ω(x)dx

• ω(x) ≥ 0 misurabile su [a, b].

• Esistono tutti i momenti µk =∫ baxkω(x)dx, k = 0, 1, . . ..

•∫ baω(x)dx > 0.

Calcolare

I(f) ≃ I(f) :=

n∑

i=1

wif(xi)

... in modo che l’ uguaglianza sia esatta fino al maggior grado possibile.Pj := polinomi di grado ≤ j Pj := polinomi monoici di grado ≤ j Prodotto scalare

(f, g) =

∫ b

a

f(x)g(x)ω(x)dx

L2[a, b] = insieme delle funzioni “a quadrato sommabile.(f, g) = 0 : f e g sono ortogonali.

Teorema. Per ogni k ∈ N , esistono polinomi pj ∈ Pk, j = 0, 1, . . . , n,tali che

(pi, pj) = δij

Questi polinomi sono:

p0(x) = 1, p−1(x) = 0,

pi+1(x) = (x− δi+1)pi(x)− γ2i+1pi−1(x), i ≥ 0

• δi+1 := (x pi, pi)/(pi, pi), i ≥ 0.

• γ2i+1 :=

0, se i = 0,

(pi, pi)/(pi−1, pi−1), i ≥ 1

67

Dato che se p ∈ Pj allora p =∑j

i=0 ci pi, risulta(p, pk) = 0 per ogni p ∈ Pk−1.

Teorema. Le radici xi, i = 1, . . . , n di pn sono

• reali e semplici,

• contenute nell’ intervallo [a, b].

Teorema. La matrice

A :=

(p0(t1) . . . p0(tn)pn−1(t1) . . . pn−1(tn)

)

e non singolare se i punti ti sono distinti.

Quindi il sistema lineare

AT

0c1. . .cn−1

= fi

ha una e una sola soluzione, ossia esiste p(x) t.c. p(x) =∑n−1

i=0 cifi(x),p(ti) = fi.

A viene chiamata Matrice di Haar edet(a) 6= 0 viene chiamata condizione di Haar.p0, . . . , pn soddisfano la condizione di Haar: si dice che formano un sistema

di Chebychev.

Teorema. Siano p0, . . . , pn polinomi ortogonali. Siano x1, . . . , xn le radicidi pn e sia w = (w1, . . . , wn)

T la soluzione del sistema lineare

(p0(x1) . . . p0(xn)pn−1(x1) . . . pn−1(xn)

)

w1

. . .wn

(p0, p0)00

(9.5)

allora i valori wi (detti pesi) sono maggiori di zero e

68

∫ baf(x)ω(x)dx =

∑ni=1wip(xi)

per ogni p ∈ P2n−1.(9.6)

Viceversa: se wi, xi soddisfano (9.6) allora pn(xi) = 0 e w = (w1, . . . , wn)soddisfa (9.5).

NON esistono 2n numeri wi, xi tali che (9.6) valga per ogni p ∈ P2n.

Esempio 9.2.1 ω(x) = 1, intervallo = [−1, 1]:polinomi di Legendreformula ricorrente:

p−1 = 0, p0 = 1

pn+1(x) =2n+ 1

n+ 1x pn(x)−

n

n+ 1pn−1(x)

definizione alternativa:

pk(x) =k!

(2k)!

dk

dxk(x2 − 1)k

Come si calcolano i wi e gli xi?

J =

δ1 γ1 0γ1 δ2 γ2 0

. . .0 γn−2 δn−1 γn−1

0 0 γn δn

(9.7)

det(J − xI) = pn(x)

pn(x) = 0 sse x ∈ λ(J)

I wi si possono calcolare a partire da V (J).

Teorema. Sia f ∈ C2n[a, b],esiste ξ ∈ [a, b] tale che

∫ b

a

f(x)ω(x)dx−n∑

i=1

wif(xi) =f (2n)(ξ)

(2n)!(pn, pn)

69

Capitolo 10

Equazioni differenzialiordinarie.

10.1 Cenni sulla risolubilita delle ODE.

Si chiama problema di Cauchy il problema di trovare una funzione y continuae derivabile nell’ intervallo I0 = [x0, x1] tale che

y′(x) = f(x, y(x)), x ∈ I0 (10.1)

y(x0) = y0 (10.2)

La condizione (10.2) si chiama condizione di Cauchy e anche condizioneiniziale.

Una funzione y che verifica la (10.1) viene chiamata un integrale dell’equazione differenziale.

La equazione (10.1) viene detta del primo ordine perche in essa inter-viene solo la derivata prima di y.

Vale il

Teorema (Cauchy-Lipschitz): se la funzione f e continua su I0 × ℜed esiste una costante reale positiva L tale che

|f(x, y)− f(x, z)| ≤ L|y − z|, (10.3)

per tutte le coppie (x, y), (x, z) ∈ I0 × ℜ, allora il problema(10.1)+(10.2) ammette una ed una sola soluzione.

Per poter risolvere numericamente un problema differenziale, non bastache la soluzione esista e sia unica, bisogna anche che essa sia ben condizio-nata. Per studiare il condizionamento del problema, dobbiamo studiare il

70

comportamento della soluzione quando i dati cambiano. A questo proposito,sia ψ(x) una funzione continua in I0, β ed ǫ numeri reali. Come problemaperturbato consideriamo il seguente:

trovare una funzione yǫ ∈ C1(I0) tale che

y′ǫ(x) = f(x, yǫ(x)) + ǫψ(x), x ∈ I0 (10.4)

yǫ(x0) = y0 + ǫβ (10.5)

Vale la

Proposizione: se valgono le condizioni del teorema di Cauchy-Lipschitz, allora il problema (10.4)+(10.5) ammette una ed una solasoluzione e per tutti gli x ∈ I0 risulta

|y(x)− yǫ(x)| ≤ |ǫβ|eL|x−x0| +

∫ x

x0

|ǫψ(s)|ds (10.6)

La relazione (10.6) implica la dipendenza continua della soluzione daidati, una caratteristica molto importante dal punto di vista della risolubilitanumerica.

Si dice che un problema e ben posto se la sua soluzione esiste, e unica edipende con continuita dai dati. Un problema differenziale e numericamenterisolubile se e non solo ben posto, ma anche ben condizionato. La relazione(10.6) mostra che il problema in esame e ben condizionato quando la quantita

maxx∈I0

eL|x−x0| (10.7)

non e “troppo grande.Consideriamo infine il problema della regolarita, ossia del grado di de-

rivabilita, della soluzione in dipendenza dalla regolarita dei dati. L’ ordi-ne di convergenza dei metodi numerici dipende infatti dalla regolarita dellasoluzione. Vale il

Teorema: sia y ∈ C1(I0) una soluzione del problema (10.1)+(10.2).Se f e derivabile p volte nel punto (x, y(x)), allora nello stesso puntola soluzione y e derivabile (p+ 1) volte e si ha

y(p+1)(x) = f (p)(x, y(x)).

71

Molte interessanti considerazioni sulle ODE e sui metodi per la lorointegrazione numerica si trovano in [EEHJ96].

Esercizio 10.1.1 Si vuole integrare numericamente il problema differenziale

y′ + y = 0, x ∈ [0, 1] (10.8)

y(0) = y0 = 1 (10.9)

Calcolare le soluzioni approssimate y(1) e y(2) che si ottengono usando rispet-tivamente gli schemi alle differenze:

(1 + h)yn+1 − yn = 0 (10.10)

yn+2 − 4yn+1 + (3− 2h)yn = 0 (10.11)

con passo h = 1/2.Calcolare le norme 1 dei vettori degli errori relativi che si commettono

utilizzando i due schemi.

Soluzione.Tramite le relazioni ricorrenti, e immediato calcolare le approssimazio-

ni richieste. L’ unico problema e il valore di y(2)1 . Poniamo y

(2)1 := y

(1)1 .

Otteniamo cosı:

x exp(-x) y1 y2

0 1.0000 1.0000 1.0000

0.5000 0.6065 0.6667 0.6667

1.0000 0.3679 0.4444 0.6667

dove y1 := y(1), y2 := y(2).La soluzione esatta del problema e y(x) = exp(−x). Le norme 1 dei

vettori degli errori relativi

∥∥e(j)∥∥

1=

2∑

i=0

∣∣∣∣∣y(xi)− y(j)

i

y(xi)

∣∣∣∣∣

sono quindi:∥∥e(1)

∥∥1≃ 0.3072,

∥∥e(2)∥∥

1≃ 0.9114.

72

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.5

1.0

1.5

2.0

2.5

3.0

h = 1/3

h = 1/5

exp(-x)

Schema instabile di ordine 2.

Come si vede, l’ errore relativo generato con lo schema del secondo ordinee piu elevato, per il fatto che la soluzione ha una componente spuria checresce rapidamente. E’ interessante vedere che cosa accade quando si usatale schema e si diminuisce il passo di integrazione. In figura si puo vedere ilcomportamento della soluzione esatta e di quelle approssimate ottenute conpassi h = 1/3 e h = 1/5: le soluzioni dell’ equazione alle differenze cresconomolto rapidamente e non convergono alla soluzione esatta.

Esercizio 10.1.2 Approssimare la soluzione del problema differenziale

y′ + y = x2, x ∈ [0, 1] (10.12)

y(0) = y0 = 3 (10.13)

utilizzando lo schema di Eulero.Analizzare la stabilita e la convergenza dello schema alle differenze che si

ottiene.

Soluzione.Lo schema di Eulero e

yn+1 = yn + hf(xn, yn).

73

In questo caso abbiamo il problema:

y′ = f(x, y) = −y + x2,

l’ equazione alle differenze risultante e

yn+1 = yn + h(−yn + x2n),

ossia, dato che xn = x0 + nh = nh:

yn+1 − (1− h)yn = n2h3.

La soluzione generale di questa equazione alle differenze e

y∗n = c(1− h)n + (hn)2 − 2(hn) + 2− h,

che, tenendo conto delle condizioni iniziali, si particolarizza in

yn = (1 + h)(1− h)n + (hn)2 − 2(hn) + 2− h.

La soluzione del problema differenziale e invece

y = e−x + x2 − 2x+ 2.

E’ facile vedere che per h → 0 e n → ∞ in modo che nh = x= cost., siha

limh→0

(1 + h)(1− h)n = e−x, limh→0

[(hn)2 − 2(hn) + 2− h

]= x2 − 2x+ 2,

ossia lo schema e convergente:

limh→0

yn = y(x).

Per quanto riguarda la stabilita, se yk e affetto da un errore ǫk, abbiamo:

yn+1 + ǫn+1 − (1− h)(yn + ǫn) = n2h3,

da cui l’ errore soddisfa alla relazione

ǫn+1 − (1− h)ǫn = 0.

Quindiǫn = ǫ0(1− h)n,

e lo schema e stabile a condizione che |1− h| ≤ 1, cioe 0 ≤ h ≤ 2.

74

0 0 0

totale

L’ errore nelle ODE

troncamentoarrotondamento

h ottimaleFigura 10.1: L’ errore nelle ODE.

75

10.2 Errori

L’ errore di troncamento tende a zero con il passo di discretizzazione, h.L’ errore di arrotondamento cresce quando h diventa piccolo.L’ errore totale, che e una composizione dei due, e minimo per un certo

valore di h, detto ottimale.

10.3 Sistemi di ODE.

Supponiamo di avere un sistema di m equazioni differenziali ordinarie nelleincognite y1(t), . . . , ym(t)

y′1 = f1(t, y1, . . . , ym). . . = . . .y′m = fm(t, y1, . . . , ym)

(10.14)

il dominio sia t ∈ (t0, T ], le condizioni iniziali

y1(t0) = y(0)1 , y2(t0) = y

(0)2 , . . . , ym(t0) = y(0)

m . (10.15)

Per risolverlo, si puø applicare a ciascuna equazione uno dei metodi stu-diati per le singole equazioni. Ad esempio, applicando il metodo di Euleroesplicito, otteniamo le relazioni ricorrenti:

y(n+1)1 = y

(n)1 + hf1(t, y

(n)1 , . . . , y

(n)m )

. . . = . . .

y(n+1)m = y

(n)m + hfm(t, y

(n)1 , . . . , y

(n)m )

(10.16)

Esempio 10.3.1 Risolvere il sistema di ODE:

y′1 = +C1y1(1− b1y1 − d2y2),y′2 = −C2y2(1− b2y2 − d1y1),

(10.17)

dove C1 = C2 = 1, b1 = b2 = 0, d1 = d2 = 1, t ∈ [0, 10], y1(0) = y2(0) = 2.Si ottiene lo schema:

y(n+1)1 = +y

(n)1 + hy

(n)1 (1− y(n)

2 ),

y(n+1)2 = −y(n)

2 + hy(n)2 (1− y(n)

1 ).

y(0)1 = 2, y

(0)2 = 2.

(10.18)

76

10.4 Runge-Kutta per sistemi di ODE.

Ecco un esempio di applicazione di uno schema del tipo Runge-Kutta (loschema di Heun) per risolvere un sistema di equazioni differenziali.

Risolvere il problema differenziale:

y′′ + 2y′ + y = x+ 1, x ∈ [0, 1] (10.19)

y(0) = 0, y′(0) = 0. (10.20)

La soluzione di questo problema e y = e−x + x− 1.L’ equazione (10.19) e equivalente al sistema di equazioni

y′ = z = f1(x, y, z)z′ = −2z − y + x+ 1 = f2(x, y, z)

(10.21)

Lo schema di Heun si scrive:

yn+1 = yn +1

2k1 +

1

2k2 (10.22)

k1 = hf(xn, yn) (10.23)

k2 = hf(xn + h, yn + k1). (10.24)

La sua ovvia generalizzazione ad un sistema di due equazioni in due incognitee

yn+1 = yn +1

2ky,1 +

1

2ky,2 (10.25)

zn+1 = zn +1

2kz,1 +

1

2kz,2 (10.26)

ky,1 = hf1(xn, yn, zn) (10.27)

kz,1 = hf2(xn, yn, zn) (10.28)

ky,2 = hf1(xn + h, yn + ky,1, zn + kz,1) (10.29)

kz,2 = hf2(xn + h, yn + ky,1, zn + kz,1) (10.30)

dove i termini ky,∗, kz,∗ indicano incrementi riguardanti y oppure z.Applicando questa generalizzazione al sistema (10.21) otteniamo le rela-

zioni

yn+1 = yn +1

2ky,1 +

1

2ky,2 (10.31)

zn+1 = zn +1

2kz,1 +

1

2kz,2 (10.32)

ky,1 = hzn, (10.33)

kz,1 = h(−2zn − yn + xn + 1) (10.34)

ky,2 = h(zn + kz,1) (10.35)

kz,2 = h(−2(zn + kz,1)− (yn + ky,1) + (xn + h) + 1) (10.36)

77

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4Soluzione esatta

Schema di Heun; h=0.1

Soluzione esatta

Risoluzione di un sistema di ODEs.

Figura 10.2: Sistema di ODE risolto con lo schema di Heun.

La figura (10.2) mostra l’ andamento della componente yn della soluzionedi questo sistema di equazioni alle differenze, usando un passo h = 0.1 eimponendo le condizioni iniziali (10.20).

10.5 Risoluzione di un BVP con shooting.

E’ dato il problema con valori al contorno (BVP)

y′′ + 2y′ = 4x+ 3ex, x ∈ [0, 1] (10.37)

y(0) = 1, y(1) = e. (10.38)

Vogliamo calcolare un’ approssimazione numerica della sua soluzione, usandoun metodo di shooting. Esso consiste nel risolvere i problemi

y′′ + 2y′ = 4x+ 3ex, x ∈ [0, 1] (10.39)

y(0) = 1, y′(0) = s (10.40)

per diversi valori del parametro s: ad ogni valore corrispondera una soluzioney∗(x; s) che in generale non soddisfa alla condizione al contorno in x = 1,ossia risulta y∗(1; s) 6= y(1). Risolvere il problema con le condizioni alcontorno (10.38) equivale a calcolare la soluzione s dell’ equazione

Φ(s) = y∗(1; s)− y(1) = 0. (10.41)

78

Si puo determinarla con qualsiasi schema per calcolare le soluzioni di un’equazione non lineare. Useremo lo schema della regula falsi

sk+1 = sk − Φ(sk)/Φ(sk)−Φ(sk−1)

sk−sk−1(10.42)

k = 1, 2, . . . , s0, s1 valori dati. (10.43)

Per integrare numericamente l’ equazione (10.37) usiamo lo schema alledifferenze che si ottiene usando le approssimazioni

y′n =∆ynh

+O(h), (10.44)

y′′n =∆2ynh2

+O(h), (10.45)

e scrivendo l’ equazione (10.37) nel punto xn = nh, ossia lo schema

yn+2 + 2(h− 1)yn+1 + (1− 2h)yn = 4h2xn + 3h2exn. (10.46)

Stimiamo ora i due valori iniziali di s. Poniamo ad esempio

s0 =y+δ

1 − y(0)

h, s1 =

y−δ1 − y(0)

h, (10.47)

h = 1/2, dove i valori y+δ1 , y−δ1 sono stati ricavati dalla relazione (10.46)

ponendo h = 1/2, y0 = y(0) = 1, y±δ2 = y(1) ± δ, essendo δ un valorearbitrario, diciamo δ = 1.

Utilizzando questi due valori iniziali e applicando lo schema (10.42) otte-niamo i risultati:

x y y0 y1 y2 yb

0.0000 1.0000 1.0000 1.0000 1.0000 1.0000

0.1000 1.0152 0.4037 0.8037 1.0458 1.0000

0.2000 1.0614 -0.0434 0.6766 1.1125 1.0300

0.3000 1.1399 -0.3639 0.6121 1.2030 1.0912

0.4000 1.2518 -0.5757 0.6051 1.3201 1.1847

0.5000 1.3987 -0.6926 0.6520 1.4662 1.3121

0.6000 1.5821 -0.7254 0.7503 1.6438 1.4747

0.7000 1.8038 -0.6822 0.8984 1.8554 1.6743

0.8000 2.0655 -0.5689 1.0956 2.1034 1.9126

0.9000 2.3696 -0.3899 1.3417 2.3901 2.1917

1.0000 2.7183 -0.1479 1.6374 2.7183 2.5137

79

0.0 0.2 0.4 0.6 0.8 1.0

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0 0.2 0.4 0.6 0.8 1.0

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

2.5

3.0

Soluzione esatta

y0

y1

Applicaz. regula falsi.

y’(0) = 0

Risoluzione di un BVP.

Figura 10.3: Applicazione del metodo di shooting con uso della regula falsi.

dove y e la soluzione esatta y = x2 − x+ ex del problema dato, yi = y(i) =yj(si) = y(xj ; si), j = 1, . . . n = 1/h sono le soluzioni dello schema(10.46), ottenute ponendo di volta in volta

y1(si) = y0 + hsi, (10.48)

dove s0 e s1 sono dati dalla (10.47) e s2 quello ottenuto applicando la regulafalsi (10.42). Notare come basta una sola iterazione con lo schema della regulafalsi per ottenere una buona approssimazione della soluzione cercata. Se sicalcola anche y(3) si vede che y(2) = y(3): non si guadagna nulla procedendooltre la prima iterazione.

yb = y(xn; 0) e la soluzione numerica che si ottiene imponendo la con-dizione iniziale esatta y′(0) = 0. E’ interessante notare come y(2) sia piuaccurata di y(xn; 0): conoscere esattamente la condizione iniziale sulla deri-vata e meno vantaggioso che conoscere il valore della funzione nell’ estremofinale dell’ intervallo.

La figura (10.3) fornisce uno schizzo delle curve che si ottengono plottandoi valori trovati.

10.6 BVP e sistemi alle DF

E’ dato il problema con valori al contorno (BVP)

y′′ + 2y′ = 4x+ 3ex, x ∈ [0, 1] (10.49)

y(0) = 1, y(1) = e, (10.50)

80

che ha la soluzioney = x2 − x+ ex. (10.51)

Si puo determinare una soluzione numerica di questo problema usandouno schema alle differenze finite. Usiamo lo schema alle differenze

yn+2 + 2(h− 1)yn+1 + (1− 2h)yk = 4h2xn + 3h2exn (10.52)

xn = nh, che si ottiene dalla equazione (10.49) impiegando note approssima-zioni dell’ operatore differenziale. Posto h = 0.25, abbiamo y0 = y(0) = 1,y4 = y(1) = e, e scrivendo le relazioni (10.52) per n = 0, 1, 2, otteniamo unsistema lineare di 3 equazioni in tre incognite

As = b (10.53)

dove

[ 2 * (h - 1) 1 0 ]

A = [ (1 - 2 * h) 2 * (h - 1) 1 ]

[ 0 (1 - 2 * h) 2 * (h - 1)]

[ -1.5000 1.0000 0 ]

A = [ 0.5000 -1.5000 1.0000 ]

[ 0 0.5000 -1.5000 ]

[ (3 * h*h - (1 - 2 * h) ) ]

b = [ (4 * h*h * h + 3 * h*h * exp(h) ) ]

[ (4 * h*h * 2*h + 3 * h*h * exp(2*h) - exp(1))]

[ -0.3125 ]

= [ 0.3033 ]

[ -2.2841 ]

Risolvendo questo sistema lineare otteniamo il risultato

x y yn

0 1.0000 1.0000

0.2500 1.0965 1.2673

0.5000 1.3987 1.5884

0.7500 1.9295 2.0522

1.0000 2.7183 2.7183

81

-1.0 -0.5 0.0 0.5 1.0

x

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

u(x;

eps)

-1.0 -0.5 0.0 0.5 1.0

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

eps=0.03125

eps=0.125

eps=0.5

eps=0

solution of 4 eps u_xx + 2 u_x = (x+1)/2

Figura 10.4: Soluzioni del problema per alcuni valori di ǫ.

10.7 Strato limite

10.7.1 Formulazione del problema

Si vuole risolvere numericamente il problema [HS93]

4ǫuxx + 2ux = (x+ 1)/2, x ∈ (−1, 1), ǫ > 0,

u(−1) = u(1) = 0(10.54)

La sua soluzione analitica e

u = −(4 ǫ− 1) e12 ǫ

− x2 ǫ

2 e1ǫ − 2

+

+x2 + (1− 4 ǫ) x+ 8 ǫ2 − 2 ǫ

4−

+(4 ǫ2 − 3 ǫ+ 1) e

1ǫ − 4 ǫ2 − ǫ

2 e1ǫ − 2

Quando ǫ → 0, la soluzione di questa equazione ha uno strato limite oboundary layer in x = −1. Il significato di questa affermazione si puo capireguardando le figure 10.4 e 10.5.

82

-1.0 -0.9 -0.8 -0.7 -0.6 -0.5

x

-0.5

-0.48

-0.46

-0.44

-0.42

-0.4

u(x;

eps)

-1.0 -0.9 -0.8 -0.7 -0.6 -0.5

-0.5

-0.48

-0.46

-0.44

-0.42

-0.4

u05

u0125

u003125

ulim=(x+1)^2/8 - 1/2

solution of 4 eps u_xx + 2 u_x = (x+1)/2

Figura 10.5: Dettaglio della figura precedente.

10.7.2 Risoluzione numerica

Definiamo una griglia in [0,1] ponendo

xi = −1 + ih, i = 0, 1, . . . , N

N intero positivo dato, h = 2/N .Lo schema alle differenze finite (FD) del secondo ordine centrato (centered

FD) e definito come

4ǫuj+1 − 2uj + uj−1

h2+2

uj+1 − uj−1

2h=xj + 1

2, j = 1, 2, . . . , N−1. (10.55)

Lo schema del primo ordine upwind (Upwind FD) come

4ǫuj+1 − 2uj + uj−1

h2+ 2

uj+1 − ujh

=xj + 1

2, j = 1, 2, . . . , N − 1. (10.56)

Quando ǫ = O(1), l’ errore di troncamento di (10.55) e O(h2), mentrequello di (10.56) e O(h).

Sia u(x) la soluzione esatta e u(x) quella approssimata. Poniamo

e(x) = u(x)− u(x), E = max0≤j≤N

|e(xj)|. (10.57)

La figura 10.6 mostra log10E in funzione di − log10 ǫ.La figura 10.7 mostra log10 |e| per alcuni valori di ǫ.

83

Figura 10.6: log10 E in funzione di − log10 ǫ.

84

Figura 10.7: log10 |e| in funzione di ǫ.

85

10.7.3 Conclusioni

Notiamo le seguenti caratteristiche dei due metodi.Differenze centrali:

• la soluzione calcolata oscilla quando ǫ→ 0;

• errore e soluzione sono limitati per N dispari e illimitati per N pari;

• E = O(h2) quando ǫ = O(1);

Differenze Upwind:

• la soluzione calcolata non oscilla quando ǫ→ 0;

• i comportamenti dell’ errore e della soluzione non dipendono dal fattoche N sia pari o dispari;

• E = O(h) quando ǫ = O(1);

86

Capitolo 11

Metodi Multigrid

11.1 Introduzione

Supponiamo di voler risolvere sul dominio Ω, costituito dal quadrato [0, 1]×[0, 1], il problema di Poisson in due dimensioni

−∆u = −∂2u(x, y)

∂x2− ∂2u(x, y)

∂y2= f(x, y), (11.1)

con condizioni al contorno di Dirichlet [QS06].Discretizziamo il dominio Ω con una griglia quadrata uniforme come quel-

la schizzata in figura 11.1, dividendo ogni lato in n parti, individuate dalleascisse e ordinate xi = ih e yj = jh, i, j = 0, . . . , n, h = 1/(n+ 1). Nel-la figura i nodi sono numerati in ordine lessicografico per colonne. Sianouij = u(ih, jh) e fij = f(ih, jh).

Approssimiamo l’ operatore di Laplace ∆ usando il classico stencil a 5punti. Otteniamo un sistema lineare di n2 equazioni in n2 incognite

4uij − ui−1,j − ui+1,j − ui−1,j−1 − ui,j+1 = h2fij , i, j = 1, . . . , n. (11.2)

Risolviamo il sistema (11.2), che in forma matriciale scriviamo A · x = b,con il metodo di rilssamento di Jacobi [QS06]

u(k+1) = −D−1(L+ U)u(k) +D−1b = u(k) +D−1r(k),

dove A = L+D+U e lo splitting di A in componente triangolare bassa (L),diagonale (D) e triangolare alta (U), r(k) = b− Au(k) e il residuo.

Se g(x) e una funzione periodica in I = [0, 1] (il nostro dominio e ilquadrato [0, 1]2, quindi la soluzione, come funzione della sola x o della sola y

87

Figura 11.1: Griglia quadrata uniforme e numerazione delle incognite.

non ci interessa al di fuori di I, possiamo assumere che sia periodica...) sottoopportune ipotesi possiamo scrivere la serie di Fourier [Ber94]

g(x) =a0

2+

+∞∑

m=1

(am cos(2πmx) + bm sin(2πmx)). (11.3)

I valori sm =√a2m + b2m rappresentano il cosiddetto spettro di potenza della

g(x): quando sm e “grande”, significa che la componente m-esima (dettam-esima armonica) della somma 11.3 e “grande” e viceversa. Notate che lacomponente m-esima ha frequenza m, quindi se sm e grande quando m≫ 1,diciamo che g(x) ha una “forte” componente in armoniche alte, una debolecomponente se sm e piccolo quando m e grande.

Poniamo n = 15, il sistema (11.2) ha N = 152 = 225 incognite. Conside-riamo l’ errore dello schema di Jacobi alla k-esima iterazione, e(k) = u−u(k).L’ iterazione di Jacobi calcola ogni nuova approssimazione u(k+1)(xi, yj) =

u(k+1)i,j usando la relazione

v(k+1)i,j =

v(k)i−1,j + v

(k)i+1,j + v

(k)i,j−1 + v

(k)i,j+1 + h2fi,j

4. (11.4)

Il valore in ogni nodo, P , viene aggiornato usando i valori calcolati al passoprecedente, associati al nodo stesso e ai nodi immediatamente vicini, ossia anord, sud, est e ovest di P . L’ aggiornamento non tiene conto dei nodi pi˘

88

20 40 60 80 100 120 1400

0.05

0.1

0.15

0.2Spettro di potenza errore alla iterazione 1

frequenza

pote

nza

20 40 60 80 100 120 1400

0.05

0.1

0.15

0.2alla iterazione 3

frequenza

pote

nza

20 40 60 80 100 120 1400

0.05

0.1

0.15

0.2alla iterazione 5

frequenza

pote

nza

20 40 60 80 100 120 1400

0.05

0.1

0.15

0.2alla iterazione 7

frequenza

pote

nza

Figura 11.2: Spettro di potenza dell’ errore in alcune iterazioni di Jacobi.

lontani da P , conseguentemente ci aspettiamo che, al crescere del numerodi iterazioni, le componenti in frequenza dell’ errore che coinvolgono puntilontani, ossia le frequenze pi˘ basse (m piccolo) , vengano ridotte meno dellecomponenti ad alta frequenza (m grande), che riguardano punti vicini nellagriglia. Questa supposizione e avvalorata dalla figura 11.2, che mostra lospettro di potenza, per le frequenze m = 1, . . . , 150 dell’ errore, alle iterazionik = 1, 3, 5, 7 (per m > 256/2 = 128 lo spettro dell’ errore ha grandezzatrascurabile). La figura mostra come per m < 50, ad esempio, le componentidello spettro dell’ errore, vengano ridotte molto meno di quelle per m > 50.

I metodi multigrid, che verrano descritti nei prossimi paragrafi, per smus-sare anche le componenti di frequenza bassa combinano diverse soluzioni ap-prossimate ottenute discretizzando il problema su pi˘ di una griglia, a diversilivelli di raffinamento.

Risolviamo il problema di Poisson, combinando il metodo di Jacobi con

89

20 40 60 80 100 120 1400

0.05

0.1

0.15

0.2Spettro di potenza errore alla iterazione 1

frequenza

pote

nza

20 40 60 80 100 120 1400

0.05

0.1

0.15

0.2alla iterazione 3

frequenza

pote

nza

20 40 60 80 100 120 1400

0.05

0.1

0.15

0.2alla iterazione 5

frequenza

pote

nza

20 40 60 80 100 120 1400

0.05

0.1

0.15

0.2alla iterazione 7

frequenza

pote

nza

Figura 11.3: Spettro di potenza dell’errore in alcune iterazioni del metodoMultigrid.

90

uno schema Multigrid. Quest’ ultimo (vedere pi˘ avanti per comprendere ilsignificato dei termini) consiste in una serie di cicli di tipo V, articolati suK = 4 livelli, ogni livello grossolano ha celle (quadrati) di lato doppio diquelle del livello immediatamente pi˘ raffinato. Ogni iterazione di livello kprevede un’ operazione di pre–smussamento (smoothing), una applicazionericorsiva del metodo e un passo di post–smussamento. La figura 11.3 mostrache tutte le componenti di frequenza m < 150 vengono egualmente ridottenelle prime 7 iterazioni.

11.2 Nested iterations

Sia dato un dominio poligonale Ω, ossia un dominio che ha una frontierafatta da segmenti. In generale un dominio “reale” avra una frontiera che euna curva arbitraria, ma noi supponiamo di poterla approssimare con unaspezzata fatta di segmenti di retta sufficientemente piccoli, senza perderemolto in accuratezza. Ricordiamo che il diametro di un insieme di punti, D,e:

diam(D) = supP,Q∈D

dist(P,Q).

Supponiamo di voler risolvere il sistema lineare

AK z = g, (11.5)

dove AK(pK × pK) e l’ approssimazione discreta di un dato operatore diffe-renziale, che caratterizza il problema da risolvere, tramite una griglia, τK , suΩ. Vogliamo utilizzare K sistemi ausiliari

Aixi = bi, i = 1, . . . , K − 1,

che sono discretizzazioni su una gerarchia di griglie dello stesso problemacontinuo, ognuna avente pi incognite, p1 < . . . < pK . Pensiamo ad esempioal dominio [0, 1]2 e alle griglie quadrate ottenute dividendo i lati in 2, 4, 8,16, ... parti.

L’ algoritmo Nested Iterations [KA03] (NI), e descritto nella tabella 11.1.L’ operatore S(A, x, b) rappresenta un solutore iterativo di sistemi lineari(rilassamento o smoothing). Le relazioni che useremo sono:

x(i+1)k = S(Ak, x

(i)k , bk), e

(i+1)k =

∣∣∣x− x(i+1)k

∣∣∣ ≃ µk

∣∣∣x− x(i)k

∣∣∣ = µke(i)k ,

dove µk e la costante di contrattivita del metodo iterativo, pensato comemetodo del punto fisso. L’ operatore Ikk−1 e un operatore di prolungamento

91

Tabella 11.1: Nested iterations.1 begin2 Choose mk ∈ N, k = 1, . . . , K;3 x := A−1

1 b1;4 for k := 1, K do5 x := Ikk−1x; x is prolonged 6 for i := 1, mk do7 x := S(Ak, x, b);8 endfor9 endfor

10 return x;11 end

(concetto spiegato meglio pi˘ avanti) che permette di interpolare una funzioneda una griglia grossolana ad una pi˘ raffinata. L’ algoritmo NI puo essereutilizzato per risolvere il sistema lineare (11.5), quando i K sistemi ausiliari,che sono discretizzazioni dello stesso problema continuo, hanno pi incognite,e risulta p1 < . . . < pK .

Se ogni metodo risolutivo S(Ak, ·, ·) ha costante di contrattivita µk, edesiste un valore µ

µk ≤ µ < 1, k = 1, . . . , K,

t.c. ponendo µk := µ il metodo NI converge, allora ha complessita ottima-le [KA03] ed e chiamato metodo multigrid. In tal caso si puo scegliere unvalore m, tale che ponendo mk = m, k = 0, . . . , K, l’ algoritmo NI converge.

Esempio 11.2.1 Un suggerimento esemplificativo. Supponiamo di appros-simare su ogni mesh τi con un metodo agli Elementi Finiti la soluzione,u(x, y), di un problema differenziale su Ω, ottenendo le approssimazioniui(x, y), i = 1, . . . , m. Supponiamo che fra i diametri, di, valgano le relazionidiami+1 = diami/2 e

u− ui = ei ≃ C diam2i , (11.6)

dove la costante C non dipende da diami. Allora abbiamo

ei ≃ 4 ei+1, ui − ui+1 ≃ 3 ei+1,

u ≃ ui+1 +ui − ui+1

3=

2 ui+1 − ui3

. (11.7)

Usando la relazione (11.7), possiamo calcolare un’ approssimazione dellasoluzione, pi˘ accurata di ui+1.

92

2

12

3

4

5 4

6

1 3

1

3 4

5

6

12

34

56

7

89101112

1314

15

16

1 7 2 8 3

9101112

4135

15 14

6

1

3 4

5

6

32

11 13 154 6

78

910 12 14

2221

2019

18 161723

242527

28 26

2930

3132

3334

3536

3738

39

4041

424344

4546

4748

4950

5152

5354

55

5657

5859

60

6162

63

64

1

716

117

218

819

2021222324252627

1228

1129

1030 9

313233343536

3713

38 4

39

14

45

6

44

42 41 40

4315

3

5

5

Figura 11.4: Una mesh (livello k = 0) e due suoi raffinamenti regolari (livellik = 1 e k = 2). La mesh pi˘ raffinata e stata ingrandita.

11.3 Griglie triangolari

Una triangolazione di Ω e una suddivisione fatta di triangoli, in cui nessunvertice di triangolo cade su un lato di un’ altro triangolo. I lati dei triangolivengono detti anche lati della triangolazione.

Supponiamo di avere un insieme di triangolazioni, τi, i = 1, . . . , K, suΩ. Diciamo che τi e la mesh a livello i. Assumiamo che queste mesh sianoraffinamenti regolari della mesh iniziale, τ1, ossia ogni triangolo di τi+1 siottiene unendo i punti medi dei lati di un triangolo di τi. Se quest’ ultima hadiametro diami, pi nodi e ti elementi, si ottiene una mesh che ha ti+1 = 4tielementi e diametro diami+1 = diami/2.

La figura 11.4 mostra un esempio di mesh triangolare e due raffinamentiregolari.

I triangoli di ogni mesh raffinata sono simili a quelli della mesh da cuiprovengono, quindi se τ1 e di Delaunay [KA03] tutte le altre lo sono.

93

Una famiglia di mesh τk, k = 1, . . ., sul dominio Ω, ognuna aventediametro diamk, e quasi uniforme [BRS02] se esiste ρ > 0 t.c.

maxdiam(BT ) : T ∈ τk ≥ ρ diamk diam(Ω), (11.8)

per ogni diamk ∈ (0, 1], dove BT e la pi˘ grande palla contenuta in T .Sia Vi lo spazio delle funzioni φ ∈ C0 lineari a tratti su τi, che sono nulle

su ∂Ω,φ = 0, su ∂Ω. (11.9)

Si noti che se n(τk) e l’ insieme dei nodi della k-esima mesh,

n(τi) ⊃ n(τj)⇒ Vi ⊃ Vj , 1 ≤ j ≤ i ≤ K. (11.10)

Nota: il fatto che ogni v ∈ Vk sia nulla su ∂Ω non e limitante, se il problemadi risolvere e di Dirichlet. Infatti, se uh → u, soluzione del problema diDirichlet omogeneo, quando u e la soluzione del problema di Dirichlet u = gsu ∂Ω, risulta

uh = uh +∑

xi∈∂Ω

g(xi)φi(x)→ u.

Quando il problema da risolvere non ha condizioni al contorno solamente diDirichlet, puo essere necessario cambiare lo spazio delle funzioni approssi-manti.

Definizione 11.3.1 Sia Nk il numero dei nodi, Pψk(i), i = 1, . . . , Nk, internia τk, ψk : 1, . . . , Nk → 1, . . . , pk. Definiamo il prodotto interno su Vk,dipendente dalla mesh (mesh–dependent inner product), come

(v, w)k := diam2k

Nk∑

i=1

v(Pψk(i))w(Pψk(i)).

La condizione (11.9) semplifica le dimostrazioni, ma i ragionamenti chefaremo possono essere estesi a spazi di funzioni che non soddisfano questacondizione. Per sviluppare alcuni esempi, definiamo anche gli spazi Wi dellefunzioni φ ∈ C0 lineari a tratti su τi. Notare che Wi ⊃ Vi. Il prodotto internosu Wk, dipendente dalla mesh, e naturalmente

(v, w)k := diam2k

pk∑

i=1

v(Pi)w(Pi).

94

Definizione 11.3.2 Sia

Ikk−1 : Vk−1 → Vk,

l’ operatore di prolungamento (coarse–to–fine operator), detto anche inter-polazione (interpolation), che facciamo coincidere con l’ operatore di immer-sione:

Ikk−1v = v, ∀v ∈ Vk−1.

Definizione 11.3.3 Sia

Ik−1k : Vk → Vk−1,

l’ operatore di restrizione o pesatura residuale (residual weighting, fine–to–coarse operator), che definiamo come il trasposto di Ikk−1, rispetto ai prodottiinterni (·, ·)k−1, (·, ·)k, ossia

(Ik−1k w, v)k−1 = (w, Ikk−1v)k = (w, v)k, ∀v ∈ Vk−1, w ∈ Vk.

Nota: quando la restrizione e il trasposto del prolungamento, e piu agevo-le provare che un metodo multigrid converge. Tuttavia sono stati propostialgoritmi in cui questo non si verifica: in questi casi non vengono date dimo-strazioni di convergenza e ci si limita ad una “validazione” dell’ algoritmo suproblemi test.

11.4 Multigrid

Lo schema multigrid usa due strategie principali [Wes92]:

• smussamento su τi: vengono smussate le componenti ad alta frequenzadell’ errore;

• correzione dell’ errore su τi−1, usando uno schema iterativo sviluppa-to su quest’ ultimo spazio: vengono smussate le componenti a bassafrequenza dell’ errore.

Descriviamo in dettaglio un algoritmo multigrid “generale”, su K griglietriangolari, τ1 ⊂ . . . ⊂ τK , implementato dalla funzione MG(k, z, g), cheviene definita ricorsivamente.

Le considerazioni che faremo valgono anche, con modifiche opportune,per griglie quadrate e rettangolari.

Innanzitutto MG(1, z0, g) := A−11 g, ossia al livello pi˘ grossolano, si

risolve esattamente il sistema lineare.

95

32

5 11 13 154 6

78

910 12 14

2221

2019

18 161723

242527

28 26

2930

3132

3334

3536

3738

39

4041

424344

4546

4748

4950

5152

5354

55

5657

5859

60

6162

63

64

1

716

117

218

819

2021222324252627

1228

1129

1030 9

313233343536

3713

38 4

39

14

45

6

44

42 41 40

4315

3

5 5

Figura 11.5: La mesh a livello 3 e i supporti di φ(13)3 e φ

(10)2 .

Oltre a due valori m1, m2 ∈ N, occorre fissare un parametro p ∈ N,che quasi sempre assume i valori 1 o 2. Se poniamo p = 1, otteniamo uncosiddetto ciclo V, se p = 2, otteniamo un ciclo W.

Definiamo l’ operatore di rilassamento o smussamento,

S(A, z, g) := z − 1

Λkr(A, z, g),

dove r(A, x, g) = g − Ax, e il residuo del sistema lineare Ax = g, e Λk eun opportuno parametro di rilassamento. Per ogni livello, k, va fissato unopportuno parametro, che e un maggiorante del raggio spettrale di Ak, ρ(Ak),ossia ad ogni livello k, valgano le relazioni ρ(Ak) ≤ Λk ≤ C/diam2

k, per un’opportuna costante C.

Per k > 1, MG(k, z0, g), detto iterazione di k-esimo livello, si articola intre passi:

1. Passo di pre–smussamento. Per 1 ≤ l ≤ m1, poniamo

zl := S(Ak, zl−1, g).

2. Passo di correzione dell’ errore. Sia

g := Ik−1k r(Ak, zm1 , g), q0 = 0;

96

Per 1 ≤ i ≤ p, siaqi := MG(k − 1, qi−1, g);

Poniamozm1+1 := zm1 + Ikk−1qp.

3. Passo di post–smussamento. Per m1 + 2 ≤ l ≤ m1 +m2 + 1,

zl := S(Ak, zl−1, g).

Notare che se m2 = 0, nessun passo di post–smussamento viene effet-tuato.

Il risultato e:MG(k, z0, g) := zm1+m2+1.

Lo schema multigrid completo per risolvere Au = f , richiede di risolveres volte l’ iterazione di k-esimo livello, che e riportata nella tabella 11.2.

La tabella 11.3 mostra il tracing dell’ esecuzione di MG(4, z3, f4), po-nendo m1 = 1, m2 = 2, p = 2, s = 3, z3 e un’ approssimazione di u sullagriglia τ3. La figura 11.6 schizza una parte del procedimento. Notare che iprolungamenti tratteggiati estendono soluzioni approssimate.

11.5 Convergenza e Costo

Sia Ω un poligono convesso. Poniamo

a(u, v) =

∫

Ω

(α∇u · ∇v + βu v) dx,

dove α, β sono funzioni smooth t.c. per qualche α0, α1, β1 ∈ R+, α0 ≤α(x) ≤ α1, 0 ≤ β(x) ≤ β1, ∀x ∈ Ω.

Sia dato il problema di Dirichlet di trovare u ∈ V :=

H1(Ω) t.c.

a(u, v) = (f, v), ∀v ∈ V.

Supponiamo di risolvere questo problema usando il multigrid, con K > 1livelli, ottenuti tramite raffinamenti regolari della mesh a livello 1, che ha t1triangoli.

Ricordiamo che ‖v‖E =√a(v, v). Valgono i risultati [BRS02]:

Teorema 11.5.1 Se m = maxm1, m2 e grande abbastanza, l’ iterazionedi k-esimo livello di un W-ciclo e una contrazione con costante contrattivaindipendente da k.

97

Tabella 11.2: Algoritmo multigrid.1 MG(k, z, g)2 begin3 if k = 14 then5 z := A−1

1 g;6 else7 z := Ikk−1z; z is prolonged 8 for s := 1, s do9 for l := 1, m1 do

10 z := S(Ak, z, g);11 endfor12 g := Ik−1

k r(Ak, z, g); r is restricted 13 q := 0;14 for i := 1, p do15 q := MG(k − 1, q, g);16 endfor17 z := z + Ikk−1q; q is prolonged 18 for l := m1 + 2, m1 +m2 + 1 do19 z := S(Ak, z, g);20 endfor21 endfor22 endif23 end

98

Tabella 11.3: Tracing di MG(4, z3, f4), s = 3.

99

p=2 $m 1$ = 1$\calP$ = prolongation

$\calR$ = restriction

$\xi$ = compute the exact solution

$\rho$ = relaxation step

$\sbar=3$$m 2$ = 2

1 2 3 4 5

16 18

17

6 7 8 9 10

to

k=4

k=1

k=2

k=3

coarserlevels

14 15 23 24 25 26

11 12 13 19 20 21 22

27

28 34 35 36 37

29

30 31 32 3340

39

Figura 11.6: Rappresentazione schematica di MG(4, z3, f4), s = 3, fino al passo 40 del tracing in tabella 11.3.

100

Teorema 11.5.2 Per qualsiasi m > 0 l’ iterazione di k-esimo livello di unV-ciclo e una contrazione con costante contrattiva indipendente da k.

Per il multigrid completo vale il teorema [BRS02]:

Teorema 11.5.3 Se (a) l’ iterazione di k-esimo livello e una contrazionecon costante contrattiva γ indipendente da k e (b) il numero di esecuzionidella iterazione di k-esimo livello, s, e grande abbastanza, allora esiste unacostante C > 0 t.c.

‖uk − uk‖E ≤ Chk|u|H2(Ω). (11.11)

dove uk e la proiezione della soluzione esatta sui nodi di τk e uk e l’ appros-simazione ottenuta con l’ algoritmo multigrid completo.

Inoltre:

Teorema 11.5.4 Il costo computazionale dell’ algoritmo multigrid completoa k livelli e O(nk), dove nk = 2t14

k.

Si noti che nk e una maggiorazione del numero di nodi in τk, essendo t1il numero di triangoli nella mesh iniziale.

Notare che la stima (11.11) e in norma dell’ energia, che riguarda essen-zialmente le derivate delle funzioni. Per la stima dell’ errore in norma L2,nelle stesse ipotesi del teorema precedente, si puo mostrare che

‖uk − uk‖L2 ≤ C ′h2k|u|H2(Ω).

11.6 Osservazioni interessanti

Ricordiamo i seguenti risultati, tratti da [BR02].Supponiamo di voler risolvere il sistema Ax = b e di aver calcolato un’

approssimazione della soluzione, v.Una misura precisa dell’ efficienza dello smussamento fornito da un me-

todo di rilassamento si puo ottenere quando la griglia di discretizzazione euniforme e l’ iterazione (detta anche relaxation sweep) percorre i punti inun dato ordine ammissibile. Ignorando i contorni ed eventuali variazioni deicoefficienti, si puo calcolare il fattore di convergenza ad ogni iterazione perogni componente di Fourier dell’ errore [Bra77, TOS00]. Definiamo allora ilfattore di smussamento, µ, come il piu grande (ovverosia il peggiore) fatto-re di convergenza per iterazione fra tutte le componenti che oscillano tropporapidamente per essere visibili sulla griglia grossolana, quindi debbono essere

101

ridotte dal rilassamento. Ad esempio consideriamo l’ operatore di Laplacediscretizzato con il classico schema delle differenze contrali a 5 punti, usandoraffinamenti regolari, che riducono i diametri di un fattore 2. Allora il metododi Gauss–Seidel (GS) che scandisce i punti di griglia in ordine lessicografico,riga per riga, fornisce µ = 0.5; scandendo in ordine red/black (detto anche ascacchiera), si ha µ = 0.25. Percio due iterazioni GS con red/black riduconodi un fattore 0.252 = 0.0625 tutte le componenti di errore invisibili per lagriglia grossolana.

Supponiamo che in un ciclo multigrid vengano eseguite m iterazioni dirilassamento (m = m1+m2, m1 e il numero di iterazioni di pre–smussamento,m2 quello delle iterazioni di post–smussamento). Ci si aspetta che questeiterazioni riducano le componenti di errore visibili nella griglia attiva, manon visibile da quelle piu grossolane, di un fattore µm, se µ e il fattore dismussamento. Quando tutte le griglie vengono attraversate, ci aspettiamoche il ciclo riduca tutte le componenti di errore almeno di un fattore µm.Teoria ed esperienza confermano che per problemi scalari ellittici regolari epiccoli valori di m, questa efficienza viene raggiunta, purche le condizionial contorno vengano rilassate correttamente e si usino appropriati operatoridi restrizione e prolungamento. Percio il fattore µ e un eccellente stimatoredella accuratezza che il metodo multigrid dovrebbe avere.

Le analisi teoriche che stimano il fattore di convergenza dei metodi mul-tigrid si applicano solo a problemi relativemente semplici e ad algortimi ar-tificiali, spesso diversi dall’ algoritmo ottimale nelle applicazioni. Le unichestime teoriche che sono realistiche e abbastanza precise, sono basate su tec-niche di analisi di Fourier locali, chiamate local mode analysis (LMA). L’analisi LMA piu semplice ed usata e quella appena descritta, che si riduce adusare il fattore di smussamento. Una stima piu elaborata si ottiene usandoun’ analisi di Fourier su piu livelli (spesso due soli), analizzando cosı sia ilrilassamento che i trasferimenti di informazioni tra griglie. Risultati detta-gliati su come calcolare questi fattori di convergenza e codici per valutarlisi trovano in [Wei01]. Nonostante l’ analisi di Fourier sia valida solo perequazioni a coefficienti costanti in domini infiniti o rettangolari, in praticale stime valgono per una piu ampia classe di problemi, percio vengono usateper sviluppare algoritmi e verificare codici, perfino per complessi probleminon lineari. Per problemi ellittici generali con coefficienti regolari a trattiin domini generali discretizzati con griglie uniformi, viene provato rigorosa-mente in [Bra89, Bra94], che per ampiezza di mesh tendente a zero, le stimedi convergenza fornite da LMA sono accurate, purche il procedimento dimultigrid sia sostanziato da algoritmi appropriati per l’ elaborazione vicino esulla frontiera. Teoria e pratica mostrano che l’ efficienza del multigrid vieneaumentata introducendo speciali passi di rilassamento in ogni zona in cui il

102

residuo e piu grande della media: vedi la adaptive relaxation rule in [BR02,Cap. 3]

In un generico algoritmo multigrid, i cicli possono essere applicati a qual-siasi approssimazione iniziale data sulla griglia piu fine. Nel Full Multigrid(FMG), ogni prima approssimazione e ottenuta usando il multigrid stesso.Gli algoritmi FMG sono piu robusti in accuratezza dei cicli multigrid di altrotipo.

Supponiamo di discretizzare una PDE su una data griglia. Definiamounita di lavoro minimale il numero di operazioni necessarie, sulla data gri-glia, per calcolare la discretizzazione piu semplice, oppure per eseguire unaiterazione del piu semplice schema di rilassamento. Un FMG accuratamentesviluppato, che usa un ciclo V o un W ad ogni livello di raffinamento ed unpaio di efficienti iterazioni di rilassamento ad ogni livello,

costa meno di 10 unita di lavoro minimali e risolve una PDE discretiz-zata con una accuratezza almeno O(h2) su una griglia di dimensioneh.

Questa efficienza ideale e stata chiamata textbook multigrid efficiency (TME).Gli algoritmi multigrid che si rifanno a griglie di discretizzazione vengono

chiamati anche metodi di multigrid geometrico (geometric multigrid). Sonospesso di difficile realizzazione e funzionano meglio con grigie uniformi.

Gli algoritmi di multigrid algebrico (algebraic multigrid, AMG) risolvonosistemi di equazioni lineari senza usare esplicitamente la geometria delle gri-glie. Sono in genere meno efficienti dei multigrid geometrici, ma di piu facileimplementazione e funzionano egualmente bene sia quando ci sono griglie benstrutturate che non. In tali metodi si parla di prolungamento e restrizione travettori, individuando sottoinsiemi, S, dell’ insieme di partenza (visto comeinsieme raffinato) delle incognite. La scelta di S si puo basare sui cosiddettischemi di compatible relaxations [BR02]. Il fattore di convergenza di questischemi e un buon stimatore della efficenza del metodo AMG che ne deriva.

11.7 Multigrid adattivo

Supponiamo che la griglia piu fine, sulla quale vogliamo risolvere il problema,τK , K > 1, in alcune zone non abbia il livello di raffinamento necessario perdescrivere il fenomeno. Vorremmo risolvere il problema su una mesh a livelloK + L, L > 0, per trattare le zone che richiedono il massimo raffinamento,ma nel resto del dominio vogliamo limitarci a raffinare a livello K.

Per semplicita supponiamo vi sia una sola zona, Z, dove bisogna raffinareoltre il livello K.

103

K

Z

domain

K−3

K−2

K−1

K+1

Figura 11.7: Esempio di raffinamento adattivo.

Applichiamo quindi il multigrid a livello K e consideriamo una griglia“fantasma” a livello K+L, sulla quale solo all’ interno della zona Z si raffinaa livello K + L, mentre nel resto del dominio si usano i prolungamenti dallivello K delle funzioni coinvolte nel multigrid.

La figura 11.7 mostra un dominio raffinato fino al livelloK e al suo internouna zona Z, raffinata a livello K + 1.

L’ algoritmo multigrid adattivo e equivalente all’ algoritmo della tabel-la 11.2, nella quale quando k > K l’ operatore Ak viene sostituito con l’ ope-ratore Ak che approssima nel dominio l’ operatore A a livello di raffinamentoK, tranne nella zona Z, dove l’ approssimazione e a livello k > K.

11.7.1 Interfacce fra sottogriglie

Supponiamo di risolvere il problema di Poisson −∆u = g su una grigliasuddivisa in piu zone, con raffinamenti diversi. Usiamo il solito schema a 5punti per approssimare l’ operatore di Laplace.

La figura 11.8 mostra una situazione in cui il dominio Ω e diviso in duezone. La prima, Ωc, in grigio nella figura, e raffinata grossolanamente tra-mite la sotto–mesh a elementi quadrati µc (c=coarse). La seconda, Ωf , eraffinata in maniera piu fine dalla sotto–griglia µf (f=fine). Supponiamo chela discretizzazione sia centrata sulle celle (cell–centered), ossia le soluzioniapprossimate vengano calcolate sui centri delle celle.

Bisogna stare attenti a come si trattano i confini tra sotto–mesh aven-ti diversi livelli di raffinamento (li chiameremo interfacce), come illustra ilragionamento in [MC96, pag. 5], che ora riformuliamo. Supponiamo di risol-vere il problema −∆cu

(c) = g(c), dove ∆c e la discretizzazione dell’ operatoredi Laplace sulla griglia grossolana, τc. Il vettore g(c) contiene i valori deltermine forzante nei nodi di τc. Risolviamo ora il problema −∆f u

(f) = g(f),dove ∆f e la discretizzazione dell’ operatore di Laplace sulla griglia fine. Per

104

vy,(i,j−1/2)

vy,(i,j+1/2)

j+1

j

j-1

vx,t

vx,b

vx,(i+1/2,j)

i i+1

i-1/4

i-1

j-1/4

j+1/4

i-3/4

A

D B

C

Figura 11.8: Due sotto–griglie a maglia quadrata. I cerchi rappresentanovalori calcolati, i quadrati valori interpolati, le “punte” valori interpolati.

105

jj+

1

j-1/4

j-1

i-1/4 i+1/4

j-3/4

i-1i i+1

vy,l vy,r

BA

C D

vx

,(i+

1/2,j

)

vx

,(i−

1/2,j

)

vy

,(i,

j+

1/2)

Figura 11.9: Analoga alla figura precedente.

106

calcolare u(f) nei nodi di τf adiacenti all’ interfaccia, servono anche valori sunodi di τc adiacenti all’ interfaccia. Supponiamo di calcolare questi valori perinterpolazione lineare ed estrapolazione lineare sui valori in τc, eventualmenteusando valori sul contorno fisico, ossia il contorno di Ω. Si scopre che questoprocedimento fornisce in generale un risultato della stessa accuratezza che siottiene usando solo la griglia piu grossolana! Infatti la soluzione u e continuasull’ interfaccia ∂Ωi = ∂Ωc ∩ ∂Ωf , ma la derivata normale ∂u/∂n non lo e.La griglia grossolana fornisce informazione a quella raffinata come se fosseuna parte di frontiera di Dirichlet, la griglia raffinata usa questa informazio-ne, ma cambia il valore della derivata normale. Se non si restituisce questavariazione della derivata alla griglia grossolana, il metodo non guadagna inaccuratezza. In altre parole, la griglia grossolana non vede gli effetti dellagriglia fine: sull’ interfaccia, dobbiamo imporre la continuita della soluzio-ne, ma anche quella della derivata normale. Queste sono le condizioni dicontinuita per l’ ellitticita (elliptic matching condition) sulle interfacce.

Per aumentare l’ accuratezza, dobbiamo costruire una discretizzazionecomposita dell’ operatore di Laplace, definita in modi differenti nelle varieregioni.

Supponiamo di avere una gerarchia di mesh uniformi, τk di diametro hk,che coprono l’ intero dominio Ω per ogni τk una gerarchia di sotto-griglie,µk,l, µk,l ⊂ τl, che ricoprono parti di τk ⊂ Ω. In ogni livello le sotto-grigliedebbono soddisfare il criterio di annidamento proprio (proper nesting), ossianessuna cella a livello l+1 rappresenta una regione adiacente a una occupatada una cella a livello l−1. In altre parole, nessuna cella di livello l e ricopertasolo parzialmente, non totalmente, da celle di livello l − 1.

Sia R l’ operatore di restrizione che mappa gli elementi geometrici di unlivello l sul livello piu grossolano l − 1.

Per semplicita supponiamo che il rapporto di raffinamento r = hl/hl+1

sia r = 2, costante fra tutti i livelli adiacenti.Poiche assumiamo che le soluzioni su mesh fini siano piu accurate, in ogni

livello distinguiamo tra regioni valide e regioni non valide. Una regione validaa livello l, µk,l,V , non e coperta da celle di griglie piu fini: µk,l,V = µl\R(µl+1).Indichiamo con µ∗

k,l i lati delle celle di livello l e con µ∗k,l,V i lati delle celle di

livello l non coperti da lati di livello l + 1. Si noti che l’ interfaccia ∂µ∗k,l+1

tra i livelli l e l+ 1 e valida a livello l+ 1, ma non a livello l. L’ operatore direstrizione viene esteso ai lati: R(µ∗

k,l+1) e l’ insieme dei lati a livello l copertida lati di livello l+1. Una variabile composita e definita sull’ unione di tuttele regioni valide, a tutti i livelli. Una variabile di livello l, u(l), e definitasu tutto µk,l, mentre una variabile composita, u(c), e definita sull’ unionedelle regioni valide su tutti i livelli. Definiamo anche campi vettoriali, chevengono perø associati ai punti medi dei lati, ossia utilizziamo le cosiddette

107

griglie sfalsate (staggered grids). Analogamente a prima definiamo vettori dilivello e vettori compositi.

Supponiamo ora di lavorare all’ interno della griglia di livello k, per cuiomettiamo il livello principale e scriviamo µl, µl,V , al posto di µk,l, µk,l,V etc.

Per approssimare l’ operatore di Laplace, nelle zone non adiacenti ainterfacce, usiamo il classico stencil a 5 punti.

Nelle zone adiacenti a una interfaccia fra i livelli l e l + 1, usiamo ladivergenza del gradiente. Ricordiamo che ∆u = ∇∇u = ∇v = ∇(vx, vy),posto (vx, vy) := ∇u. Associamo ai punti medi dei lati i valori di vx e vy.

Sia u = u(c) l’ approssimazione composita della soluzione esatta, u; sia u(l)

l’ approssimazione di livello l. Definiamo la discretizzazione della divergenzaa livello l, Dl, con la formula tipo MAC projection [Min96] (vedi figura 11.8)

(∇l v)ij = (Dl(vx, vy))ij =

vx,(i+1/2,j) − vx,(i−1/2,j)

hl+vy,(i,j+1/2) − vy,(i,j−1/2)

hl, (11.12)

dove

vx,(i+1/2,j) =ui+1,j − ui,j

hl, vx,(i−1/2,j) =

ui,j − ui−1,j

hl,

vy,(i,j+1/2) =ui,j+1 − ui,j

hl, vy,(i,j−1/2) =

ui,j − ui,j−1

hl.

La discretizzazione della divergenza a livello l viene definita ignorando i livellipiu fini e usando la discretizzazione standard (11.12). Definiamo anche unadivergenza composita, D(l,l+1), che opera tra due livelli, l e l+1. Assumiamo

che il campo vettoriale v(l) = (v(l)x , v

(l)y ) possa essere esteso a tutto µ∗

l , l’insieme di tutti i lati di celle in µl, inclusa l’ interfaccia ∂µ∗

l+1 = µ∗l+1 ∩

µ∗l . La divergenza composita su µl puo essere calcolata tramite Dl, con una

correzione dovuta al livello piu fine l+1. L’ operatore composito di divergenzapuo essere definito come

(D(l,l+1)v)ij = (Dlv(l))ij + (Dl,Rδv(l+1))ij ,

dove

δv(l+1) = 〈v(l+1)〉 − v(l) = I ll+1v(l+1) − v(l), su R(∂µ∗

l+1).

L’ operatore 〈·〉 rappresenta il valor medio. Per calcolare D(l,l+1) in modoefficiente, abbiamo introdotto un registro di gradiente, δv(l+1), definito suR(∂µ∗

l+1), che memorizza la differenza sull’ interfaccia tra i valori della re-strizione di v(l+1) a livello l e i valori a livello l. Il registro δv(l+1) appartiene

108

al livello l+1 perche rappresenta informazione su ∂µ∗l+1, perø ha spaziatura e

indici del livello l. Abbiamo introdotto anche la divergenza di riflusso, Dl,R,che e lo stencil Dl, applicato ai vettori sull’ interfaccia con il livello l + 1. Alivello l, Dl,R puo essere definito come

(Dl,Rδv(l+1))ij =1

hl

∑

p∈Q

s · (δv(l+1))p,

dove Q e l’ insieme dei lati che sono anche interfacce con il livello l + 1.Il fattore s vale s := +1 se p e un lato alto della cella (i, j) e s := −1 sep e un lato basso. Si noti che Dl,R influenza l’ insieme di celle di livello limmediatamente adiacenti all’ interfaccia.

Usando un volume di controllo possiamo esplicitare D(l,l+1). Ad esempionel caso della figura 11.8, fissiamo il nostro volume di controllo attorno allacella, della griglia grossolana che ha centro nel nodo (i, j). Abbiamo:

(D(l,l+1)v)ij =v

(l)x,(i+1/2,j) − 〈v

(l+1)x,(i−1/2,j)〉

hl+v

(l)y,(i,j+1/2) − v

(l)y,(i,j−1/2)

hl, (11.13)

dove

〈v(l+1)x,(i−1/2,j)〉 =

vx,t + vx,b2

,

vx,t =u

(l,I)i−1/4,j+1/4 − u

(l+1)i−3/4,j+1/4

hl+1

, vx,b =u

(l,I)i−1/4,j−1/4 − u

(l+1)i−3/4,j−1/4

hl+1

,

i valori u(l,I)i−1/4,j+1/4, u

(l,I)i−1/4,j−1/4, si ottengono interpolando sull’ interfaccia

∂µl+1, con un metodo usato ad esempio in [Min96], ossia

• prima si calcolano i valori di u nei punti (i, j + 1/4), (i, j − 1/4), perinterpolazione quadratica in direzione parallela all’ interfaccia, sui nodidella griglia grossolana.

• Poi si effettua un’ altra interpolazione quadratica in direzione normaleall’ interfaccia per stimare u nei punti (i−1/4, j+1/4), (i−1/4, j−1/4).

Usare questi valori interpolati al bordo, equivale a imporre la continuita del-la derivata normale all’ interfaccia. Per quanto riguarda i nodi adiacentiall’ interfaccia che stanno nella sotto–mesh piu fine, applichiamo lo stencilstandard, usando i valori interpolati che ci sono serviti per calcolare le ap-prossimazioni nei nodi adiacenti all’ interfaccia, che stanno nella sotto–meshpiu grossolana Ad esempio, riferendosi alla figura 11.8, per calcolare il Lapla-ciano nel nodo (i− 3/4, j− 1/4), applichiamo lo stencil standard ai valori in

109

(i−3/4, j+1/4), (i−3/4, j−3/4), (i−5/4, j−1/4), che sono nella sotto–meshpiu fina e quindi disponibili, assieme al valore in (i − 1/4, j − 1/4) (puntoD nella figura), che e stato interpolato per calcolare il Laplaciano nel nodo(i, j) Per quanto riguarda i centri di celle che stanno nella zona raffinata,adiacenti all’ interfaccia, una volta calcolati i valori che servono nella rela-zione (11.13), possiamo utilizzarli per valutare il Laplaciano con lo stencilstandard anche nei punti suddetti.

Ancora un esempio: nel caso della figura 11.9, fissiamo il nostro volumedi controllo attorno alla cella della griglia grossolana che ha centro nel nodo(i, j). Abbiamo:

(D(l,l+1)v)ij =v

(l)y,(i,j+1/2) − 〈v

(l+1)y,(i,j−1/2)〉

hl+v

(l)x,(i+1/2,j) − v

(l)x,(i−1/2,j)

hl,

dove

〈v(l+1)y,(i,j−1/2)〉 =

vy,l + vy,r2

,

vy,r =u

(l,I)i+1/4,j−1/4 − u

(l+1)i+1/4,j−3/4

hl+1

, vy,l =u

(l,I)i−1/4,j−1/4 − u

(l+1)i−1/4,j−3/4

hl+1

,

i valori u(l,I)i+1/4,j−1/4, u

(l,I)i−1/4,j−1/4, si ottengono in maniera analoga a prima,

interpolando sull’ interfaccia ∂µl+1.Siccome il Laplaciano e un operatore che coinvolge derivate seconde, la sua

approssimazione alle differenze centrali richiede la divisione per h2. Se unaquantita interpolata che viene utilizzata nell’ approssimazione sull’ interfacciae O(hp), l’ errore sull’ interfaccia e O(hp−2). L’ interpolazione quadraticaha errore di ordine p = 3, che implica un’ accuratezza O(h) sull’ interfaccia.Nonostante l’ accuratezza nelle altre parti del dominio sia piu elevata, O(h2),dato che l’ interfaccia e un insieme di co-dimensione 1, possiamo accontentarcidi un O(h) su di essa e avere comunque globalmente un O(h2).

Dato che vogliamo usare l’ interpolazione quadratica ovunque sia possibi-le per collegare griglie fini e grossolane, dobbiamo usare stencils differenti incasi speciali come interfacce ad angolo (vedi figura 11.10). Se uno dei puntidello stencil cade in una zona raffinata grossolanamente, lo shiftiamo in mo-do da usare solo celle in µl−1\R(µl). Se non esiste uno stencil “grossolano”parallelo all’ interfaccia, abbassiamo l’ ordine di interpolazione e usiamo lecelle della mesh grossolana che abbiamo. Usando punti di mesh grossolanee fini nelle interpolazioni e usando gli stessi gradienti per entrambi i tipi digriglie, otteniamo il collegamento necessario per risolvere il problema delleinterfacce. Nel caso della figura 11.10, fissiamo il nostro volume di control-lo attorno alla cella della griglia grossolana che ha centro nel nodo (i, j) e

110

vy,(i,j+1/2)

j+1

j

j-1

vx,t

vx,b

vx,(i+1/2,j)

i i+1

i-1/4

i-1

j-1/4

j+1/4

i-3/4

A

D B

C

E F

G

vy

,r

vy

,l

i+1/4

Figura 11.10: Due sotto–griglie a maglia quadrata, con interfaccia angolare.

111

otteniamo

(D(l,l+1)v)ij =v

(l)x,(i+1/2,j) − 〈v

(l+1)x,(i−1/2,j)〉

hl+v

(l)y,(i,j+1/2) − v

(l)y,(i,j−1/2)

hl,

dove

〈v(l+1)x,(i−1/2,j)〉 =

vt + vb2

,

vx,t =u

(l,I)i−1/4,j+1/4 − u

(l+1)i−3/4,j+1/4

hl+1

, vx,b =u

(l,I)i−1/4,j−1/4 − u

(l+1)i−3/4,j−1/4

hl+1

,

i valori u(l,I)i−1/4,j+1/4, u

(l,I)i−1/4,j−1/4, si ottengono in maniera analoga a prima,

interpolando sull’ interfaccia ∂µl+1.La discretizzazione del gradiente, G, e associata ai centri dei lati delle

celle. Il gradiente composito G(l,l+1), e definito su tutti i lati validi nel dominiomultilivello. Sui lati che non sono di interfaccia

(G(l,l+1),xu)i+1/2,j =ui+1,j − ui,j

hl, (G(l,l+1),yu)i,j+1/2 =

ui,j+1 − ui,jhl

.

(11.14)Per calcolare G(l,l+1) su un’ interfaccia, interpoliamo i valori di u sia sul livellogrossolano, che su quello fine.

L’ operatore Gl viene definito estendendo a tutti i lati in µ∗k,l,V l’ ope-

ratore G(l,l+1), che e definito solo su ∂µ∗k,l,V . Lontano da ∂µ∗

k,l,V , usiamo lostencil (11.14) per griglia interna, mentre sulle interfacce usiamo la proce-dura di interpolazione descritta, per calcolare i valori nelle celle virtuali cheservono.

Sia ancora u la discretizzazione composita della soluzione. La discretizza-zione del Laplaciano composito e definita come la divergenza del gradiente:

L(l,l+1)u = D(l,l+1)G(l,l+1)u, (11.15)

Nelle zone interne, la relazione si riduce al classico operatore a 5 punti,

Llu(l) = DlGlu(l). (11.16)

Sulle interfacce, l’ interpolazione fornisce i valori “fantasma” necessari. Sullato grossolano di un’ interfaccia, la relazione (11.15) diventa

L(l,l+1)u = Llu(l) +Dl,R(δGl+1u(l+1)), (11.17)

δGl+1u(l+1) = 〈Gl+1u

(l+1)〉 − Glu(l). (11.18)

112

A

C B

ii-1 i+1

j

j-1

j+1

i-3/4 i-1/4

j-1/4

j+1/4

µl+1 µl

vx,b

vx,t

Figura 11.11: Due sotto–griglie a maglia quadrata, in una meshvertex-centered. La cella che ha centro nel nodo (i, j) e stata evidenziata.

Mesh vertex–centered

Consideriamo ora il caso in cui i valori del potenziale non siano associati aicentri delle celle, ma ai nodi della mesh, ossia abbiamo una rappresentazionedetta vertex–centered.

La figura 11.11 e analoga alla figura 11.8, ma schematizza una situazionevertex–centered. Supponiamo che un nodo P = (xi, yj), giaccia su un’ inter-faccia che nella mesh di livello k, τk, divide una sottomesh di livello l, µk,l, dauna piu raffinata di livello l + 1, µk,l+1. Se cio accade, P puo essere pensatosia come nodo di una, come dell’ altra sottomesh. I valori associati a questonodo possono essere quindi di pertinenza delle due mesh. Ad esempio, ui,j

113

puo indicare sia il potenziale nel nodo P , calcolato pensando P ∈ µk,l, che il

potenziale calcolato per P ∈ µk,l+1. Il simbolo u(k)i,j rappresenta la variabile

di livello k, che riguarda l’intera mesh τk, quindi per indicare la variabilea livello k, valuata in P ∈ µk,l, useremo il simbolo u

(k,l)i,j ; quando invece si

assuma P ∈ µk,l+1, scriviamo u(k,l+1)i,j . I valori della variabile composita u

restano sempre indipendenti dalle sottomesh.Per quanto riguarda i nodi di interfaccia pensati nella mesh grossolana,

P ∈ µk,l, applichiamo una strategia analoga a quella usata per valori cell-centered. Con riferimento alla figura 11.11, fissiamo il nostro volume dicontrollo attorno alla cella della griglia grossolana che ha centro nel nodo(i, j). Poniamo:

(D(l,l+1)v)i,j =u

(k,l)i+1,j − 2u

(k,l)i,j − u(k,l+1)

i−1,j

h2l

+u

(k,l+1)i,j+1 − 2u

(k,l)i,j − u(k,l+1)

i,j−1

h2l

. (11.19)

Notare che, laddove possibile, usiamo i valori della sottomesh fine.Per quanto riguarda i nodi di interfaccia pensati nella mesh raffinata,

P ∈ µk,l+1, distinguiamo due tipi di nodi:

• nodi, chiamiamoli di tipo G, che stanno anche nella sotto-griglia gros-solana (es. il nodo (i, j) nella figura 11.11);

• nodi di tipo F , che appartengono solo alla mesh fine (es. il nodo (i, j−1/2) nella figura 11.11);

Per i nodi di tipo G, prendiamo ad esempio G = (i, j), poniamo:

(D(l,l+1)v)i,j =u

(k,I)i+1/2,j − 2u

(k,l+1)i,j − u(k,l+1)

i−1/2,j

h2l+1

+u

(k,l+1)i,j+1/2 − 2u

(k,l+1)i,j − u(k,l+1)

i,j−1/2

h2l+1

,

(11.20)

dove u(k,I)i+1/2,j si ottiene per interpolazione quadratica perpendicolare all’ in-

terfaccia.Per i nodi di tipo F , ad esempio F = (i, j − 1/2) poniamo

(D(l,l+1)v)i,j =u

(k,I)i+1/2,j−1/2 − 2u

(k,l+1)i,j−1/2 − u

(k,l+1)i−1/2,j−1/2

h2l+1

+

u(k,l+1)i,j − 2u

(k,l+1)i,j−1/2 − u

(k,l+1)i,j−1

h2l+1

, (11.21)

dove u(k,I)i+1/2,j−1/2 e il valore ottenuto

114

Ghost region

1 2 3 4 5 6 7 8 10 11

2

3

4

5

6

7

8

9

10

11

9

1

Figura 11.12: Mesh test, vertex centered.

• prima calcolando il valore di u nel punto (i + 1, j − 1/2), per interpo-lazione quadratica in direzione parallela all’ interfaccia.

• Poi effettuando un’ altra interpolazione quadratica in direzione normaleall’ interfaccia per stimare u nel punto (i+ 1/2, j − 1/2).

Test

L’ implementazione di questi schemi sulle interfacce e un’ operazione com-plessa.

Per verificare che sia corretta, possiamo usare il seguente test: suppo-niamo di voler risolvere il problema (11.1) nel dominio Ω = [−1, 1]2. Lasoluzione test, u(x, y) sia una Gaussiana. Prendiamo un rettangolo, Q, con-tenuto nella zona in cui u(x, y) ≫ 0 (vedi figura 11.14). Per ogni mesh, τk,di livello l e diametro hk, costruiamo una sotto–mesh, µk,l, di diametro hk/2,che discretizza Q. Verifichamo che l’ errore del metodo Multigrid + AMRe, soddisfi le stime asintotiche, ossia sia e = O(h), nelle zone di interfaccia,e = O(h2) altrove.

115

Ghost region

1234567891011121314151617

1 3 7 9 11 132 4

56 8 10 12 14

1516

17

Figura 11.13: Raffinamento della mesh precedente.

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1−0.8

−0.6−0.4

−0.2 0

0.2 0.4

0.6 0.8

1

Figura 11.14: Problema test per le interfacce.

116

11.7.2 Adattamento della griglia

Come possiamo controllare l’ errore locale e raffinare la mesh solo nelle zonein cui l’ errore e “grande”? Non conoscendo la soluzione esatta, non possiamostimare l’ errore, quindi abbiamo bisogno di stimatori a priori dell’ errorelocale.

Lo illustriamo con un esempio. Supponiamo di voler risolvere il solitoproblema di Poisson

−∇ ∇u = f,

su un dominio poligonale Ω ⊂ R2. Discretizziamo il dominio con una meshtriangolare τk.

Sia µ un lato in τk, Tµ l’ unione dei due triangoli che hanno µ in comune.Consideriamo l’ indicatore locale di errore (local error indicator) [BRS02]

εµ(uk)2 =

∑

T⊂Tµ

h2T‖f +∇ ∇uk‖2L2(T ) + (11.22)

hµ‖[n ∇uk]n‖L2(µ), (11.23)

dovehT = |T |1/2, hµ = |µ|,

l’ operatore [·]n rappresenta il salto lungo la normale. Si noti che [n v]n eindipendente dall’ orientazione della normale.

Si puo provare che [BRS02]

|eh|H1(Ω) ≤ γ

(∑

µ

εµ(uk)2

)1/2

, (11.24)

dove γ e una costante che dipende solo dall’ errore di interpolazione.Si puo anche considerare una versione dell’ errore orientata ai triangoli,

invece che ai lati, ossia

εT (uk)2 = h2

T‖f +∇ ∇uk‖2L2(T ) + (11.25)∑

µ⊂∂T

hµ‖[n ∇uk]n‖L2(µ). (11.26)

Vale una relazione tipo (11.24).Un’ altra stima dell’ errore locale per il problema di Poisson e proposta

ed analizzata in [D96]. Se abbiamo due livelli, k e k + 1, possiamo stimareil gradiente dell’ errore ‖∇ek‖, usando i salti delle derivate normali sui latidella mesh. Nella versione orientata ai lati, essa si basa sulla misura di errore

η2µ = |µ| · ‖[∂nuk]µ‖20; µ,

117

dove µ e un qualsiasi lato interno alla mesh di discretizzazione. L’ operatore[∂n·]µ calcola il salto della derivata normale su µ. La norma ‖·‖0; µ e la norma

L2 sul lato µ. Equivalente e la misura basata sui triangoli

η2T =

1

2

∑

µ⊂∂T\∂Ω

|µ| · ‖[∂nuk]µ‖20; µ.

Bisogna raffinare la mesh laddove l’ indicatore di errore locale, ε, e “gran-de”. Tuttavia non e sicuro che dove ε e piccolo, anche l’ errore sia piccolo:influenze a distanza possono rendere grande l’ errore anche dove ε e piccolo.

Come stabilire se ε e grande? Data una tolleranza θ, possiamo decideread esempio di raffinare l’ elemento se

εT > θ · uk(T ),

dove uk(T ) e il valore (o un valore medio) di uk in T .Se εT e grande, effettuiamo la suddivsione regolare di T in 4 triangoli ad

esso simili.Un approccio semplice che vale per problemi privi di singolarita e quello

proposto ad esempio in [JFH+98]. Se u e la soluzione approssimata delproblema, per un elemento T della griglia, si definisce il local truncationerror

ǫT = |T | |∇u||u| .

Per ogni zona, Z = ∪Mi=1Ti, del dominio (definita in base al problema), sidefinisce l’ average truncation error

ǫZ =1

M

M∑

i=1

ǫTi.

La zona Z viene raffinata se ǫZ > θ, essendo θ una tolleranza definita in baseal problema.

118

Appendice A

Progetto di un contenitore

A.1 Il problema

La legge di stato di un gas perfetto e data dalla eq.:

pV = nRT (A.1)

dove p e la pressione, V e il volume, n e il numero di grammi-molecole, T ela temperatura assoluta ed R e la costante universale dei gas.L’eq. (A.1) e valida per pressioni e temperature comprese in un intervallolimitato ed inoltre alcuni gas la soddisfano meglio di altri.Per i gas reali l’eq. (A.1) e stata corretta da van der Waals nella relazione:

(p+

n2a

V 2

)(V

n− b)

= RT (A.2)

dove le costanti a e b sono determinate empiricamente per ogni specifico gas.Volendo progettare dei contenitori per l’anidride carbonica(CO2) e l’ossigeno(O2) per diverse combinazioni dei valori di p e T , e necessario calcolareaccuratamente con la (A.2) il volume

v = V/n (A.3)

occupato da una grammo-molecola.Sono assegnati i seguenti dati:

119

R = 0.082054 l· atm/(K· mol)

a = 3.592 atm l2/mol2

per CO2

b = 0.04267 l/mol

a = 1.360 atm l2/mol2

per O2

b = 0.03183 l/mol

Le pressioni e le temperature di progetto sono:

p(atm) 1 1 1 10 10 10 100 100 100T(K) 300 500 700 300 500 700 300 500 700

A.2 Compiti

Scrivere un programma che, per una assegnata tabella di valori di p e T eper ogni specifico gas, calcola con la (A.1) il valore v0 iniziale e trova il valorefinale vfin risolvendo la (A.2) con lo schema di Newton-Raphson. Stampareper i due gas assegnati e per ogni valore di p e T di progetto, v0, vfin ed ilnumero di iterazioni richieste per ottenere vfin.

120

Appendice B

Modelli di crescita

B.1 Introduzione

I modelli di crescita di popolazioni assumono che la loro densita, p, siaproporzionale alla popolazione esistente, ossia

dp

dt= kp, (B.1)

dove k non dipende dalla concentrazione di nutrimento disponibile. Anchequando il nutrimento non scarseggia, la crescita in ambiente chiuso vienelimitata dalle sostanze di rifiuto che gli esseri viventi producono. Questiprodotti inibiscono la crescita quando la densita raggiunge un valore massimopm. In questa situazione, la (B.1) si modifica nella

dp

dt= Kp(pm − p). (B.2)

Prendiamo come esempio una popolazione di cellule in un liquido di coltura.Se la dimensione di p e [p] = cellule/litro, e il tempo e misurato in giorni, Kha dimensione [K] = (litri/cellule)/giorno. Posto p(0) = p0, la soluzione di(B.2) e:

p(t) =pm

1 + (pm/p0 − 1) exp(−Kpmt). (B.3)

Questa funzione viene chiamata modello di crescita logistica, o semplicementelogistica [Wik08]. Il suo andamento e esemplificato in figura B.1.

Come esempio di applicazione del modello logistico, consideriamo il pro-blema di stimare l’evoluzione del mercato di un certo tipo di elaboratori.All’istante t = 0 sono stati venduti pochi elaboratori, diciamo p0 = 100. Enoto che dopo t = t = 60 settimane, ne sono stati venduti 25.000 e che K

121

0

10000

20000

30000

40000

50000

60000

70000

0 20 40 60 80 100

dens

ity

time

(63200.0 /(1.0 + (63200.0 / 10.0 - 1.0) * exp(-2.0e-6 * 63200.0 * x)))63200.0

Figura B.1: Curva logistica.

vale 2 × 10−6. Si vuole estrapolare il numero di elaboratori venduto, pT ,dopo T = 90 settimane. Poniamo [p] = numero elaboratori, misuriamo iltempo in settimane, allora [K] = 1/(elaboratori · settimane). Sostituendo leinformazioni disponibili per t = t nella eq. (B.3), si ottiene una equazionenon lineare, che va risolta nella variabile pm. Sostituendo pm nella (B.3), sipuo infine calcolare pT = p(T ).

B.2 Compiti

Calcolare la quantita pm,

• usando lo schema di Newton-Raphson,

• usando lo schema della secante variabile,

• usando lo schema della secante fissa.

Siask = |(pk+1 − pk)/pk+1|

lo scarto relativo alla k-esima iterazione. Arrestare le iterazioni quando sk <τ , τ = 10−2, 10−4, 10−6. Calcolare i corrispondenti valori di p(T ).

122

Tabulare il numero di iterazioni effettuate, le approssimazioni intermediee gli scarti relativi.

Tabulare e graficare le soluzioni numeriche ottenute, in maniera oppor-tuna (non bisogna produrre tabelle piu lunghe di una pagina, ne graficiilleggibili).

123

Appendice C

Flusso di calore

C.1 Il problema

In molti campi dell’ingegneria si pone il problema di determinare la tempe-ratura in corpi solidi. In problemi stazionari e per solidi omogenei ed isotropila temperatura T e governata dall’equazione di Laplace:

∂2T

∂x2+∂2T

∂y2= 0

Le derivate seconde possono approssimarsi con schemi alle differenze finite(si veda il numero 9.10 del testo Metodi Numerici di G. Gambolati). Su unagriglia regolare come quella della figura C.1.

i,j i+1,j

i,j-1

i,j+1

i-1,j

1

2

3

4

C

5

6

7

8

9

∆

∆

x

y

T=0

T=0

T=100

T=100 o

o

o

o

C C

C

Figura C.1: Stencil per l’ approssimazione delle derivate parziali ediscretizzazione di una piastra.

L’ approssimazione alle differenze finite delle derivate seconde sul nodo i,

124

j e:∂2T

∂x2=Ti+1,j − 2Ti,j + Ti−1,j

∆x2

∂2T

∂y2=Ti,j+1 − 2Ti,j + Ti,j−1

∆y2

N.B.: ∆x2 = (∆x)2,∆y2 = (∆y)2.Assumendo ∆x = ∆y (reticolo regolare quadrato) l’equazione di Laplace

sul nodo i, j e approssimata da:

Ti+1,j + Ti−1,j + Ti,j+1 + Ti,j−1 − 4Ti,j = 0

Se i nodi del reticolo sono n, ne risultera un sistema lineare di n equazioninelle temperature nodali.

Si consideri ora la piastra di figura C.1.I lati della piastra sono mantenuti alle temperature ivi indicate. Le tem-

perature sui 9 nodi numerati si ottengono risolvendo il sistema:

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

T1

T2

T3

T4

T5

T6

T7

T7

T8

T9

=

−100−100−200

00−100

00−100

C.2 Compiti

Costruire un programma che risolva il sistema assegnato prima con unoschema numerico diretto e poi con lo schema di sovrarilassamento:

(ωL+D)xk+1 = [(1− ω)D − ωU ]xk + ωb

Partendo da ω = 1 e finendo con ω = 1.5 con passo ∆ω = 0.05 si risolva ilsistema per i diversi valori di ω. Si fissi la tolleranza in uscita ε come segue:

| xk+1 − xk |< ε

125

con ε ≤ 10−2, essendo | . | la norma euclidea.Riportando in grafico il numero di iterazioni al variare di ω si stimi un

valore approssimato di ωopt.

Facoltativo: Osservando che la matrice gode di proprieta A ed e coe-rentemente ordinata si calcoli esattamente ωopt con la formula:

ωopt =2

1 +√

1− µ21

=2

1 +√

1− λ1,s

essendo µ1 l’autovalore massimo della matrice di Jacobi del sistema assegnatoe λ1,s quello della matrice di iterazione di Seidel.

Si calcoli λ1,s applicando il metodo delle potenze ai vettori differenza:

dk+1 = xk+1 − xk

per ottenere:

λ1,s = limk→∞| dk+1 || dk |

Si riporti in grafico semilogaritmico la norma euclidea | dk | in funzionedi k. Si calcoli quindi, utilizzando il segmento rettilineo del grafico, il fattoredi convergenza e lo si confronti col valore trovato di λ1,s.

126

Appendice D

Calcolo di angoli

D.1 La Cefalometria

L’ introduzione della radiografia cefalometrica1, avvenuta nel 1934 da partedi Hofrath in Germania e di Broadbent negli Stati Uniti, ha fornito unostrumento di ricerca e clinico per lo studio delle malocclusioni (problemi dimasticazione) e delle disarmonie scheletriche loro connesse.

In origine il loro scopo era la ricerca di un modello di crescita, ma fuben presto chiaro che si poteva utilizzarle per analizzare le proporzioni den-tofacciali e chiarire le basi anatomiche delle malocclusioni. Queste analisi,realizzate in vari modelli dalle singole scuole, tendono ad offrire al clinico lapossibilita di interpretare le relazioni reciproche tra le piu importanti unitafunzionali ed estetiche del volto. Si basano sull’analisi di valori lineari edangolari individuati da reperti anatomici standard.

Per trovare letteratura sull’argomento, si puo effettuare una ricerca suWEB con parole chiave: Cefalometria Tweed, oppureCefalometria Ricketts, oppure Cefalometria Steiner.

Ad oggi la ricerca si sta dirigendo verso analisi in 3D del cranio, ma l’inda-gine standard (compatibile con i costi e i tempi) e quella della Teleradiografiadel cranio in proiezione Latero-Laterale, standardizzata come segue: La di-stanza della sorgente radiogena dal piano mediosagittale del paziente deveessere di 60” e la pellicola deve essere a 15” da questo piano.

127

Figura D.1: Radiografia esemplificativa.

128

D.2 Descrizione del problema

La figura D.1 fornisce un esempio di Teleradiografia. Sono evidenziati 6punti. A partire da rette che passano per questi punti (esse sono le proiezionilaterali di piani ideali, noi le chiameremo per semplicita assi), si costruisceun triangolo, detto triangolo di Tweed.

Si usano tre assi:

• Asse di Francoforte, che va dal punto piu alto del meato acustico esterno(chiamato Porion, nella figura abbreviato in “Po”) al punto inferioredel margine orbitale (chiamato Orbitale, abbreviato in “Or”).

• Asse Mandibolare, che va dal punto piu basso posteriore del corpo dellamandibola (Gnation, “Gn”) al punto piu basso e avanzato del corpodella mandibola (Menton, “Me”).

• Asse Incisivo, costituito dall’asse dell’incisivo inferiore.

I tre angoli formati dagli assi sono detti

• FMA: angolo tra asse di Francoforte e asse Mandibolare;

• FMIA: tra asse di Francoforte e asse Incisivo;

• IMPA: tra asse Mandibolare e asse Incisivo.

Si prendono tutti positivi; la loro somma, esendo gli angoli interni di untriangolo, deve essere π = 180o.

D.3 Approccio risolutivo

Vogliamo risolvere il problema generale di calcolare l’angolo sotteso da duevettori, individuati fornendo due coppie ordinate di punti, ad esempio (Me,Gn) e (Or, Po) nella figura D.1.

Consideriamo i due vettori v1 = P2 − P1, v2 = P4 − P3 (nell’esempioprecedente, sarebbe P1 = Me, P2 = Gn, P3 = Or, P4 = Po).

Poniamo P(0)1 = P1, P

(0)2 = P3, P

(0)i = (x

(0)i , y

(0)i ).

Al variare di α1 e α2 in R, le due rette, r1 e r2, sono date dai punti [Lyo95]:

r1 = P(0)1 + α1v1, r2 = P

(0)2 + α2v2,

1Con la collaborazione del dott. I. Gazzola

129

oppure equivalentemente dai punti (xi, yi), tali che

v(2)i (xi − x(0)

i )− v(1)i (yi − y(0)

i ) = 0, i = 1, 2.

Il loro punto di intersezione e la soluzione (x, y) del sistema

v

(2)1 (x− x(0)

1 )− v(1)1 (y − y(0)

1 ) = 0

v(2)2 (x− x(0)

2 )− v(1)2 (y − y(0)

2 ) = 0

ossia v

(2)1 x− v

(1)1 y = v

(2)1 x

(0)1 − v

(1)1 y

(0)1

v(2)2 x− v

(1)2 y = v

(2)2 x

(0)2 − v

(1)2 y

(0)2

(D.1)

Il coseno dell’angolo α formato dalle due rette e:

c = cos(α) = cos(r1r2) =v1 · v2

‖v1‖2‖v2‖2,

quindi l’angolo (in radianti) formato dalle due rette e:

α = arccos(c),

ossia in gradi

αo =360α

2π.

D.4 Compiti

Scrivere un programma in MATLAB che

• visualizza la radiografia disponibile nella pagina WEB$CN/Esercitazioni/ortodon/Rx.jpg;

• consente di acquisire due coppie ordinate di punti (P1, P2), (P3, P4);

• disegna i vettori P2 − P1, P4 − P3, sull’immagine, segnando il punto incui si incontrano le due rette cui appartengono;

• calcola l’ angolo compreso e ne riporta il valore sull’immagine, in gradie radianti.

Costruire un esempio in cui il sistema (D.1) e mal condizionato e discutereil caso.

130

Appendice E

Schemi di rilassamento

E.1 Introduzione

Supponiamo di voler risolvere un sistema lineare le cui incognite rappresen-tano valori associati ai nodi di una griglia uniforme (vedi figura E.1).

I metodi iterativi noti come metodi di rilassamento sulle coordinate (coor-dinate relaxation), per la risoluzione di tali sistemi lineari calcolano ogni

nuova componente della soluzione, x(k+1)m , associata al nodo (i, j), utilizzan-

do solo le componenti relative ai nodi adiacenti (i− 1, j), (i+1, j), (i, j− 1),(i, j + 1).

Il miglioramento ad ogni iterazione e dunque locale e per ridurre l’errorea livello globale occorre “propagare” la riduzione su tutto il vettore.

Cerchiamo di approfondire questa affermazione, studiando il compor-tamento di uno di questi metodi, quello di Gauss–Seidel [QS06], che perrisolvere il sistema

Ax = b, (L+D + U)x = b, (L+D)x = −Ux + b

usa la formulazione

(L+D)x(k+1) = −Ux(k) + b,

dalla quale si ottiene

x(k+1) = −(L+D)−1Ux(k) + (L+D)−1b,

ossia

x(k+1)i =

1

aii

(bi −

i−1∑

j=1

aijx(k+1)j −

n∑

j=i+1

aijx(k)j

)

131

x

y

i

j(i, j)

Figura E.1: Griglia uniforme quadrata.

E.2 Compiti

• Costruire una matrice strettamente diagonalmente dominante, di ordinen = 2m > 64.

• Applicare il metodo di Gauss–Seidel per risolvere il sistema, calcolandoad ogni iterazione il vettore errore, e(k), e la sua trasformata discreta diFourier, F (k). Arrestare le iterazioni quando il residuo relativo e minoredi τ = 10−4. L’iterazione finale sia la K-esima.

• Disegnare

– le componenti degli errori e(k), k = 0, 2, 4, 8, ..., k < K;

– le componenti di F (k), k = 0, 2, 4, 8, ..., k < K;

– i valori di F(k)i , i = 1, n/2, n al variare di k.

Discutere i grafici ottenuti, determinando se gli andamenti delle compo-nenti in frequenza delle armoniche piu basse dell’ errore decrescono con k,oppure no.

132

Appendice F

Valutazione di sistemi

F.1 Descrizione del problema

Consideriamo la rete tandem con feedback [Bal99] illustrata in Fig. F.1,con processo di arrivo geometrico di parametro λ, nodi a servizio geometricodi parametro p1 e p2. Un utente che esce dal nodo 2 lascia il sistema conprobabilita p. Con probabilita 1−p ritorna al nodo 1. Supponiamo che ad undato istante vi siano n1 utenti nel nodo 1 e n2 utenti nel nodo 2. La coppia(n1, n2) rappresenta lo stato del sistema. L’insieme E = (n1, n2)t.c.ni ≥0 e lo spazio degli stati. La funzione π(n1, n2) = (π1, π2, . . . , πn)

T e ladistribuzione stazionaria di stato [Bal99]1.

La porzione iniziale del diagramma degli stati del processo Markovianoassociato e illustrata in Fig. F.2.

Ipotizziamo di avere un numero massimo n di utenti ammessi nel sistema(arrivi con perdita). Il diagramma degli stati e allora finito. Le probabilitadi transizione sono riassunte nella tabella F.1.

La matrice Q delle probabilita di transizione del processo ha ordine

1Nota: nel testo della Prof. Balsamo, i vettori sono vettori riga

Figura F.1: Rete tandem con feedback, a tempo discreto.

133

Figura F.2: Diagramma degli stati. NOTA: p1 = µ1, p2 = µ2.

da a P Descrizione(n1, n2) (n1 + 1, n2) λ Arrivo di un elemento(n1, n2) (n1 − 1, n2 + 1) p1 Passaggio interno al sistema(n1, n2) (n1, n2 − 1) p2p Abbandono del sistema(n1, n2) (n1 + 1, n2 − 1) p2(1− p) Uscita e rientro

Tabella F.1: Probabilita di transizione.

134

stato (0,0) (1,0) (0,1) (2,0) (1,1) (0,2) (3,0) (2,1)num. 1 2 3 4 5 6 7 8stato (1,2) (0,3) (4,0) (3,1) (2,2) (1,3) (0,4)num. 9 10 11 12 13 14 15

Tabella F.2: Numerazione degli stati.

M = (n + 1)(n + 2)/2. Essa gode della proprieta Q · 1 = 1, essendo1 = (1, 1, . . . , 1)T . La soluzione stazionaria e la soluzione positiva non nulladel sistema [Bal99]

QTπ = 0, Q = Q− I, (F.1)

soggetta alla condizione di normalizzazione πT1 = 1. Si puø dimostrare cheil rango di Q e M − 1, quindi esiste una sola soluzione π al problema.

Ordiniamo i nodi del grafo in modo che il nodo corrispondente allo stato(i, j), sia l’ m-esimo della numerazione, m = j+1+k(k+1)/2, k = i+j ≤ n.La tabella F.2 mostra la numerazione degli stati quando n = 4.

Poniamo a = p1, b = p2p, c = p2(1 − p). Quando n = 4, la matriceQ(15× 15) e:

γ1 λ 0 0 0 0 0 0 0 0 0 0 0 0 00 γ2 a λ 0 0 0 0 0 0 0 0 0 0 0b c γ3 0 λ 0 0 0 0 0 0 0 0 0 00 0 0 γ4 a 0 λ 0 0 0 0 0 0 0 00 b 0 c γ5 a 0 λ 0 0 0 0 0 0 00 0 b 0 c γ6 0 0 λ 0 0 0 0 0 00 0 0 0 0 0 γ7 a 0 0 λ 0 0 0 00 0 0 b 0 0 c γ8 a 0 0 λ 0 0 00 0 0 0 b 0 0 c γ9 a 0 0 λ 0 00 0 0 0 0 b 0 0 c γ10 0 0 0 λ 00 0 0 0 0 0 0 0 0 0 γ11 a 0 0 00 0 0 0 0 0 b 0 0 0 c γ12 a 0 00 0 0 0 0 0 0 b 0 0 0 c γ13 a 00 0 0 0 0 0 0 0 b 0 0 0 c γ14 a0 0 0 0 0 0 0 0 0 b 0 0 0 c γ15

(F.2)dove γi = 1−∑M

j=1,i6=j Qij .

135

F.2 Compiti

Risolvere il sistema (F.1), per λ = 0.1, p1 = 0.1, p2/p1 = 0.2, 0.4, 0.8, 1.0,2.0, 4.0, 6.0, 8.0, p = 1/2.

Definiamo il condizionamento della matrice Q come κ(Q) = σ1/σM−1,dove σi, i = 1, . . . ,M , sono i valori singolari di Q. Studiare e graficare ilcondizionamento del sistema al variare dei parametri.

Al variare dei parametri, calcolare le soluzioni π e

• disegnare l’andamento di π(0, 0) (probabilita di permanenza nello statoinattivo);

• disegnare l’andamento della probabilita di funzionamento a regime mas-simo:

πmax =∑

(n1,n2) t.c. n1+n2=n

π(n1, n2)

136

Appendice G

Information Retrieval

G.1 Formulazione del problema

Vogliamo effettuare delle ricerche in una collezione di d documenti. Dataun’interrogazione, vogliamo estrarre i documenti rilevanti per quell’interro-gazione, ossia che contengono informazioni pertinenti l’interrogazione stessa.La tabella G.1 riporta un esempio, tratto da [BDJ99], in cui i “documenti”sono titoli di articoli.

All’interno dei documenti individuiamo i termini “significativi”, che chia-meremo semplicemente termini. La tabella G.2 individua i termini salienti.Abbiamo d = 5, t = 6. Costruiamo la matrice A, che nella riga i-esimae colonna j-esima riporta il numero di occorrenze dell’i-esimo termine nelj-esimo documento.

A =

1 0 0 1 01 0 1 1 11 0 0 1 00 0 0 1 00 1 0 1 10 0 0 1 0

Normalizziamola, ossia dividiamo ogni colonna c per ‖c‖2. Otteniamo cosıla matrice A.

A =

0.5774 0 0 0.4082 00.5774 0 1.0000 0.4082 0.70710.5774 0 0 0.4082 0

0 0 0 0.4082 00 1.0000 0 0.4082 0.70710 0 0 0.4082 0

137

Abbiamo ottenuto un modello vettoriale [SM83] della collezione di docu-menti. Il documento j-esimo della collezione e rappresentato da un vettoredj di t elementi. L’intera collezione e rappresentata da una matrice At×d(t ≫ d) di rango r, le cui colonne rappresentano documenti e le righe ter-mini. Il generico elemento Aij e funzione della frequenza con cui l’i-esimotermine compare nel j-esimo documento. Il valore determina l’importanzadel termine nel rappresentare la semantica del documento. Un’interrogazioneviene rappresentata da un vettore q di t elementi.

Per determinare la similarita fra l’interrogazione q e il j-esimo documentodj della collezione, misuriamo il coseno dell’angolo compreso fra i due vettori

cosj(q) =djT q

‖dj‖2‖q‖2=

(Aej)T q

‖Aej‖2‖q‖2, (G.1)

dove ej e il j-esimo vettore coordinato. Quanto pi˘ grande e cosj(q), tanto piui due vettori sono paralleli, ossia il documento e pertinente all’interrogazione,o rilevante.

G.1.1 Latent Semantic Indexing

Utilizziamo la tecnica chiamata Latent Semantic Indexing (LSI) [BDJ99],per effettuare, data un’interrogazione q, un’efficiente ricerca dei documentipertinenti, basata sulla decomposizione a valori singolari (vedi par. 5.3),A = UΣV T . Sia U = [u1 . . . ut], V = [v1 . . . vd], Σ = diag(σ1, . . . , σd),σ1 ≥ . . . ≥ σr > σr+1 = . . . = σd = 0. Poniamo Σk = diag(σ1, . . . , σk),Uk = [u1 . . . uk], Vk = [v1 . . . vk]. Calcoliamo la matrice ridotta di rango k,k < r,

Ak = UkΣkVkT . (G.2)

Ak contiene solo k ≤ r ≤ d componenti linearmente indipendenti di A, quindise k < r e uno spazio di dimensione ridotta. Sostituiamo nella formula (G.1)la matrice A con la sua approssimazione di rango k. Essendo le matrici U eV ortonormali, otteniamo

cosj(q) =(Akej)

T q

‖Akej‖2‖q‖2=sjT (Uk

T q)

‖sj‖2‖q‖2, (G.3)

sj = ΣkVkT ej. Il costo per calcolare cosj(q) usando la (G.3) e inferiore a

quello necessario per valutare la (G.1). Il vantaggio non si ferma qui: pocheAk cattura la maggior parte della struttura semantica dalla base di dati, usarela (G.3) riduce i problemi di sinonimia che affliggono i metodi di ricerca basatisul modello vettoriale [BDJ99].

138

G.1.2 Valutazione delle prestazioni

Classifichiamo i documenti della collezione come riportato nella tabella G.3.Due parametri particolarmente importanti nella valutazione delle presta-

zioni dei metodi per la ricerca ed il recupero delle informazioni [SM83] sonoil richiamo R e la precisione P . Sono definiti nel seguente modo:

R =RELRET

RELRET +RELNRET, (G.4)

P =RELRET

RELRET +NRELRET. (G.5)

Il richiamo misura la capacita di recuperare tutti i documenti rilevanti, men-tre la precisione misura la capacita di recuperare solo documenti rilevanti.A valori maggiori corrisponde una maggiore efficienza.

G.2 Obiettivi

Realizzare un programma Matlab per determinare la precisione P in funzionedi R e k. Verificare il comportamento riportato nella Tabella G.4.

Come base di sperimentazione si utilizzi la collezione campione e le inter-rogazioni delle quali e data la rappresentazione vettoriale all’ URL$CN/CalcoloNumerico/Esercitazioni/LSI/.

G.3 Parte facoltativa

Studiare l’andamento

• di P in funzione di k per almeno due diversi valori di R.

• di P in funzione di R per almeno tre diversi valori di k.

Scegliere opportunamente i valori non esplicitamente fissati in questotesto e tracciare grafici riassuntivi.

139

j Documento1 How to Bake Bread Without Recipes2 The Classic Art of Viennese Pastry3 Numerical Recipes: The Art of Scientific Computing3 Breads, Pastries, Pies and Cakes: Quantity Baking Recipes5 Pastry: A Book of Best French Recipes

Tabella G.1: Collezione di documenti

i Termine1 bak(e,ing)2 recipes3 bread(s)3 cakes5 pastr(y,ies)6 pies

Tabella G.2: Termini considerati

Rilevanti Non rilevanti(RELevant) (Not RELevant)

Recuperati RELRET NRELRET(RETrieved)

Non recuperati RELNRET NRELNRET(Not RETrieved)

Tabella G.3: Classificazione dei documenti

140

i R P k1 0.50 0.1579 822 0.75 0.2000 823 1.00 0.0610 824 0.50 0.1429 305 0.75 0.1667 306 1.00 0.1429 307 0.50 0.2143 158 0.75 0.1600 159 1.00 0.1351 1510 0.50 0.1875 811 0.75 0.2000 812 1.00 0.1220 8

Tabella G.4: Comportamento rilevato.

141

Appendice H

Costo computazionale

H.1 Introduzione

Si vuole stimare il costo computazionale di due algoritmi numerici per risol-vere un sistema lineare Ax = b, di ordine n.

• Sia G lo schema di eliminazione di Gauss, senza pivoting.

• Sia J il metodo di Jacobi.

Consideriamo la matrice tridiagonale A che ha gli elementi tutti nulli,tranne

Aij =

8, se i = j,−1, |i− j| = 1

, i, j = 1, . . . , n.

Poniamo A = A + p, p = 10−5. Chiamiamo p parametro di perturbazione.Sia s(n) uno schema risolutivo applicato alla risoluzione di un sistema di

ordine n. Indichiamo con C[s(n)] il suo costo computazionale, ossia il nume-ro di operazioni floating point necessarie per risolvere il sistema. PoniamoC[⊕] = C[⊗] = C[⊘] = 1, dove ⊕, ⊗, ⊘ sono la somma, moltiplicazione edivisione floating point.

Abbiamo

• C[G(n)] = O(2n3/3) [BCM88], poiche occorrono n3/3 moltiplicazionie n3/3 addizioni per eseguire l’algoritmo.

• C[J(n)] = O(2m · n2).

La seconda stima si ottiene usando i risultati del paragrafo 6.4. Bisognacalcolare l’errore relativo, eG, che si commette con il metodo di Gauss. SiaeG ≃ 10−k, allora m = [−k/ log(0.25)] + 1 e C[J(n)] = O(2m · n2). Nota:[x]=parte intera di x.

142

H.2 Compiti

Scrivere un programma MATLAB che, utilizzando l’ opzione profile e lafunzione cputime calcola il numero di operazioni floating point eseguite daivari metodi e il tempo speso per completarne l’esecuzione su una o piu mac-chine. Disegnare i grafici che riportano l’ andamento dei tempi e del numerodi operazioni, assieme ai grafici dei costi computazionali teorici. Confrontarei tempi di CPU per l’esecuzione degli schemi di Gauss e Jacobi al variare din, 4 ≤ n ≤ 512.

Calcolare il punto di inversione, ossia l’ordine oltre il quale Jacobi diventacomputazionalmente piu conveniente di Gauss.

H.3 Approfondimenti

Il costo delle operazioni floating point non e lo stesso: moltiplicazioni e divi-sioni costano piu delle somme. Scrivendo un opportuno programma e usandol’ opzione profile e possibile stimare il costo delle singole operazioni. Po-nendo il costo della somma C⊕ = 1, siano C⊗ e C⊘ i costi delle altre dueoperazioni. Ricalcolare i costi teorici degli schemi e graficarli assieme ai costidi CPU effettivamente valutati con MATLAB.

Ricalcolare il punto di inversione, e confrontarlo con il valore precedente.

H.4 Memorizzazione in forma sparsa

La matrice A ha in percentuale pochi elementi non nulli. Memorizzarla informa compatta, usando l’opzione sparse di MATLAB. Supponiamo di vo-ler risolvere il sistema Ax = b, utilizzando il metodo di Jacobi, anche seper una matrice tridiagonale simmetrica il metodo piu indicato e l’algorit-mo di Thomas. Eseguendo i prodotti Az in forma sparsa, stimare il costocomputazionale del nuovo schema e confrontarlo con i costi effettivi.

H.5 Facoltativo

Determinare l’ accuratezza della stima del costo di Jacobi al variare delparametro di perturbazione p.

143

Appendice I

Deformazione di progetto

I.1 Introduzione

Sono stati costruiti alcuni alberi per un’ imbarcazione con leghe sperimentalidi alluminio. Sono state eseguite delle prove per definire la relazione tra losforzo σ (forza per unita di area trasversale dell’albero) applicato al materialee allungamento relativo ε = ∆l/l, essendo l la lunghezza dell’albero e ∆l ilsuo allungamento sotto l’azione dello sforzo σ. Le relazioni sperimentali traσ ed ε sono riportate nella tabella I.1.

Si vuole determinare la variazione di lunghezza ∆lp di ogni albero (cheassumiamo tutti alti 10 m) per uno sforzo

σp = 735 kg/cm2

che e quello massimo di progetto calcolato per una determinata sollecitazionesotto l’azione del vento.

Per risolvere questo problema occorre interpolare o approssimare conqualche funzione i dati della tabella, quindi sostituire σp di progetto nella

σ(kg/cm2) ε(cm/cm)ǫ1 ǫ2 ǫ3

1.8000e+02 5.0000e-04 1.0000e-03 1.2000e-035.2000e+02 1.3000e-03 1.8000e-03 2.0000e-037.2000e+02 2.0000e-03 2.5000e-03 3.7000e-037.5000e+02 4.5000e-03 2.6197e-03 3.8197e-038.0000e+02 6.0000e-03 2.8276e-03 4.0276e-031.0000e+03 8.5000e-03 3.7655e-03 5.9655e-03

Tabella I.1: Relazioni sperimentali tra σ ed ε.

144

funzione interpolante e approssimante per ottenere εp e quindi calcolare ∆lpcome:

∆lp = εpl

I.2 Compiti

Scrivere un programma che fornisca i valori interpolati εp, corrispondentia σp = 735 kg/cm2, con polinomi di grado crescente dal 1 fino al 5. Siconfrontino tra loro i 5 valori per εp cosı ottenuti.

Si riportino in grafico gli andamenti dei polinomi di 5 grado nell’intervallo

180 ≤ σ ≤ 1000 kg/cm2.

Si rappresentino anche i valori sperimentali di σ ed ε e si analizzi il diversocomportamento tra il profilo sperimentale e quello del polinomio che passaper tutti i dati osservati.

Ne scaturira l’osservazione che l’interpolazione e poco adatta per questoproblema anche se i valori trovati per εp con i polinomi interpolati di ordinecrescente sono abbastanza stabili.

Si provveda allora a sostituire le interpolazioni con approssimazioni, ades. regressioni lineari ai minimi quadrati [QS06]:

ε = a0 + a1 σ

Si scriva una subroutine per il calcolo dei coefficienti a0 ed a1 e si stiminoi valori di εp, confrontandoli coi valori ricavati per interpolazione.

Il problema con la regressione lineare e che essa fornisce valori non rea-listici nell’origine, vale a dire una deformazione percentuale diversa da zerocon sforzo nullo, mentre in un materiale elastico per σ = 0 deve essere ε = 0.Una soluzione alternativa che evita l’inconveniente appena lamentato e quel-la di costruire regressioni lineari sui logaritmi (in base opportuna, ad esempio10) dei dati, cioe:

log ε = b0 + b1 log σ

La relazione esplicita per ε e allora:

ε = 10b0 · σb1

da cui si ricava ε(0) = 0.Si approssimino i logaritmi dei dati sperimentali con una retta ai minimi

quadrati calcolando b0 e b1. Si calcolino i valori εp ottenuti con questa ultimaapprossimazione e li si confronti coi precedenti.

145

Infine si riportino in grafico i dati sperimentali, i polinomi interpolatori di5 grado, le rette ai minimi quadrati e le curve ai minimi quadrati sui valorilogaritmici. Si tenti una discussione critica dei risultati e quindi dei 3 diversialgoritmi numerici che li hanno generati.

146

Appendice J

Grafica 2D

J.1 Formulazione del problema

Vogliamo rappresentare semplici figure geometriche nel piano e possibili tra-sformazioni affini su di esse, ossia trasformazioni che conservano gli allinea-menti.

Limitiamoci a rappresentare poligoni convessi. Ricordiamo che

Definizione J.1.1 Un insieme I e convesso se, dati due suoi punti qualsiasi,il segmento di retta che li congiunge e tutto contenuto in I.

Usando un approccio ormai classico in questo campo, rappresentiamo ipunti in coordinate omogenee [FvDFH90]. Ogni punto P del piano e unvettore di tre componenti: P = (x, y,W )T . Se W = 1, le coordinate sidicono normalizzate. Almeno una delle tre coordinate deve essere diversa dazero: (0, 0, 0)T non e un punto. Due terne v1 e v2 rappresentano lo stessopunto P se e solo se v1 = cv2, 0 6= c ∈ R. I punti P = (x, y,W )T tali cheW = 0 vengono chiamati punti all’infinito, o direzioni. Li chiameremo anchevettori–direzione.

Una trasformazione affine A e individuata da una matrice non singolare

A =

m11 m12 m13

m21 m22 m23

0 0 1

(J.1)

Considereremo solo i seguenti tre tipi di trasformazioni.

• Traslazioni. Sono spostamenti rigidi dell’origine del sistema di coordi-nate, che lasciano invariata l’orientazione degli assi. Traslare il puntoP = (x, y,W )T delle quantita tx e ty rispettivamente lungo gli assi

147

x e y, ottenendo cosı il punto P ′ = (x′, y′,W ′)T , equivale a eseguirel’operazione

P ′ =

x′

y′

W ′

=

1 0 tx0 1 ty0 0 1

xyW

= TP. (J.2)

T viene detta matrice di traslazione.

• Rotazioni. Ruotare una figura di θ radianti rispetto all’ origine degliassi equivale a trasformare ogni punto P nel punto

P ′ =

x′

y′

W ′

=

cos θ − sin θ 0sin θ cos θ 0

0 0 1

xyW

= R(θ)P. (J.3)

R(θ) e la matrice di rotazione di un angolo θ.

• Scalature. Scalare una figura di un fattore sx lungo x e di un fattoresy lungo y equivale a trasformare ogni punto P nel punto

P ′ =

x′

y′

W ′

=

sx 0 00 sy 00 0 1

xyW

= SP. (J.4)

S e la matrice di scalatura.

J.2 Obiettivi

Implementare delle funzioni per

• disegnare curve lineari a tratti aperte e chiuse;

• determinare l’intersezione di due segmenti;

• determinare le intersezioni di due curve lineari a tratti;

• disegnare la spline cubica naturale che passa per un insieme di puntiassegnati Pi = (xi, yi,Wi)

T , i = 1, . . . , n (senza utilizzare la funzioneMatlab che la calcola direttamente). Se xi ≥ xi+1, per qualche i, noncalcolarla e segnalare la situazione.

148

J.3 Rappresentazioni

• Rappresentare un segmento s in forma parametrica

s = P0 + tv,

dove P0 e un estremo del segmento, v il vettore che individua la suadirezione e la sua lunghezza, 0 ≤ t ≤ 1 il parametro.

• Individuare una spezzata di n lati, P, fornendo l’elenco dei suoi ver-tici, percorsi in senso antiorario, P = P1, . . . , Pn. In un poligono(spezzata chiusa), il primo punto e l’ultimo coincidono.

J.4 Compiti

Esercizio J.4.1 Utilizzando le funzioni messe a punto, disegnare il quadratoQ di lato unitario (e il riferimento cartesiano), centrato in P = (0, 0, 1)T , coni lati paralleli agli assi.

• Ruotarlo di θ = π/6.

• Traslarlo ulteriormente di tx = 1, ty = −1.

• Scalarlo di sx = 2, sy = 0.5.

Sia Q′ il nuovo poligono ottenuto con queste trasformazioni.

• Disegnare i segmenti che uniscono i punti di intersezione di Q e Q′.

Siano dati i punti della tabella T , disponibile all’ URL$CN/Esercitazioni/graf2d.crd.

Esercizio J.4.2 Disegnare le seguenti coppie di poligoni e determinarne glieventuali punti di intersezione.

• P1 = T1, T2, T3, T4, T1, P2 = T6, T7, T8, T9, T6.

• P3 = T1, T2, T3, T4, T1, P4 = T6, T7, T8, T6.

Esercizio J.4.3 Dati i punti Ti, i = 11, . . . , 31, disegnare la spline cubicanaturale che li unisce. Disegnare anche il polinomio interpolatore a tratti ela retta ai minimi quadrati sugli scarti verticali.

149

Nota

Particolare attenzione va data al condizionamento dei vari sistemi lineari chevengono risolti.

J.5 Parte facoltativa

Implementare la seguente strategia. Ogni vertice che entra in una nuovafigura viene inserito, se non e gia presente, nella tabella T . Le spezzate sonoindividuate dagli indici dei loro vertici, numerati secondo la tabella T .

Modificare le routines grafiche usando questo diverso approccio.

150

Appendice K

Analisi di carico

K.1 Il problema

La sezione trasversale di una imbarcazione e riprodotta in figura K.1.La forza che il vento esercita sulla vela, espressa in kg per m di altezza

dell’albero, varia con l’elevazione z come riportata nella figura. L’espressioneanalitica di tale pressione lineare e:

p(z) = 3003z

5 + 3zexp(−2z/10) kg/m

Assumendo che la cima dell’albero si sposti di un valore trascurabile sottol’azione del vento, si calcoli la forza F risultante della pressione del vento ela distanza verticale da O della corrispondente retta di applicazione.

La forza F e data dall’integrale:

F =

∫ 10

0

3003z

5 + 3zexp(−2z/10)dz

mentre d si ottiene dal rapporto:

d =

∫ 10

0z p(z)dz

∫ 10

0p(z)dz

=

∫ 10

0z300 3z

5+3zexp(−2z/10)dz

F

K.2 Compiti

Scrivere un programma che calcoli F e d utilizzando sia la formula dei trapeziche quella di Cavalieri-Simpson. Dividere l’intervallo di integrazione un nu-mero crescente n di sottointervalli a partire da n = 2 per finire con n = 128.

151

a

F

d

1 m

vento

10 m

z=10 m

z=0 mT

Figura K.1: Sezione trasversale di una imbarcazione.

Tabulare, al crescere di n, i valori ottenuti per F e per d coi due algoritmi.Commentare i risultati.

152

Appendice L

Moto di un giocattolo

L.1 Formulazione del problema

Un’ auto giocattolo si muove su un piano illimitato a velocita costante. Adun assegnato istante l’ asse anteriore inizia a girare a destra ad un tassocostante. L’ asse posteriore dell’ auto e fisso, quello anteriore puo girareindefinitamente.

Determinare il cammino percorso dall’ automobilina.

L.2 Procedimento risolutivo.

Formalizzando la soluzione come in [AM88], sia (x, y) il centro dell’ asseanteriore dell’ auto e s la lunghezza del cammino percorso dal centro stesso(vedi figura L.1). La distanza l tra gli assi sia 1. Sia ξ l’ angolo formatodall’ auto con il centro degli assi di riferimento e sia θ l’ angolo che l’ asseanteriore forma con l’ asse principale dell’ auto.

Supponiamo che il guidatore inizi a girare il volante quando x0 = 0, y0 = 0e che sia ξ0 = θ0 = 0.

Allora il cammino dell’ automobilina, ricordando che sia ds/dt che dθ/dtsono costanti, e governato dalle equazioni (vedi [AM88]):

dx/ds = cos(ξ + θ)dy/ds = sin(ξ + θ)dξ/ds = sin(θ)dθ/ds = k

(L.1)

k rappresenta la velocita con cui il guidatore gira il volante. Il sistemadi equazioni (L.1) assieme alla condizioni iniziali permette di calcolare latraiettoria dell’ automobilina al variare di k.

153

Figura L.1: Rappresentazione del problema.

154

L.3 Obiettivi

In figura L.2 sono riprodotte traiettorie per alcuni valori di k. L’ obiettivo equello di ricalcolare queste traiettorie integrando le equazioni (L.1).

L.4 Suggerimenti

Ecco alcuni modi per provare ad integrare numericamente le equazioni.

• Approssimare le 4 equazioni con lo schema alle differenze

du/ds = ∆u/h+O(h). (L.2)

• Generalizzare a sistemi di 4 equazioni gli schemi di Runge-Kutta pro-posti in [Gam94].

• Si puo provare ad utilizzare anche un pacchetto di manipolazione sim-bolica (es: Mathematica, Macsyma) o numerica (MATLAB), oppureuna libreria numerica (IMSL) per integrare le equazioni.

• Risolvere analiticamente le ultime due equazioni. Sostituire i risultatinelle prime due. Il sistema si riduce a:

dx/ds = cos(α(s))dy/ds = sin(α(s))

, α(s) =1

k(1− cos(ks)) + ks. (L.3)

Esso puø essere integrato ad esempio con lo schema dei trapezi. Con-frontare le soluzioni con quelle trovate con gli altri schemi.

L.5 Esempio di risoluzione

Applicando lo schema alle differenze di Eulero (L.2) al sistema otteniamo leequazioni:

xi+1 = xi + h cos(ξi + θi)yi+1 = yi + h sin(ξi + θi)ξi+1 = ξi + h sin(θi)θi+1 = θi + hk

(L.4)

Usando il pacchetto di elaborazione numerica MATLAB (vedi [The06]),si puo facilmente calcolare una soluzione, usando un programma come ilseguente:

155

Figura L.2: Traiettorie ottenute scegliendo alcuni valori di k. Le linee con-tinue rappresentano le traiettorie dell’ asse anteriore, quelle tratteggiate letraiettorie dell’ asse posteriore.

156

0 1 2 3 4 5 6 7

-7

-6

-5

-4

-3

-2

-1

0

0 1 2 3 4 5 6 7

-7

-6

-5

-4

-3

-2

-1

0

MATLAB + Finite Differences

k = -0.2

Figura L.3: Traiettoria ottenuta con lo schema di Eulero.

% risoluzione numerica del problema moto di un’auto, da

% SIAM Review 37(2) 1995

k = -0.2

sfin=20

nmax=100

h = sfin/nmax

x = 0*(1:(nmax+1));

y = x;

xi = x;

th = x;

i=1

s0 = 0

s(i) = 0;

x(i) = 0;

xi(i) = 0;

th(i) = 0;

for i=1:nmax

s(i+1) = s0 + i*h;

x(i+1) = x(i) + h*cos(xi(i) + th(i));

y(i+1) = y(i) + h*sin(xi(i) + th(i));

157

xi(i+1) = xi(i) + h*sin(th(i));

th(i+1) = th(i) + h*k;

end

plot(x,y)

[x’ y’]

%

In figura L.3 si vede la traiettoria calcolata per k = −0.2, ponendo h =0.2. Confrontandola con la corrispondente figura L.2 si vede che l’ accordo esoddisfacente.

158

Appendice M

Onde di traffico

M.1 Introduzione

Vogliamo modellizzare lo scorrimento delle auto in una corsia di un’ auto-strada1. Per semplificare il problema. supponiamo che le auto non possanosuperarsi e che non vi siano caselli nel tratto che si studia2

M.2 Descrizione e compiti

Seguiamo la linea di ragionamento proposta in [Hab98].Supponiamo di concentrare ogni auto nel suo baricentro, riducendola cosı

ad un punto. Sia u(x, t) il campo di velocita, ossia il valore u(x, t) e la velocitaassociata al punto x, all’ istante t. Se nel punto, in quell’ istante non si trovaalcuna auto, u(x, t) = 0.

Sia ρ(x, t) la densita del traffico, una funzione liscia3 in x e t.L’ ipotesi che si possano usare funzioni continue per studiare il problema,

va sotto il nome di ipotesi del continuo.Indichiamo con q il flusso di traffico. Il valore q(x, t) e la quantita di auto

che passa nell’ unita di tempo per il punto x, all’ istante t. Si puø vedere chevale la relazione

q(x, t) = ρ(x, t)u(x, t). (M.1)

Il principio di conservazione, afferma che la densita ρ e il campo di ve-locita, u(ρ), pensato come funzione della densita, sono legati dalla relazio-

1Con il contributo di A. Ceccato.2Queste assunzioni riducono la complessita del problema, ma ci allontanano anche

molto dalla situazione reale!3Per capire come si possa introdurre una funzione continua ρ(x, t), a partire da una

distribuzione discreta di auto puntiformi, si veda ancora [Hab98], par. 58.

159

ne [Hab98]

∂ρ

∂t+∂q

∂x= (M.2)

∂ρ

∂t+∂(ρ · u(ρ))

∂x= 0. (M.3)

Usando la regola della catena, dalla (M.2), otteniamo la relazione

∂ρ

∂t+dq(ρ)

dρ

∂ρ

∂x= 0. (M.4)

Ancora usando la regola della catena, si puø scrivere

dρ

dt=∂ρ

∂t+dx

dt

∂ρ

∂x. (M.5)

Confrontando le equazioni (M.4) e (M.5), possiamo dire che la densita ecostante nel tempo, ossia

dρ

dt= 0, (M.6)

sedx

dt=dq(ρ)

dρ≡ q′(ρ). (M.7)

Questa relazione si puø interpretare dicendo che se un osservatore, che ha laposizione x, si muove con velocita dx/dt data dalla relazione (M.7), alloraegli osserva una densita di traffico costante.

La densita tende a propagarsi in forma di onde, si parla di onde di densita.La relazione (M.7) ci dice che le onde di densita si propagano alla velocitaq′(ρ), che viene chiamata velocita di densita locale. Dato un punto x0, nota ladensita all’ istante 0, ρ(x0, 0) = ρ0, partendo dall’ istante t = 0 e integrandola funzione q′(ρ), si ottiene una traiettoria

x∗0(t) = x0 +

∫ t

0

q′(ρ)dt, (M.8)

lungo la quale la densita e costante. Lungo questa traiettoria, inoltre, lanostra equazione alle derivate parziali (M.4) si riduce all’ equazione alle deri-vate ordinarie (M.6). Le curve x∗0(t) vengono chiamate curve caratteristichedella (M.4).

Fissato x0 e noto ρ0, dato che ρ e costante lungo le linee caratteristiche,la densita all’ istante t > 0, nel punto x∗0(t), e ρ(x∗0, t) = ρ0.

160

Supponiamo di conoscere la densita iniziale, ρ0(x) = ρ(x, 0), all’istantet = 0, in ogni punto x. Possiamo calcolare una soluzione approssimata delproblema differenziale

∂ρ∂t

+ dq(ρ)dρ

∂ρ∂x

= 0, x ∈ [a, b], t ∈ [0, T ]

ρ(x, 0) = ρ0(x),(M.9)

usando il seguente metodo.Sia umax il limite di velocita sulla strada, che supponiamo nessuno superi!

Sia ρmax la densita massima raggiungibile, quella alla quale le auto si trovanoparafango contro parafango (in realta la densita massima corrisponde ad autoseparate tra loro da un piccolo spazio).

Per semplificare la trattazione, supponiamo che il campo di velocita di-penda dalla densita in modo lineare, esattamente

u(ρ) = umax

(1− ρ

ρmax

). (M.10)

Essendo q = ρu,

dq

dρ=

d (ρ umax(1− ρ/ρmax))dρ

= (M.11)

umax

(1− 2ρ

ρmax

). (M.12)

Poiche (x∗i (t), t) sono i punti di una curva caratteristica, ρ(x, t) e costantesu questa curva, quindi la relazione (M.8) diventa

x∗i (t) = xi(0) + umax

(1− 2ρ(xi(0))

ρmax

)t. (M.13)

Le curve caratteristiche sono dunque delle rette.Ripetendo questo ragionamento per ogni xi, otteniamo una soluzione del

problema (M.9), sui punti di N caratteristiche. Approssimando questi valoricon una funzione ρ(x, t) abbiamo una soluzione ρ del problema (M.9). Adesempio possiamo interpolare linearmente in t.

Questo schema per risolvere equazioni differenziali iperboliche [TS64] edetto metodo delle caratteristiche.

Consideriamo la situazione in cui nell’ intervallo [−1, 1] vi e una barriera,che rallenta il traffico e lo blocca, perciø non vi sono auto nella regione x > 1,mentre nella regione x < −1 le auto sono a densita massima, parafango–contro–parafango. Cosa accade se all’ istante t = 0 la barriera viene rimossa?

161

0

50

100

150

200

250

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

densita‘ (auto/km)

posizione (km)

Figura M.1: Caso in cui viene rimossa una barriera.

Modellizziamo questa situazione con il problema (M.9), in cui poniamo

ρ(x, 0) =

ρmax, se x < −1,(1− x)ρmax/2, se − 1 ≤ x ≤ 1,

0, altrimenti.. (M.14)

Qui e nel seguito, poniamo ρmax = 250 auto/Km, umax = 100 Km/h; misu-riamo le posizioni in Kilometri.

Esercizio M.2.1 Risolvere e disegnare i grafici delle caratteristiche a partiredai punti xi = −2 + i∆x, ∆x = 4/M , i = 0, . . . ,M = 40, nell’intervallo0 ≤ t ≤ T = 0.027 (0.027 ore ≃ 100 sec.). Ricavare e disegnare la densita ρdel problema, nel rettangolo R = [−2, 2]× [0, T ].

Esaminiamo il caso in cui una regione di bassa densita e seguita da unadi alta densita. In questo caso, si forma un’ onda di compressione, che sipropaga all’indietro sull’ asse delle x. Modellizziamo questa situazione con ilproblema (M.9), in cui poniamo

ρ(x, 0) = |(1 + π/2 + arctan(20 ∗ x))ρmax/(1 + π)|. (M.15)

In questa situazione, vi sono delle zone in cui le caratteristiche si inter-secano. In queste zone il metodo delle caratteristiche non funziona. Essopredice valori diversi di ρ(x, t), a seconda della caratteristica per (x, t) che siconsidera. In questi punti bisogna sviluppare un’analisi pi˘ sofisticata.

162

0

50

100

150

200

250

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

densita‘ (auto/km)

posizione (km)

Figura M.2: Caso di congestione davanti alla colonna.

Esercizio M.2.2 Risolvere e disegnare i grafici delle caratteristiche neglistessi punti dell’esercizio precedente. Ricavare e disegnare la densita ρ delproblema negli stessi punti dell’esercizio precedente. Non tracciare le carat-teristiche oltre i punti in cui si intersecano. Non calcolare i valori di ρ nellezone in cui le caratteristiche si intersecano.

Quando lungo una strada il traffico diventa piu denso, le caratteristichesi intersecano e la teoria predice contemporaneamente due diversi valori perla densita. Il metodo non e applicabile.

Abbiamo assunto che il campo di velocita e la densita del traffico fosse-ro funzioni continue (ipotesi del continuo). Rimovendo questa assunzione,possiamo trattare il caso in cui le caratteristiche si intersecano.

Ipotizziamo quindi che la densita di traffico e il campo di velocita sianodiscontinui in un certo punto, xs, dello spazio e che tale discontinuita si possapropagare nel tempo secondo una data legge xs(t). La discontinuita vienechiamata shock e il punto xs viene detto punto di shock.

Consideriamo una regione x1 < x < x2, tale che x1 < xs(t) < x2. Ilnumero di auto contenute e dato da:

N(t) =

∫ x2

x1

ρ(x, t) dx (M.16)

Questo integrale continua ad essere ben definito nonostante la discon-tinuita. La variazione della quantita di auto, N(t), presente nella regione,

163

0

50

100

150

200

250

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

densita‘ (auto/km)

posizione (km)

Figura M.3: Densita discontinua.

dalla legge di conservazione delle auto e data da:

dN(t)

dt=

d

dt

∫ x2(t)

x1(t)

ρ(x, t) dx = q(x1, t)− q(x2, t). (M.17)

Dividiamo l’intervallo in due parti∫ x2

x1

ρ dx =

∫ xs

x1

ρ dx+

∫ x2

xs

ρ dx. (M.18)

Dato che xs dipende dal tempo, anche x1 e x2 non sono costanti. Le leggedi derivazione di integrali con estremi dipendenti dal tempo e:

d

dt

∫ β(t)

α(t)

f(x, t) dx =dβ

dtf(β, t)− dα

dtf(α, t) +

∫ β(t)

α(t)

∂f(x, t)

∂tdx. (M.19)

Questa formula non vale se vi e una discontinuita nell’intervallo di inte-grazione, ecco perche lo dividiamo in due parti. Applicando la formula aidue integrali otteniamo:

d

dt

∫ xs

x1

ρ dx =dx−sdt

ρ(x−s , t)−dx1

dtρ(x1, t) +

∫ xs

x1

∂ρ(x, t)

∂tdx; (M.20)

d

dt

∫ x2

xs

ρ dx =dx2

dtρ(x2, t)−

dx+s

dtρ(x+

s , t) +

∫ x2

xs

∂ρ(x, t)

∂tdx; (M.21)

164

Introduciamo l’operatore di salto (jump operator)

[z(xs)] = [z] = z(x−s )− z(x+s ). (M.22)

La (M.19) diventa allora

d

dt

∫ x2

x1

ρdt = [ρ(xs)]dxsdt

+

dx2

dtρ(x2, t) −

dx1

dtρ(x1, t) +

∫ xs

x1

∂ρ

∂tdx +

∫ x2

xs

∂ρ

∂tdx = q(x1, t)− q(x2, t).

Se x1 e x2 sono vicini allo shock, il contributo dei due integrali e trascu-rabile, inoltre (dx2/dt)ρ(x2, t) ≃ (dx1/dt)ρ(x1, t), quindi (dx2/dt)ρ(x2, t) −(dx1/dt)ρ(x1, t) ≃ 0. In definitiva si ha:

[ρ]dxsdt

= [q]. (M.23)

Nei punti di discontinuita questa condizione di shock sostituisce l’equa-zione differenziale che e valida negli altri casi.

Per esempio, consideriamo un flusso di traffico che viaggia con densita ρ0

e che viene improvvisamente bloccato da un semaforo rosso situato nel puntox = 0. Ipotizziamo che le auto siano in grado di fermarsi istantaneamente.In x = 0 il traffico e bloccato (u = 0) e la densita e massima in ogni istantedi tempo successivo a quello iniziale (vedi figura M.3). Il modello cosı pro-posto presenta caratteristiche che si intersecano. Introduciamo uno shock:le densita prima e dopo lo shock sono rispettivamente la densita iniziale ρ0

e la densita massima ρmax. Nel punto di shock deve essere soddisfatta lacondizione:

dxsdt

=[q]

[ρ](M.24)

La posizione iniziale dello shock e nota e costituisce una condizione delproblema: xs(0) = 0.

Si ha quindi:[q]

[ρ]=u(ρmax)ρmax − u(ρ0)ρ0

ρmax − ρ0

(M.25)

Ricordando che u(ρmax) = 0, l’eq. (M.24) diventa

dxsdt

=−u(ρ0)ρ0

ρmax − ρ0< 0 (M.26)

165

La velocita dello shock e dunque costante e negativa e risolvendo l’equa-zione differenziale otteniamo

xs =−u(ρ0)ρ0

ρmax − ρ0

· t.

Esercizio M.2.3 Considerare la densita:

ρ(x, 0) =

ρmax

3, x < 0

ρmax, x > 0.(M.27)

Risolvere e disegnare i grafici delle caratteristiche, ricavare e disegnare ladensita ρ negli stessi punti dell’esercizio precedente.

166

Appendice N

Analisi della dinamica

N.1 Introduzione

E’ dato il sistema di equazioni differenziali:

y′1 = −[(−π y1 y4 + π y2 y3 − y1) K + 2 y1D] /2y′2 = −[(π y2 y4 + π y1 y3 − y2) K + 2 y2D] /2y′3 = 2 π y1 y2K − 4 y3Dy′4 = (π y2

2 − π y21) K − 4 y4D

(N.1)

N.2 Compiti

(a) Analizzarne la dinamica rispetto al parametro K, posto D = 1.

(b) Date le condizioni iniziali

y1(0) = ǫ, y2(0) = ǫ, y3(0) = 0, y4(0) = 0 (N.2)

e postoT = 20, D = 1, K = 4, ǫ = 10−6,

calcolare le soluzioni del sistema che si ottengono usando almeno duediversi metodi alle differenze finite, fra i quali:

• Runge-Kutta del secondo ordine;

• rappresentazione con differenze in avanti delle derivate e “conge-lamento” dei secondi membri agli istanti in cui le quantita sononote.

Tabulare e graficare le soluzioni numeriche ottenute.

167

Appendice O

Sistemi Hamiltoniani

O.1 Introduzione

La risoluzione di problemi riguardanti la dinamica 1 di sistemi meccaniciinizio con Galilei(1638) e con i Principia di Newton(1687). Newton ridussela descrizione del movimento di un punto-massa ad un sistema di equazionidifferenziali.

Lo studio del movimento di sistemi complessi (corpi rigidi o collegatitramite corde o aste) comportava serie difficolta, fino a quando Lagrangetrovo un modo elegante di trattarlo [Gan75]. Supponiamo di avere un sistemameccanico a d gradi di liberta, detto sistema Lagrangiano, la cui posizione edescritta dalle quantita q = (q1, ..., qd), chiamate coordinate generalizzate.

Lagrange usa le quantita

T = T (q, q), U = U(q), (O.1)

che sono rispettivamente l’energia cinetica e l’ energia potenziale del sistema.La funzione

L = T − U (O.2)

viene chiamata il Lagragiano del sistema. Le coordinate q1(t), ..., qd(t) obbe-discono alle equazioni differenziali:

d

dt(∂L

∂q) =

∂L

∂q, (O.3)

che vengono dette le equazioni Lagragiane del sistema. La risoluzione diqueste equazioni permette di descrivere il movimento di un tale sistema, apartire da condizioni iniziali date.

1Sezione sviluppata con la collaborazione di E. Mion.

168

m

lφ

Figura O.1: Pendolo dell’esempio 1.

Esempio O.1.1 Supponiamo di avere un pendolo (Fig.O.1) di massa m,collegato ad un’asta rigida senza massa di lunghezza l. E’ un sistema ad unsolo grado di liberta. La sua energia cinetica e

T =1

2ml2 φ2 (O.4)

mentre l’energia potenziale e:

U = mgz = mg l · cosφ. (O.5)

Il Langragiano del problema e:

L = T − U =1

2ml2 φ2 −mg l · cosφ. (O.6)

Hamilton nel 1934 semplifico la struttura delle equazioni di Lagrange etrasformo il problema in una forma che possiede una interessante simmetria.Egli

• introdusse le variabili di Poisson, ricavando i momenti coniugati ogeneralizzati

pk =∂L

∂qk(q, q), k = 1, . . . , d. (O.7)

169

• considero una nuova funzione, chiamata Hamiltoniana

H := pT q − L(q, q). (O.8)

E’ possibile dimostrare che le equazioni di Lagrange sono equivalenti alleequazioni di Hamilton

pk = −∂H∂qk

(p, q), qk =∂H

∂pk(p, q), k = 1, . . . , d. (O.9)

Esempio O.1.2 Riprendendo l’esempio 1 possiamo scrivere:

pφ =∂L

∂φ= ml2φ (O.10)

e risolvendo rispetto a φ otteniamo

φ =pφml2

. (O.11)

L’ hamiltoniana del problema e:

H(pφ, φ) =p2φ

2ml2−mgl · cosφ (O.12)

Ponendo p = pφ e q = φ otteniamo

H(p, q) =p2

2ml2−mgl · cos q. (O.13)

Possiamo dunque scrivere le equazioni del moto

p = −∂H∂q

(p, q) = −mg l sin q (O.14)

q =∂H

∂p(p, q) =

p

m l2. (O.15)

O.1.1 Compiti

Supponendo che la la massa m del pendolo sia 1 Kg, che la lunghezza dell’astasia l = 1 m e che l’ accelerazione di gravita agente sul punto sia g = 1 m/s2,determinare il moto del pendolo utilizzando i metodi di Eulero Esplicito e

170

Eulero Implicito 2. Calcolare il moto per i valori h = 0.2 s e (p0, q0) = (0, 0.5).Disegnare due grafici che mostrano l’evoluzione temporale del sistema nellecoordinate p e q, 0 ≤ t ≤ 100 s. Il piano (p, q) viene chiamato spazio dellefasi (phase space). Disegnare anche i grafici che mostrano l’andamento di φin funzione di t.

O.2 Sistemi Simplettici

I sistemi Hamiltoniani hanno una importantissima proprieta, chiamata sim-pletticita del flusso.

Per illustrarne il significato, consideriamo un parallelogramma che giacesu un piano. Il parallelogramma e individuato da due vettori

ξ = (ξp, ξq)T , η = (ηp, ηq)T , (O.16)

ed e formato dai punti

P = tξ + sη t.c. 0 ≤ s, t ≤ 1. (O.17)

Nel caso di dimensione d = 1, definiamo l’area orientata, Ao del paralle-logramma

Ao(P ) = det

(ξp ηp

ξq ηq

)= ξpηq − ξqηp. (O.18)

Quando d > 1, consideriamo la somma delle aree orientate della proie-zione di P nelle coordinate (pi, qi)

ω(ξ, η) :=

d∑

i=1

det

(ξpi ηpiξqi ηqi

)=

d∑

i=1

(ξpi ηqi − ξqi ηpi ). (O.19)

Questa e una mappa bilineare che agisce sui vettori di R2d. In notazionematriciale la mappa ha la forma

ω(ξ, η) = ξTJη, J = (0 I−I 0

), (O.20)

dove I e la matrice identica.

2Ricordiamo che i due schemi sono

yn+1 = yn + hf(yn) (Eulero Esplicito)

yn+1 = yn + hf(yn+1) (Eulero Implicito)

Per applicare quest’ultimo metodo al problema del pendolo, bisogna risolvere ad ogni passoun’ equazione non lineare. Usate lo schema di Newton–Raphson per farlo.

171

Definizione O.2.1 Una funzione lineare A : R2d → R2d e chiamata sim-plettica se per ogni ξ, η ǫ R2d:

ATJA = J, (O.21)

ossia ω(Aξ,Aη) = ω(ξ, η).

Nel caso d = 1, possiamo dire che l’area del parallelogramma si preserva.In generale, la simpletticta garantisce che la somma delle aree orientate delleproiezioni di P su (pi, qi) e uguale a quella del parallelogramma trasformatoA(P ).

Dato un sistema del tipo:

u = a(u, v) (O.22)

v = b(u, v), (O.23)

possiamo calcolarne una soluzione numerica tramite lo schema

un+1 = un + h a(un+1, vn), (O.24)

vn+1 = vn + h b(un+1, vn). (O.25)

Questo schema viene chiamato Metodo Simplettico di Eulero [HLW02]. Siprova che questo schema preserva la simpletticita, mentre Eulero Esplicito edEulero Implicito non la preservano. Questo comporta che la dinamica a lun-go termine del sistema modellizzato non viene correttamente approssimatausando gli ultimi due schemi.

O.2.1 Compiti

Risolvere l’esercizio del pendolo utilizzando il Metodo Simpettico di Eulero.Utilizzare i valori h = 0.3 s, p0 = 0 e q0 = 0.7, 1.4, 2.1. Disegnare ungrafico che mostra l’ evoluzione temporale del sistema nelle coordinate p e q,0 ≤ t ≤ 100 s. Disegnare anche i grafici che mostrano l’andamento di φ infunzione di t. Confrontare i nuovi grafici con quelli ottenuti precedentemente.Qual e piu corretto?

172

Appendice P

Test dello Shock Tube

P.1 Descrizione del problema

Supponiamo di avere un tubo (vedi figura P.1) diviso in due parti uguali dauna membrana1. Nella meta sinistra e contenuto un fluido l, caratterizzatoda alta densita, dl, e pressione, pl. Nella meta destra, invece, vi e fluido r, didensita dr ≪ dl, e pressione pr ≪ pl.

Al tempo t0 = 0, la membrana viene rimossa ed il fluido l inizia a dif-fondere entro il fluido r, originando uno shock che si propaga verso la partedestra del tubo, mentre un’onda, chiamata onda di rarefazione, si propagain direzione opposta, alla velocita locale del suono.

Si vengono a creare cinque distinte regioni (vedi figura P.2):

• il fluido imperturbato l;

• l’onda di rarefazione;

• una regione di densita e pressione costante (interfaccia);

• lo shock;

• il fluido imperturbato r.

1Sezione sviluppata con la collaborazione di C. Gheller, R. Favaretto, A. Missio, M.Guidolin.

173

0

Left zone Right zone

dL ≫ dR

pL, dL pR, dR

pL ≫ pR

Figura P.1: Rappresentazione del tubo.

−2 −1.5 −1 −0.5 0 0.5 1 1.5 20.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

posizione

dens

ità

Shock Tube − Densità

ZONA SINISTRA

ONDA DIRAREF.

POSTRAREF.

POSTSHOCK

ZONADX

SHOCK

TESTA

CODA

Figura P.2: Regioni che si vengono a creare dopo la rimozione dellamembrana.

174

P.2 Sistema fluidodinamico ideale

L’evoluzione di un sistema fluidodinamico ideale in una dimensione vienedescritta dalle seguenti equazioni idrodinamiche [CF48, LL87]:

∂ρ

∂t+∂[ρv]

∂x= 0, (P.1)

∂ρv

∂t+∂[ρv2 + p]

∂x= 0, (P.2)

∂E

∂t+∂[(E + p)v]

∂x= 0, (P.3)

e dall’equazione di stato:

E =p

γ − 1+

1

2ρv2, (P.4)

dove E indica l’energia, ρ la densita, p la pressione e γ e la costante adiaba-tica; γ = 5/3 per un fluido barionico non relativistico.

Ponendo w1 = ρ, w2 = ρ v, w3 = E, e f1 = ρ v = w2, f2 = ρv2 + p,f3 = (E + p)v, le equazioni si possono riscrivere

∂wk∂t

+∂fk∂x

= 0, k = 1, 2, 3. (P.5)

Per determinare l’ evoluzione del sistema, bisogna quindi risolvere leequazioni (P.5).

P.3 Metodi di integrazione numerica

Descriveremo due metodi numerici, lo schema di Lax e quello di Lax-Wendroff,che useremo per risolvere il test dello shock tube.

Entrambi gli schemi fanno uso di metodi di approssimazione alle differenzefinite.

Sia xi la coordinata dell’i-esimo punto della griglia.Per determinare l’evoluzione temporale di wk possiamo considerarne lo

sviluppo in serie di Taylor rispetto al tempo di riferimento t:

wk(xi, t) = wk(xi, t) +

[∂wk(xi, t)

∂t

]

t=t

(t− t) + (P.6)

+1

2

[∂2wk(xi, t)

∂t2

]

t=t

(t− t)2 + ... (P.7)

175

Siano: tn l’n-esimo istante temporale; t = tn; ∆t = tn+1 − tn.Scrivendo la (P.7) al tempo tn+1 otteniamo:

wk(xi, tn+1) = wk(xi, tn) +

[∂wk(xi, t)

∂t

]

t=tn

(t− tn) + (P.8)

+1

2

[∂2wk(xi, t)

∂t2

]

t=tn

(t− tn)2 + ... (P.9)

Usando la (P.5), otteniamo

wk(xi, tn+1) = wk(xi, tn)−[∂fk(x, t)

∂x

]

x=xi

∆t− (P.10)

+1

2

[∂

∂t

[∂fk(x, t)

∂x

]

x−xi

]

t=tn

∆t2 + ... (P.11)

Trascurando i termini di ordine superiore al primo e approssimando laderivata con le differenze centrali, otteniamo

wk(xi, tn+1) = wk(xi, tn)−fk(xi+1, tn)− fk(xi−1, tn)

2∆x∆t (P.12)

P.4 Integrazione temporale

La scelta di un ∆t opportuno si rivela estremamente delicata dato che, se egrande si ha instabilita numerica, se e troppo piccolo, il costo computazionalediventa alto e si verificano altri problemi numerici.

Consideriamo una griglia unidimensionale caratterizzata da N punti spa-ziali, definita sull’intervallo [0, L], xi = i∆x, ∆x = L/N , i = 0, . . . , N .

Se |v| e la velocita con cui si propaga la perturbazione, per garantire sta-bilita numerica, il passo temporale deve soddisfare alla cosiddetta condizionedi Courant

∆t ≤ C∆x

|v| . (P.13)

dove 0 < C ≤ 1 e una costante, chiamata coefficiente di Courant.In questo modo durante il tempo ∆t il fenomeno non si puo propagare piu

di ∆x nello spazio, quindi la sua variazione e abbastanza limitata. Poniamo

α =∆x

|vm|, |vm| = max

i|vi|. (P.14)

Poiche possiamo avere velocita piccole nella simulazione, fissiamo unvalore β > 0 che dipende dal problema in esame, e poniamo

∆t =

Cα, se α < β,Cβ, α ≥ β

(P.15)

176

P.4.1 Metodo di Lax

Il metodo di Lax e uno schema alle differenze finite del primo ordine. Apartire dalla eq. (P.12) si pone

wk(xi, tn+1) =wk(xi+1, tn) + wk(xi−1, tn)

2− (P.16)

+fk(xi+1, tn)− fk(xi−1, tn)

2∆x∆t (P.17)

In letteratura si prova che pregi di questo metodo sono la limitata spesacomputazionale e l’assenza di fenomeni come oscillazioni non fisiche o laperdita di positivita di grandezze che fisicamente sono positive.

P.4.2 Il metodo di Lax-Wendroff

Questo e uno schema alle differenze finite del secondo ordine. L’eq. (P.7) sipuo riscrivere

wk(tn+1) = wk(tn) +

[∂wk∂t

]

t=t

∆t+1

2

[∂2wk∂t2

]

t=t

∆t2 (P.18)

Usando l’equazione (P.5) possiamo scrivere:

∂2wk∂t2

= −∂∂t

∂fk∂x

= −∂∂t

(∂x

∂t

∂fk∂x

). (P.19)

Ma dall’equazione (P.5) si ha anche che:

∂x

∂t= − ∂fk

∂wk. (P.20)

Sostituendo nella (P.19) otteniamo

∂2wk∂t2

=∂

∂x

(∂fk∂wk

∂fk∂x

). (P.21)

Possiamo, quindi, riscrivere la (P.18) come

wk(tn+1) = wk(tn)−∂fk∂t

∆t+1

2

∂

∂x

(∂fk∂wk

∂fk∂x

)∆t2. (P.22)

E stato dimostrato che tale approssimazione puo essere sostituita da uno sche-ma a due step, detto schema di Lax-Wendroff-Richtmyer [RM67]. Poniamo

w(i)k,n := wk(xi, tn).

177

Quantita Dimensioni (SI) Parte sinistra Parte destradensita M/L 1.0 0.125velocita L/T 0.0 0.0pressione M/T 2 2/3 2/30

Tabella P.1: Condizioni iniziali.

1. Si calcolano i valori delle grandezze ausiliarie w(i)k al tempo tn+1/2 =

tn + ∆tn/2, ∆tn = tn+1 − tn,

w(i+1)k,n+1/2 =

1

2

[w

(i+1)k,n + w

(i)k,n

]− ∆t

2∆x

[f

(i+1)k,n − f (i)

k,n

]. (P.23)

A partire dai w(i)k,n+1/2, si calcolano i valori f

(i)k,n+1/2. Ad esempio, f1 =

ρ v, quindif

(i)1,n+1/2 = ρ

(i)n+1/2 v

(i)n+1/2, i = 0, . . . , N.

2. All’ istante tn+1 abbiamo

w(i)k,n+1 = w

(i)k,n −

∆t

∆x

(f

(i+1)k,n+1/2 − f

(i)k,n+1/2

). (P.24)

Pregio di questo metodo e la maggior accuratezza rispetto al metodo diLax. In particolare, e molto contenuto il fenomeno della diffusione numerica,anche se, in prossimita delle zone in cui compaiono forti gradienti, compaionodelle oscillazioni non fisiche. La spesa rispetto allo schema di Lax e maggiore,sia come occupazione di memoria (si devono mantenere in memoria le gran-dezze al tempo tn+1/2), sia come costo computazionale (l’algoritmo richiededue steps ad ogni passo temporale).

P.5 Compiti

Implementare gli schemi di Lax e Lax–Wendroff, confrontando i risultati conquelli analitici, forniti dal codice Matlab shock.m, reperibile nella cartella$CN/Esercitazioni/Shocktube.

Porre γ = 5/3, ∆t = 0.05. Le condizioni iniziali sono riportate nellatabella P.1. L’energia interna iniziale a destra e a sinistra si puø calcolareusando l’ equazione (P.4). Effettuare i confronti negli istanti temporali t =1, 2, 4, 8, 16, ponendo alternativamente C = 0.9, 0.8, 0.6. Il passo temporalesoddisfa la condizione di Courant?

Quanto vale α?

178

Bibliografia

[AM88] J. C. Alexander e J. H. Maddocks. On the maneuvering ofvehicles. SIAM J. Appl. Math., 48:38–51, 1988.

[Atk78] K. E. Atkinson. An introduction to Numerical Analysis. JohnWiley & Sons, New York, 1978.

[Bal99] S. Balsamo. Sistemi di elaborazione dell’ informazione: va-lutazione delle prestazioni di sistemi. Dispense per il cor-so di Sistemi di Elaborazione dell’ Informazione, Universitadi Venezia, Corso di Laurea in Informatica. Disponibile all’URL http://www.dsi.unive.it/∼balsamo/dispense.html,ultimo accesso: 17 sttembre 2007, 1999.

[BCM88] D. Bini, M. Capovani, e O. Menchi. Metodi Numerici per l’Algebra Lineare. Zanichelli, Bologna, 1988.

[BDJ99] M. W. Berry, Z. Drmac, e E. R. Jessup. Matrices, vector spacesand information retrieval. SIAM Review, 41(2):335–362, 1999.

[Ber94] M. Bertsch. Istituzioni di Matematica. Bollati Boringhieri,Torino, 1994.

[BR02] A. Brandt e D. Ron. Multigrid solvers and multileveloptimization strategies. Manuscript, 2002.

[Bra77] A. Brandt. Multi-level adaptive solutions to boundary-valueproblems. Math. Comp., 31(138):333–390, April 1977.

[Bra89] A. Brandt. Rigorous local mode analysis of multigrid. In J. Man-del, S. F. McCormick, J. E. Dendy, C. Farhat, G. Lonsdale, S. V.Parter, J. W. Ruge, e K. Stuben, (a cura di), Proceedings ofthe Fourth Copper Mountain Conference on Multigrid Methods,pagina 438, Philadelphia, PA, USA, 1989. SIAM.

179

[Bra94] A. Brandt. Rigorous quantitative analysis of multigrid, I: Con-stant coefficients two-level cycle with L2-norm. SIAM J. Num.Anal., 31(6):1695–1730, December 1994.

[BRS02] S. C. Brenner e L. Ridgway Scott. The Mathematical Theory ofFinite Element Methods. Springer Verlag, new York, ii edizione,2002.

[Car91] G. Cariolaro. La teoria unificata dei segnali. UTET, Torino,1991.

[CF48] R. Courant e K. O. Friedrichs. Supersonic Flow and ShockWaves. Interscience, New York, 1948.

[CLRS05] T. H. Cormen, C. E. Leiserson, R. L. Rivest, e C. Stein.Introduzione agli algoritmi e strutture dati. McGraw-Hill, 2005.

[Com91] V. Comincioli. Analisi Numerica. McGraw-Hill Italia, Milano,1991.

[D96] W. Dorfler. A convergent adaptive algorithm for Poisson’s equa-tion. SIAM Journal of Numerical Analysis, 33(3):1106–1124,1996.

[Dem97] J. W. Demmel. Applied Numerical Linear Algebra. SIAM,Philadelphia, PA, 1997.

[EEHJ96] K. Eriksson, D. Estep, P. Hansbo, e C. Johnson. ComputationalDifferential Equations. Cambridge University Press, CambridgeMA, 1996.

[FvDFH90] J. D. Foley, A. van Dam, S. K. Feiner, e J. F. Hughes. ComputerGraphics. Principles and Practice. Addison-Wesley, Reading,MA, ii edizione, 1990.

[Gam94] G. Gambolati. Lezioni di Metodi Numerici per Ingegneria eScienze Applicate. Cortina, Padova, 1994.

[Gan75] F. R. Gantmacher. Lectures in Analytical Mechanics. MIRPublishers, Moscow, 1975.

[Gol91] D. Goldberg. What every computer scientist should know aboutfloating-point arithmetic. ACM Computing Surveys, 23(1):5–48,1991.

180

[Gv89] G. Golub e C. F. van Loan. Matrix Computation. The JohnsHopkins University Press, Baltimore, ii edizione, 1989.

[Hab98] R. Haberman. Mathematical Models. SIAM, Philadelphia,PA, 1998. Unabridged republication of Prentice-Hall book,Englewood Cliffs, NJ, 1977.

[HLW02] E. Hairer, C. Lubich, e G. Wanner. Geometric NumericalIntegration. Springer-Verlag, Berlin, 2002.

[HS93] W. Huang e D. M. Sloan. A new pseudospectral method withupwind features. IMA J. Num. Anal., 13:413–430, 1993.

[JFH+98] P. J. Jessee, W. A. Fiveland, L. H. Howell, P. Colella, e R. B.Pember. An adaptive mesh refinement algorithm for the ra-diative transport equation. Journal of Computational Physics,139:380–398, 1998.

[KA03] P. Knaber e L. Angermann. Numerical Methods for Elliptic andParabolic Partial Differential Equations. Springer, New York,2003.

[LL87] L. D. Landau e E. M. Lifschitz. Fluid Mechanics. Oxford Univ.Press & Butterworth–Heinemann, Oxford, 1987.

[Lyo95] L. Lyons. All you wanted to know about mathematics, volume 1.Cambridge University Press, Cambridge, 1995.

[MC96] D. Martin e K. Cartwright. Solving poisson’s equation usingadaptive mesh refinement. available on internet, october 1996.

[Min96] M. L. Minion. A projection method for locally refined grids.Journal of Computational Physics, 127:158–178, 1996.

[Mv96] J. D. Murray e W. vanRyper. Encyclopedia of Graphics FileFormats. O’Reilly and Associates, Inc., Sebastopol CA, secondaedizione, 1996.

[Nie03] Y. Nievergelt. Scalar fused multiply-add instructions producefloating-point matrix arithmetic provably accurate to the pe-nultimate digit. ACM Transactions on Mathematical Software,29(1):27–48, March 2003.

[Par80] B. N. Parlett. The Symmetric Eigenvalue Problem. Prentice-Hall, Englewood Cliffs NJ, USA, 1980.

181

[PH05] D. A. Patterson e J. L. Hennessy. Computer organization anddesign. The hardware/software interface. Elsevier, 2005.

[PZ92] G. Pini e G. Zilli. Esercizi di Calcolo Numerico. Imprimitur,Padova, 1992.

[QS06] A. Quarteroni e F. Saleri. Introduzione al Calcolo Scientifico.Springer Verlag Italia, iii edizione, 2006.

[RM67] R. D. Richtmyer e K. W. Morton. Difference Methods for Initial-Value Problems. Interscience, New York, ii edizione, 1967.

[SB80] J. Stoer e R. Bulirsch. Introduction to Numerical Analysis.Springer-Verlag, New York, 1980.

[SM83] G. Salton e M. J. McGill. Introduction to Modern InformationRetrieval. McGraw-Hill, New York, 1983.

[SP02] F. Sartoretto e M. Putti. Introduzione alla Programmazioneper elaborazioni numeriche. Edizioni Libreria Progetto, Padova,2002.

[Sta04] W. Stallings. Architettura e organizzazione dei calcolatori.Pearson, Milano, sesta edizione, 2004.

[The06] The MathWorks Inc., Natick, MA. MATLAB. Com-plete documentation., 2006. Reperibile all’ URL:http://www.mathworks.com/access/helpdesk/help-

/techdoc/matlab.html, ultimo accesso: 1 settembre2006.

[TOS00] U. Trottenberg, C. W. Oosterlee, e A. Schuller. Multigrid.Academic Press, London, 2000.

[TS64] A. N. Tychonov e A. A. Samarsky. Partial Differential Equationsof Mathematical Physics. Holden-Day, San Francisco, 1964.

[Wei01] R. Weinands. Extended Local Fourier Analysis for Mul-tigrid: Optimal Smoothing, Coarse Grid Correction andPreconditioning. PhD thesis, University of Koln, 2001.

[Wes92] P. Wesseling. An Introduction to Multigrid Methods. John Wileyand Sons Ltd, New York, 1992. Corrected Reprint. Philadelphia:R. T. Edwards Inc., 2004.

182

[Wik08] Wikipedia. Logistic function. Reperibile all’ URL http://en.

wikipedia.org/wiki/Logistic_function, Ultimo accesso: 8marzo 2008.

183