Metodi numerici in Matlab

1. Metodi numericiin Matlab Gabriella Puppo

2. Argomenti trattati • Algebra lineare • Interpolazione • Studio della matrice del filo elastico • Risoluzione di alcuni problemi di elasticita’

3. Algebra lineare • Risoluzione di sistemi • Risoluzione di sistemi nel senso dei minimi quadrati • Calcolo del numero di condizionamento

4. Risoluzione di sistemi lineari con Matlab Matlab risolve il sistema lineare Ax=b con il comando: x=A\b. Se A è quadrata, questo comando implica i seguenti passi: - Calcola la fattorizzazione LU con pivoting sulle colonne - Risolvi i due sistemi triangolari Se A è rettangolare (sistemi sovradeterminati), Matlab calcola la soluzione nel senso dei minimi quadrati, cioè: - Calcola la fattorizzazione QR con pivoting sulle colonne - Risolve il sistema triangolare superiore Se il condizionamento di A è elevato, stampa un messaggio di warning, ma calcola ugualmente la soluzione.

5. Esempio >> a=[2 -1 0; -1 2 -1; 0 -1 2]; >> b=ones(3,1); >> x=a\b x = 1.5000 2.0000 1.5000

6. Attenzione! Matlab calcola una soluzione anche quando il sistema non ammette soluzione. Nell’esempio seguente, a è singolare: >> a=[1 2 3; 4 5 6; 7 8 9]; >> b=[1; 1; 0]; >> x=a\b; Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.541976e-018. >> norm(a*x-b) ans = 5.0990 Non è chiaro il modo in cui Matlab calcola la soluzione in questo caso.

7. Fattorizzazione LU Matlab calcola la fattorizzazione LU di una matrice con il comando: [l,u] = lu(a) oppure con il comando: [l,u,p]=lu(a) Nel primo caso, l contiene anche gli scambi di riga effettuati. Se la fattorizzazione è PA = LU, allora l=P-1 L. Nel secondo caso, l=L, u=U, p=P.

8. Esempio Notare che l contiene la matrice triangolare L con le righe scambiate, cioè l=P-1 L >> a=[1 2 3; 4 5 6; 7 8 9]; >> [l,u]=lu(a) l = 0.1429 1.0000 0 0.5714 0.5000 1.0000 1.0000 0 0 u = 7.0000 8.0000 9.0000 0 0.8571 1.7143 0 0 0.0000

9. Per forzare la soluzione di un sistema tramite fattorizzazione LU, devo quindi dare i comandi: >> [l,u]=lu(a); >> y=l\b; >> x=u\y;

10. Fattorizzazione QR Matlab esegue la fattorizzazione A=QR mediante il comando: [q,r] = qr(a) Per forzare la soluzione di un sistema lineare mediante fattorizzazione QR devo quindi dare i comandi: >> [q,r]=qr(a); >> x=r\q'*b;

11. Soluzione di sistemi sovradeterminati Un sistema lineare con una matrice M per N, con M > N, e’ un sistema sovradeterminato, che si puo- risolvere nel senso dei minimi quadrati, usando il comando / Esempio. Risolvere il sistema Ax=b, con: >> a=[1 1 1; 1 2 4; 1 -1 1;1 -2 4] ; >> b=[1; 1; 1; 2];

12. La soluzione e’: >> x=a\b x = 0.8333 -0.2000 0.1667 Questa non e’ la soluzione del sistema in senso classico, infatti: >> a*x-b ans = -0.2000 0.1000 0.2000 -0.1000 Quindi il residuo e’ diverso da zero. La soluzione trovata e’ il vettore che minimizza la norma del residuo.

13. Numero di condizionamento Matlab stima il numero di condizionamento di una matrice mediante il comando cond. Precisamente: cond(a,p) calcola il numero di condizionamento della matrice a in norma p. Se p non c’e’, si assume che p=2. Esercizio: calcolare il numero di condizionamento della matrice del filo elastico, per N=10 e N=100. Inoltre, disegnare un grafico con il numero di condizionamento in funzione di N.

14. Matrice del filo elastico % Questo script disegna un grafico contenente % il numero di condizionamento della matrice % di rigidita' del filo elastico in funzione di % N, con N=1:200. nmax=200; for n=1:nmax a=tridiag(n); c(n)=cond(a); end plot(c) title('Numero di condizionamento')

15. Si ottiene questo andamento:

16. Interpolazione polinomiale • Interpolazione polinomiale “esatta” • Interpolazione polinomiale nel senso dei minimi quadrati • Studio dell-andamaneto di un insieme di dati

17. Interpolazione polinomiale Sono assegnati dei dati (x i, f i) con i=1,…,M. Si vuole trovare un polinomio P N di grado N tale che: P N (x i) = f i per i = 1,…,M . Ci interessano due casi: 1) N = M - 1 (interpolazione polinomiale globale) 2) N < M - 1 (interpolazione nel senso dei minimi quadrati) Matlab gestisce entrambi i casi con la function polyfit

18. function polyfit La function polyfit calcola i coefficienti del polinomio che interpola dei dati, precisamente: >> p=polyfit(x,f,n) Interpola i dati contenuti nei vettori x ed f con un polinomio di grado n. Il vettore p contiene i coefficienti del polinomio interpolatore, con la convenzione: p(x) =p(1)*x^n +p(2)*x^(n-1) + … +p(n+1). Se n = length(x)-1, si ha interpolazione “esatta”. Se n < length(x)-1, si ha interpolazione nel senso dei minimi quadrati

19. function polyval La function polyval valuta un polinomio su un insieme di valori x: >> px = polyval(p,x) Il vettore px contiene i valori del polinomio con coefficienti p, calcolato sui punti x.

20. Esempio Interpolare la funzione F(x) = exp(-x 2) * sin(3x) con due polinomi, usando 5 e 10 punti equispaziati sull’intervallo [0 2]

21. Calcola i coefficienti dei due polinomi, p4 e p9 % Interpola la funzione f(x)=exp(-x^2)*sin(3x) % usando un polinomio di grado 4 e un polinomio di grado 9 % f=inline('exp(-x.^2) .* sin(3*x)'); x1=linspace(0,2,5); % griglia di interpolazione f1=f(x1); p4=polyfit(x1,f1,4); %trova i coefficenti del polinomio x2=linspace(0,2,10); % griglia di interpolazione f2=f(x2); p9=polyfit(x2,f2,9); %trova i coefficenti del polinomio

22. Disegna il grafico dei due polinomi, sulla stessa figura % Disegna il grafico dei due polinomi e della funzione f xg=linspace(0,2,201); fg=f(xg); % calcola il polinomio sulle ascisse del grafico p4_xg=polyval(p4,xg); p9_xg=polyval(p9,xg); plot(xg,fg,xg,p4_xg,xg,p9_xg) legend('f','p4','p9')

23. Andamento del numero di condizionamento della matrice del filo elastico, in funzione di n. Vorrei calcolare per quale valore di m cond(A h (n)) ~ n m Secondo la stima teorica, dovrei trovare m=2. Provo quindi ad interpolare i valori del numero di condizionamento con un polinomio di grado 2, e a sovrapporre le due curve, cioe’ il numero di condizionamento, ed il polinomio, per vedere se il polinomio rappresenta bene l’andamento della curva.

24. Questo e’ l’andamento trovato:

25. E questo e’ lo script: % Interpola i dati sul condizionamento con un polinomio % di grado 2 nmax=100; for n=1:nmax c(n)=cond( tridiag(n) ); end ni=1:nmax; p=polyfit(ni,c,2); % Valuta il polinomio sui dati ni pn=polyval(p,ni); plot(ni,c,'go'); hold on; plot(ni,pn,'r','Linewidth',2)

26. … ma se non sapessi m a priori, come potrei determinarlo? Esercizio. Calcolare diversi polinomi di interpolazione dei dati sul condizionamento e osservare i coefficienti dei polinomi ottenuti.

27. Esercizio. Paragonare il numero di condizionamento della matrice del filo elastico per le seguenti condizioni al contorno: 1) Dirichlet: u(0) = u(1) = 0 2) Misto: u’(0) = 1, u(1) = 0.

28. Alcuni problemi di elasticita’ • Carico uniforme • Come passare una function come argomento ad un’altra function • Carico puntiforme

29. Risolviamo ora il problema del filo elastico nel caso seguente: -  u = f in (0,1) u(0) = u(1) = 0 nel caso in cui f (x) = -1. Consideriamo una griglia uniforme, con N punti interni. Per il vettore di carico, prendiamo: b i ~ h f(x i )

30. function [u,x]=filo_unif(n) % Risolve il problema del filo elastico % per condizioni al bordo omogenee di % Dirichlet, con carico uniforme % e N punti interni a=tridiag(n); h=1/(n+1); % spaziatura di griglia b=-h^2*ones(n,1); % ho portato h a sinistra u=a\b; % Ora aggiungo le condizioni al contorno u(2:n+1)=u; u(1)=0; u(n+2)=0; % Calcola il vettore delle ascisse for i=0:n+1; x(i+1) = i*h; end

31. Scriviamo ora una function che ci permetta di calcolare la posizione di equilibrio del filo, per carichi f da definire dall’esterno. Ho bisogno di passare il nome della function come argomento Questa operazione in Matlab e’ piuttosto delicata. Cerchiamo di ricordare come si procede.

32. Passare una function come argomento Molte funzioni “lavorano” su altre funzioni, predefinite dall’utente. Supponiamo di voler scrivere una function che stampa il grafico di una funzione sull’intervallo [a,b]. La function richiesta deve avere in input il nome della function di cui vogliamo il grafico e l’intervallo [a,b]. La sua intestazione sarà quindi: function plot_funz(funz,a,b) La difficoltà nasce dal fatto che per Matlab (prima della versione 6) il parametro “funz” deve essere una stringa di caratteri, contenente il nome della function.

33. Supponiamo quindi di aver creato un M-file f.m che contiene la function di cui vogliamo il grafico: function f=f(x) % Calcola la funzione exp(x)*cos(4*x) f=exp(x).*cos(4*x); Poiché il nome della function deve essere passato come stringa di caratteri, la chiamata a plot_funz per disegnare il grafico di f(x) sull’intervallo [0,5] sarà: >> plot_funz('f',0,5) Oppure, in Matlab 6, si puo’ usare anche uno handle (@): >> plot_funz(@f,0,5)

34. Function feval All’interno della function plot_funz, è necessario valutare la funzione passata come argomento sui valori x delle ascisse. Per far questo devo usare la function feval: f=feval(funz,x); Questa istruzione trasforma la stringa funz , associando a funz il file funz.m . Poi, esegue il comando funz(x). Possiamo ora dare il listato della function plot_funz

35. Listato di plot_funz function plot_funz(funz,a,b) % PLOT_FUNZ(FUNZ,A,B): disegna il grafico della funzione % FUNZ sull'intervallo [a,b], usando 201 punti. % FUNZ deve essere una stringa di caratteri h=(b-a)/200; x=a:h:b; f=feval(funz,x); %Valuta funz nei valori x plot(x,f)

36. Inline functions Per costruire funzioni semplici, posso usare l’istruzione inline, che genera funzioni di una riga. Consideriamo il file f.m: function f=f(x) % Calcola la funzione exp(x)*cos(4*x) f=exp(x).*cos(4*x); Un modo più semplice di generare la stessa funzione è: >> f=inline('exp(x).*cos(4*x)') f = Inline function: f(x) = exp(x).*cos(4*x)

37. Attenzione! L’oggetto f creato dall’istruzione inline è una stringa di caratteri. Quindi se voglio disegnare il grafico di f usando la function plot_funz, ho due possibilità: 1) -creo un file f.m, che contiene la funzione desiderata. - chiamo plot_funz con il comando: >> plot_funz('f',0,5) 2) -creo f come inline function e poi chiamo plot_funz: >> f=inline('exp(x).*cos(4*x)'); >> plot_funz(f,0,5)

38. Ora possiamo costruire una function che risolva il problema del filo elastico, con un carico assegnato dall’esterno. Per risolvere un problema specifico, dovremo creare una function come file.m o come oggetto inline che contenga l’espressione del carico

39. Ho ottenuto il seguente listato: function [x,u]=filo_dir(f,n) % Risolve il problema del filo elastico su (0,1) % per condizioni al bordo omogenee di % Dirichlet, con carico definito dalla funzione % F e N punti interni a=tridiag(n); h=1/(n+1); % spaziatura di griglia % calcola la posizione dei nodi interni: for i=1:n; x(i) = i*h; end b=h^2*feval(f,x); b=b'; % ho portato h a sinistra u=a\b; % Ora aggiungo le condizioni al contorno u(2:n+1)=u; u(1)=0; u(n+2)=0; x(2:n+1)=x; x(1)=0; x(n+2)=1;

40. Esempio. Risolviamo il problema del filo elastico soggetto ad un carico puntiforme, cioe’: f(x) = - delta( x-0.5 )

41. function che descrive un carico puntiforme: function f=delta(x) % F=DELTA(X): Costruisce un vettore riga, che contiene % un carico a forma di delta, quindi F(N/2)=-N, dove % N e' la lunghezza di X, cioe' F(N/2) = -N su un intervallo % di ampiezza 1/N n=length(x); f=zeros(1,n); if mod(n,2)==0 % controlla se n e' pari f(n/2)=-n/2; f(n/2 +1)= -n/2; else f((n+1)/2)=-n; end

42. Risolviamo il problema per diversi valori di N: [x,u]=filo_dir(@delta,9); plot(x,u,'r','Linewidth',2) hold on [x,u]=filo_dir(@delta,99); plot(x,u,'g','Linewidth',2) [x,u]=filo_dir(@delta,999); plot(x,u,'c','Linewidth',2)

43. Si ottiene il seguente grafico: