Metodi numericiin Matlab

Gabriella Puppo

Argomenti trattati

• Algebra lineare

• Interpolazione

• Studio della matrice del filo elastico

• Risoluzione di alcuni problemi di elasticita’

Algebra lineare

• Risoluzione di sistemi

• Risoluzione di sistemi nel senso dei minimi quadrati

• Calcolo del numero di condizionamento

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 lasoluzione 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 diwarning, ma calcola ugualmente la soluzione.

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

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 questocaso.

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.

Esempio

>> 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 0u = 7.0000 8.0000 9.0000 0 0.8571 1.7143 0 0 0.0000

Notare che l contiene lamatrice triangolare Lcon le righe scambiate,cioè l=P-1 L

Per forzare la soluzione di un sistema tramite fattorizzazioneLU, devo quindi dare i comandi:

>> [l,u]=lu(a);>> y=l\b;>> x=u\y;

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;

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];

La soluzione e’:

>> x=a\bx = 0.8333 -0.2000 0.1667

Questa non e’ la soluzione del sistema in senso classico, infatti:

>> a*x-bans = -0.2000 0.1000 0.2000 -0.1000

Quindi il residuo e’ diverso da zero.La soluzione trovata e’ il vettore cheminimizza la norma del residuo.

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.

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);endplot(c)title('Numero di condizionamento')

Si ottiene questo andamento:

Interpolazione polinomiale

• Interpolazione polinomiale “esatta”

• Interpolazione polinomiale nel senso dei minimi quadrati

• Studio dell-andamaneto di un insieme di dati

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

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 gradon. 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

function polyvalLa function polyval valuta un polinomio su un insieme di valori x:

>> px = polyval(p,x)

Il vettore px contiene i valori del polinomio con coefficientip, calcolato sui punti x.

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

% 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 interpolazionef1=f(x1);p4=polyfit(x1,f1,4); %trova i coefficenti del polinomio x2=linspace(0,2,10); % griglia di interpolazionef2=f(x2);p9=polyfit(x2,f2,9); %trova i coefficenti del polinomio

Calcola i coefficienti dei due polinomi, p4 e p9

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

Disegna il grafico dei due polinomi, sulla stessa figura

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 condizionamentocon 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.

Questo e’ l’andamento trovato:

% Interpola i dati sul condizionamento con un polinomio% di grado 2nmax=100;for n=1:nmax c(n)=cond( tridiag(n) );endni=1:nmax;p=polyfit(ni,c,2);% Valuta il polinomio sui dati nipn=polyval(p,ni);plot(ni,c,'go'); hold on; plot(ni,pn,'r','Linewidth',2)

E questo e’ lo script:

… 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.

Esercizio.

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

Alcuni problemi di elasticita’

• Carico uniforme

• Come passare una function come argomento ad un’altra function

• Carico puntiforme

Risolviamo ora il problema del filo elastico nel casoseguente:

- 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 )

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 internia=tridiag(n);h=1/(n+1); % spaziatura di grigliab=-h^2*ones(n,1); % ho portato h a sinistrau=a\b;% Ora aggiungo le condizioni al contornou(2:n+1)=u; u(1)=0; u(n+2)=0;% Calcola il vettore delle ascissefor i=0:n+1; x(i+1) = i*h; end

Scriviamo ora una function che ci permetta di calcolare laposizione di equilibrio del filo, per carichi f da definiredall’esterno.

Ho bisogno di passare il nome della function come argomentoQuesta operazione in Matlab e’ piuttosto delicata.Cerchiamo di ricordare come si procede.

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.

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)

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 funzil file funz.m . Poi, esegue il comando funz(x).Possiamo ora dare il listato della function plot_funz

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 caratterih=(b-a)/200;x=a:h:b;f=feval(funz,x); %Valuta funz nei valori xplot(x,f)

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)

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)

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

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 internia=tridiag(n);h=1/(n+1); % spaziatura di griglia% calcola la posizione dei nodi interni:for i=1:n; x(i) = i*h; endb=h^2*feval(f,x); b=b'; % ho portato h a sinistrau=a\b;% Ora aggiungo le condizioni al contornou(2:n+1)=u; u(1)=0; u(n+2)=0;x(2:n+1)=x; x(1)=0; x(n+2)=1;

Ho ottenuto il seguente listato:

Esempio.

Risolviamo il problema del filo elastico soggetto ad un caricopuntiforme, cioe’:

f(x) = - delta( x-0.5 )

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/Nn=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

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)

Si ottiene il seguente grafico: