MATLAB

1. MATLAB

2. Annalisa Pascarella pascarel@dima.unige.it www.aula.dimet.unige.it

3. …di cosa parliamo oggi… • Sistemi lineari Ax=b • A quadrata singolare e non • A m x n • Fattorizzazione LU • Cenni sulla stabilità • M-file • Esercizi

4. matrice dei coefficienti matrice dei termine noti Sistemi Lineari Teorema di Rouchè-Capelli • rango(A) = rango(A|b) <=> il sistema è risolubile • rango(A) = rango(A|b) = n => unica soluzione • rango(A) = rango(A|b) = p < n => ∞n-p soluzioni

5. Risoluzione di sistemi Lineari Ax = b • Sostituzione • Cramer • SE A è quadrata E invertibile • Algoritmo di Gauss • operazioni elementari sulle righe • pivotizzazione (parziale o totale) -> importante quando si implementa l’algoritmo al calcolatore

6. il simbolo nn è quello della divisione!! Sistemi Lineari non singolari • la soluzione è calcolata mediante l’algoritmo Gaussiano con pivot parziale • tempo richiesto minore del calcolo dell’inversa x = A\b x = inv(A)*b

7. Esempio A = [1 1 1; 1 1 -1; 1 -1 1]; b = [3 2 2]’; det(A) x = A\b ci assicuriamo che il det di A sia diverso da zero

8. A quadrata singolare A = [3 4 -1; 5 2 3; 0 1 -1]; b = [14 14 2]’; rank(A) rank([A b]) per vedere se il sistema è risolubile confrontiamo il rango di A con quello della matrice completa (A|b) questo sistema è risolubile in quanto rango(A)=rango(A|b)=2 => ∞ soluzioni

9. Il comando rref • Per studiare e risolvere un sistema qualunque si deve ridurre la matrice completa (A|b) • Si usa il comandorref(reduced row echelon form ) rref([A b]) rrefmovie([A b])

10. A non quadrata A = [1 0 1; 0 -1 0]; b = [0 1]’; rank(A) rank([A b]) rref([A b]) questo sistema è risolubile in quanto rango(A)=rango(A|b)=2 => ∞ soluzioni

11. …riassumendo… • A quadrata NON singolare • Matlab ci fornisce la soluzione • A quadrata singolare o A m x n • rref ci restituisce la matrix ridotta • le soluzioni le dobbiamo scrivere noi a partire dalla matrix ridotta ottenuta x = A\b rank(A) rank([A b]) rref([A b])

12. Esercizi • Studiare e risolvere, eventualmente, i seguenti sistemi lineari: • Lanciare per l’ultima matrice il comando rrefmovie([A b])

13. Fattorizzazione LU dove: • U è la matrice triangolare superiore ottenuta da A mediante l’algoritmo di Gauss con pivotizzazione parziale • L è una matrice quadrata invertibile e “a meno di permutazioni delle righe” è una matrice triangolare inferiore con tutti 1 sulla diagonale

14. LU in MAtlab A=[1 0 1 1 3 2 1 -3 -8] det(A) [L U] = lu(A)

15. Risoluzione di un sistema con LU sistema triangolare In Matlab [L U] = lu(A); y = L\b; x = U\y Il tempo complessivo richiesto dai tre comandi equivale a quello richiesto dall’algoritmo Gaussiano

16. Quando conviene usare LU? Quando si hanno più sistemi con la stessa matrice dei coefficienti A conviene decomporre una sola volta A in LU e risolvere i sistemi lineari con i comandi precedenti => risparmio di tempo [L U] = lu(A); y = L\bi; i=1,…,k x = U\y;

17. Esercizio • Dopo essersi accertati che il sistema ammette un’unica soluzione trovare la soluzione con Gauss e mediante la decomposizione LU

18. Stabilità • Un problema è ben condizionato se le perturbazioni dei dati nn influenzano “eccessivamente” i risultati • Nel caso di un algoritmo per indicare il suo comportamento rispetto alla propagazione degli errori si usa il termine stabilità • Un algoritmo è stabile se un errore “nn si propaga tanto durante il calcolo”

19. …concretamente… A=[10 7 7 8; 7 10 5 9; 7 5 5 6 ; 8 9 6 10]; b=[32 31 23 33]'; det(A) x=A\b

20. …perturbando i dati… b=[32 31 23 33]'; db = [-0.1 0.1 0.1 -0.1]’ x=A\(b+db) Sistema mal condizionato!!!

21. Esercizio • Generare la matrix di Hilbert col comando • Creare un vettore colonna xesatto di lunghezza n=10 formato dagli elementi 1,2,…,10 • Creare la matrix dei coefficienti b = A*xesatto • Risolvere il sistema lineare Ay = b • y coincide con xesatto? A = hilb(10)

22. M-file • E’ possibile scrivere degli script in Matlab • cliccando su new • File -> New -> M-file

23. Editor • Dopo aver scritto i comandi • si salva il file cliccando su save • si assegna un nome al file (es. prova_lu.m)

24. Richiamare lo script • si definiscono A e b • si richiama lo script prova_lu

25. Esercizi • Provare a scrivere uno script per risolvere un sistema triangolare con il metodo di sostituzione all’indietro • Studiare il seguente sistema

26. Esercizio • Sia H la matrice di Hankel 7x7 del vettore v = (7,6,…,1) (si genera col comando hankel(v) ) • Costruire una matrix A t.c. • le prime n-1 righe e n-1 colonne siano tratte da H • l’ultima colonna sia la successione 6,5,…,1 • l’ultima riga sia la successione 3*7-1,3*7-4,…1 • Risolvere se possibile il sistema lineare Ax = b, dove • b=(0 1 2 2 2 1 1)T • Scrivere le soluzioni