Metodi numerici
in 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 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.
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
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.
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
0
u=
7.0000 8.0000 9.0000
0 0.8571 1.7143
0
0 0.0000
Notare che l contiene la
matrice triangolare L
con le righe scambiate,
cioè l=P-1 L
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;
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\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.
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);
end
plot(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 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
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.
Esempio
Interpolare la funzione F(x) = exp(-x 2) * sin(3x) con due polinomi,
usando 5 e 10 punti equispaziati sull’intervallo [0 2]
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
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')
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.
Questo e’ l’andamento trovato:
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)
… 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) = 0
2) 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 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 )
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
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.
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 funz
il 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 caratteri
h=(b-a)/200;
x=a:h:b;
f=feval(funz,x); %Valuta funz nei valori x
plot(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
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;
Esempio.
Risolviamo il problema del filo elastico soggetto ad un carico
puntiforme, 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/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
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: