Calcolo numerico

Calcolo numerico

1

Indice

1. Il problema generale ..................................................................................................................................... 4

2. Metodi numerici per sistemi lineari .............................................................................................................. 5

2.1 Fattorizzazione LU ............................................................................................................................................... 7

2.2 Metodi iterativi ................................................................................................................................................... 9

2.3 Stabilità di sistemi lineari .................................................................................................................................. 15

3. Problema agli autovalori ............................................................................................................................. 18

3.1 Metodo delle potenze ...................................................................................................................................... 18

3.2 Metodo delle potenze inverse ......................................................................................................................... 20

3.3 Metodo delle potenze inverse con shift .......................................................................................................... 21

4. Determinazione delle radici di equazioni non lineari ................................................................................... 22

4.1 Metodo di bisezione ......................................................................................................................................... 22

4.2 Metodo delle corde .......................................................................................................................................... 24

4.3 Metodo delle secanti ........................................................................................................................................ 24

4.4 Metodo di Newton............................................................................................................................................ 25

4.5 Problema di punto fisso .................................................................................................................................... 25

4.5 Criteri d’arresto ................................................................................................................................................. 27

4.6 Sistemi di equazioni non lineari ....................................................................................................................... 28

5. Approssimazione di funzioni e dati ............................................................................................................. 30

5.1 Polinomio di interpolazione di Lagrange ......................................................................................................... 30

5.2 Interpolazione di funzioni ................................................................................................................................. 31

5.3 Interpolazione composita ................................................................................................................................. 34

5.4 Spline cubiche ................................................................................................................................................... 35

5.5 Approssimazione ai minimi quadrati ............................................................................................................... 37

6. Calcolo dell’integrale di ........................................................................................................................... 39

6.1 Calcolo approssimato dell’integrale di ......................................................................................................... 39

6.2 Calcolo approssimato dell’integrale di con aree trapezoidali ..................................................................... 40

+

6.3 Calcolo approssimato dell’integrale di a punti ................................................................................. 40

6.4 Formula di Cavalieri-Simpson ........................................................................................................................... 42

6.5 Andamento dell’errore di quadratura ............................................................................................................. 42

6.6 Grado di esattezza e ordine di accuratezza ..................................................................................................... 44

6.7 Formule di quadratura di Gauss ....................................................................................................................... 44

7. Equazioni differenziali ordinarie .................................................................................................................. 46

7.1 Modello differenziale ........................................................................................................................................ 46

7.2 Problema di Cauchy .......................................................................................................................................... 46

7.3 Metodi numerici per equazioni differenziali ordinarie ................................................................................... 48

7.4 Metodo di Eulero esplicito ............................................................................................................................... 48

7.5 Metodo di Eulero implicito ............................................................................................................................... 49

7.6 Metodo di Crank-Nicholson ............................................................................................................................. 50

2

7.7 Andamento dell’errore di approssimazione .................................................................................................... 50

7.8 Assoluta stabilità ............................................................................................................................................... 53

7.9 Problema modello di assoluta stabilità ............................................................................................................ 54

7.10 Valutazione dell’assoluta stabilità .................................................................................................................. 54

7.11 Metodi di Runge-Kutta ................................................................................................................................... 56

7.12 Metodo di Heun .............................................................................................................................................. 58

7.13 Metodo di Runge-Kutta di ordine 4 esplicito ................................................................................................ 58

3

1. Il problema generale

Il punto di partenza è un problema fisico da risolvere, che viene convertito in un problema

matematico mediante la descrizione del problema stesso secondo equazioni, semplificando il

problema, comportando l’introduzione dell’errore di modello. Nel momento in cui la soluzione

analitica del problema matematico non sia ottenibile bisogna ricorrere al problema numerico,

in grado di fornire una soluzione approssimata e, quindi, che introduce un ulteriore errore,

l’errore numerico. La soluzione numerica deve essere il più prossima possibile a quella reale

(comunque incognita).

Ad esempio, nel calcolo di area sottesa da una curva, l’integrale può essere calcolato

numericamente al calcolatore utilizzando il metodo dei rettangoli; con il numero di rettangoli

→ ∞ l’errore tende a zero e la soluzione numerica tende a quella esatta. L’errore numerico,

quindi, è un errore controllabile: aumentando il numero di intervalli N è possibile far decrescere

l’errore. errore

L’ottenimento finale della soluzione numerica comporta un ulteriore errore, detto

computazionale , che dipende dall’aritmetica finita con cui lavora il calcolatore: i numeri sono

rappresentabili con una quantità finita di cifre, fatto che introduce un errore nel caso di numeri

irrazionali, come il π, che nel caso di una somma si possono amplificare, portando ad una

soluzione errata. L’errore computazionale, infatti, non è controllabile e può pregiudicare il

funzionamento di un algoritmo tecnicamente corretto.

4

2. Metodi numerici per sistemi lineari

Dati una matrice ∈ ℝ e un vettore ∈ ℝ si vuole determinare il vettore ∈ ℝ tale che

= (1)

Il problema (1) deve essere ben posto, ovvero la soluzione deve esistere unica. Ciò è vero nel

momento in cui vale almeno una delle condizioni seguenti, tutte equivalenti fra loro.

i () =

ii det() ≠ 0

−1

iii ∃

La matrice A deve avere rango massimo e deve essere invertibile. A questo punto, la soluzione

analitica è data dalla (2). −1

= (2)

−1

Il calcolo della matrice inversa , tuttavia, è un problema complesso che crea errori

computazionali non controllabili: la soluzione numerica non può essere ottenuta mediante

soluzione dell’equazione (2).

metodo di Cramer

Il fornisce un diverso algoritmo per il calcolo della soluzione analitica, grazie

alla regola dei determinanti, tale per cui )

det(

= (3)

det ()

ove è la matrice a cui viene sostituita l’-esima riga con il vettore termine noto . Questo

algoritmo, tuttavia, impiega 3( − 1)! operazioni elementari per calcolare la soluzione e,

dunque, al crescere di il tempo di calcolo aumenta in modo proibitivo: un calcolatore in grado

9

di eseguire 10 operazioni elementari al secondo, ovvero da 1 gigaFLOPS , con = 15 è in grado

1

di calcolare la soluzione in circa 12 ore. Anche questo metodo, dunque, non è adatto alla

risoluzione del problema numerico.

Caso particolare: si considera una matrice triangolare inferiore L (5), tale che

= (4)

Per cui deve valere det() ≠ 0 ⇔ ≠ 0 = , in modo che esista l’inversa di L.

1 Per FLOPS si intende FLoating point Operations Per Second, ovvero il numero di operazioni a virgola

mobile eseguite in un secondo 5

0 0 0 0

11

0 0 0

21 22

= ⋮ ⋮ ⋱ 0 0 (5)

⋮ ⋮ ⋮ ⋱ 0

[ ⋮ ⋮ ⋮ ⋮ ]

Nel caso di matrice triangolare inferiore, tutti i termini al di sopra della diagonale sono nulli.

L’equazione (4), quindi, è facilmente risolvibile, infatti

1

= → = (6)

11 1 1 1

11

1 ( )

+ = → = − (7)

21 1 22 2 2 2 2 21 1

22

−1

1

= − ∑

( ) (8)

=1

sostituzione in avanti

L’algoritmo di permette di ottenere le componenti del vettore soluzione.

for

Al calcolatore, l’algoritmo può essere implementato con un semplice ciclo . Questo algoritmo,

2

inoltre, impiega solamente 0.5 ∙ operazioni e dunque permette di ottenere una soluzione in

tempi modesti, anche con elevato.

Il caso di matrice triangolare superiore è analogo e risolvibile mediante il duale algoritmo di

sostituzione all’indietro . −1

1

() = − ∑ (, )()

(() )

(, ) =0

−

1

() = − ∑ (, − )( − )

(() )

(, ) =0

La generica matrice A, dunque, può essere ricondotta al prodotto fra una matrice triangolare

inferiore L e una matrice triangolare superiore U.

= (9)

Le matrici L e U devono essere invertibili. A questo punto si deve risolvere il sistema di seguito,

equivalente alla (1) =

{ (10)

=

2

La risoluzione del sistema costa circa flops, ma deve essere preceduta dalla fattorizzazione

3

LU (equazione 9), che costa circa operazioni.

6

Si noti che per la matrice L deve valere la condizione = 1, ∀ = 1, . . , per avere un’ottima

allocazione di memoria. Per sistemi ad ordine molto elevato, infatti, serve una gran quantità

di memoria disponibile ad allocazione per risolvere il problema; l’allocazione viene ottimizzata

sovrascrivendo alla matrice A gli ingressi delle matrici L e U: per fare ciò è necessario che la

diagonale di L non contenga informazioni, in quanto è disponibile spazio per salvare i dati di una

sola matrice, sostituendo la diagonale di A.

2.1 Fattorizzazione LU

Data la matrice ∈ ℝ , la matrice triangolare superiore viene costruita in modo univoco e

con passaggi. Ottenuta , ∃! tale che = , ovvero è garantita l’unicità di entrambe le

matrici.

Il metodo di fattorizzazione LU sfrutta il metodo di eliminazione di Gauss, che si avvale di tre

semplici operazioni applicabili alle righe della matrice , senza alterare il problema, ovvero

swap

moltiplicazione di una riga per uno scalare, operazioni di somma algebrica tra le righe e di

righe.

Esempio 1. Fattorizzazione LU della seguente matrice quadrata di ordine = 3.

1 1 3 1 1 3 1 1 3

= → →

[ ] [ ] [ ]

2 3 5 0 −1 −1 0 −1 −1 (11)

=−2 =−

7 8 9 0 −1 −12 0 0 −11

=−7

Opportuni passaggi permettono di ottenere una matrice triangolare superiore con = 3

passaggi, ovvero si è ricavata la matrice 1 1 3

= [ ]

0 −1 −1 (12)

0 0 −11

L’algoritmo, per funzionare in casi più generali, deve prevedere il pivoting: se al -esimo passo

dell’algoritmo il termine diagonale è nullo, è necessario scambiare la -esima riga con una riga

successiva. L’algoritmo più efficiente prevede lo scambio fra la -esima riga e la riga avente

l’elemento pivotale a modulo massimo, in modo da contenere gli errori computazionali.

L’algoritmo così modificato implementa il metodo di fattorizzazione LU con pivoting.

Dal punto di vista matriciale, scambiare di posizione due righe prevede l’introduzione della

matrice di permutazione , tale che = (13)

Nella risoluzione del problema (1), MATLAB calcola tre matrici, L, U e P, resituite come output.

[L, U, P] = lu(A)

Nel caso in cui non avvengano permutazioni la matrice P sarà uguale alla matrice identità.

7

= ⇔ = (14)

Come si vede dalla (14), la matrice P tiene conto del pivoting nel termine noto.

Teorema 1. Condizione necessaria e sufficiente a che esista unica la matrice .

Data la matrice ∈ ℝ , esistono uniche le matrici , ∶ = se e solo se tutte le

sottomatrici principali di A di ordine 1,2, … , sono non singolari, ovvero hanno determinante

diverso da zero.

Il Teorema 1 è di difficile dimostrabilità, per cui nella pratica si considerano delle condizioni

sufficienti, applicabili a classi particolari di matrici.

Teorema 2. Per matrici a predominanza diagonale per righe o per colonne e per matrici

simmetriche e definite positive, la matrice esiste ed è unica.

Definizione 1. La matrice A è simmetrica se = . +

Definizione 2. La matrice A è definita positiva se ∀ ∈ ℝ , ≠ 0, vale ∈ ℝ .

Definizione 3. La matrice A è simmetrica e definita positiva se i suoi autovalori sono tutti reali e

+

()

maggiori di zero, ovvero se ∈ ℝ , ∀ = 1, … , .

Definizione 4. La matrice A si dice a predominanza diagonale stretta per colonne nel momento

| |

in cui vale > | |, ≠ , ∀ = 1, … , .

Definizione 5. La matrice A si dice a predominanza diagonale stretta per righe, analogamente, se

| |

vale > | |, ≠ , ∀ = 1, … , .

3

Osservazione 1. La fattorizzazione LU è un algoritmo che costa circa 2 / 3 operazioni. Nel caso

in cui sia simmetria e definita positiva (sdp) è possibile utilizzare la fattorizzazione di

3

Cholesky, algoritmo risolubile utilizzando la metà dei flops, /3.

= (15)

Osservazione 2. La matrice non viene memorizzata nel calcolatore per intero, ma vengono

inseriti solamente i valori non nulli e i relativi indici di riga e colonna, in modo da ottimizzare la

scrittura della memoria allocabile.

Definizione 6. La matrice A di ordine è sparsa se ha circa (o più) elementi non nulli.

In MATLAB è possibile determinare gli elementi non nulli di una matrice usando la funzione

.

nnz(A)

Nel caso in cui la matrice sia un matrice sparsa, la fattorizzazione non è conveniente, in quanto

aumenta la memoria da allocare per il problema. In questi casi, quindi, i metodi diretti non sono

fill-in

una buona scelta: si ha il fenomeno di , per cui la matrice fattorizzata non è più sparsa. Con

problemi ad ordine molto grande, ovvero per dell’ordine di 10 o maggiore, avere una matrice

6

sparsa è una condizione fondamentale per garantire la corretta soluzione del problema: è

necessario calcolare la soluzione utilizzando un metodo senza fattorizzazione.

8

2.2 Metodi iterativi

I metodi iterativi implementano algoritmi di risoluzione di problemi = in cui è