Appunti di Calcolo Numerico

Appunti di Calcolo Numerico

Università degli Studi di Trento

Dipartimento di Ingegneria Civile, Ambientale e Meccanica

Matteo Franzoi

Ultimo aggiornamento: 8 marzo 2019

I seguenti appunti sono stati scritti durante il corso di Calcolo Numerico

tenuto dal Professor Vincenzo Casulli durante l’anno accademico 2017/2018.

Gli appunti possono contenere errori; nel caso il lettore ne riscontrasse è

invitato a comunicarli.

Indice

I Algebra lineare 1

Teorema 1.1 1

1 Metodo di eliminazione di Gauss 2

1.1 Metodo di eliminazione di Gauss . . . . . . . . . . . . . . . . 3

2 Sistemi tridiagonali 6

Teorema 1.2 7

Lemma 1 7

Lemma 2 7

2.1 Algoritmo di Thomas . . . . . . . . . . . . . . . . . . . . . . 10

Teorema 1.3 12

3 Sistemi non lineari 13

3.1 Metodo di Newton . . . . . . . . . . . . . . . . . . . . . . . . 14

−

3.2 Caso particolare: F (x) = A x b (sistema lineare) . . . . . . 20

3.2.1 Metodo di Jacobi . . . . . . . . . . . . . . . . . . . . . 20

3.2.2 Metodo di Gauss - Seidel o prima variante del metodo

di Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.2.3 Metodo S.O.R o seconda variante del metodo di Jacobi 21

4 Sistemi debolmente non lineari 23

Teorema 1.4 23

5 Metodi del gradiente 24

5.1 Metodo del gradiente semplice . . . . . . . . . . . . . . . . . 26

5.2 Metodo del gradiente coniugato . . . . . . . . . . . . . . . . . 27

6 Calcolo del raggio spettrale di una matrice 33

6.1 Metodo delle potenze . . . . . . . . . . . . . . . . . . . . . . . 33

II Approssimazione di funzioni 35

7 Funzioni discrete 35

8 Interpolazione lineare a tratti 35

3

9 Interpolazione parabolica a tratti 36

9.1 Interpolazione cubica a tratti . . . . . . . . . . . . . . . . . . 37

10 Interpolazione di Lagrange 37

10.1 Interpolazione parabolica a tratti . . . . . . . . . . . . . . . . 37

10.2 Interpolazione cubica a tratti . . . . . . . . . . . . . . . . . . 38

10.3 Generalizzando per polinomi di ordine k . . . . . . . . . . . . 39

11 Spline cubiche 39

1

11.1 Spline cubiche di classe C . . . . . . . . . . . . . . . . . . . 39

2

11.2 Spline cubiche di classe C . . . . . . . . . . . . . . . . . . . 41

12 Spline naturali 41

13 Minimi quadrati 43

III Approssimazioni di integrali 44

14 Metodo dei trapezi 45

Teorema 3.1 46

Teorema 3.2 47

15 Metodo di Simpson 48

Teorema 3.4 49

Teorema 3.5 50

16 Metodo di Romberg 50

16.1 Caso particolare . . . . . . . . . . . . . . . . . . . . . . . . . 52

17 Differenza finita in avanti 52

18 Differenza finita all’indietro 53

19 Differenza finita centrata 53

IV Problema ai valori iniziali 55

20 Metodo di Eulero 56

Lemma 4.1 57

4

Lemma 4.2 57

20.1 Errore nel metodo di Eulero . . . . . . . . . . . . . . . . . . . 58

Teorema 4.1 58

Teorema 4.2 59

20.2 Errori di arrotondamento nel metodo di Eulero . . . . . . . . 59

Teorema 4.3 59

21 Metodo di Runge – Kutta 60

21.1 I metodi di Runge – Kutta . . . . . . . . . . . . . . . . . . . 62

22 Serie di Taylor 64

22.1 Primo ordine . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

22.2 Secondo ordine . . . . . . . . . . . . . . . . . . . . . . . . . . 65

23 Instabilità 65

23.1 Condizioni di stabilità . . . . . . . . . . . . . . . . . . . . . . 66

23.2 Equazione alle differenze finite lineare a coefficienti costanti

del secondo ordine . . . . . . . . . . . . . . . . . . . . . . . . 67

24 Metodo di Eulero esplicito 70

25 Metodo di Eulero implicito 71

25.1 Sistema di equazioni lineari . . . . . . . . . . . . . . . . . . . 71

26 Soluzioni periodiche 73

Teorema 4.4 73

26.1 Metodo della secante . . . . . . . . . . . . . . . . . . . . . . . 74

V Problema ai limiti per equazioni differenziali ordinarie 75

27 Metodo alle differenze centrate 77

Teorema 5.1 80

Teorema 5.2 81

28 Metodo Upwind 83

28.1 Metodo Upwind generalizzato . . . . . . . . . . . . . . . . . . 85

Teorema 5.3 86

Teorema 5.4 86

5

29 Problemi debolmente non lineari (per problemi al contorno) 87

29.1 Differenze centrate . . . . . . . . . . . . . . . . . . . . . . . . 87

29.2 Metodo di Newton . . . . . . . . . . . . . . . . . . . . . . . . 89

29.3 Metodo Upwind . . . . . . . . . . . . . . . . . . . . . . . . . . 90

29.4 Convergenza delle differenze finite per problemi ai limiti . . . 91

Lemma 5.1 91

29.5 Stima dell’errore per il metodo alle differenze finite centrate

per problemi debolmente non lineari . . . . . . . . . . . . . . 93

Teorema 5.5 93

29.6 Convergenza per le differenze finite centrate per sistemi de-

bolmente non lineari . . . . . . . . . . . . . . . . . . . . . . . 95

Teorema 5.6 95

29.7 Metodo upwind per sistemi debolmente non lineari . . . . . . 95

Lemma 5.2 95

29.8 Stima dell’errore per il metodo upwind . . . . . . . . . . . . . 96

Teorema 5.7 96

29.9 Convergenza del metodo upwind . . . . . . . . . . . . . . . . 96

Teorema 5.8 96

29.10Condizioni al contorno . . . . . . . . . . . . . . . . . . . . . . 97

29.10.1 Condizioni al contorno miste . . . . . . . . . . . . . . 98

29.10.2 Reticolo . . . . . . . . . . . . . . . . . . . . . . . . . . 100

30 Metodo ibrido 101

30.1 Convergenza . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

31 Considerazioni 104

32 Metodo agli elementi finiti 105

33 Autovalori e autofunzioni 109

◦

VI Equazioni differenziali a derivate parziali del II ordi-

ne 111

34 Equazione ellittica 111

35 Equazione parabolica 111

36 Equazione iperbolica 112

Introduzione

È impossibile fare calcoli precisi con numeri reali (irrazionali) o nu-

meri razionali; questo perché hanno infinite cifre decimali che sarebbero

impossibili da scrivere tutte.

Per questo motivo si usa l’approssimazione.

Parte I

Algebra lineare

Si consideri un sistema lineare di n equazioni in n incognite

 · · ·

a x + a x2 + a x3 + + a x = b

11 1 12 13 1n n 1







 · · ·

a x + a x2 + a x3 + + a x = b

 21 2 22 23 2n n 2





 · · ·

a x + a x2 + a x3 + + a x = b ∀

a , b noti i, j = 1, 2, . . . , n

31 1 32 33 3n n 3 ij i

..



 .









 · · ·

a x + a x2 + a x3 + + a x = b

 n1 1 n2 n3 nn n n

che si può riscrivere nella forma matriciale associata

Ax = b (0.1)

dove A è la matrice dei coefficienti, x è il vettore delle incognite e b è il

vettore dei termini noti.

     

a a a ... a x b

11 12 13 1n 1 1

a a a ... a x b

21 22 23 2n 2 2

     

     

a a a ... a x b

A = , x = , b =

31 32 33 3n 3 3

     

     

.. .. ..

     

. . .

     

b

a a a . . . a x

n1 n2 n3 nn n n

Teorema 1.1. Il sistema A x = b ammette una e una sola soluzione se e

6

solo se det A = 0.

Nel caso in cui det A = 0 il sistema può avere infinite soluzioni oppure

nessuna.

Metodo di Cramer. det(A )

i ∀i

x = = 1, 2, . . . , n (0.2)

i det(A)

dove A è la matrice che si ottiene sostituendo alla i-esima colonna di A il

i

vettore dei termini noti b. 1

Dalla (0.2) si può notare la necessità che il denominatore sia non nullo.

Il metodo di Cramer risulta però dispendioso dal punto di vista del

calcolo per il fatto che calcolare n determinanti non è un processo facile.

Per un sistema a 30 equazioni in 30 incognite, infatti, il tempo per ri-

solvere il sistema con un calcolatore molto potente, attraverso il metodo di

15

Cramer, sarebbe ”infinito” (10 secoli).

1 Metodo di eliminazione di Gauss

Definizione. ∀j

1. la matrice A è triangolare inferiore se e solo se a = 0 > i;

ij ∀j

2. la matrice A è triangolare superiore se e solo se a = 0 < i;

ij

∀|i −

3. la matrice A è tridiagonale se e solo se a = 0 j| > 1 e

ij

6 ∀|i − ≤

a = 0 j| 1;

ij nj=1,

P |a |

|a | ≥

4. la matrice A è diagonale dominante se e solo se ij

ii j6 = i

∀i = 1, 2, . . . , n; T n

≥ ∀x ∈

5. la matrice A è definita positiva se e solo se x A x 0 e

R

T

risulta x A x = 0 nel solo caso in cui x = 0̄.

Matrice triangolare inferiore

Per una matrice triangolare inferiore il sistema (0.1) risulta



a x = b

11 1 1









a x + a x = b

 21 1 22 2 2







a x + a x2 + a x3 = b

31 1 32 33 3

..



 .









 · · ·

a x + a x2 + a x3 + + a x = b

 n1 1 n2 n3 nn n n

6

da cui, imponendo a = 0, si ricava

11 b

1

x =

1 a 11

Nel caso in cui risultasse a = 0 il sistema sarebbe singolare.

11

Proseguendo con la seconda equazione

−

b a x

2 21 1

x =

2 a 22

2

6

se e solo se a = 0.

22

È quindi facilmente risolvibile partendo dalla prima equazione e arrivan-

do alla n-esima equazione.

Matrice triangolare superiore

In una matrice triangolare superiore il sistema (0.1) risulta



a x = b

nn n n







a x + a x = b

 n−1,n n n−1,n−1 n−1 n−1

 ..

.







 · · ·

a x + a x + + a x + a x + a x = b

 1n n 1,n−1 n−1 13 3 12 2 11 1 1

Si procede in modo analogo alla matrice triangolare inferiore partendo,

invece, dalla n-esima equazione e finendo con la prima.

Si cerca, quindi, di semplificare i casi in modo che risulti una matrice

A triangolare superiore/inferiore. Questa semplificazione avviene tramite il

metodo di eliminazione di Gauss.

1.1 Metodo di eliminazione di Gauss

Il metodo di eliminazione di Gauss consiste in due diversi processi:

(i) eliminazione;

(ii) sostituzione.

Triangolare superiore

Per trasformare un sistema lineare del tipo (0.1) in un sistema triangolare

superiore, per prima cosa, si consideri non trasformato il sistema con apice

”(1)” tale che

 (1) (1) (1) (1) (1)

· · ·

a x + a x + a x + + a x = b

1 2 3 n

11 12 13 1n 1





 (1) (1) (1) (1) (1)

 · · ·

a x + a x + a x + + a x = b

 1 2 3 n

21 22 23 2n 2





 (1) (1) (1) (1) (1)

· · ·

a x + a x + a x + + a x = b

1 2 3 n

31 32 33 3n 3

..



 .







 (1) (1) (1) (1) (1)

 · · ·

a x + a x + a x + + a x = b

 nn n

1 2 3 n

n1 n2 n3 a

6 in modo

Si assuma a = 0 e si moltiplichi la prima riga per il fattore 21

11 a

11

tale da avere delle righe linearmente dipendenti.

3

Considerando il sistema con apice ”(2)” quello trasformato una volta, si

può scrivere  (1)

a (1)

(2) (1) 21

 − a

a = a



 12

22 22 (1)

 a



 11



 (1)

 a

 (2) (1) (1)

21

 −

a = a a



 23 23 13

(1)

a 11

..



 .







 (1)

 a

 (2) (1) (1)

 21

−

a = a a



 2n 2n 1n

(1)

 a

 11

o più in generale (1)

 a (1)

(2) (1) 21

−

 a

a = a

 1j

2j 2j

 (1)

 a

 11



 (1) 6 ∀j

a = 0, = 2, 3, . . . , n

11

(1)

 a



 (2) (1) (1)

21

−



b = b b

 2 2 1

 (1)

 a

11

Cosı̀ facendo si è eliminato la variabile x dalla seconda equazione, ottenendo

1

 (1) (1) (1) (1) (1)

· · ·

a x + a x + a x + + a x = b

1 2 3 n

11 12 13 1n 1





 (2) (2) (2) (2)

 · · ·

0 + a x + a x + + a x = b

 2 3 n

22 23 2n 2





 (2) (2) (2) (2)

· · ·

0 + a x + a x + + a x = b

2 3 n

32 33 3n 1

...











 (2) (2) (2) (2)

 · · ·

0 + a x + a x + + a x = b

 nn n

2 3 n

n2 n3

Per eliminare x da tutte le equazioni

1 (1)

 a (1)

(2) (1) i1 ∀j

−

 a = 2, 3, . . . , n

a = a

 1j

ij ij

 (1)

 a

 11



 ∀i = 2, 3, . . . , n

(1)

 a



 (2) (1) (1)

i1

−



b = b b

 1

i i

 (1)

 a 11

Volendo eliminare anche x , x , . . . x , si generalizza ancora

2 3 n

(k)

 a

(k+1) (k) (k)

ik

− ∀j



a = a a = k + 1, . . . , n

 ij ij kj

 (k)

 a

 kk



 ∀k −

= 1, 2, . . . , n 1 (1.1)

(k)

 a



 (k+1) (k) (k)

ik

 − ∀i

b = b b = k + 1, . . . , n

 i i k

 (k)

 a kk (k) 6

assumendo sempre il pivot a = 0.

kk 4

Esempio

Si consideri il seguente sistema lineare

(

x = 2

2

x = 1

1

La matrice dei coefficienti risulta

0 1

A = 1 0 −1 6

e il suo determinante è non nullo det(A) = = 0.

Risulta però impossibile applicare il metodo di Gauss perché non è

6

soddisfatta la condizione che a = 0.

11

È però possibile scambiare la prima e la seconda riga, in modo tale che la

matrice A risulti

1 0

A = 0 1

Scambiando le righe, il determinante non varia; in questo modo, però, è

possibile applicare il metodo di eliminazione di Gauss.

Tornando al caso generale, gli elementi della matrice A possono essere

numeri razionali o irrazionali e quindi soggetti ad errori.

(k) (k)

a a (k)

Se il fattore < 1 allora nell’operazione a l’errore del termine

ik ik kj

(k) (k)

a a

kk kk

(k)

a viene smorzato.

kj (k)

a (k)

Al contrario, se > 1 l’errore di a viene amplificato.

ik kj

(k)

a

kk (k)

Si ricerca quindi il valore di a maggiore della k-esima colonna in modo

kk

(k)

a

da avere il rapporto più piccolo possibile. Questo processo è detto

ik

(k)

a

kk

pivoting.

Finito il processo di eliminazione si passa alla sostituzione.

L’ultima equazione avrà come soluzione

(n)

b

n

x =

n (n)

a nn

mentre, la penultima penultima equazione, che è del tipo

(n−1) (n−1) (n−1)

a x + a x = b

n−1 n

n−1,n−1 n−1,n n−1

5

avrà come soluzione (n−1) (n−1)

−

b a x n

n−1 n−1,n

x =

n−1 (n−1)

a

n−1,n−1

Generalizzando si ottiene

(i) (i)

nj=i+1

P

−

b a x j

i ij ∀i − −

= n 1, n 2, . . . , 1 (1.2)

x =

i (i)

a ii

Si applichi il metodo di eliminazione di Gauss in un sistema di 30 equazioni

in 30 incognite. Il suddetto sistema è composto da

n−1

n−1 X

X −

− − (n k)

(n k)(n k + 1) + k=1

k=1

operazioni. Il primo termine si riferisce al numero di moltiplicazioni nel-

la fase di eliminazione, mentre il secondo termine si riferisce alla fase di

sostituzione.

In totale si hanno 3 2

n n 5

−

N OP = + n

3 2 6

Per n = 30 risulta N OP = 9425 moltiplicazioni.

30

Il tempo di calcolo è 9425 −6

' ·

T = 5 10 s

9

·

2 10

che è un valore decisamente minore rispetto al calcolo con il metodo di

Cramer.

2 Sistemi tridiagonali

Si consideri il sistema A x = b tale che



a x + a x = b

11 1 12 2 1









a x + a x + a x = b

 21 1 22 2 23 3 2







a x + a x + a x = b (2.1)

32 2 33 3 34 4 3

.. .. ..



 . . .











a x + a x = b

 n,n−1 n nn n n

6

Il sistema (2.1) ha come matrice dei coefficienti

 

a a 0

11 12

a a a

21 22 23

 

 

a a a

32 33 34

 

 

.. .. ..

 

. . .

 

0 a a

n,n−1 nn

Teorema 1.2. ”Se il sistema (2.1) è tridiagonale a diagonale dominante e

∀i

a < 0 = 1, 2, . . . , n (elementi sulla diagonale principale) e a > 0

ii j,j+1

(elementi sulla diagonale superiore), a > 0 (elementi sulla diagonale

j+1,j

∀j −

inferiore) = 1, 2, . . . , n 1 allora la soluzione del sistema (2.1) esiste ed

è unica.” · · · ∀i

Lemma 1. ”Se x = a x + a x + + a x dove a > 0 = 1, 2, . . . , k

1 1 2 2 i

k k

k

P a = 1 allora

e i

i=1 ≤ ≤

min(x , x , . . . , x ) x max(x , x , . . . , x )

1 2 1 1

k k

Dimostrazione. · · · ≤ · · ·

x = a x + a x + + a x a max x + a max x + + a max x

1 1 2 2 1 i 2 i i

k k k

1≤i≤k 1≤i≤k 1≤i≤k

e quindi k

X a ( max x ) = max(x , x , . . . , x )

i i 1 2 n

1≤i≤k

i=1

Lo stesso procedimento vale anche per il minimo.

NB. Si hanno le uguaglianze tra x e il minimo/massimo se e solo se

· · ·

x = x = = x = x

1 2 k · · · ∀i

Lemma 2. ”Se x = a x + a x + + a x dove a > 0 = 1, 2, . . . , k

1 1 2 2 i

k k

ki=1