Calcolo numerico

T

y = Q b

x = Ry

y=Q’b;

x=solupper(R,y);

2 Tabella Complessitá

2.1 Algoritmo di Thomas

Utile per le matrici tridiagonali

Complessitá Computazionale: O(n) u = d

1 1

f or i = 2, ..., n

c

l = i

i u

i−1

−

u = d l b

i i i i−1 3

3 Matrice di Vandermonde

La matrice di Vandermonde é sempre malcondizionata, essa viene rappre-

sentata dal polinomio di Lagrange.

4 Interpolazione

4.1 Polinomio di Lagrange

Dati n+1 punti, il polinomio di Lagrange é di grado n ed é unico.

Lagrange = y L (x) + y L (x) + ... + y L (x) = Σy L

0 0 1 1 n n i i

(x−x )(x−x )...(x−x )

n

0 1 −

L (x) = al numeratore manca (x x )

i i

−x −x −x

(x )(x )...(x )

n

0 1

i i i ω(x) n−1

f (ξ)

Errore di Interpolazione R(x) = (n+1)!

− − −

con ω(x) = (x x )(x x )...(x x )

0 1 n

4

5 Estrapolazione

5.1 Polinomi di Chebyshev

∗

T = cos(n arccos(x))

n

Si tratta di polinomi di grado crescente:

n n−1

• Il coefficiente di x vale 2 in T (x)

n

• −1 ≤ ≤

T (x) 1 (perché si tratta di una funzione coseno)

n

• Sono funzioni pari per n pari e dispari per n dispari

I polinomi di Chevyshev sono linearmente indipendenti.

cos(2i+1)π

b−a a+b

Zeri del polinomio: x = ( ) +

)(

i 2 2

2(n+1)

5.1.1 Nodi di Chebyshev ∗ |ω(x)|

Questi nodi rendono l’interpolazione con l’errore minimo. ω = max =

x∈[a,b]

||ω(x)||

∞

5.2 Polinomio di Newton

Il polinomio di interpolazione di Newton si costruice con la tabella delle

differenze divise, gli elementi della diagonale sono i coefficienti del polinomio

− − −

p (x) = f [x ] + f [x x ](x x ) + f [x x ] + f [x x ](x x )(x x ) + ...

n 0 0 1 0 0 1 0 1 0 1 (n(n+1)

Utilizzando lo schema di Horner generalizzato la complessitá é di O 2

5.3 Interpolazione di Birkoff Hermite

É una generalizzazione dei polinomi di Taylor, Lagrange e Newton, esso

mantiene la struttura di f(x) meglio degli altri algoritmi di interpolazione.

Si impongono condizioni sui valori che deve assumere un polinomio

i punti prefissati e condizioni sui valori delle derivate.

5

6 Metodo dei Grafi

Fattori di Amplificazione y

x ±

Somma: ε =

x±y ε ε

x±y x±y

Prodotto: ε = ε + ε

xy x y

−

Divisione: ε = ε ε

xy x y

= αε

Potenza: ε α x

x

√ 1

n

Radice: x = ε x

n

7 Metodi di convergenza

Condizione necessaria sufficiente per una matrice convergente:

ρ(A) < 1

Condizione necessaria ma non sufficiente: 1

|det(G)| |tr(G)|

< 1 e < n Velocitá di convergenza: R (G)

∞

7.1 Metodo di Iacobi ⇒

Se il raggio spettrale ρ < 1 Metodo di Iacobi converge

−logρ(J)

Velocitá di convergenza= −

Teorema: Sia A di ordine n simmetrica e sia 2D A definita positiva.

Allora A é definita positiva se e solo se il metodo di Iacobi converge.

7.2 Metodo di Gauss - Seidel

−1

−

G = (D L) U

Se A é:

• strettamente diagonale dominante

• irriducibilmente diagonale dominante ||G|| ||J||

allora sia Iacobi che Gauss-Seidel convergono e <

∞ ∞

Teorema: Sia A quadrata di ordine n simmetrica con a > 0 allora a é

ii

definita positiva, se e solo se, il metodo di Gauss-Seidel é convergente.

6

7.3 SOR

Tecnica di rilassamento o estrapolazione:

• ω = 1: fornisce il metodo di partenza

• ω > 1: sovrarilassamento

• ω < 1: sottorilassamento Teorema: Sia A simmetrica con a > 0 e

ii

sia 0 < ω < 2 allora A é definita positiva se e solo se il metodo SOR

converge.

7.4 Convergenza Globale di Newton

f (x )

− k

Iterazione di Newton: x m

k f (x )

k+1

Newton converge globalmente se

• f (a) < 0 e f (b) < 0

0

• 6 ∈

f (x) = 0 [a, b]

00 00

• ≤ ≥

f (x) 0 oppure f (x) 0

0 0

• |f ≤ − |f ≤ −

(b)| (b a)|f (b)| oppure (a)| (b a)|f (a)|

7.5 Altre considerazioni

Se A é tridiagonale con elementi diagonali non nulli, valgono le seguenti

asserzioni:

• −µ

se µ é autovalore della matrice di Iacobi J anche é autovalore di J

• se µ é autovalore di J e vale la relazione:

2 2 2

−

(λ + ω 1) = ω λµ allora λ é autovalore di G(ω)

• se λ é un autovalore non nullo di G(ω) e vale la relazione

2 2 2

−

(λ + ω 1) = ω λµ allora µ é un autovalore di J

8 Info Generali

Problema di malcondizionamento: permutando di poco i dati iniziali si

trovano soluzioni diverse r=b-Aw

Se r tende a 0 allora probabilmente w é soluzione attendibile, infatti:

7

||x−w|| ||r||

−1

||A||||A ||

Err = <=

relativo ||x|| ||b|| ||r||

−1

||A||||A || ||A||

k(A) = oppure k(A) = ||b|| ≥

In una norma naturale (per matrici quadrate), vale che k(A) 1

−1

||I|| ||AA || ≤

1 = = k(A)

se k(A) 1 é malcondizionata

'

se k(A) 1 é bencondizionata

Il numero di condizionamento di una matrice esprime quanto una matrice

é vicina alla singolaritá

Piú in generale:

max||Ax||

k(A) = min||Ax||

Tuttavia il calcolo dell’inversa é molto complesso, si occorre dunque ad

una stima del numero di condizione:

||x−w||

Sottostima: ||r|| 8

9 Codici Utili

% es 1

clear all

close all

for n=0:2:16

s=1+sum(1./[2:2:n])-log(n);

fprintf(’n=%g appr=%g\n’,n, s);

end

% es 2

close all

clear all

n=10;

A=hilb(n);

sol=[1:2:2*n-1]’;

%%%%%%%%%%%%%

b=A*sol;

[R,p] = chol(A);

if p~=0

error(’A non definita positiva’);

else

y=sollower(R’,b);

x=solupper(R,y);

fprintf(’Soluzione calcolata con operatore Matlab\n’)

xx=A\b;

for i=1:n

fprintf(’x=%g xx=%g\n’,x(i),xx(i));

end

fprintf(’\n’);

r=b-A*x;

resnorm=norm(r,inf)/norm(b,inf);

fprintf(’norm(r)/norm(b)=%g\n’,resnorm);

ee=sollower(R’,r);

ee=solupper(R,ee);

KA=norm(ee,inf)/norm(r,inf)*norm(A,inf);

fprintf(’ stima del numero di condizione =%g\n’, KA);

fprintf(’stima dell’’errore relativo=%g\n’, KA*resnorm);

end 9

%es 3

clear all

close all

A=[1 0 1

0 1 1

2 -1 2];

b=[1 0 1]’;

[L,R,P]=gauss2(A);

y=sollower(L,b(P));

x=solupper(R,y);

% es 4

close all

clear all

n=30;

D=10+1:2:10+2*n-1;

D=D’;

D_1=-1+1./[2:n+1];

D_1=D_1’;

D1=-sqrt(1:2:2*n-2);

D1=[0;D1’];

A=spdiags([D_1, D,D1],[-1 0 1],n,n);

xex=ones(n,1);

b=A*xex;

x=sparse(n,1);

maxit=100;

tol=1e-6;

[ x, k ] =g_s(A,b,x,maxit,tol);

err=norm(x-xex,inf);

fprintf(’errore=%g\n’,err);

fprintf(’numero iterazioni=%g\n’,k);

fprintf(’residuo normalizzato=%g\n’,norm(b-A*x,inf)/norm(b,inf));

10