Laboratorio Matlab - €¦ · MATLAB pu`o essere utilizzato in modo diretto per calcolare semplici...

Laboratorio Matlab

Matematica I:Calcolo differenziale, Algebra lineare, Probabilita e Statistica

G.Aletti, G.Naldi, L.Pareschi

Introduzione a MATLAB

Introduzione a MATLAB

Questo capitolo introduttivo si propone lo scopo di fornire una breve guida a MATLAB.Per maggiori approfondimenti, oltre ai laboratori MATLAB presentati nel seguito, ri-mandiamo a testi piu specifici 1.Per iniziare a lavorare con MATLAB in ambiente Windows solitamente occorre selezio-nare l’icona di MATLAB nel desktop o nella barra comandi, mentre se si usa il sistemaoperativo Unix (o qualche derivato come Linux) e necessario scrivere matlab, seguito dainvio, in una finestra al prompt dei comandi. Una volta avviato, MATLAB vi offre di-verse finestre; fisseremo la nostra attenzione sulla finestra dei comandi caratterizzata dalprompt >> se state utilizzando la versione Professional Edition o dal prompt EDU>> nelcaso abbiate a disposizione la versione Student Edition. Per comodita nel seguito utiliz-zeremo sempre >> per indicare il prompt della riga comandi della finestra principale diMATLAB.

Prime operazioni

MATLAB puo essere utilizzato in modo diretto per calcolare semplici espressioni mate-matiche:

>> 8 + 3 - 1ans =

10>> 1/sqrt(2)ans =

1.4142

La funzione sqrt e una funzione matematica elementare di MATLAB che calcola la radicequadrata del suo argomento. La risposta e del tipo ans=..., dove ans e una variabilegenerata automaticamente da MATLAB quando un’espressione non e assegnata a unavariabile definita dall’utente. Le variabili definite dall’utente possono essere utilizzate inaltre espressioni matematiche

>> a = 8 + 3 - 1a =

10>> b = 1/sqrt(2)b =

1.4142>> c = a + bc =

11.4142

1Si veda oltre alla Reference Guide di MATLAB, il testo MATLAB: concetti e progetti, ed.Apogeo 2002, degli stessi autori.

Introduzione a MATLAB 3

Per evitare la visualizzazione del risultato basta concludere l’espressione con un punto evirgola ;. Per visualizzare il contenuto di una variabile e sufficiente scriverne il nome/ieventualmente separati con una virgola: ,

>> a,b,ca =

10b =

1.4142c =

11.4142

Piu espressioni o piu istruzioni MATLAB possono convivere sulla stessa riga, basta chesiano separate da un punto e virgola (se vogliamo sopprimere la visualizzazione di qualcherisultato) oppure da una virgola (se vogliamo visualizzare il risultato). Per esempio

> a=2+3/5; b=sqrt(a), c=a*bb =

1.6125c =

4.1924

I nomi delle variabili possono essere lunghi un massimo di 31 caratteri con la distinzionetra lettere maiuscole e minuscole, per esempio le variabili a e A sono distinte. La primalettera di una variabile deve essere un carattere alfabetico (a-z, A-Z), mentre dalla se-conda lettera in avanti possiamo utilizzare un qualsiasi carattere alfanumerico incluso ilsimbolo underscore “ ”.

Variabile Significato

ans valore ultima operazione eseguita non assegnata a una variabilei,j unita immaginaria,

√−1pi approssimazione di π, 3.14159265...eps precisione macchinarealmax massimo numero macchina positivo rappresentabilerealmin minimo numero macchina positivo rappresentabileInf ∞, ossia un numero maggiore di realmaxNaN Not a Number, tipicamente il risultato di un’espressione 0/0computer tipo di computerversion versione di MATLAB

Tabella 2.4 Le principali variabili predefinite in MATLAB.

Alcune variabili sono predefinite in MATLAB per la memorizzazione di valori speciali.Parte di queste variabili sono riportate nella Tabella 2.4. Il significato delle variabilieps, realmax, realmin, Inf e NaN e strettamente legato all’aritmetica e all’insieme deinumeri macchina utilizzato.I principali operatori disponibili in MATLAB sono indicati nella Tabella 2.5. Vi sono gliusuali operatori aritmetici e altri operatori preceduti da un punto, vedremo in seguito ilsignificato delle operazioni corrispondenti.I numeri in virgola mobile (la rappresentazione che un calcolatore fa dei numeri reali)non sono l’unico tipo di dato numerico ammesso. Un altro tipo di dato utile in varieapplicazioni e presente in MATLAB sono i numeri complessi e le operazioni con questi.Un numero complesso z, in forma algebrica z = Re z + iIm z (con Re z, Im z partereale e parte immaginaria, due valori reali), puo essere scritto in modo simile anche inMATLAB:

4 Introduzione a MATLAB

Operatore Significato

+ addizione- sottrazione* moltiplicazione/ divisioneˆ elevamento a potenza.* moltiplicazione termine a termine per vettori./ divisione termine a termine per vettori.ˆ elevamento a potenza termine a termine per vettori

Tabella 2.5 Le principali operazioni definite in MATLAB.

>> a=3+4i;>> b=2.6+3.3i;>> a+bans =

5.6000 + 7.3000i>> a*bans =

-5.4000 +20.3000i

Ricordiamo che l’unita immaginaria, rappresentata dalle variabili i e j, e tale che i2 =j2 = −1. Le forme a=3+4i, a=3+4*i, a=3+4j, a=3+4*j sono accettate e sono equivalenti.Oltre alle operazioni di base, molte delle funzioni comunemente disponibili per unacalcolatrice scientifica sono presenti in MATLAB. Per esempio

>> y=cos(pi/4)+sin(pi/4)y =

1.4142

utilizza le funzioni trigonometriche coseno cos e seno sin, con argomento π/4, sommatali valori e infine assegna il risultato alla variabile y.Si noti che il risultato dell’espressione cos (π/4)+sin (π/4) e visualizzato utilizzando solo

quattro cifre decimali. Questa e l’impostazione del formato standard di MATLAB. Epossibile modificarla tramite il comando format. Per esempio

>> format long>> cos(pi/4)+sin(pi/4)ans =

1.41421356237309>> format short

abilita prima il formato a 14 cifre decimali (format long), calcola il risultato, poi riattivail formato standard a 4 cifre decimali (format short)2.Per numeri molto piccoli o molto grandi e comoda la notazione esponenziale scientifica incui si esprime un valore numerico come un numero (maggiore di 1 e minore di 10, dettomantissa), moltiplicato per una opportuna potenza di 10, detta esponente. Per esempioil numero 31415000000 sara rappresentato come 3.1415e+010 che significa 3. 1415×1010,mentre il numero 0.000001234 come 1.2340e-006, ossia 1. 234× 10−6. Alcune funzionimatematiche elementari predefinite in MATLAB, sono riportate nella Tabella 2.6, peruna lista piu ampia si provi il comando help elfun.

2La modifica della visualizzazione di un risultato tramite format non ha nulla a che vederecon l’effettiva precisione con cui MATLAB effettua il calcolo. Tutte le operazioni in MATLABsono eseguite in doppia precisione, questo vuol dire avere almeno 15 cifre decimali significative.


Funzione Significato

sin senocos cosenoasin arcosenoacos arcocosenotan tangenteatan arcotangenteexp esponenzialelog logaritmo naturale, base elog2 logaritmo in base 2log10 logaritmo in base 10sqrt radice quadrataabs valore assolutosign la funzione segno

Tabella 2.6 Alcune funzioni matematiche predefinite in MATLAB.

Vettori e matrici

Se le operazioni viste nella sezione precedente sono comuni a molte calcolatrici tascabilinon altrettanto si puo dire della capacita di effettuare calcoli con matrici e vettori. In taleambito MATLAB esprime al meglio le sue potenzialita (non a caso l’acronimo MATLABsta per MATrix LABoratory). Dall’algebra lineare sappiamo come i vettori possono essereconsiderati casi particolari di matrici. Per esempio un vettore di dimensione 8 nella formariga o colonna e, rispettivamente, una matrice di tipo 1 × 8 (una riga e otto colonne) o8×1 (otto righe e una colonna). Questo e esattamente il modo in cui MATLAB definiscei vettori. Non solo, per MATLAB anche le quantita scalari sono semplici matrici di tipo1× 1. Per memorizzare il vettore x = (1,3,5,7,9,11,13,15)T nella forma riga in MATLABpossiamo scrivere

>> x = [1 3 5 7 9 11 13 15]x =

1 3 5 7 9 11 13 15

Le parentesi quadre delimitano gli elementi del vettore, mentre gli spazi bianchi le singolecomponenti del vettore riga. Il vettore colonna puo invece essere ottenuto come segue

>> x = [1;3;5;7;9;11;13;15]x =

13579

111315

Le componenti in questo caso sono delimitate da un punto e virgola che serve per indicarela fine della riga. Per accedere al valore di una componente del vettore (riga o colonna)basta scrivere

>> x(4)ans =

7


Se cambiamo il valore di una componente del vettore dobbiamo ricordarci di usare ilpunto e virgola per evitare la visualizzazione dell’intero vettore

>> x(4)=1;

Per passare da vettori riga a vettori colonna si utilizza il simbolo di apostrofo, per esempio

>> x = [1 3 5 7 9 11 13 15]’x =

13579

111315

Dal punto di vista dell’algebra lineare l’operazione che si effettua e la trasposizione.La funzione length in MATLAB consente di determinare la lunghezza di un vettore

>> length(x)ans =

8

Se ora vogliamo inserire la matrice 2× 4

A =

8>: 1 3 5 7

9 11 13 15

9>; ,

possiamo scrivere

>> A = [1 3 5 7 ; 9 11 13 15]A =

1 3 5 79 11 13 15

Anche in questo caso gli spazi separano gli elementi per colonna e il punto e virgola separale righe. Potremo accedere agli elementi della matrice in maniera naturale utilizzando icorrispondenti indici

>> A=(1,3)ans=

5

Se vogliamo quindi cambiare l’elemento A(1,3) basta scrivere

>> A=(1,3)=4;

La funzione MATLAB che ci consente di determinare le dimensioni di una matrice e size

>> size(A)ans =

2 4

Se vogliamo controllare le variabili che attualmente MATLAB ha memorizzato nell’areadi lavoro, Workspace, possiamo usare i comandi who e whos. Il primo comando, who,elenca i nomi delle variabili attualmente presenti mentre whos aggiunge informazioniriguardanti il tipo di variabile e l’occupazione di memoria,

>> whoYour variables are:

A x


>> whosName Size Bytes Class

A 2x4 64 double arrayx 1x8 64 double array

Grand total is 16 elements using 128 bytes

Come si puo notare una variabile reale occupa 8 byte di memoria. Per rimuovere levariabili dalla memoria lavoro di MATLAB si utilizza l’istruzione clear

>> clear x A

Un operatore di fondamentale importanza per costruire vettori equispaziati e per operarecon insiemi di indici e la notazione due punti . Il vettore x utilizzato in precedenza, x =[1 3 5 7 9 11 13 15], potrebbe essere generato come

>> x = 1:2:15;

Altri esempi sono i seguenti

>> y=1:9y =

1 2 3 4 5 6 7 8 9>> z=10:-2:2z =

10 8 6 4 2

Se il valore iniziale e maggiore di quello finale e il passo e positivo viene creato un vettorevuoto

>> x=7:1x =

[]

Si noti che il vettore vuoto (tabella vuota) si denota in MATLAB con []. Il vettorevuoto ha lunghezza zero, dimensione 0×0. La notazione due punti “funziona” anche convalori non interi, per esempio

>> x=0:0.2:1x =

0 0.2000 0.4000 0.6000 0.8000 1.0000

Nel caso in cui il passo non sia intero puo risultare difficile l’uso dell’operatore duepunti. Per esempio se vogliamo creare un vettore con un numero prefissato di puntiequispaziati all’interno di un dato intervallo. In questo caso e preferibile utilizzare ilcomando linspace. Per esempio

>> a = 0; b=1; n=8;>> x = linspace(a,b,n)x =

Columns 1 through 70 0.1429 0.2857 0.4286 0.5714 0.7143 0.8571

Column 81.0000

restituisce un vettore riga x di lunghezza n con la proprieta che l’elemento di indice ivale x(i) = a + (i− 1) · (b− a)/(n− 1). Un uso particolarmente efficace della notazionedue punti si ha nella gestione di indici di vettori e matrici. Per esempio utilizzare i duepunti “:” singolarmente significa selezionare ogni elemento di una data riga o colonna


>> A = [1 2 3; 4 5 6; 7 8 9];>> A(:,1)ans =

147

>> A(3,:)ans =

7 8 9

Naturalmente e possibile specificare un intervallo di indici ed estrarre cosı parti di vettori(o matrici)

>> x=0:0.2:1.0;>> x(2:5)ans =

0.2000 0.4000 0.6000 0.8000>> A = [1 2 3 ; 4 5 6 ; 7 8 9];>> A(:,1:2)ans =

1 24 57 8

o parti di righe e colonne in matrici e sottomatrici

>> A = [1 2 3; 4 5 6; 7 8 9];>> A(3,1:2)ans =

7 8>> A(2:3,1:2)ans =

4 57 8

La notazione due punti puo essere usata anche per assegnare in modo rapido nuovi valoria righe e colonne di matrici

>> A = [1 2 3; 4 5 6; 7 8 9];>> A(2,:)=2:2:6A =

1 2 32 4 67 8 9

Esistono molte funzioni MATLAB che consentono di costruire particolari matrici e vet-tori, alcune di queste sono riportate per convenienza nella Tabella 2.7. Queste funzioniMATLAB possono essere utilizzate con un diverso numero di parametri. Per esempioabbiamo

>> B=zeros(2,3)B =

0 0 00 0 0

>> Z = zeros(2)Z =

0 00 0

Analogamente per altre funzioni quali ones, rand, eye. Attenzione, nel caso la matriceabbia elementi complessi la trasposizione non si limita a scambiare righe con colonne, maapplica l’operazione di coniugazione a ogni componente


Funzione Azione

linspace vettore riga di elementi equispaziatilogspace vettore riga di elementi equispaziati in scala logaritmicazeros matrice contenente solo elementi uguali a zeroones matrice contenente solo elementi uguali a unorand matrice contenente numeri pseudo-casuali in [0,1]eye matrice identitadiag matrice diagonale

Tabella 2.7 Alcune funzioni predefinite per la costruzione di particolari matrici inMATLAB.

>> A=[1+3i 2+5i ; 6+7i 1-4i]A =

1.0000 + 3.0000i 2.0000 + 5.0000i6.0000 + 7.0000i 1.0000 - 4.0000i

>> A’ans =

1.0000 - 3.0000i 6.0000 - 7.0000i2.0000 - 5.0000i 1.0000 + 4.0000i

La funzione eye genera invece la matrice identita

>> I3 = eye(3)I3 =

1 0 00 1 00 0 1

Operazioni tra matrici

La maggior parte delle funzioni predefinite in MATLAB accettano come argomenti arraya piu indici. Per esempio

>> n = 7;>> x = linspace(0,pi,n);>> c = cos(x);>> [x’ c’]ans =

0 1.00000.5236 0.86601.0472 0.50001.5708 0.00002.0944 -0.50002.6180 -0.86603.1416 -1.0000

L’istruzione c=cos(x) applicata a un vettore x restituisce un vettore c di uguali dimen-sioni e tipo con la proprieta che l’elemento di indice i vale c(i) = cos(x(i)). L’istruzione[x’ c’ s’] crea una matrice le cui colonne sono i vettori trasposti x’, c’ e s’. Uno deivantaggi principali dell’uso di funzioni vettoriali si ha in termini di velocita di esecuzionedelle operazioni rispetto alla corrispondente esecuzione in forma scalare.

Consideriamo alcune operazioni che coinvolgono array a piu indici e che agiscono si-multaneamente su tutte le componenti dell’array stesso (si veda la Tabella 2.5). Vediamosubito degli esempi,


>> x = [1 3 5 7 9];>> y = [50 40 20 30 10];>> 4*xans =

4 12 20 28 36>> x+yans =

51 43 25 37 19>> y-xans =

49 37 15 23 1

Le tre operazioni appartengono all’algebra lineare e sono la moltiplicazione di un vettoreper uno scalare, la somma e la differenza tra vettori. In maniera del tutto analoga leprecedenti operazioni agiscono su matrici

>> A = [ 1 2 3 ; 4 5 6 ; 7 8 9];>> B = [ 1 1 1 ; 2 2 2 ; 3 3 3];>> C = A+BC =

2 3 46 7 8

10 11 12>> 2*Bans =

2 2 24 4 46 6 6

Al contrario le operazioni puntuali sono tipiche di MATLAB e non hanno necessariamenteun corrispondente dal punto di vista dell’algebra lineare. Per esempio

>> x.*yans =

50 120 100 210 90>> y./xans =

50.0000 13.3333 4.0000 4.2857 1.1111>> x.^2ans =

1 9 25 49 81

L’istruzione x.*y effettua la moltiplicazione puntuale tra vettori e fornisce un vettore conla proprieta che ogni sua componente e uguale al prodotto delle corrispondenti componen-ti dei singoli vettori x e y. Allo stesso modo agisce la divisione puntuale y./x restituendoun vettore le cui componenti sono il risultato della divisione delle corrispondenti com-ponenti di y per quelle di x. Infine l’espressione x. 2 eleva ogni componente del vettorex al quadrato. Per vettori x, y con le stesse dimensioni e anche possibile l’operazionepuntuale y. x, ogni componente del vettore y e elevato alla potenza determinata dallacorrispondente componente del vettore x

>> y.^xans =

1.0e+010 *0.0000 0.0000 0.0003 2.1870 0.1000

Le stesse operazioni possono essere applicate nel caso di vettori colonna o piu in generale

nel caso di matrici. E essenziale che gli operandi siano dello stesso tipo e abbiano lestesse dimensioni. Le uniche eccezioni a questa regola si verificano quando le precedentioperazioni vengono applicate coinvolgendo una matrice e una costante. In tal caso MAT-LAB considerera la costante come un vettore di pari dimensioni avente tutte componenticostanti


>> x = 1:7;>> x+2ans =

3 4 5 6 7 8 9>> 2-xans =

1 0 -1 -2 -3 -4 -5>> 2./xans =

2.0000 1.0000 0.6667 0.5000 0.4000 0.3333 0.2857

Consideriamo ora il problema del calcolo di una tabella di valori di una funzione nonpredefinita in MATLAB, per esempio la funzione razionale

f(x) =120− 60x + 12x2 − x3

120 + 60x + 12x2 + x3

sull’intervallo [0,1]. Per costruire una tabella di 5 valori possiamo utilizzare la sequenzadi istruzioni

>> x = linspace(0,1,5);>> y = (120-60*x+12*x.^2-x.^3)./(120+60*x+12*x.^2+x.^3);>> [x’ y’]ans =

0 1.00000.2500 0.77880.5000 0.60650.7500 0.47241.0000 0.3679

Utilizzando le operazioni vettoriali, i valori della funzione in corrispondenza della suc-cessione di punti definita dal vettore x sono creati con una singola istruzione vettoriale eassegnati al vettore y.

Funzioni avanzate

MATLAB e molto piu di un semplice software in grado di eseguire calcoli numerici(anche se in forma molto sofisticata come abbiamo visto nei paragrafi precedenti). ConMATLAB e possibile realizzare grafici di funzioni anche in piu dimensioni e realizzareveri e propri programmi. Non solo, tramite il Symbolic Toolbox con MATLAB e possibileanche eseguire calcoli simbolici.

Cenni di grafica

Uno dei punti di forza di MATLAB e la sua capacita grafica che consente di rappresentaredati memorizzati in vettori e matrici in molti modi differenti. Grafici bidimensionalivengono utilizzati per visualizzare la variazione di una variabile rispetto a un’altra, peresempio il grafico di una funzione y = f(x). Se vogliamo visualizzare in forma grafica lafunzione sin(x) sull’intervallo [0,2π], possiamo per esempio scrivere

>> ezplot(’sin(x)’,[0,2*pi]);

Ma come ha fatto MATLAB a realizzare il grafico della funzione? Se procediamo in mododiverso possiamo capire quale automatismo sia insito nella funzione ezplot. Costruiamoprima una successione di valori della funzione del tipo


>> n = 31;>> x = linspace(0,2*pi,n);>> y = sin(x);

dove in questo caso abbiamo scelto 31 punti equispaziati. Per realizzare il grafico bastaaggiungere l’istruzione

>> plot(x,y)

In questo modo otteniamo una rappresentazione grafica della tabella di valori ottenu-ta semplicemente raccordando con segmenti di retta nel piano cartesiano (x,y) i vertici(x(i),y(i)) in modo ordinato al variare di i da 1 a 31. La scelta della scala di visualiz-zazione e automatica. Se vogliamo evidenziare i vertici della poligonale cosı costruitapossiamo utilizzare l’istruzione

>> plot(x,y,’-o’)

che disegna una linea con un circoletto in corrispondenza dei vertici. Non cambia lafunzione ma i parametri in ingresso, in questo caso i parametri aggiunti specificano lostile della linea del grafico.

Colore Simbolo Linea

y giallo . punto - linea continuam magenta o circoletto : linea punteggiatac ciano x per -. linea puntor rosso + piu - - linea tratteggiatag verde * asteriscob blu s quadratinow bianco d diamantek nero v triangolo

Tabella 2.8 Principali opzioni del comando plot in MATLAB.

Esistono molte possibilita per personalizzare un grafico come quello creato nell’esem-pio precedente, nella Tabella 2.8 riportiamo le opzioni principali. Possiamo inoltrecommentare il grafico aggiungendo delle scritte sugli assi coordinati e mettendo un titolo

>> title(’Grafico della funzione sin(x)’)>> xlabel(’asse x’)>> ylabel(’asse y’)

Se vogliamo produrre un risultato visivamente privo di “spigoli”, sara sufficiente aumen-tare il numero di punti n in modo tale che i segmenti di raccordo siano cosı piccoli dafornire l’impressione di una curva continua (vedi la Figura 2.5).

E possibile cambiare scala nei grafici tramite la funzione axis che permette di stabilirela grandezza della finestra in cui vedere il grafico. Per esempio

>> x = linspace(0,2*pi);>> y = sin(x);>> plot(x,y);>> a = axisa =

0 7 -1 1>> axis([a(1) pi a(3) a(4)]);


0 1 2 3 4 5 6−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

asse x

asse

y

Grafico della funzione sin(x)

0 1 2 3 4 5 6−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1Grafico della funzione sin(x)

asse

y

asse x

Figura 2.5 Rappresentazione grafica dei valori della funzione sin(x) in [0,2π] con 31punti e 101 punti.

per ottenere solo la parte di grafico sull’intervallo [0,π]. La funzione axis([xmin xmaxymin ymax]) impone che l’intervallo di valori di x sia compreso tra xmin e xmax e quellodi y tra ymin e ymax. Il comando axis usato da solo fornisce un vettore contenete i valori

attuali della finestra grafica. E possibile inoltre sovrapporre grafici di piu funzioni nellastessa finestra grafica. Per esempio se vogliamo realizzare i grafici delle funzioni seno ecoseno sull’intervallo [0,2π] nella stessa finestra grafica dobbiamo costruire i vettori divalori per le due funzioni:

>> x = linspace(0,2*pi);>> y1 = cos(x);>> y2 = sin(x);

Per disegnarne i grafici sovrapposti possiamo scrivere (utilizzando due tratti di lineadiversi)

>> plot(x,y1,’-’)>> hold on>> plot(x,y2,’--’)>> hold off

Il comando hold on mantiene l’attuale contenuto della finestra grafica in modo tale chei grafici successivi si sovrappongano nell’attuale finestra grafica. Il comando hold offritorna all’impostazione originale.Lo stesso risultato puo essere ottenuto in forma piu compatta tramite la singola funzioneplot per grafici multipli

>> plot(x,y1,x,y2)

che utilizza in modo automatico linee di colore differente per i diversi grafici. Infine sevogliamo realizzare i due grafici evidenziando per esempio i punti di intersezione delledue curve possiamo scrivere

>> plot(x,y1,’-’,x,y2,’--’,[1 5]*pi/4,[1 -1]/sqrt(2),’o’);

Per aggiungere una legenda all’interno del grafico possiamo infine utilizzare la funzionelegend

>> legend(’Coseno’, ’Seno’, ’Intersezione’)


0 1 2 3 4 5 6 7−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1Coseno Seno Intersezione

Figura 2.6 Sovrapposizione dei grafici delle funzioni seno e coseno con una legenda.

Nella Figura 2.6 si mostra il risultato. Per aggiungere la griglia grid on, per toglierlagrid off. Come per la funzione hold, scrivere solo grid permette di passare da on aoff e viceversa.Nella Tabella 2.9 riportiamo le principali funzioni che consentono di aggiungere commentia un grafico. Ulteriori funzioni possono essere sperimentate utilizzando il comando helpgraph2d.

Funzione Significato

axis prescrive i valori minimi e massimi sugli assi x e ytitle inserisce un titolo nel graficoxlabel inserisce un nome per l’asse xylabel inserisce un nome per l’asse ygrid inserisce una griglia sugli assi x e ylegend inserisce una legenda per identificare i diversi graficitext inserisce una stringa di testo in una posizione specificata

Tabella 2.9 Funzioni per perfezionare i grafici in MATLAB.

M-files

Le successioni di comandi viste negli esempi precedenti possono essere memorizzate di-

rettamente in un file di testo. In questo modo si crea un M-file MATLAB. E possibilecreare un M-file tramite un editor di testi qualunque. L’ambiente di sviluppo MATLABfornisce un editor orientato alla scrittura di M-file. Per attivare tale editor e sufficienteselezionare il menu File nella finestra principale di MATLAB e poi New e quindi M-file. In alternativa, una rapida scorciatoia consiste nell’utilizzare il comando edit dallafinestra principale di MATLAB.Supponiamo di aver scritto uno script M-file di nome esempio.m all’interno del qualeabbiamo posto in successione i comandi MATLAB per generare il grafico della funzionesin(x), ossia


n = 31;x = linspace(0,2*pi,n);y = sin(x);plot(x,y,’-o’)

Per salvare il nostro script e assegnargli il nome prestabilito, e sufficiente selezionare ilmenu Save as... all’interno del menu File dalla finestra dell’editor di MATLAB.Digitando il nome assegnato allo script nella finestra principale di MATLAB, scrivendoquindi esempio, si richiama e si esegue questo M-file, ottenendo lo stesso risultato che sisarebbe avuto scrivendo uno a uno i comandi elencati nello script. Per una minima ma-nipolazione dei file su disco, MATLAB mette a disposizione alcuni comandi, i principalisono riassunti nella Tabella 2.10.

MATLAB MS-DOS UNIX

dir dir lstype type catdelete del rmcd cd cdpwd chdir pwdmkdir mkdir mkdircopyfile copy cp

Tabella 2.10 Alcuni comandi per la gestione del file system.

Il calcolo simbolico

Molti linguaggi di programmazione consentono di effettuare calcoli numerici analoghi aquelli visti finora. Tali linguaggi dispongono di librerie numeriche opportune che consen-tono di eseguire in modo semplificato operazioni molto complesse. Al contrario calcolisimbolici (ossia calcoli analoghi a quelli che noi eseguiamo con carta e penna quandosviluppiamo il quadrato di un binomio (a + b)2 = a2 + b2 + 2ab) possono essere effettuatida pochi programmi specifici quali per esempio Maple, Mathematica, MATLAB, Scilab eOctave. In particolare MATLAB puo eseguire calcoli simbolici attraverso il Symbolic Ma-thematics toolbox che utilizza le routine di calcolo simbolico di Maple. Il calcolo simboliconon comporta le approssimazioni introdotte dal calcolo numerico. Il prezzo da pagare perottenere questo ricade sulla velocita d’esecuzione e sulla necessita di particolari strutturedati. In generale i risultati numerici possono essere ottenuti molto piu velocemente tra-mite un calcolo numerico invece di una valutazione numerica di un calcolo simbolico. Laparola chiave principale e syms che consente di dichiarare quali variabili sono simboliche.L’uso risulta poi estremamente intuitivo

>> syms a b>> expand((a+b)^2)ans =a^2+2*a*b+b^2

Abbiamo usato il comando expand che forza MATLAB ad “espandere” l’espressionesimbolica. Ritroviamo la nota formula per il calcolo di radici di equazioni di secondogrado tramite

>> syms a b x>> solve(’a*x^2+b*x+c=0’)


ans =[ 1/2/a*(-b+(b^2-4*a*c)^(1/2))][ 1/2/a*(-b-(b^2-4*a*c)^(1/2))]

Nel corso del presente testo si e fatto largo uso del calcolo simbolico in quanto consentesostanzialmente di risolvere la maggior parte dei problemi analitici che vengono affrontati.Per esempio il calcolo di limiti, di integrali e di derivate. Si considerino i seguenti esempi(rinviamo ai rispettivi capitoli per una trattazione piu approfondita)

>> syms x>> limit(1/x,x,0,’right’)ans =inf>> limit(1/x,x,0,’left’)ans =-inf>> limit(1/x,x,inf)ans =0

>> syms a b>> int(x^2,a,b)ans =1/3*b^3-1/3*a^3>> diff(exp(-x^2))ans =-2*x*exp(-x^2)

0Preliminari

0.1 Numeri macchina

Una approssimazione inevitabile che si incontra in MATLAB o in ambienti di calcolosimili e dovuta alla rappresentazione dei numeri reali in un calcolatore. In particolareabbiamo a disposizione solo un numero finito di cifre. I numeri macchina sono rappresen-tati usualmente considerando l’idea della notazione scientifica dove un numero e espressocome il prodotto di una quantita per una opportuna potenza di una assegnata base. Ilnumero macchina x e quindi nella forma

x = s× � d1β−1 + d2β

−2 + . . . + dtβ−t�× βp,

dove β e la base (per esempio β = 2), l’insieme delle t cifre d1d2. . . dt si chiama mantissae l’esponente p e un intero. Ogni cifra soddisfa 0 ≤ dk ≤ β − 1, k = 1,. . . ,t, e l’espo-nente p ∈ [L,U ]. Il numero di cifre t e fissato per ogni insieme di numeri macchina chestiamo considerando. Per vedere il piu grande numero macchina e il piu piccolo numeromacchina,

>> realmax , realminans =

1.7977e+308ans =

2.2251e-308

In generale dato un numero reale z dobbiamo trovare un suo rappresentante tra i numerimacchina, lo indicheremo con fl(z). A meno che z sia gia un numero macchina si incorrein un errore. Per z 6= 0 e arrotondando al numero macchina piu vicino l’errore relativo(ovvero l’errore z−fl(z) valutato rispetto a una quantita di riferimento, noi considereremolo stesso z) che si commette e��z − fl(z)

z

�� ≤ eps =1

2β1−t,

dove la quantita eps e detta precisione macchina. La precisione macchina e legataall’approssimazione ottenibile con l’insieme dei numeri macchina a disposizione. InMATLAB

>> epsans =

2.2204e-016

e si puo verificare un’interessante proprieta

>> 1+eps > 1ans =

18 Preliminari

1>> 1+eps/2 > 1ans =

0

ovvero eps e il piu piccolo numero “sentito” dall’aritmetica dei numeri macchina (perquantita dell’ordine di 1). Ricordiamo che il risultato di un confronto, o di una espressionelogica, e diverso da zero (uguale a 1 solitamente) se e vero, e nullo se e falso, per esempio

>> 3>=5ans =0>>3<7ans =

1

Riassumendo, bisogna prestare molta attenzione quando si effettuano operazioni di cal-colo tramite un calcolatore o una semplice calcolatrice tascabile. Infatti sara inevitabilela presenza di errori nel risultato dovuta all’approssimazione numerica.

0.2 Calcoli simbolici

Va precisato che MATLAB e in grado di eseguire anche calcoli esatti, trattando i numericome simboli. Per fare cio abbiamo bisogno del Symbolic Toolbox. Abbiamo che il numero

>> t=0.1t =

0.1000

puo essere rappresentato nella sua forma razionale tramite

>> sym(t)ans =

1/10

Utilizzando l’opzione ’e’ otteniamo la forma razionale del numero piu la differenza trala rappresentazione razionale teorica e l’attuale rappresentazione in termini di numerimacchina

>> sym(t,’e’)ans =

1/10+eps/40

Utilizzando l’opzione ’d’ visualizziamo la sua rappresentazione decimale con 32 cifre

>> sym(t,’d’)ans =

.10000000000000000555111512312578

Interessante infine e anche l’uso dell’istruzione format rat che visualizza l’approssima-zione tramite razionali di un qualunque numero inserito. Per esempio

>> format rat>> sqrt(2)ans =

1393/985

1Funzioni

1.1 Calcolo di polinomi

Tra le funzioni elementari le funzioni polinomiali sono tra quelle piu utilizzate in varialgoritmi numerici. In MATLAB un polinomio e rappresentato come un vettore checontiene i coefficienti ordinati dal coefficiente del monomio di grado massimo al monomiodi grado minore. Per esempio la funzione p(x) = x3 − 2x− 5 si puo rappresentare comep=[1 -2 -5] (se il polinomio ha alcuni coefficienti nulli essi vanno sempre inseriti nelvettore per non alterare il grado del polinomio). Per calcolare una funzione polinomialein corrispondenza di un valore x occorre utilizzare la funzione polyval. Per esempio con

>> p=[1 -3 -5];>> x=3;>> y=polyval(p,x)y =

-5

Si ricorda che l’uso del punto e virgola esegue un comando senza visualizzare il risultato.Se vogliamo realizzare il grafico del polinomio, dobbiamo innanzitutto creare una

griglia di valori. Un calcolatore infatti per realizzare il grafico di una qualunque funzionef(x) utilizza in generale una griglia equispaziata di valori {xk} nei quali valuta la funzioneyk = f(xk), individuando cosı un insieme di punti del piano caratterizzati dalle coppie(xk,yk). Tali punti sono poi raccordati graficamente in modo da creare l’effetto di unalinea continua.Una serie di 200 valori nell’intervallo [−2,5] e le corrispondenti valutazioni del polinomiovengono costruiti in MATLAB tramite

>> x=linspace(-2,5,200);>> y=polyval(p,x);

Per disegnare, per punti, il grafico in tale intervallo basta utilizzare plot(x,y) (si vedala Figura 1.1).Per il calcolo delle radici di un polinomio, continuando l’ultimo esempio1

>> roots(p)ans =

4.1926-1.1926

1Occorre fare attenzione. Non e detto che tutte le radici di un polinomio appartenganoall’insieme dei numeri reali. In generale la radice di un polinomio appartiene all’insieme deinumeri complessi (insieme C estensione dell’insieme R che permette di avere un ambiente adattoper le radici di tutti i polinomi).

20 Funzioni

−2 −1 0 1 2 3 4 5−8

−6

−4

−2

0

2

4

6

Figura 1.1 Grafico di p(x) = x2 − 3x− 5 con plot.

Per moltiplicare due polinomi si puo utilizzare conv mentre per dividere deconv,

>> p=[1 -2 5]; q=[1 2 3];>> c=conv(p,q)c =

1 0 4 4 15>> [a,r]=deconv(c,p)a =

1 2 3r =

0 0 0 0 0

dove nell’ultima espressione r rappresenta il resto della divisione di c con p.

1.2 Altre funzioni elementari

Molte altre funzioni elementari sono elencate in Appendice B. Ulteriori funzioni si posso-no ottenere componendo le funzioni elementari, per esempio la funzione lx=log(abs(x))calcola il logaritmo del valore assoluto di x.Abbiamo detto che per disegnare per punti il grafico di una data funzione occorre prepa-rare l’insieme dei punti x in cui calcolare la funzione in esame e poi considerare istruzioniper la grafica. In alternativa e possibile utilizzare fplot, o ezplot per esempio

>> ezplot(’sin(2*x)+cos(x)’,[-2*pi 2*pi -2 2]);

disegna il grafico della funzione f(x) = sin 2x + cos x, per x ∈ [−2π, 2π] e considerandouna finestra grafica con ordinate y comprese tra −2 e 2 (si veda la Figura 1.2). Puorisultare utile definire un oggetto che contenga una funzione in modo tale da poterlautilizzare in piu contesti, nell’esempio precedente si poteva procedere utilizzando inlinenella forma

>> f=inline(’sin(2*x)+cos(x)’);>> ezplot(f,[-2*pi 2*pi -1 1]);

ed ottenere lo stesso risultato.Nel caso si esplorino regioni non appartenenti al dominio delle funzioni in uso vi sara unsegnale di avvertimento, warning

Funzioni 21

−6 −4 −2 0 2 4 6−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

x

sin(2 x)+cos(x)

Figura 1.2 Grafico di f(x) = sin 2x + cos x con ezplot.

>> log(0)Warning: Log of zero. ans =

-Inf>> 1/sin(0)Warning: Divide by zero. ans =

Inf

Il valore predefinito Inf equivale a +∞ e potrebbe essere utilizzato formalmente neicalcoli, per esempio 1/Inf fornisce il valore zero. Nel caso di indeterminazione il valoreassociato e NaN (Not-a-Number), per esempio l’espressione sin(0)/0 fornisce questo valorecome esito. Esistono anche altre costanti predefinite, per esempio le cifre disponibili diπ sone memorizzate in pi, lo stesso valore puo essere ottenuto in molti modi

>> format long>> 4*atan(1)ans =

3.14159265358979

Possiamo verificare qualche proprieta del logaritmo e delle funzioni trigometrichecon qualche esempio

>>l1= log(2.3 * 5.4) , l2=log(2.3)+log(5.4)l1 =

2.5193l2 =

2.5193>> 10^(log10(100))ans =

100>> acos(1)ans =

0>> cos(2*pi/3)^2+sin(2*pi/3)^2ans =

1.0000

Per le funzioni trigonometriche gli argomenti devono essere espresse in radianti. Nel casosi abbiano i gradi occorre convertirli in radianti, per esempio per il seno di 30◦

>> sin(30*pi/180)ans =

22 Funzioni

0.5000

1.3 Funzioni e calcolo simbolico

Nell’ambiente di calcolo simbolico possiamo avvalerci di un “calcolatore funzionale” ov-vero possiamo lavorare in modo analogo a quanto potremmo fare con carta e penna.Per esempio ritroviamo

>> syms s t>> expand(cos(s+t))ans = cos(s)*cos(t)-sin(s)*sin(t)

oppure

>> syms a x y>> expand(a^(x+y))ans = a^x*a^y

Se consideriamo invece le operazioni tra funzioni osserviamo che tramite il comandosym(’f(x)’) possiamo creare funzioni astratte. Per esempio

>> f=sym(’x^2-1’) ; g=sym(’exp(x)’);>> f+gans = x^2-1+exp(x)>> f*gans = (x^2-1)*exp(x)>> f/gans = (x^2-1)/exp(x)

Per valutare una funzione astratta in un punto possiamo utilizzare subs e scrivere con ledefinizione precedenti

>> syms x>> subs(f,x,3)ans =

8>> subs(g,x,2)ans =

7.3891

Per la composizione di funzioni, utilizzando le due funzioni f e g definite sopra si usa ilcomando compose

>> compose(g,f)ans = exp(x^2-1)>> compose(f,g)ans = exp(x)^2-1

mentre per determinare la funzione inversa il comando finverse

>> f=sym(’(exp(x)-exp(-x))/(exp(x)+exp(-x))’);>> finverse(f)ans = log(1/(-1+x)*(-(-1+x)*(1+x))^(1/2))

2Limiti e continuita

2.1 Calcolo di limiti

Affrontare il calcolo di un limite dal punto di vista numerico puo portare a brutte sor-prese dovute all’approssimazione del mondo dei numeri reali con i numeri macchina. Peresempio non possiamo affrontare il calcolo del limite

limn→+∞

�1 +

1

n

�n

,

semplicemente aumentando il valore di n sperando di avvicinarci al valore limite, ilnumero e. Infatti, se il valore n e troppo grande segue che 1/n e minore della precisionemacchina e non viene piu “sentito” nelle operazioni, quindi (1+1/n) = 1. Da quel valoredi n in avanti la somma (1+1/n) resta uguale a uno e il limite numerico sembra attestarsiproprio su 1. Per esempio

>> n=10^16;>> (1+1/n)^nans =

1

Fenomeni simili possono avvenire anche per altre “anomalie” dovute alla necessita diricorrere ai numeri macchina e alle operazioni tra questi, per esempio si possono perderecifre significative o avere errori detti di cancellazione.Dal punto di vista del calcolo simbolico abbiamo a disposizione la funzione limit, qualcheesempio e piu chiaro di qualsiasi descrizione formale

>> syms x t ;>> limit(sin(x)/x,x,0)ans =1

>> limit((x+1)/(x-1),x,1,’right’)ans =inf

>> limit((x+1)/(x-1),x,1,’left’)ans =-inf

>> limit((3*x^2+1)/(-5*x^2+x+2),x,Inf)ans =-3/5

>> limit((1+1/t)^t,t,Inf)ans =

exp(1)>> limit((t-2)/(t^2-4),t,2)

24 Limiti e continuita

ans =1/4

2.2 Sistemi dinamici discreti

Un interessante applicazione del concetto di limite si ha nel caso di simulazione di semplicimodelli rappresentabili come successioni numeriche. In particolare possiamo pensare amodelli a tempo discreto ovverosia a successioni {xk}k∈N dove l’indice k indica un certoistante, si immagini una specie di cronometro che segni il tempo (secondi, mesi, anni)k = 0, 1, 2,. . ., mentre xk indica una certa quantita al tempo k. Tipico e il caso in cui ilvalore xk+1 e legato a uno o piu valori passati. Ci limitiamo al solo caso

xk+1 = g(xk),

ovverosia a successioni definite per ricorrenza (utilizzando una differente terminologia,abbiamo un sistema dinamico discreto).Un esempio di carattere ecologico riguarda la dinamica di una popolazione di salmonied e dovuto al biologo W. Ricker. L’evoluzione della popolazione (studiata da Ricker sulfiume Fraser in Canada) puo essere scritta nel modo seguente

xk+1 = xk er(1−xk),

dove xk indica il numero di salmoni (normalizzato rispetto a una popolazione di riferimen-to), r e un parametro del modello legato alla modalita di riproduzione e di spostamentodella popolazione di salmoni, l’indice k conta il tempo (di solito in anni). La domandasul modello essenzialmente e: come cambia la dinamica della popolazione al variare delparametro r e del numero inziale di salmoni x0? Abbiamo estinzione, xk → 0, oppurecrescita illimitata, xk → +∞, o un valore stazionario, xk → L > 0, o altro ancora?Per esempio, la popolazione puo diventare periodica? La risposta a queste domandepotrebbe essere una prima indicazione sulla gestione della popolazione, per esempio suquale quantita di salmoni sia permesso pescare in un certo tempo.Non analizzeremo dal punto di vista teorico il modello ma alcuni esperimenti numericisi possono facilmente fare. Di seguito scriviamo le istruzioni MATLAB per generare unvettore x con N valori, per un fissato valore del parametro r e un dato iniziale x0.

>> N=50;>> x=zeros(N,1);>> x0=0.5; x(1)=x0; xk=x0; r=1; %parametri>> for k = 2:N, xk = xk * exp(r * (1-xk)); x(k) = xk; end

La struttura for...end rappresenta un ciclo, ossia permette di ripetere un certo numerodi volte un blocco di istruzioni nelle quali la variabile k parte da 2 e arriva a N conincrementi di uno. Le restanti istruzioni dovrebbero essere chiare, non ci soffermeremonel descriverle essendo chiaro il contesto. Il blocco di istruzioni precedenti potrebbeconvenientemente essere memorizzato in un M-file. Nella Figura 2.1 mostriamo alcuniesempi di simulazione per differenti valori del parametro r ma per il medesimo datoiniziale x0=0.5.Ogni grafico e ottenuto tramite le istruzioni precedenti e utilizzando poi l’istruzioneplot(x,’-o’). Come si vede al variare di r abbiamo differenti comportamenti: da un va-lore limite positivo, r=1, ad un comportamento estremamente disordinato, r=3, passandoper dinamiche di tipo periodico.

2.3 Metodi iterativi

Metodi iterativi, ovverosia metodi che generano successioni {xk}k∈N con opportune pro-prieta, sono utilizzati in MATLAB per la ricerca degli zeri di equazioni. Vengono genera-lizzate procedure simili a quella adottata nel Metodo di bisezione per la dimostrazione del

Limiti e continuita 25

0 20 40 600

0.5

1

1.5

2

0 20 40 600

0.5

1

1.5

2

0 20 40 600

0.5

1

1.5

2

0 20 40 600

0.5

1

1.5

2

r=1 r=2

r=2.5 r=3

Figura 2.1 Alcuni esempi di simulazione del Modello di Ricker.

Teorema degli zeri. Vediamo un esempio tratto dalla Meccanica Celeste. Nel XVII secoloKeplero ha descritto la cinematica dei pianeti in moto intorno al Sole, riassumendone lecaratteristiche in tre leggi

• un pianeta descrive un orbita ellittica intorno al Sole, il Sole occupa uno dei duefuochi di tale ellisse;

• la linea che congiunge il pianeta al Sole descrive uguali aree in tempi uguali;

• esiste una correlazione tra il periodo T dell’orbita di ciascun pianeta e il raggiomedio Rm dell’orbita, R3

m = k T 2, dove k e una costante positiva.

La seconda legge riguardante le aree e descritta nella Figura 2.2, i punti A, B descrivonola posizione del pianeta in due tempi successivi. Il medesimo intervallo di tempo vieneimpiegato dal pianeta per passare da C a D, le aree A1 e A2 sono uguali. Per disegnarel’ellisse con semiassi a e b basta descrivere l’equazione in forma parametrica x(t) = a cos t,y(t) = b sin t, con t ∈ [0,2π]. Dalle leggi di Keplero si deduce la relazione tra il tempo t el’angolo al centro ϑ descritto dal pianeta,

2πt = T (ϑ− e sin (ϑ)),

dove T e il periodo di rivoluzione, e l’eccentricita dell’orbita ellittica, e ∈ (0,1). Pertrovare l’angolo ϑ in corrispondenza di un certo tempo t occorre risolvere l’equazionescritta sopra. La funzione MATLAB utile allo scopo e la funzione fzero. Per esempioconsideriamo un pianeta con T = 219 giorni ed eccentricita approssimativamente ugualea 0. 5 (il periodo corrisponde al periodo del pianeta Venere). Volendo conoscere l’angolocorrispondente al tempo t = 20 (giorni) possiamo operare come segue,

>> f=inline(’2*pi*20-219*theta+0.5*219*sin(theta)’);>> fzero(f,0.1)Zero found in the interval: [-0.924, 1.124].ans =

0.9925

26 Limiti e continuita

Sole

A1

A2

A

B

C

D

Figura 2.2 Illustrazione della seconda legge di Keplero.

Il valore 0.1 e il valore iniziale, essendo angoli esprime un valore in radianti, da cuiiniziare la ricerca dello zero. Non abbiamo fatto un’analisi teorica per l’unicita dellostesso zero, l’esistenza discende direttamente dal Teorema dei valori intermedi. Variandoil tempo t e possibile ricostruire la legge ϑ = ϑ(t), almeno per punti.

3Derivata di una funzione

3.1 Calcolo di derivate

Anche per le derivate il calcolo numerico tramite opportuni rapporti incrementali e unaquestione piuttosto delicata (sempre per via della necessita di basarsi sull’insieme deinumeri macchina e sul fatto che possiamo incappare in errori di cancellazione o altro).Memorizzati i valori di una funzione in un vettore y la funzione diff permette di calco-lare le differenze tra due valori successivi del vettore. Assegnato il passo h il rapportoincrementale approssimato si puo quindi valutare come ri=diff(y)/h. Il vettore ri hauna componente in meno rispetto al vettore y. Con le seguenti istruzioni (a parte lescritte) si ottiene il disegno nella Figura 3.1.

>> N=120;>> x=linspace(0,2*pi,N);>> h=x(2)-x(1);>> y=sin(x);>> ri=diff(y)/h;>> subplot(211),plot(x,y);>> subplot(212),plot(x(1:N-1),ri);

Abbiamo approssimato con un rapporto incrementale la derivata della funzione senonell’intervallo [0,2 ∗ pi]. L’istruzione subplot consente di visualizzare grafici multiplinella stessa finestra.

Per visualizzare il grafico dell’errore assoluto, differenza tra valore approssimato evalore vero in modulo, possiamo scrivere

>> xd=x(1:N-1);>> plot(xd,abs(ri-cos(xd)));

Non bisogna essere troppo ottimisti pensando che basti un passo h molto piccolo. Vedia-mo un esempio per una approssimazione della derivata della funzione f(x) = sin 10x inx = 1 con un passo piccolo, il valore vero e 10 cos 10 ≈ −8. 3907.

>> h=eps/4;>> (sin(10*(1+h))-sin(10))/hans =

0

Nell’ambito del calcolo simbolico la stessa funzione diff permette di calcolare la funzionederivata. Qualche esempio

>> f1=sym(’sin(10*x)’);>> f2=sym(’exp(x) * cos(2*x)’);>> f3=sym(’x^2+1/x’);>> diff(f1)

28 Derivata di una funzione

0 1 2 3 4 5 6 7−1

−0.5

0

0.5

1

0 1 2 3 4 5 6 7−1

−0.5

0

0.5

1

sin(x)

diff(sin(x))/h

Figura 3.1 Approssimazione di sin′(x).

ans = 10*cos(10*x)>> diff(f2)ans = exp(x)*cos(2*x)-2*exp(x)*sin(2*x)>> diff(f3)ans = 2*x-1/x^2

3.2 Ricerca di massimi e minimi

Per la ricerca dei massimi e dei minimi di una funzione e possibile utilizzare la funzionefmin o fminbnd. Tali funzioni calcolano il punto di minimo (globale) in un certo intervallodi una funzione specificata. Si noti che per il calcolo del massimo di una funzione f bastaconsiderare il problema del minimo della funzione −f . Per esempio con

>> fminbnd(’cos’,3,4)ans =

3.1416

calcoliamo il punto di minimo della funzione coseno nell’intervallo [3,4] ovverosia unaapprossimazione di π. I metodi per la ricerca del punto di minimo costruiscono unasuccessione di valori che sotto opportune ipotesi tende verso il punto di minimo. Con

>> fminbnd(’cos’,3,4,optimset(’TolX’,1e-12,’Display’,’iter’))Func-count x f(x) Procedure

1 3.38197 -0.971249 initial2 3.61803 -0.888633 golden3 3.23607 -0.995541 golden4 3.13571 -0.999983 parabolic5 3.1413 -1 parabolic6 3.14159 -1 parabolic7 3.14159 -1 parabolic8 3.14159 -1 parabolic9 3.14159 -1 parabolic

Optimization terminated successfully: the current x satisfies thetermination criteria using OPTIONS.TolX of 1.000000e-012 ans =

3.1416

Derivata di una funzione 29

vengono visualizzati i passi compiuti per calcolare una approssimazione di π con 12 cifresignificative esatte. La colonna Procedure indica il tipo di algoritmo selezionato passodopo passo, ovvero il tipo di strategia per la selezione del prossimo punto della succes-sione approssimante.

3.3 Problemi di ottimizzazione

Affrontiamo un problema di ottimizzazione. Supponiamo di dover gestire il traffico perun tratto di strada, per esempio un tratto di galleria. Tasso di flusso delle macchine edefinito come il numero di mezzi che passano in un posto fissato nell’unita di tempo. Ilflusso dei mezzi dipendera da vari fattori tra cui:(i) velocita del traffico;

(ii) distanza tra due mezzi consecutivi;

(iii) lunghezza dei mezzi.Consideriamo una velocita media v e una media delle distanze tra mezzo e mezzo denotatacon d, inoltre sia b la lunghezza media dei mezzi. Se ci mettiamo in un dato punto diosservazione, quando passa un mezzo il tempo medio che dobbiamo aspettare per vederneun altro e uguale a

T =d + b

v.

Il tasso r di flusso e uguale a

r =1

T=

v

d + b.

Assumiamo come parametri v e d, ovviamente la velocita v influisce sulla distanza d.Consideriamo la distanza minima d (per evitare incidenti) che dividiamo in due parti conrelazioni empiriche legate alla velocita (i valori sono comunque nell’unita di misura sceltaper le distanze). Abbiamo la distanza pensata dal conduttore,

dp = |v|,e la distanza utile per la fermata

df = |v|2/20.Posto d = dp + df , il tasso di flusso e dato da

r =v

(|v|+ |v|2/20 + b),

visto che la velocita sara considerata positiva possiamo togliere il modulo nei termini aldenominatore. Il grafico di r e mostrato nella Figura 3.2. Vogliamo trovare il punto dimassimo ovvero il massimo flusso ottenibile e valutare la corrispondente velocita mediada suggerire. Con una lunghezza media di b = 4 (metri) abbiamo,

>> rm=inline(’-v / (v+v^2/20+4)’);>> fminbnd(rm,0,100)ans =

8.9443

dove si considera la funzione rm = −r perche cerchiamo il massimo con un algoritmo cheseleziona solo punti di minimo. Evidentemente la risposta, se pensata come chilometriall’ora, non sembra molto realistica. Possiamo rilassare i vincoli e considerare la distanza

d = dp + α · df ,

con α ∈ (0,1) parametro. Possiamo utilizzare il calcolo simbolico per trovare il puntostazionario, che sara anche il punto di massimo per come e l’andamento del grafico delflusso.

30 Derivata di una funzione

0 10 20 30 40 50 60 70 80 90 1000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

velocità

r = mezzi /tempo

Figura 3.2 Andamento flusso del traffico.

>> syms v a b x ;>> ra=sym(’v / (v+a*v^2/20+b)’);>> d=diff(ra,v)d = 1/(v+1/20*a*v^2+b)-v/(v+1/20*a*v^2+b)^2*(1+1/10*a*v)>> vm=solve(d,’v’)vm = [ 2/a*5^(1/2)*(a*b)^(1/2)] [ -2/a*5^(1/2)*(a*b)^(1/2)]

La funzione solve risolve la corrispondente equazione, e a indica il parametro α. Tra ledue soluzioni selezioniamo solo quella positiva,

vm =p

(20b/α).

Per un fissato b e possibile generare una tabella con velocita per il flusso massimo e ladistanza minima permessa.

4Grafici di funzioni e approssimazioni

4.1 Grafici di funzioni

Per disegnare il grafico di una funzione abbiamo gia visto che possiamo utilizzare lafunzione ezplot. Con le seguenti istruzioni

>> subplot(221),ezplot(’cos(x)*exp(-x)’);>> subplot(222),ezplot(’cos(x)*exp(-x)’,[-10*pi 4*pi]);>> subplot(223),ezplot(’cos(x)*exp(-x)’,[-3*pi 3*pi]);>> subplot(224),ezplot(’cos(x)*exp(-x)’,[-pi pi]);

si ottengono i grafici di Figura 4.1 che si riferiscono alla stessa funzione f(x) = cos xe−x

ma con domini differenti. Bisogna fare attenzione a come collezionare le varie informazio-

−6 −4 −2 0 2 4−40

−20

0

20

40

60

x

cos (x) exp(−x)

−30 −20 −10 0 10

−3

−2

−1

0

1

2x 10

9

x

−5 0 5−200

−100

0

100

200

300

x−2 0 2

−15

−10

−5

0

x

Figura 4.1 Grafici della funzione cos xe−x in vari intervalli.

ni, i quattro grafici nella figura non sembrano della stessa funzione se non parzialmente.Si possono comunque analizzare anche i grafici delle derivate successive avvalendosi dellafunzione diff e ancora di ezplot.Proviamo a ripercorrere i passi suggeriti per lo studio di una funzione con l’ausilio del

32 Grafici di funzioni e approssimazioni

calcolo simbolico. Consideriamo la funzione

f(x) =x

2x + 1e−x.

Il dominio della funzione e determinato dai valori x tali che il denominatore non si annulli,(2x + 1) 6= 0 ⇒ x 6= −1/2. La funzione non sembra presentare particolari simmetrie.Per le intersezioni con gli assi coordinati, f(0) = 0, il grafico di f passa per l’origine(0,0). Valutiamo i limiti per x che tende agli estremi degli intervalli che compongono ildominio. Il limite sinistro e il limite destro per x → −1/2 valgono rispettivamente

>> syms x;>> f=sym(exp(-x)*x/(2*x+1));>> limit(f,x,-1/2,’left’)ans = inf>> limit(f,x,-1/2,’right’)ans = -inf

dove inf significa ∞ mentre -inf significa −∞. La retta x = −1/2 risulta essere unasintoto verticale. Per i limiti per x → ±∞ abbiamo invece:

>> limit(f,x,-inf,’right’)ans = inf>> limit(f,x,inf,’left’)ans =

0

quindi y = 0 risulta essere un asintoto orizzontale. Per la derivata prima abbiamo

>> fd= simplify(diff(f))fd = -exp(-x)*(2*x^2+x-1)/(2*x+1)^2

dove l’uso della funzione simplify effettua eventuali semplificazioni presenti dopo il cal-colo della derivata. Tramite la funzione solve possiamo calcolare gli zeri della derivata,che risultano essere gli zeri del numeratore e quindi tramite root gli seri del polinomio(2x2 + x− 1)

>> solve(fd)ans = [ -1] [ 1/2]

Il segno sara determinato dal segno del polinomio −(2x2 + x− 1), da cui

f ′(x) > 0 x ∈ (−1,1/2), f ′(x) x ∈ (−∞,− 1) ∪ (1/2, +∞).

Dal segno della derivata prima possiamo valutare gli intervalli in cui la funzione f cresceo decresce. Avremo inoltre che x = −1 risulta essere un punto di minimo locale mentreil punto x = 1/2 e un punto di massimo locale. Per stabilire se i punti stazionari sonopunti di massimo o di minimo possiamo anche utilizzare la derivata seconda

>> fd2=simplify(diff(fd))fd2 = exp(-x)*(4*x^3+4*x^2-3*x-6)/(2*x+1)^3>> x=-1; subs(fd2)ans =

8.1548>>x=0.5; subs(fd2)ans =

-0.4549

Per valutare il segno della derivata seconda occorre valutare il comportamento dellafunzione razionale che vi compare. Posto N(x) = (4x3 +4x2−3x−6), D(x) = (2x+1)3,per il numeratore N abbiamo

Grafici di funzioni e approssimazioni 33

>> roots([4 4 -3 -6])ans =

-1.0279 + 0.6034i-1.0279 - 0.6034i1.0558

quindi essendoci un’unica radice reale N(x) > 0 per x > 1. 0558, N(x) < 0 per x <1. 0558. Per il denominatore D abbiamo semplicemente D(x) < 0 per x < −1/2, mentreD(x) > 0 per x > −1/2. Riassumendo N(x)/D(x) > 0 per x ∈ (−∞,−1/2)∪ (1. 0558,+∞), mentre N(x)/D(x) < 0 per x ∈ (−1/2,1. 0558). Nel primo caso, essendo e−x > 0 perogni x, la derivata seconda e positiva e la funzione f risulta essere convessa, nel secondocaso la funzione risulta essere concava. I punti con ascissa x = −1/2 e x = 1. 0558risultano essere punti di flesso. Il grafico di f e schematizzato nella Figura 4.2.

−3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2

−6

−4

−2

0

2

4

6

8

10

12

14

x

exp(−x) x/(2 x+1)

Figura 4.2 Grafico qualitativo della funzione f(x) = e−xx/(2x + 1).

4.2 Approssimare dati e funzioni

Il polinomio di Taylor rappresenta solo una delle svariate possibilita a disposizione. Aven-do a disposizione dei dati, che rappresentano solo alcuni valori di una funzione oppureche rappresentano misure di un fenomeno di cui stiamo cercando un modello, possiamoconsiderare una approssimazione con un polinomio di un certo grado. A grandi lineeabbiamo a disposizione due strade

• pretendere che il polinomio passi esattamente per tutti i punti;

• lasciare che il polinomio indichi solo l’andamento qualitativo e parzialmente quan-titativo dei dati.

Nel primo caso si parla di interpolazione e il polinomio deve essere di grado opportuno(il grado deve essere uguale al numero di punti meno uno), inoltre tutte le ascisse deipunti devono essere differenti per punti diversi altrimenti non stiamo parlando di unafunzione p(x) (in corrispondenza di una ascissa non possiamo pretendere che il polinomiop assuma due valori differenti). Nel secondo caso si parla di fitting e le ricadute dal puntodi vista statistico saranno considerate nel seguito di questo libro.


Consideriamo, per esempio, i dati corrispondenti ai censimenti della popolazione degliStati Uniti dal 1790 al 1990. I dati sono accessibili dall’ambiente MATLAB perche fannoparte di insiemi di dati utilizzabili come test. Per leggerli e memorizzarli

>> load census ;

A questo punto vengono creati due vettori, pop e cdate: il primo contiene i dati deicensimenti, il secondo gli anni corrispondenti. In tutto abbiamo 21 dati rappresentatiin Figura 4.3. Possiamo trovare, con la funzione polyfit, il polinomio interpolatoree alcuni polinomi che tentano il fitting. Sperimentiamo alcune possibilita e, per ognimodello polinomiale scelto, facciamo la “previsione” per il 2000,

>> pi=polyfit(cdate,pop,20);prev(1)=polyval(pi,2000);>> pi=polyfit(cdate,pop,1);prev(2)=polyval(pi,2000);>> pi=polyfit(cdate,pop,2);prev(3)=polyval(pi,2000);>> pi=polyfit(cdate,pop,5);prev(4)=polyval(pi,2000);>> prevprev =

143.3750 219.4586 274.6221 281.0268

1750 1800 1850 1900 1950 20000

50

100

150

200

250

Censimenti

anni

Figura 4.3 Dati censimento USA.

Come si puo constatare le previsioni, indipendentemente dall’unita di misura (qui sonomilioni di persone), sono molto differenti per i vari polinomi (di grado 20, 1, 2, 5). Trovareil modello corretto o fare previsioni attendibili e un compito piuttosto delicato.Per ottenere valori intermedi rispetto ad una tabella data possiamo utilizzare la funzioneinterp1. Per esempio generando alcuni valori con la funzione humps possiamo provarevarie opzioni ed ottenere il grafico di Figura 4.4,

>> x=linspace(0,1,10); y=humps(x); xi=linspace(0,1,50);>> y1=interp1(x, y, xi, ’nearest’);>> y2=interp1(x, y, xi, ’linear’);>> y3=interp1(x, y, xi, ’spline’);>> y4=interp1(x, y, xi, ’cubic’);

Le etichette nearest, linear, spline, cubic stanno ad indicare il metodo con cui si eriempita la tabella partendo dai 10 valori iniziali. I primi due consistono nello scegliere ilvalore piu vicino o una media tra i valori piu vicini. Gli ultimi due metodi consistono in un

Grafici di funzioni e approssimazioni 35

0 0.2 0.4 0.6 0.8 10

20

40

60

80

100

0 0.2 0.4 0.6 0.8 10

20

40

60

80

100

0 0.2 0.4 0.6 0.8 10

20

40

60

80

100

0 0.2 0.4 0.6 0.8 10

20

40

60

80

100

nearest linear

spline cubic

Figura 4.4 Varie possibilita per interpolare una tabella.

approccio interessante che consiste nel costruire un certo numero di polinomi consecutivi,in questo caso di grado 3, che si raccordano con una certa regolarita. Questa idea estata pienamente sviluppata nell’ambito della progettazione assista dal calcolatore e neldisegno industriale. Le stesse idee si possono adottare in piu dimensioni e, per esempio,generando una prima tabella con la funzione MATLAB peaks, abbiamo

>> [x,y,z]=peaks(10);>> [x1,y1]=meshgrid(linspace(-3,3,40),linspace(-3,3,40));>> z1=interp2(x, y, z,x1,y1,’nearest’);>> z2=interp2(x, y, z,x1,y1, ’linear’);>> z3=interp2(x, y, z,x1,y1,’spline’);>> z4=interp1(x, y, z,x1,y1,’cubic’);

e si possono ottenere i grafici nella Figura 4.5. Nella figura ogni singolo grafico e ottenuto

−3−2

−10

12

3

−3

−2

−1

0

1

2

3−8

−6

−4

−2

0

2

4

6

8

−3−2

−10

12

3

−3

−2

−1

0

1

2

3−6

−4

−2

0

2

4

6

8

−3−2

−10

12

3

−3

−2

−1

0

1

2

3−8

−6

−4

−2

0

2

4

6

8

−3−2

−10

12

3

−3

−2

−1

0

1

2

3−8

−6

−4

−2

0

2

4

6

8

Figura 4.5 Interpolazione di una tabella bidimensionale.


tramite la funzione mesh, per esempio mesh(x1,y1,z1).

5Integrali, aree, primitive

5.1 Calcolo numerico di integrali

Anche nell’ambito dell’integrazione possiamo sfruttare le potenzialita del calcolo simbo-lico oppure utilizzare approssimazione numeriche. Quest’ultime si basano soprattuttosulla possibilita di considerare l’integrale come limite di una opportuna somma finita. Inparticolare si approssima un integrale definito come somma pesataZ b

a

f(x)dx =X

wkf(xk),

dove la scelta dei punti xk e dei pesi wk > 0 caratterizza le varie formule di integrazione(formule di quadratura numerica). Per esempio

>> f=inline(’exp(-x).*sin(x)’);>> quad(f,0,pi)ans =

0.5216>> f=inline(’sqrt(1+cos(x))’);>> quad8(f,0,1)ans =

1.3560

dove le funzioni per l’integrazione numerica sono quad e quad8 (la seconda fornisce usual-mente una migliore approssimazione anche se e piu costosa).Tramite il calcolo simbolico siamo in grado di calcolare integrali in modo esatto. Peresempio avremo

>> syms x>> int(x^2)ans = 1/3*x^3>> int(log(x),0,1)ans =

-1>> int(x*exp(-x))ans = -x*exp(-x)-exp(-x)>> int(cos(x),0,pi/2)ans =

1

38 Integrali, aree, primitive

5.2 Periodi e lunghezze

Vogliamo determinare il periodo di un pendolo semplice. Questo periodo e legato all’in-tegrale ellittico di prima specie definito come

K =

Z π/2

0

dϑp1− x2 sin2(ϑ)

,

dove x2 < 1 e un parametro reale. Per esempio con x2 = 0. 4 possiamo valutare

>> f=inline(’1./sqrt(1-0.4*sin(t))’,’t’);>> quad8(f,0,pi/2)ans =

1.8374

dove il risultato deve essere interpretato in base all’unita di misura del tempo stabilite.Una interessante applicazione dell’integrale definito consiste nel calcolo della lunghezzadi una curva espressa in forma parametrica. In generale per una curva nello spazioparametrizzata tramite le funzioni, supposte derivabili,

x(t), y(t), z(t), t ∈ [a,b],

il calcolo della lunghezza della curva stessa misurata dal punto (x(a),y(a),z(a)) al punto(x(b),y(b),z(b)) vale Z b

a

p(x′(t))2 + (y′(t))2 + (z′(t))2 dt.

Per esempio consideriamo per la curva

x(t) = sin(2t), y(t) = cos(3t), z(t) = t, t ∈ [0,2π].

Un grafico nello spazio tridimensionale della curva si ottiene tramite plot3 nella forma

>> t=linspace(0,2*pi);>> plot3(sin(2*t),cos(3*t),t);

ed e mostrato nella Figura 5.1. Per la valutazione dell’integrale

−1

−0.5

0

0.5

1

−1

−0.5

0

0.5

10

1

2

3

4

5

6

7 t=b

t=a

Figura 5.1 Qual’e la lunghezza di questa curva?.

Integrali, aree, primitive 39

Z 2π

0

q4 cos2(2t) + 9 sin2(3t) + 1 dt,

possiamo invece operare come segue.

>> curva=inline(’sqrt(4*cos(2*t).^2 +9 * sin(3*t).^2+1)’,’t’);>> lung=quad8(curva,0,2*pi)lung =

16.6635

Se la funzione e molto oscillante oppure se presenta delle singolarita, per esempio derivatainfinita, le funzioni di calcolo numerico possono incontrare delle difficolta che segnalanocon una sequenza di avvertimenti,

>> f=inline(’sin(1./x)’);>> quad8(f,0.01,1)Warning: Recursion level limit reached in quad8. Singularitylikely.

In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123

...

10 20 30 40 50 60 70 80 900

5

10

15

20

25

30

35

40

45

50

V

P

W

Ps

Pi

Figura 5.2 Calcolo di pressione media effettiva.

5.3 Integrali partendo da una tabella

Se abbiamo a disposizione una tabella di valori e non una funzione con la sua espressioneanalitica possiamo utilizzare formule che approssimano la tabella di valori con funzionilineari a tratti, funzione trapz. Per esempio da un diagramma sperimentale di volumeV e pressione P vogliamo ricavare l’area racchiusa in un ciclo (si pensi ad un ciclo di unmotore per esempio),

W =

ZPdV.

40 Integrali, aree, primitive

Alcuni dati sono rappresentati nella Figura 5.2, dove Ps e Pi indicano la pressionesuperiore e inferiore, V e il volume.I dati sono memorizzati prelevandoli dall’esterno, ovvero utilizzando la funzione load. Inogni caso abbiamo

>> V=[14 16.25 18.2 20.7 23 32.1 41.2 51 59 68 77 81.6 86.2];>> PI=[15.1 11 9.6 8 7.2 4.5 3.0 2 1.5 1 1.5 1.6 3.5];>> PS=[16.2 47 42.1 37.2 33 21 14.6 11.2 8.7 7.2 5. 4.5 3.5];

e il calcolo

>> load dati PS PI V>> trapz(PS-PI)ans =

181.1500

dove ovviamente, dovendo calcolare un’area compresa tra due curve, abbiamo fatto ladifferenza tra i valori di pressione.

6Integrali impropri, serie ed equazioni

differenziali

6.1 Serie di Taylor

Tramite le possibilita del calcolo simbolico possiamo valutare le somme parziali della seriedi Taylor di una data funzione. Per esempio con

>> taylortool(’atan(x)’)

si ottiene il disegno mostrato nella Figura 6.1. Non avendo specificato altro, MATLAB

−6 −4 −2 0 2 4 6

−3

−2

−1

0

1

2

3

Taylor Series Approximation

TN

(x) = x−1/3 x3+1/5 x5−1/7 x7

Figura 6.1 Un esempio di formula di Taylor.

considera la somma parziale fino a n = 7. Altre funzioni interessanti sono taylor (peril calcolo della serie di Taylor ma senza rappresentazione grafica) e symsum che permettela somma di una serie in modo formale. Per esempio

42 Integrali impropri, serie ed equazioni differenziali

>> s=sym(’exp(-x)’);>> taylor(s)ans = 1-x+1/2*x^2-1/6*x^3+1/24*x^4-1/120*x^5

Tra i molti modi per il calcolo di π vi sono anche le serie di potenze. La scelta del tipo dirappresentazione influisce notevolmente sulla velocita di convergenza ovvero sul numerodi termini necessari per aver un certo numero di cifre esatte. Per esempio la serie diLeibniz 1−1/3+1/5−1/7. . . converge a π/4 ma in modo estremamente lento. Abbiamoinvece

arctan(1/2) + arctan(1/3) = arctan(1) = π/4,

e, grazie alla rappresentazione in serie di Taylor della funzione arcotangente, si ottengonodelle serie efficienti dal punto di vista del calcolo.

>> t=sym(’atan(x)’);>> t1=taylor(t,10);>> p=4*(subs(t1,’x’,’1/2’)+subs(t1,’x’,’1/3’))p = 498668825/158723712>> sym(498668825/158723712,’d’)ans = 3.1417411974336890700953972554998

6.2 Un primo modello di crescita

Grande attenzione e dedicata alla soluzione di equazioni differenziali sia dal punto nu-merico che dal punto di vista simbolico. Le equazioni differenziali costituiscono infattiuno strumento importante nella modellizzazione matematica. Per l’aspetto simbolico epossibile utilizzare dsolve, per esempio

>> dsolve(’Dx = -a*x’)ans = exp(-a*t)*C1>> x = dsolve(’Dx = -a*x’,’x(0) = 1’)x = exp(-a*t)

Dal punto di vista numerico esistono diverse funzioni perche vi sono differenti difficoltanumeriche per l’approssimazione di equazioni differenziali (o sistemi di equazioni diffe-renziali). Consideriamo, per esempio, la simulazione del seguente modello di crescitatumorale

u′(t) =λu

α

�1−

�u

ϑ

�α�,

dove u(t) rappresenta il numero di celle tumorali presenti al tempo t, λ e un fattoredi crescita, α e ϑ sono dei parametri determinati dalle caratteristiche del tumore chestiamo analizzando. L’equazione scritta si puo analizzare dal punto di vista analitico,vogliamo fare una simulazione numerica. Modelli simili sono importanti, una volta chesi siano determinati i parametri, per simulare il comportamento del tumore nel caso didifferenti terapie. Consideriamo il caso α = 0. 5, λ = 0. 01, ϑ = 5 e u(0) = 100. Stiamonormalizzando le quantita in gioco, valori realistici, per esempio per la popolazione in-iziale, potrebbero essere dell’ordine di 1011. Dobbiamo definire una funzione che ricevain ingresso il tempo t e il valore u e ritorni il valore della derivata prime u′

>> f=inline(’0.02*u*(1-(u/50)^0.5)’,’t’,’u’)f =

Inline function:f(t,u) = 0.02*u*(1-(u/50)^0.5)

Con

>> [t,u]=ode45(f,[0 100],100);>> plot(t,u);

Integrali impropri, serie ed equazioni differenziali 43

si ottiene uno dei grafici nella Figura 6.2. La funzione ode45 necessita della funzionecon cui descrivere l’equazione differenziale, il vettore [0 100] indica l’intervallo di tempot ∈ [0,100] in cui valutare la soluzione, l’ultimo argomento denota la condizione iniziale.

0 10 20 30 40 50 60 70 80 90 10060

70

80

90

100

110

120

130

θ = 50

θ = 150

θ = 100

Figura 6.2 Un esempio di formula di simulazione di un modello di crescita tumorale.

6.3 Sistemi e modelli

Abbiamo gia accennato nel paragrafo sulle equazioni differenziali che un’equazione delsecondo ordine puo essere trasformata in un sistema di equazioni lineari. Per esempio,posto u′ = y, l’equazione u′′(t) + 4u′(t) + 13u(t) = 0 puo essere riscritta come(

u′(t) = y(t)

y′(t) = −4y(t)− 13u(t).

Risolviamo il problema di Cauchy per questa equazione con u(0) = 1, u′(0) = y(0) = 1.Abbiamo

>> f=inline(’[0 1 ; -13 -4]*u’,’t’,’u’);>> [t,u]=ode45(f,[0 5],[1 1]);>> plot(t,u);

e il risultato e mostrato in Figura 6.3. Infine un esempio storico, il modello preda preda-tore di Lotka-Volterra. Abbiamo due specie in competizione, u1 e u2, la cui dinamica egovernata dal sistema (

u′1(t) = ku1(t)− au1(t)u2(t)

u′2(t) = −ru1(t) + bu1(t)u2(t)

dove a,b,k,r sono costanti positive. La variabile u1 rappresenta le prede mentre la variabileu2 i predatori. I termini presenti nel sistema riguardano l’evolversi delle due popolazioniin base ai propri tassi di accrescimento e agli eventuali incontri preda-predatore. InMATLAB e gia disponibile una funzione lotka per la simulazione di questo sistema. Iparametri scelti sono k = 1, r = 1, a = 0. 01, b = 0. 02. Un esempio di simulazione emostrato in Figura 6.4 e ottenuto tramite

44 Integrali impropri, serie ed equazioni differenziali

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−2

−1.5

−1

−0.5

0

0.5

1

1.5

u

y

Figura 6.3 Soluzione di un’equazione del secondo ordine tramite un sistema di equazionidifferenziali.

>> [t,u]=ode45(’lotka’,[0 20],[80 100]);>> plot(t,u);

dove siamo partiti con 80 prede e 100 predatori. Come si puo vedere nella figura abbiamoun comportamento di tipo periodico.

0 2 4 6 8 10 12 14 16 18 2020

40

60

80

100

120

140

160

u2

u1

Figura 6.4 Un esempio di simulazione di un sistema preda-predatore.

7Vettori e matrici algebrici

7.1 Operazioni con matrici e vettori

Una delle caratteristiche principali di MATLAB e la sua abilita nell’operare con matrici.In Appendice B abbiamo descritto alcune di queste possibilita. Nel seguito vedremoalcuni semplici esempi. Inseriamo in MATLAB le seguenti matrici A e B

>> A = [3 5 7 ; 9 11 15; -2 1 5]A =

3 5 79 11 15

-2 1 5>> B = [1 3 ; 5 7 ; -1 0]B =

1 35 7

-1 0

Possiamo calcolare direttamente il prodotto 2B, AB e A3−A2 + 3I tramite le istruzioni

>> 2*Bans =

2 610 14-2 0

>> A*Bans =

21 4449 104-2 1

>>A^3-A^2+3*eye(3)ans =

514 1101 19591215 2596 4599-12 51 148

Se definiamo la matrice C data dalle prime due colonne di A possiamo sommarla a B

>> C = A(:,1:2)C =

3 59 11

-2 1

46 Vettori e matrici algebrici

>> C+Bans =

4 814 18-3 1

Dati ora i vettori colonna x e y

>> x = [1; -1; 3 ; 7]x =

1-137

>> y = [2; 4; -3; 9]y =

24

-39

Possiamo calcolare il prodotto scalare xT y ed il prodotto yT x tramite

>> x’*yans =

52>> x*y’ans =

2 4 -3 9-2 -4 3 -96 12 -9 27

14 28 -21 63

7.2 Alcune funzioni avanzate

Ricordiamo come tramite il Symbolic Toolbox sia possibile utilizzare anche simboli alposto di numeri. Ad esempio abbiamo

>> syms a>> a*Aans = [ 3*a, 5*a, 7*a] [ 9*a, 11*a, 15*a] [ -2*a, a, 5*a]>> syms x1 x2 x3>> x=[x1; x2; x3]x = [ x1] [ x2] [ x3]>> A*xans = [ 3*x1+5*x2+7*x3] [ 9*x1+11*x2+15*x3] [ -2*x1+x2+5*x3]

Inoltre il calcolo dell’inversa puo essere effettuato tramite l’istruzione MATLAB inv. Lebasi teoriche dell’algoritmo utilizzato saranno presentate nel prossimo capitolo. Abbiamoquindi

>> A=[1 2; 3 4];>> A1=inv(A)A1 =

-2.0000 1.00001.5000 -0.5000

>> A1*Aans =

1.0000 0

Vettori e matrici algebrici 47

0.0000 1.0000>> A*A1ans =

1.0000 00.0000 1.0000

7.3 Matrici e immagini

Una delle applicazioni principali delle matrici la si deve alle immagini digitali. Un’imma-gine digitale infatti non e altro che una matrice (o piu matrici) nella quale ogni elementorappresenta un’informazione sulla luminosita (o sul colore) della stessa. Ogni elementodella matrice viene detto pixel. In MATLAB una immagine monocromatica di m × npixels e memorizzata in una matrice con le medesime dimensioni, mentre una immaginea colori e memorizzata in una matrice multidimensionale m×n×3 (ossia tramite tre ma-trici ognuna delle quali contiene informazioni relativamente ad una delle tre componentidi colore fondamentali: rosso, verde e blu). Vediamo un semplice esempio di immaginemonocromatica. La successione di istruzioni

>> colormap(’gray’);>> I=eye(10);>> image(255*I);

realizza tramite il comando image che associa ad ogni elemento della matrice un colorel’immagine nella Figura 7.1. L’istruzione colormap consente di associare all’immagine

1 2 3 4 5 6 7 8 9 10

1

2

3

4

5

6

7

8

9

10

Figura 7.1 Matrice identita e immagini.

la mappa di colori delle scale di grigio nella quale lo 0 corrisponde al nero ed il 255 albianco. L’operazione di prodotto per uno scalare quindi corrisponde ad una variazionedella luminosita dell’immagine (si provi con image(50*I)).

Possiamo infine visualizzare come l’azione di trasposizione corrisponda ad una riflessionedell’immagine rispetto alla diagonale principale della matrice tramite le istruzioni

>> colormap(’gray’);>> A=zeros(7);>> A(3,4)=1; A(4,3:5)=1; A(5,2:6)=1;>> image(255*A);>> figure;>> image(255*A’);

48 Vettori e matrici algebrici

1 2 3 4 5 6 7

1

2

3

4

5

6

7

1 2 3 4 5 6 7

1

2

3

4

5

6

7

Figura 7.2 Trasposizione e immagini.

Il risultato e riportato nella Figura 7.2. Si provi help image per esplorare ulteriorifunzioni MATLAB per la gestione delle immagini.

8Sistemi lineari

8.1 Sistemi lineari

La risoluzione di sistemi lineari puo essere affrontata in MATLAB in modo estremamenteefficiente (originariamente MATLAB e stato progettato proprio per svolgere questo com-pito). Non deve stupire quindi il fatto che MATLAB riservi addirittura uno specificooperatore, indicato con ’\’ e detto backslash, per la risoluzione di un sistema. Vediamoun esempio

>> A=[1 3 2; -1 1 2; 2 4 -2];>> b=[3;-3;10];>> x=A\bx =

21

-1

L’uso come si vede e immediato ma bisogna fare attenzione che non stiamo “dividendo”A per b ma stiamo risolvendo il sistema lineare Ax = b. Abbiamo gia osservato comeil metodo di eliminazione di Gauss rappresenti un algoritmo applicabile direttamente alcalcolatore. L’operatore backslash nel caso di matrici quadrate si basa su questo algo-ritmo. Uno degli aspetti importanti nell’ambito della risoluzione di un sistema lineare el’efficienza di un metodo risolutivo, ossia quanto paghiamo in termini di costo computa-zionale per ottenere il risultato. Un metodo per misurare il costo computazionale e quellodi valutare il numero di operazioni necessarie al calcolatore per ottenere la soluzione. Perl’eliminazione di Gauss applicata ad una matrice quadrata di dimensione n il numero dioperazioni e n3/3.Nel caso di sistemi inconsistenti MATLAB risolve comunque il problema approssimandola soluzione nel senso dei minimi quadrati. Parleremo di questo nel prossimo capitolo,ma e importante ricordare da subito che l’operatore backslash fornisce sempre ’qualcosa’che non e’ detto sia il risultato cercato. Per esempio

>> A=[1 3 2; -1 1 2; 2 4 -2; 1 3 2];>> b=[3;-3;10; -3];>> x=A\bx =

0.12501.3750-2.1250

>> A*xans =

-0.0000-3.000010.0000

50 Sistemi lineari

-0.0000

La soluzione x di fatto non e quindi una soluzione. Infatti il sistema in esame risulta essereinconsistente e MATLAB automaticamente ne ha calcolato una soluzione nel senso deiminimi quadrati. La verifica dell’inconsistenza puo essere effettuata tramite il Teoremadi Rouche-Capelli e la funzione rank di MATLAB. Utilizzando i valori di A e b definitiin precedenza avremo infatti

>> rank(A)ans =

3>> rank([A b])ans =

4

Nel caso di sistemi che ammettono infinite soluzioni l’operatore backslash si comporta inmodo diverso a seconda delle situazioni. Se la matrice e quadrata ci segnala il fatto che ildeterminante e nullo e non calcola nessuna soluzione, se invece la matrice e rettangolareci segnala un’eventuale deficienza del rango rispetto al rango massimo possibile e poicalcola una sola soluzione del sistema.

>> A=[1 3 2; -1 1 2; 2 6 4];>> b=[3;-3;6];>> x1=A\bWarning: Matrix is singular to working precision. x1 =

InfInfInf

>> A=[1 3 2; -1 1 2; 2 6 4; 4 12 8];>> b=[3;-3;6;12];>> x1=A\bWarning: Rank deficient, rank = 2 tol = 1.2243e-014. x1 =

3.00000.0000

0

La forma generale della soluzione puo essere calcolata sfruttando il Teorema di struttura(continuando l’esempio precedente)

>> x0=null(A,’r’)x0 =

1-11

>> A*x0ans =

0000

Ogni vettore del tipo x = x1 + λx0 e dunque soluzione del nostro problema. La funzionenull applicata ad A calcola una soluzione del sistema Ax = 0. L’opzione ’r’ forza MAT-LAB a fornire una soluzione in formato razionale.Nel caso di matrice quadrata, dopo averne eventualmente verificato il rango, sara suffi-ciente risolvere il sistema caratterizzato dalle equazioni l.i.

>> A=[1 3 2; -1 1 2; 2 6 4];>> b=[3;-3;6];>> rank(A)ans =

2

Sistemi lineari 51

>> rank([A b])ans =

2>> x1=A(1:2,:)\b(1:2)x1 =

3.00000.0000

0

e procedere come nel caso precedente.Tramite il calcolo simbolico possiamo risolvere anche sistemi nei quali la soluzione dipendeda uno o piu parametri. Per esempio

>> syms x alpha>> A=[alpha 1 1; 1 alpha 1; 1 1 alpha];>> b=[1; 2; 3];>> A\bx = [ (alpha-4)/(alpha-2+alpha^2)] [ 2/(alpha+2)][ 3/(alpha-2+alpha^2)*alpha]

La soluzione quindi e unica per α 6= −2,1. Possiamo esplorare i casi esclusi tramite ilcomando subs

> A1=subs(A,alpha,-2)A1 =

-2 1 11 -2 11 1 -2

>> A1\bWarning: System is inconsistent. Solution does not exist.> In C:\MATLAB6p1\toolbox\symbolic\@sym\mldivide.m at line 29ans = [ Inf] [ Inf] [ Inf]>> A2=subs(A,alpha,1)A2 =

1 1 11 1 11 1 1

>> A2\bWarning: System is inconsistent. Solution does not exist.> In C:\MATLAB6p1\toolbox\symbolic\@sym\mldivide.m at line 29ans = [ Inf] [ Inf] [ Inf]

Quindi che per α = −2,1 il sistema e inconsistente.

8.2 Determinante e inversa

Il calcolo del determinante puo essere effettuato tramite il comando det. Per esempio

>> A=[0 1 3; 2 4 1; 1 -1 0];>> det(A)ans =

-17

Ritroviamo le regole delle diagonali per matrici 2×2 e pe rmatrici 3×3 tramite il calcolosimbolico

>> syms a b c d>> A=[a b; c d];>> det(A)ans = a*d-b*c

52 Sistemi lineari

>> syms a1 a2 a3 b1 b2 b3 c1 c2 c3>> A=[a1 a2 a3; b1 b2 b3; c1 c2 c3];>> det(A)ans = a1*b2*c3-a1*b3*c2-b1*a2*c3+b1*a3*c2+c1*a2*b3-c1*a3*b2

Infine possiamo calcolare l’inversa tramite la funzione inv. Avremo per esempio

>> A=[2 1 1; 4 1 0; -2 2 1];>> A1=inv(A)A1 =

0.1250 0.1250 -0.1250-0.5000 0.5000 0.50001.2500 -0.7500 -0.2500

Concludiamo osservando che a volte puo essere conveniente forzare MATLAB alla visua-lizzazione tramite razionali utilizzando format rat. L’inversa precedente risulta

>> format rat>> A=[2 1 1; 4 1 0; -2 2 1];>> A1=inv(A)A1 =

1/8 1/8 -1/8-1/2 1/2 1/25/4 -3/4 -1/4

8.3 Interpolazione

La risoluzione di sistemi lineari e un problema che interviene in numerose applicazioni.Per esempio nel caso dell’interpolazione polinomiale che abbiamo gia affrontato nel labo-ratorio Matlab del Capitolo 4. Consideriamo una funzione f incognita di cui conosciamosolo una serie di valori yi = f(xi) in corrispondenza dei punti xi, i = 1, . . . ,n. Vogliamoricostruire la funzione f utilizzando una sua approssimazione tramite un polinomio. Unesempio elementare e dato dal modo in cui MATLAB realizza il grafico della funzionetramite plot. La funzione continua visualizzata da MATLAB raccordando tutti i punticon una retta e detta spline lineare e rappresenta un polinomio a tratti di grado 1. Nelcaso generale cerchiamo un polinomio p(x) tale che p(xi) = yi per ogni i = 1, . . . ,n. Sipuo dimostrare che esiste unico un polinomio di grado n− 1 che soddisfa alle condizioniprecedenti. Infatti se

p(x) = a0 + a1x + . . . + an−1xn−1,

dalle condizioni p(xi) = yi, i = 1, . . . ,n otteniamo il sistema lineare di n equazioni nellen incognite a0, a1, . . ., an−18>>>>>>>>>>>>>>:

1 x1 x21 . . . xn−1

1

1 x2 x22 . . . xn−1

2

......

... . . ....

1 xn x2n . . . xn−1

n

9>>>>>>>>>>>>>>;8>>>>>>>>>>>>:

a0

a1

...

an−1

9>>>>>>>>>>>>; =

8>>>>>>>>>>>>:y1

y2

...

yn

9>>>>>>>>>>>>; .

La matrice del sistema, detta matrice di Vandermonde, ha determinante diverso da zerose i punti xi sono tutti distinti ed in tal caso esistono unici a0, a1, . . ., an−1, ossia esisteun unico polinomio interpolante. Vediamo un esempio in MATLAB.

>> x=[-1 0 2 5]’;>> y=[3 0 4 1]’;>> V=[x.^0 x x.^2 x.^3];>> a=V\y;>> xi=linspace(-1,5);

Sistemi lineari 53

>> yi=polyval(a(4:-1:1),xi);>> plot(x,y,’-o’,xi,yi,’-’);

I risultati sono riportati nella Figura 8.1.

−1 0 1 2 3 4 5−1

0

1

2

3

4

5

6

7

spline lineare

polinomiointerpolante

Figura 8.1 Interpolazione lineare a tratti e polinomiale.

9Autovalori e autovettori

9.1 Calcolo di autovalori

La funzione principale per il calcolo degli autovalori in MATLAB e la funzione eig. L’usoe abbastanza immediato

>> A=[4 -3 1; 2 -1 1; 0 0 2];>> eig(A)ans =

2.00001.00002.0000

Possiamo inoltre calcolare gli autovettori tramite

>> [V,D]=eig(A)V =

0.8321 0.7071 0.22760.5547 0.7071 0.4411

0 0 0.8681D =

2.0000 0 00 1.0000 00 0 2.0000

La matrice V e formata da autovettori unitari ed e tale che AV = V D.Tramite il calcolo simbolico possiamo per esempio calcolare la forma degli autovettoriper una generica matrice di dimensione 2

>> syms a b c d>> A=[a b; c d];>> eig(A)ans = [ 1/2*a+1/2*d+1/2*(a^2-2*a*d+d^2+4*b*c)^(1/2)] [1/2*a+1/2*d-1/2*(a^2-2*a*d+d^2+4*b*c)^(1/2)]

9.2 Fenomeni di risonanza

Una dei risvolti fisici di maggior rilievo degli autovalori lo si ha nell’ambito dello studiodelle frequenze di risonanza. Consideriamo un semplice modello di edificio a n piani rap-presentato da un sistema di n masse uguali sovrapposte. Supponiamo che ogni piano siacollegato ai piani adiacenti in modo elastico. Assumiamo che ogni massa possa scivolaresenza frizione e che la tendenza a conservare la configurazione iniziale sia proporzionale

56 Autovalori e autovettori

allo spostamento orizzontale del piano. Questo po essere modellizzato tramite l’azionesu ogni piano di due forze elastiche indotte dai piani adiacenti del tipo F = −kx dove xindica lo spostamento relativo del piano.Se indichiamo con xi(t) la coordinata della massa che rappresenta il piano i-esimo avremoche su ogni piano interno agiscono le due forze

F(1)i = −k(xi − xi−1), F

(2)i = −k(xi − xi+1), i = 2, . . . ,n− 1.

Sul primo piano agira la forza −kx1−k(x1−x2) mentre sull’ultimo solo la forza −k(xn−xn−1). L’applicazione della legge di Newton F = ma = mx′′ fornisce le equazioni delmoto in forma matriciale

Mx′′ = Kx,

dove x = (x1, . . . ,xn)T , M = mI con I matrice identita e la matrice di massa n× n e Ke la matrice di stiffness

K =

8>>>>>>>>>>>>>>>>>:−2k k 0 . . . 0

k −2k k . . . 0

0. . .

. . .. . . 0

0 . . . k −2k k

0 . . . 0 k −k

9>>>>>>>>>>>>>>>>>; .

A questo punto simuliamo l’azione di un terremoto tramite un’oscillazione orizzontale delterreno avente la legge E cos(ωt) con ampiezza E e pulsazione ω. Derivando due volteotteniamo la forza inerziale F = ma = mEω2 cos(ωt) che agisce su ogni piano dell’edificio.Il nostro modello finale sara quindi descritto dal seguente sistema di equazioni differenzialiordinarie non omogeneo del secondo ordine

x′′ = Ax + Eω2 cos(ωt)e,

con A = M−1K ed e = (1, . . . ,1).Gli autovalori della matrice A hanno un ruolo chiave nella risoluzione del precedentesistema. Infatti tale soluzione ha la forma x(t) = xp(t)+xc(t) dove xp(t) e una soluzioneparticolare del problema e xc(t) la soluzione del corrispondente sistema omogeneo. Lasoluzione del corrispondente sistema omogeneo puo essere effettuata tramite il calcolodegli autovalori λi della matrice A e dei rispettivi autovettori vi, mentre una soluzioneparticolare e data da xp(t) = C cos(ωt) dove C e il vettore soluzione del sistema lineare(A+ω2I)C = −E. Quest’ultimo sistema ha soluzione se e solo se la matrice A+ω2I e nonsingolare, ossia se ω2 6= −λi per ogni autovalore di A. Nei sistemi meccanici tipicamentegli autovalori di A sono tutti negativi ed i valori ωi =

√−λi sono detti frequenze naturalidel sistema. Nel caso in cui ω sia una delle frequenze naturali il sistema risulta singolareed abbiamo il fenomeno cosiddetto della risonanza.

Al fine di studiare il problema numericamente con MATLAB vediamo di calcolarele frequenze naturali del sistema. Supponiamo di avere un edificio a 10 piani e che lamassa di ogni piano sia di 20 tonnellate. La forza elastica che tiene i piani legati tra lorosia caratterizzata da k = 100 tonnellate/metro. La matrice A assume quindi la forma

A =

8>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>:

−10 5 0 0 0 0 0 0 0 0

5 −10 5 0 0 0 0 0 0 0

0 5 −10 5 0 0 0 0 0 0

0 0 5 −10 5 0 0 0 0 0

0 0 0 5 −10 5 0 0 0 0

0 0 0 0 5 −10 5 0 0 0

0 0 0 0 0 5 −10 5 0 0

0 0 0 0 0 0 5 −10 5 0

0 0 0 0 0 0 0 5 −10 5

0 0 0 0 0 0 0 0 5 −5

9>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>;

Autovalori e autovettori 57

Possiamo calcolare gli autovalori λi della matrice A, le relative frequenze naturali ωi edi periodi 2π/ωi

>> A=-10*eye(10,10)+5*diag(ones(1,9),-1)+5*diag(ones(1,9),1);>> A(10,10)=-5;>> lambda=eig(A);>> omega=sqrt(-lambda);>> [lambda omega 2*pi/omega]ans=

-5.0000 2.2361 2.8099-16.2349 4.0293 1.5594-13.6534 3.6951 1.7004-18.2624 4.2735 1.4703-19.5557 4.4222 1.4208-10.7473 3.2783 1.9166-7.7748 2.7883 2.2534-2.6695 1.6339 3.8456-0.9903 0.9951 6.3138-0.1117 0.3342 18.8005

Il modello quindi avra un comportamento critico se la frequenza di oscillazione del terrenorisulta prossima ad una delle frequenze naturali del sistema. Nella Figura 9.1 mostriamogli spostamenti nel tempo dell’ultimo piano dell’edificio nel caso di una frequenza naturalee di una frequenza non naturale.

0 10 20 30 40 50 60 70 80 90 100−20

−15

−10

−5

0

5

10

15

20

tempo

x n(t)

0 10 20 30 40 50 60 70 80 90 100−20

−15

−10

−5

0

5

10

15

20

x n(t)

tempo

Figura 9.1 Vibrazioni di risonanza all’ultimo piano per 2π/ω = 3. 8456 (sinistra) e casonon risonante per 2π/ω = 5 (destra) in funzione del tempo nel modello digrattacielo a 10 piani.

9.3 Minimi quadrati e approssimazione di funzioni

Innanzitutto osserviamo che il calcolo delle norme di vettore e di matrice e implementatoin MATLAB tramite la funzione norm. Per esempio

>> x=[1 -2 3 4 -8];>> norm(x,1)ans =

18

58 Autovalori e autovettori

>> norm(x,inf)ans =

8>> norm(x,2)ans =

9.6954>> A=[1 2 4; -1 7 3; -11 2 0];>> norm(A,1)ans =

13>> norm(A,inf)ans =

13>> norm(A,2)ans =

11.5506

Consideriamo un problema analogo a quello dell’interpolazione polinomiale visto nel la-boratorio Matlab dei Capitoli 4 e 8. Data una funzione f di cui sono noti una serie divalori yi = f(xi) nei nodi xi, i = 1, . . . ,n cerchiamo il polinomio di grado fissato chemeglio approssima la funzione nel senso dei minimi quadrati. In altre parole dato

p(x) = a0 + a1x + a2x2 + . . . + amxm,

imponiamo le condizioni p(xi) = yi, i = 1, . . . ,n. Sappiamo che se m = n − 1 ed inodi sono distinti abbiamo un sistema lineare con una unica soluzione. Se m < n − 1 ilsistema e sovradeterminato e non ammette soluzione, quindi ne cerchiamo una nel sensodei minimi quadrati. Possiamo riutilizzare i dati del Capitolo 8 con n = 4 e cercare laretta (m = 1) e la parabola (m = 2) dei minimi quadrati. Possiamo scrivere

>> x=[-1 0 2 5]’;>> y=[3 0 4 1]’;>> V1=[x.^0 x];>> a1=V1\y;>> V2=[x.^0 x x.^2];>> a2=V2\y;>> xi=linspace(-1,5);>> y1i=polyval(a1(2:-1:1),xi);>> y2i=polyval(a2(3:-1:1),xi);>> plot(x,y,’o’,xi,y1i,xi,y2i);

I risultati sono riportati nella Figura 9.2.

Autovalori e autovettori 59

−1 0 1 2 3 4 5−1

0

1

2

3

4

5

6

7

m=1

m=2

m=3

Figura 9.2 Retta dei minimi quadrati, parabola dei minimi quadrati e polinomiointerpolante.

10Vettori geometrici, rette e piani

10.1 Trasformazioni lineari

Abbiamo gia parlato di immagini digitali e della loro rappresentazione tramite matrici.Vediamo a questo punto come sia possibile trasformare un immagine utilizzando oppor-tune trasformazioni lineari o combinazioni di queste. In MATLAB e possibile disegnare a“mano libera” tramite il mouse una figura per punti. Per fare cio utilizziamo il comandoginput nella forma

>>[x,y]=ginput;>> plot(x,y,’o’);

La “raccolta di punti” terminera non appena avremo premuto il tasto “Invio”. In alter-nativa per un migliore inserimento di un numero prefissato di n punti con visualizzazionedei dati si possono usare le istruzioni

>> n=10;>> axis([-1 1 -1 1])>> hold on>> for i=1:n [x(i),y(i)]=ginput(1); plot(x(i),y(i),’o’); end

A questo punto tramite le matrici bidimensionali che caratterizzano le trasformazionilineari di R2 in R2 possiamo agire geometricamente sul disegno. Definiamo quindi lematrici di riflessione rispetto l’asse x, rispetto la retta y = x e di rotazione tramite icomandi

>> A1=[1 0; 0 -1];>> A2=[0 1; 1 0];>> theta=pi/2;>> c=cos(theta); s=sin(theta);>> A3=[c -s; s c];

Vediamo l’azione delle trasformazioni lineari sul nostro disegno

>> D=[x’;y’];>> D1=A1*D;>> D2=A2*D;>> D3=A3*D;>> x1=D1(1,:); y1=D1(2,:);>> x2=D2(1,:); y2=D2(2,:);>> x3=D3(1,:); y3=D3(2,:);>> plot(x,y,’o-’,x1,y1,’o-’);>> figure>> plot(x,y,’o-’,x2,y2,’o-’);>> figure

62 Vettori geometrici, rette e piani

>> plot(x,y,’o-’,x3,y3,’o-’);

L’uso del comando figure consente l’apertura di piu finestre grafiche. Nella Figura 10.1riportiamo i risultati ottenuti.

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Figura 10.1 Riflessioni rispetto all’asse delle x e alla retta y = x e rotazione di π/2 insenso antiorario di un disegno.

Possiamo chiaramente effettuare le stesse operazioni anche su figure geometriche bendefinite. Per esempio possiamo agire su un ellisse definendo la matrice di punti D come

>> t = linspace(0,2*pi);>> x=1+3*cos(t); y=1+2*sin(t);>> D = [x;y];

Nella Figura 10.2 sono illustrate le diverse azioni delle trasformazioni.

−4 −3 −2 −1 0 1 2 3 4 5 6−4

−3

−2

−1

0

1

2

3

4

−4 −3 −2 −1 0 1 2 3 4 5 6−4

−3

−2

−1

0

1

2

3

4

−4 −3 −2 −1 0 1 2 3 4−4

−3

−2

−1

0

1

2

3

4

Figura 10.2 Riflessioni rispetto all’asse delle x e alla retta y = x e rotazione di π/2 insenso antiorario di un ellisse.

10.2 Grafici in piu dimensioni

Vediamo come realizzare con MATLAB grafici di funzioni di piu variabili. Va precisatoche MATLAB fornisce numerose possibilita grafiche e che in questa sede ci limiteremoad analizzare le principali. In modo analogo a quanto fatto in una dimensione bisognera

Vettori geometrici, rette e piani 63

innanzitutto definire una matrice di punti in corrispondenza dei quali valutare la funzione.Per esempio consideriamo la funzione

f(x,y) = (x− y) sin(x2 + y2),

nell’intervallo [−2,2]× [−2,2]. Tramite i comandi

>> m=40;n=40;>> x=linspace(-10,10,m);>> y=linspace(-10,10,n);>> [X,Y]=meshgrid(x,y);

generiamo due matrici X e Y di m × n punti tali che l’elemento xij = (xj ,xj) e yij =(yi,yi).Per creare il grafico della superficie e sufficiente aggiungere le istruzioni

>> Z=Z=(X-Y).*sin(X.^2+Y.^2);>> mesh(X,Y,Z);

La matrice Z avra la proprieta

Z(i,j) = f(X(i,j),Y (i,j)), i = 1, . . . ,m j = 1, . . . ,n.

Il risultato ottenuto e riportato in Figura 10.3 assieme al grafico delle curve di livello.Le curve di livello si possono ottenere tramite la funzione contour(X,Y,Z). La rappre-

−2

−1

0

1

2

−2

−1

0

1

2−4

−3

−2

−1

0

1

2

3

4

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Figura 10.3 Grafico di superficie e curve di livello di f(x,y) = (x− y) sin(x2 + y2).

sentazione grafica realizzata da mesh e l’analogo in piu dimensioni del comando plot. Inquesto caso si raccordano i punti a quattro a quattro con “pezze” piane (parti di piano chepassano per i quattro punti). Utilizzando invece la funzione surf le varie pezze vengonocolorate opportunamente. In generale il grafico di superfici puo essere personalizzato indiversi modi. Le opzioni piu comuni sono date dai comandi axis, title, xlabel, ylabele zlabel per il cui uso, analogo al caso bidimensionale, si rimanda all’help in linea.Nel caso di campi vettoriali possiamo utilizzare le funzioni quiver in R2 e quiver3 inR3. Per esempio con le istruzioni

>> x=linspace(-2,2,20);>> y=linspace(-2,2,20);>> [X,Y]=meshgrid(x,y);>> px=X.*exp(-X.^2-Y.^2);>> py=-Y.*exp(-X.^2-Y.^2);>> quiver(X,Y,px,py)

otteniamo il grafico della Figura 10.4.

64 Vettori geometrici, rette e piani

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Figura 10.4 Campo vettoriale di F(x,y) = (x exp(−x2 − y2),− y exp(−x2 − y2)).

11Probabilita

11.1 Variabili aleatorie discrete

Sia X una v.a. discreta che puo assumere i valori 3,4,8 con probabilita 1/4,1/2,1/4rispettivamente. Inseriamoli come vettori

>> X = [3;4;8]X =

348

>> P=[1/4;1/2;1/4]P =

0.25000.50000.2500

Controlliamo che P possa essere una legge di probabilita. Possiamo verificare che non cisono valori negativi tramite il comando find

>> find( (P<0) )ans =

Empty matrix: 0-by-1

Inoltre la somma degli elementi di P deve dare 1

>> sum(P)ans =

1

Possiamo realizzare il grafico a barre della distribuzione di probabilita con il comandobar nella forma

>> bar(X,P)

Sempre con il comando find possiamo calcolare P ({X > 3}). Infatti, troviamo gli“indici buoni” di X per cui X > 3

>> ind_buo = find( ( X > 3 ) )ind_buo =

23

e quindi sommiamo le probabilita corrispondenti

66 Probabilita

3 4 80

0.1

0.2

0.3

0.4

0.5

Figura 11.1 Uso del comando bar per visualizzare la distribuzione di probabilita.

>> PR= sum ( P(ind_buo) )PR =

0.7500

Calcoliamo ora il valore atteso di X. Poiche E(X) =P

i xipxi , in forma matriciale

si trattera del prodotto scalare xT px, ossia

>> EX = X’ * PEX =

4.7500

Grazie al Teorema del valore atteso di g(X), per calcolare ad esempio E(sin(X)),possiamo scrivere

>> EsinX = ( sin(X) )’ * PEsinX =

-0.0958

Calcoliamo la varianza di X

>> VAR_X = ( ( X - EX ).^2 )’ * PVAR_X =

3.6875

e attraverso la formula della varianza V ar(X) = E(X2)− (E(X))2

>> Var_X = ( X.^2 )’ * P - ( EX ).^2Var_X =

3.6875

Avremmo anche potuto usare la funzione MATLAB var nella forma var(X,P) ed ottenerelo stesso risultato.Ricordiamo come tramite il Symbolic Toolbox sia possibile utilizzare anche simboli alposto di numeri.

>> syms x1 x2 x3 real>> syms p1 p2 positive>> x=[x1; x2; x3]x =[ x1][ x2][ x3]

>> p=[p1; p2; 1-(p1+p2)]p =

Probabilita 67

[ p1][ p2][ 1-p1-p2]

Abbiamo inserito p3 = 1 − (p1 + p2) in modo da garantire che p1 + p2 + p3 = 1. Ora,controlliamo la formula della varianza

Var Xz }| {E(X2)| {z }EX2

−(E(X)| {z }EX

)2

=

VAR Xz }| {E((X − E(X)2) = V ar(X)

>> EX = x’ * p; %calcola EX>> EX2= (x.^2)’ * p; %calcola EX2>> Var_X = EX2 - (EX).^2; %calcola Var_X>> VAR_X=( (x-EX).^2 )’ *p; %calcola VAR_X>> simplify(VAR_X - Var_X) %calcola-semplifica i conti VAR_X-Var_Xans = 0

11.2 Variabili aleatorie doppie X e Y

Analizziamo ora la coppia di v.a. discrete (X,Y ), dove X puo assumere i valori 3,4,8 eY puo assumere i valori −1,0,3,4. Inseriamoli come vettori

>> X = [ 3 ; 4 ; 8 ];>> Y = [ -1 ; 0 ; 3 ; 4 ];

La legge congiunta verra assegnata mediante la matrice

P = pij = p(X,Y )(xi,yj), i = 1,2,3, j = 1,2,3,4.

Inseriamo la matrice P

>> P= 1/10*ones(3,4)P =

0.1000 0.1000 0.1000 0.10000.1000 0.1000 0.1000 0.10000.1000 0.1000 0.1000 0.1000

Queste due matrici rappresenta una legge congiunta. Infatti i suoi elementi sono nonnegativi ma

Pi

Pj p(X,Y )(xi,yj) > 1

>> find( (P<0) )ans =

Empty matrix: 0-by-1>> sum( sum( P ) )ans =

1.2000

Inseriamo la “buona matrice” P

>> P(3,3:4)= 0P =

0.1000 0.1000 0.1000 0.10000.1000 0.1000 0.1000 0.10000.1000 0.1000 0 0

Questa e una legge congiunta, poiche adesso

68 Probabilita

>> sum( sum ( P ) )ans =

1

Ricordando che il comando sum somma per colonne, troviamo le distribuzioni marginali

>> P_Y = ( sum ( P ) )’ %distribuzione marginale di YP_Y =

0.30000.30000.20000.2000

>> P_X = ( sum ( P’ ) )’ %distribuzione marginale di XP_X =

0.40000.40000.2000

Verifichiamo se X e Y sono indipendenti, ossia

>> P_X * ( P_Y )’ %matrice con fenomeni X e Y indipendentians =

0.1200 0.1200 0.0800 0.08000.1200 0.1200 0.0800 0.08000.0600 0.0600 0.0400 0.0400

Per il criterio di indipendenza, i due eventi X e Y non sono indipendenti, essendoquest’ultima matrice diversa da P.

Costruiamo ora le due matrici XX e YY tramite il comando meshgrid

>> [XX,YY]=meshgrid(X,Y)XX =

3 4 83 4 83 4 83 4 8

YY =-1 -1 -10 0 03 3 34 4 4

Queste matrici servono per il calcolo di una funzione g(X,Y ). Se vogliamo avere adesempio la matrice Sij = sin(xi + 3yj), bastera scrivere

>> S = sin( XX + 3 * YY )S =

0 0.8415 -0.95890.1411 -0.7568 0.9894-0.5366 0.4202 -0.96140.6503 -0.2879 0.9129

Calcoliamo, il valore atteso di g(X,Y )

>> sum( sum( S’ .* P ) )ans =

0.0502

Calcoliamo ora E(X + Y ) e mostriamo che e uguale a E(X) + E(Y )

>> sum( sum( ( XX + YY ) .* P ) )ans =

5.5000

Probabilita 69

>> EX = sum( X’ * P_X )EX =

4.4000>> EY = sum( Y’ * P_Y )EY =

1.1000>> EX + EYans =

5.5000

11.3 Catene di Markov

Nella dinamica delle catene di Markov abbiamo una matrice P , detta matrice di transi-zione, i cui elementi rappresentano delle probabilita di transizione tra diversi stati. Piuprecisamente pij e la probabilita condizionata che il sistema si trovi “domani” nello statoj essendo “oggi” nello stato i. La somma di ogni riga di P e uguale ad 1. Indicato con p0

il vettore iniziale di probabilita degli stati si e interessati a sapere cosa succede al variaredi n al vettore pn definito come

pn+1 = Ppn = P · · ·P| {z }n volte

p0 = P np0, n ≥ 0.

La variabile n caratterizza una successione temporale discreta. In sostanza lo stato delsistema al passo n + 1 dipende solo dallo stato del sistema al passo n. Tale proprieta esoddisfatta da quasi tutti i sistemi nell’ambito della fisica e della scienza dei calcolatori1.Si noti che anche P n e una matrice di transizione avente righe a somma 1, un suo genericoelemento di indici i e j rappresenta quindi la probabilita che il sistema si trovi dopo npassi nello stato j trovandosi nello stato i all’istante iniziale. Si e interessati a saperecosa succede per n crescente. Cosa possiamo dire sul “comportamento” della matrice P n

e di pn? Si puo dimostrare che matrici del tipo di P hanno un autovalore uguale a 1 (siveda il Capitolo 9).Vediamo un esempio numerico con l’ausilio di MATLAB. Supponiamo che gli stati rap-presentino delle stanze di un labirinto. Un topolino inizialmente in una stanza si muovein questo labirinto in maniera casuale attraverso le altre stanze collegate tra loro. Sup-poniamo inoltre che se una stanza ha m collegamenti il topolino possa scegliere ognunadi esse con uguale probabilita 1/m. Una matrice che rappresenta una possibile dinamicadei cammini del topolino nel caso di 9 stanze e la seguente matrice di transizione P

P =

0BBBBBBBBBBB@

0 1/2 0 0 0 1/2 0 0 01/3 0 1/3 0 1/3 0 0 0 00 1/2 0 0 0 1/2 0 0 0

1/2 0 0 0 1/2 0 0 0 00 1/4 0 1/4 0 1/4 0 1/4 00 0 1/3 0 1/3 0 0 0 1/30 0 0 0 0 0 0 1 00 0 0 0 1/2 0 1/2 0 00 0 1/2 0 0 1/2 0 0 0

1CCCCCCCCCCCAin cui ogni elemento pij della matrice stessa rappresenta la probabilita di andare dallastanza i alla stanza j (pij = 0 significa che le stanze non sono collegate). Supponiamoche il topolino inizialmente si trovi nella prima stanza, ossia p0 = {1,0,0,0,0,0,0,0,0}.Calcoliamo il valore assoluto degli autovalori di P

1Nella realta la dimensione di questi sistemi e cosı elevata da rendere l’analisi impossibile, percui spesso si preferiscono modelli alternativi.

70 Probabilita

>> [V,D]=eig(P);>> abs(D)ans =

0.93841.00000.68550.74300.28640.28640.29540.00000.1499

A parte l’autovalore uguale ad 1, indichiamo con λ∗ il piu grande modulo dei restantiautovalori. Si ottiene λ∗ = 0. 9384 < 1. Ricordando il Capitolo 9, da PV = V D segueP = V DV −1 e quindi P n = V DnV −1. Dunque poiche λ∗ < 1 per n → ∞ si conservasolo l’elemento diagonale uguale ad 1 nella matrice D. In altre parole limn→∞Dn =diag{0,1,0,0,0,0,0,0,0}. Quindi se λ∗ < 1 la matrice di transizione converge ad unamatrice stazionaria P∞. In tal caso esiste un unico vettore di probabilita stazionariap∞ = P∞p0 tale che pn = P n p0 → p∞ per ogni vettore di probabilita degli stati inizialep0. Otteniamo

>> p0=[1 0 0 0 0 0 0 0 0];>> Dinf=diag([0 1 0 0 0 0 0 0 0]);>> Pinf=real(V*Dinf*inv(V));>> Pinf(1,:) % La matrice ha righe tutte ugualians =

0.0755 0.1572 0.1468 0.0461 0.1845 0.1887 0.0461 0.0922 0.0629>> pinf=Pinf*p0;>> pinf(1) % Il vettore ha componenti tutte ugualians =

0.0755

Essendo la matrice P∞ reale, l’estrazione della parte reale tramite real e di fatto super-flua. Tale comando migliora pero la visualizzazione dei risultati finali (si provi senza).Indipendentemente dallo stato iniziale la stanza piu probabile in cui trovare il topolino ela numero 6 con probabilita 0. 1887, mentre la probabilita che il topolino sia ancora nellaprima stanza e 0. 0755.Nella Figura 11.2 e visualizzato il comportamento delle componenti del vettore pn pern = 1, . . . ,50.

0 5 10 15 20 25 30 35 40 45 500

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figura 11.2 Comportamento delle componenti del vettore pn per n = 1, . . . ,50.

12Modelli probabilistici

12.1 Generazione di leggi discrete

Generiamo una matrice di v.a. con legge di Bernoulli si parametro 1/4

>> A =(rand(3,4)>.75)A =

1 0 0 00 1 1 10 1 0 0

La funzione rand(m,n) genera una matrice m× n di numeri pseudo-casuali indipendenticompresi tra 0 e 1. Un computer non puo produrre “caso” ma solo pseudo-caso. Infatti ilcalcolatore genera numeri con algoritmi numerici deterministici in grado di superare unaserie di test statistici che conferiscono a tali numeri una “apparente casualita”. Torneremosu questo aspetto nel laboratorio Matlab del prossimo capitolo. Il comando rand simulauna v.a. continua con risultati tra 0 e 1 e con densita 1 costante. Il risultato che si avradando piu volte il comando sara quindi diverso e indipendente dal precedente.L’operatore logico (A>B) restituisce 1 quando il predicato e vero, 0 se e falso. Nel nostrocaso la probabilita che esca 1 e la probabilita che una legge uniforme su [0,1] sia maggioredi 0. 75, ossia 1/4. In generale (rand(m,n)>(1-p)) restituisce quindi una matrice m× ndi v.a. di Bernoulli di parametro p. Se sommo le colonne della matrice, ottengo un vettoredi n colonne distribuite come una binomiale di parametri m e p. Nell’esempio

>> B=sum(A)B =

1 2 1 1

abbiamo che B e un vettore di 4 elementi, ciascuno dei quali ha legge binomiale Xi ∼B(3,1/4).Per generare un vettore X di 4 elementi, ciascuno dei quali ha legge di Poisson diparametro 3, utilizzeremo l’approssimazione binomiale X ≈ B(1000,3/1000)

>> X=sum( (rand(1000,4) > ( 1 - 3/1000 ) ) )X =

3 1 2 3

Digitiamo ora

>> X= floor(log(rand(1,1))/log(1-1/4))+1X =

6

72 Modelli probabilistici

dove abbiamo usato la funzione floor che restituisce la parte intera dell’argomento.Avremo che X = k se e solo se

k − 1 ≤ log(rand(1,1))

log(1− 1/4)< k

ossia se e solo se (si ricordi che 3/4 < 1),

(3

4)

k−1

≥ rand(1,1) > (3

4)

k

.

Poiche rand(1,1) genera una v.a. uniforme U , la probabilita con cui X = k accade e

P ({X = k}) = P({(3

4)

k−1

≥ U > (3

4)

k}) = (3

4)

k−1

− (3

4)

k

=1

4(3

4)

k−1

ossia la legge di X e una geometrica di parametro 1/4. Per generare G vettore 1 × 4 div.a. indipendenti con legge geometrica di parametro 1/25

>> X = floor(log(rand(1,4))/log(1-1/25))+1X =

8 29 5 14

12.2 Generazione di leggi continue

Digitiamo ora

>> X = -log(rand(1,1))/2X =

0.2558

Avremo che X > t se e solo se

− log(rand(1,1))

2> t

ossia se,rand(1,1) < e−2t.

La probabilita con cui X > t accade e quindi e−2t ossia la legge di X e una esponenziale diparametro 2. Per generare E vettore 1×4 di v.a. indipendenti esponenziali di parametro1/25

>> E = -log(rand(1,4))/(1/25)E =

31.7329 16.0619 49.2324 14.3825

Per generare una v.a. di legge gaussiana possiamo operare in due modi:

1. utilizzare la funzione randn(m,n) che genera una matrice m× n di numeri pseudo-casuali indipendenti distribuiti come una normale standard (con parametri 0 e 1).Avremo dunque che

>> A =randn(1,4)A =

-0.5883 2.1832 -0.1364 0.1139

e un vettore di 4 elementi distribuiti come una normale standard. A questo punto,per generare un vettore N di 4 elementi distribuiti come una normale di media 3 evarianza 5

Modelli probabilistici 73

>> N = 3 + sqrt(5)*AN =

1.6845 7.8818 2.6950 3.2548

2. utilizzare il Teorema Limite Centrale: M(s)n ≈ N(0,1). La standardizzata della

media di n leggi bernoulliane di parametro 1/2 e

M (s)n =

X1+···+Xnn

− 1/2p1/(4n)

=2√n

(X1 + · · ·+ Xn − n/2).

Utilizziamo n = 1000. Avremo dunque che

>> X = 2 / sqrt(1000) * ( sum((rand(1000,1) > 1/2 )) - 1000/2 )X =

0.8854

e distribuito come una normale standard (per generare un vettore di m elemen-ti, avremmo scritto rand(1000,m) invece di rand(1000,1)). A questo punto, pergenerare N distribuito come una normale di media 3 e varianza 5

>> N = 3 + sqrt(5)*XN =

4.9799

12.3 La funzione Φ

MATLAB calcola l’approssimazione numerica di erf(x) =R x

02/√

π e−s2ds con la funzione

erf(x). Attraverso questa funzione vogliamo calcolare Φ(x) = P ({X ≤ x}) (la funzionedelle tavole di una distribuzione normale standard). Abbiamo che

Φ(x) =R x

−∞ fX(t) dt =R 0

−∞ fX(t) dt +R x

0fX(t) dt

=1

2+

Z x

0

1√2π

e− t2

2 dt(t=s/

√2)

=1

2+

Z x/√

2

0

1√π

e−s2ds

=1

2+

1

2

Z x/√

2

0

2√π

e−s2ds =

1

2+

1

2erf

�x√2

�Possiamo quindi far disegnare il grafico la funzione Φ(x) (si veda la Figura 12.1) o calcolareP ({X ≤ 1. 25}), scrivendo

>> x= [-4:.01:4]; %ascisse per il grafico>> y= 1/2 + 1/2 * erf( x/sqrt(2) ); %ordinate per il grafico>> plot(x,y) %disegna grafico>>>> 1/2 + 1/2 * erf( -1.25/sqrt(2) ) %calcola P({ X <= 1.25 })ans =

0.1056

74 Modelli probabilistici

−4 −2 0 2 40

0.2

0.4

0.6

0.8

1

Figura 12.1 Grafico della funzione Φ.

13Campionamento statistico e descrizione dei

dati

13.1 Indici descrittivi e istogrammi

Carichiamo in un vettore A i dati dell’Esempio 13.7 e calcoliamo qualche indice descrittivotramite i comandi mean, var, std e median.

>> A = [1.67 1.61 1.84 1.73 1.6 1.79 1.8 1.84 1.86 ...1.71 1.77 1.8 1.9 1.94 1.77 1.85 1.84 1.88 1.7 1.91];>> mean(A) %mediaans =

1.7905>> var(A) %varianzaans =

0.0092>> std(A) %dev. standard ( =sqrt(var(A)) )ans =

0.0960>> median(A) %medianaans =

1.8000

La funzione MATLAB hist(x,N) consente di realizzare l’istogramma di frequenza di xsuddividendo l’intervallo [min(X), max(X)] in N sottointervalli (barre) di uguale ampiez-za e calcolando quanti valori di x cadono in un dato intervallo.

>> hist(A,8)

MATLAB ha una sua formula per calcolare gli estremi delle barre. Per selezionare gliestremi delle barre, occorre utilizzare la funzione histc(X,Y), dove X e il vettore dei datie Y e il vettore con gli estremi delle barre. La funzione histc non disegna il grafico, marestituisce l’altezza delle barre i cui estremi sono nel vettore Y .

>> Y = [1.60:0.05:2.00];>> alt= histc(A,Y)alt =

2 2 2 3 5 3 3 0 0>> bar(Y,alt,’histc’)

In quest’ultima riga si e usata la funzione bar (che disegna un bar-chart) con l’opzione’hist’ (che indica che in Y vi sono gli estremi delle barre e non i punti medi).Costruiamo l’istogramma di densita della Tabella 13.5. Prima di proseguire, notiamo che

76 Campionamento statistico e descrizione dei dati

1.6 1.7 1.8 1.9 20

1

2

3

4

5

1.6 1.7 1.8 1.9 2

1

2

3

4

5

Figura 13.1 Istogrammi dei dati dell’Esempio 13.7, ottenuti con barre assegnate daMATLAB e imponendo gli estremi delle barre.

hist restituisce i dati per il grafico con una barra in piu dell’ampiezza pari all’ultimaassegnata.1

>> Y=[1.60 1.72 1.81 1.87 1.95 2];>> num_oss= histc(A,Y);>> ser_stat= num_oss / 20;>> G=diag(ones(5,1),1)-diag(ones(6,1));>> G(6,:)=G(5,:);>> amp_int= (G*Y’)’;>> alt_bar= (ser_stat./amp_int)alt_bar =

2.0833 3.3333 4.1667 2.5000 0 0>> bar(Y,alt_bar,’histc’)>> sum(alt_bar.*amp_int) %controllo area=1ans =

1

13.2 Retta dei minimi quadrati

Inseriamo i dati dell’Esempio 13.8 e calcoliamo MXn , MY

n , SX e SY .

>> Eta=[ 36; 48; 51; 54; 57; 60];>> Alt=[ 86; 90; 91; 93; 94; 95];>> MnX= mean(Eta)MnX =

51>> MnY= mean(Alt)MnY =

91.5000>> S_X= std(Eta)S_X =

8.4853>> S_Y= std(Alt)S_Y =

3.2711

1Questo perche MATLAB assegna i valori all’intervallo di sinistra e quindi si riserva una barraper i valori uguali al massimo valore che abbiamo inserito.

Campionamento statistico e descrizione dei dati 77

>> M=corrcoef([ Eta Alt ])M =

1.0000 0.99440.9944 1.0000

La matrice M di correlazione va letta nel seguente modo

M =

8>>: r(X,X) r(X,Y )

r(Y,X) r(Y,Y )

9>>;Calcoliamo quindi a e b per la retta dei minimi quadrati. Abbiamo che

>> a=M(1,2)*S_Y/S_Xa =

0.3833>> b= MnY - a * MnXb =

71.9500

e quindi disegniamo la Figura 13.3.

>> plot(Eta,Alt,’.’) %scatterplot>> hold on>> xx=[35 65]; %per una retta bastano due punti...>> plot(xx, a*xx + b) %disegna la retta dei minimi quadrati

13.3 Distribuzione statistica

Vediamo ora una simulazione per mostrare il Teorema 13.4. Simuliamo quindi molti dati(n = 10000) di una distribuzione gaussiana standard e vediamo se F ∗(2) approssimabene P ({X ≤ 2}) = Φ(2) ≈ 0. 9772.

>> X=randn(10000,1);>> F_2= sum( ( X <=2 ) )/10000F_2 =

0.9763

L’operatore logico ( X <=2 ) restituisce 1 se gli elementi sono minori o uguali a 2, 0altrimenti. La loro somma sum conta quindi quanti elementi nelle colonne di X sonominori o uguali a 2. Ripetendo l’esperimento con 4 simulazioni.

>> X=randn(10000,4);>> F_2= sum( ( X <=2 ) )/10000F_2 =

0.9787 0.9766 0.9765 0.9755>> Phi_2=1/2 + 1/2 * erf( 2/sqrt(2) ); %Phi(2) (vedere Capitolo 12)>> F_2/Phi_2 %rapporto valore stimato/valore veroans =

1.0015 0.9993 0.9992 0.9982

13.4 Generatori di numeri casuali

Il problema della generazione di sequenze di numeri pseudo-casuali distribuiti uniforme-mente al calcolatore e un problema estremamente complesso e affascinante. Ricordiamoche nella distribuzione uniforme ogni possibile numero in un determinato intervallo eugualmente probabile (ad esempio il lancio di un dado genera una sequenza di numeri

78 Campionamento statistico e descrizione dei dati

casuali distribuiti uniformemente in [1,6]). Ci limiteremo in questo breve paragrafo adaccennare alcune idee di base. Supponiamo di volere generare un numero casuale di 4cifre, ossia un numero tra 0000 e 9999. L’dea del metodo middle-square, ideato da Johnvon Neumann nel 1946, e la seguente. Si parte da un valore iniziale, detto seme a c cifre,per esempio 1234, lo si eleva al quadrato ottenendo un numero con il doppio 2c di cifre,nel nostro caso 01522756, di queste teniamo solamente le c cifre di mezzo, ossia 5227.Ripetendo il procedimento otteniamo 27321529 e quindi 3215 e cosı via. In MATLABper esempio

>> n=10;>> seme=1234;>> x=zeros(1,n); x(1)=seme;>> for i=1:n-1, st=num2str(x(i)^2,’%08d’); x(i+1)=str2double(st(3:6)); end>> xx =

Columns 1 through 71234 5227 3215 3362 3030 1809 2724

Columns 8 through 104201 6484 422

Ogni nuovo numero nella successione e determinato univocamente dal suo predecessore.La successione generata quindi non potra essere casuale ma avra solo il carattere di ap-parente casualita. Dal punto di vista MATLAB abbiamo usato la funzione num2str checonsente di convertire un numero in una stringa, e la funzione str2double che esegue ilprocedimento inverso. Ogni successione di numeri generati da questo algoritmo inizieraa ripetersi prima o poi. Il numero di numeri della sequenza prima che intervenga unaripetizione e detto periodo della sequenza. La lunghezza di tale periodo puo essere consi-derata una misura della bonta del generatore di numeri pseudo-casuali. Sfortunatamentequesto metodo puo degenerare in sequenze con periodi molto brevi, ad esempio a partiredal valore 0 la sequenza ha sempre periodo 1, oppure partendo con 43 e numeri a 2 cifreotteniamo la sequenza 43, 84, 05, 02, 00, 00, . . ..Nel 1948 per la prima volta venne proposto un generatore di numeri casuali distribuitiuniformemente detto generatore lineare congruenziale o piu brevemente LCG. Il metodoLCG, analogamente al metodo middle-square, ha bisogno di un seme per inizializzare lasequenza di numeri secondo la seguente regola deterministica

xn+1 = (axn + c)mod m, n ≥ 0,

dove a, c ed m sono opportuni numeri interi costanti. La notazione “mod m”, ossiamodulo m, significa che axn + c viene diviso per m e xn+1 posto uguale al resto interodella divisione. Quindi xn+1 assume valori interi tra 0,1,2, . . . ,m− 1.Il problema della scelta dei migliori valori per a, c ed m e il punto cruciale del metodo.Cosa si intende pero per migliori? Un aspetto importante e la lunghezza del periodo dallaquale si vede subito che m dovra essere molto grande. Un secondo aspetto consiste nelgarantirsi che per un dato m i valori di a e c siano tali che la successione abbia periodomassimo. Una delle scelte piu popolari2 e

m = 231 − 1, a = 75, c = 0.

Da un punto di vista pratico i generatori di numeri casuali restituiscono invece di xn+1

il suo valore diviso per m, questo al fine di evitare numeri troppo grandi. Per esempio ilmetodo LCG in MATLAB potra essere realizzato come

>> n=10000;>> m=2147483647; a=16807; c=0;>> seme = 1234;>> x=zeros(1,n); x(1)=seme;>> for i=1:n-1, x(i+1)=mod(a*x(i)+c,m); end>> x=x./m;

2Il metodo LCG con questa scelta era alla base dell’algoritmo utilizzato da MATLAB per lagenerazione di numeri casuali distribuiti uniformemente fino alla versione 4.0.

Campionamento statistico e descrizione dei dati 79

dove abbiamo utilizzato la funzione MATLAB mod dall’ovvio significato.Il valore massimo generato sara quindi 0,99999999953434 ed il minimo 0,00000000046566.Dopo avere generato una tabella di 100 valori pseudo-casuali compresi tra 0 ed 1 possiamoverificarne l’uniforme distribuzione tramite

>> N=50;>> hist(x,N);

Nel caso della distribuzione uniforme tale valore teorico dovrebbe essere esattamentepari a n/N ossia 10000/50 = 200 nell’esempio specifico. Nella figura 13.2 confrontiamo ilmetodo LCG con il generatore rand di MATLAB.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

50

100

150

200

250LCG

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

50

100

150

200

250rand

Figura 13.2 Istogramma di frequenza dei generatori uniformi di numeri casuali.

Una richiesta piu importante nel valutare la bonta di un generatore uniforme di numeripseudo-casuali e l’assenza di correlazione tra i numeri generati dall’algoritmo. In altreparole non deve emergere nessuna relazione funzionale tra xn e xn+1 per n > 0. Questaproprieta puo essere verificata graficamente realizzando il grafico bidimensionale dei punti(xn,xn+j) per j > 0. In tale grafico non dovranno comparire linee, forme o altre struttureregolari. In figura 13.3 riportiamo il risultato per j = 1 con 1000 punti ottenuto con ilgeneratore LCG con le scelte “ottimali” di a, b e c e modificato scegliendo come valorim = 31, a = 13 e c = 0.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figura 13.3 Assenza di correlazione (sinistra) e forte correlazione (destra) nel graficodei punti (xn,xn+1) in un generatore uniforme di numeri pseudo-casuali.

14Metodi statistici inferenziali

14.1 Intervalli di confidenza

Costruiamo un intervallo di confidenza per la media e per la varianza di un numero moltoalto (n = 1000) di simulazioni di una X ∼ U(0,1). Le nostre simulazioni vengono raccoltenel vettore A.

>> A=rand(1000,1);

Avremo, per esempio, che le osservazioni tra la 45-esima e la 50-esima saranno

>> A(45:50)’ans =

0.7266 0.8288 0.3637 0.0351 0.6318 0.4423

Disegniamo l’istogramma con i dati che abbiamo

>> B=histc(A,[0:1/10:1]);>> bar([0:1/10:1]+1/20,B,’histc’)

ottenendo, a meno della casualita, un istogramma simile a quello nella Figura 14.1.

0 0.2 0.4 0.6 0.8 10

20

40

60

80

100

120

140

Figura 14.1 Istogramma delle mille simulazioni da distribuzione uniforme.

Costruiamo un intervallo di confidenza con affidabilita 95% per il valore atteso delladistribuzione simulata, che sappiamo essere 1/2.

82 Metodi statistici inferenziali

>> MA=mean(A);>> SA=std(A);>> ES=MA-1.96*SA/sqrt(1000) %estremo sinistro intervalloES =

0.4942>> ED=MA+1.96*SA/sqrt(1000) %estremo destro intervalloED =

0.5295

Ora, costruiamo 20 simulazioni di intervalli di confidenza e confrontiamoli sul grafico colvalore vero del parametro, ossia 1/2, come mostrato nella Figura 14.2.

>> A=rand(1000,20);>> MA=mean(A);>> SA=std(A);>> ES=MA-1.96*SA/sqrt(1000);>> ED=MA+1.96*SA/sqrt(1000);>> alt=[1:20];>> figure>> hold on>> for i=1:20, plot([ES(i) ED(i)],[i i]), end>> plot([1/2 1/2],[0 21],’--’), hold off

0.45 0.5 0.550

5

10

15

20

Figura 14.2 Venti intervalli di confidenza a confronto col parametro vero.

Il valore vero non appartiene a tutti gli intervalli, infatti dalla definizione di intervallo diconfidenza

P (ϑ ∈ [a,b]) = βsegue che prima di effettuare le simulazioni, avevamo

P (1/2 ∈ [a,b]) = 0.95.

Qual e la probabilita che tutte e 20 le simulazioni contengano il valore vero? Dalla leggebinomiale, se Y ∼ B(20,0.95), P ({Y = 20}) ≈ 0.36.

14.2 Test d’ipotesi

Utilizziamo gli stessi dati per un test di adattamento alla legge uniforme, ossia controllia-mo se il simulatore di MATLAB e “uniforme”, come dichiara. Dividiamo l’intervallo [0,1]

Metodi statistici inferenziali 83

in 10 classi, ognuna di ampiezza 1/10. Avremo quindi N = 10, p0k = 1/10 (distribuzione

uniforme), n = 1000. Per ogni intervallo, mi aspetto 100 simulazioni, come si verificasemplicemente dalla relazione n0

k = n ∗ p0k.

>> VO=histc(A,[0:1/10:1]); %numero di osservazioni in ogni intervallo>> VO=VO(1:10,:);>> VA=100*ones(10,20);

Il nostro test, a priori, rigettera l’ipotesi di adattamento se D > d1−α, con D ∼ χ2N−1.

Dalle tavole, d0.99 ∼ 21. 666. Impostiamo quindi i 20 test (abbiamo in memoria 20simulazioni).

>> D=sum((VO-VA).^2./VA) %distanza pivotale per i test Chi QuadroD =

Columns 1 through 72.4600 8.4200 2.3200 4.8600 2.3200 4.9600 6.6000

Columns 8 through 143.7400 6.0400 13.5000 10.7000 13.2600 7.9800 11.2600

Columns 15 through 205.0800 10.1400 3.7800 5.0200 5.3000 2.8000

Per questa simulazione, in nessun caso possiamo rifiutare l’ipotesi di adattamento, comesi puo piu facilmente osservare con l’operatore logico >

>> (D > 21.666)ans =

Columns 1 through 120 0 0 0 0 0 0 0 0 0 0 0

Columns 13 through 200 0 0 0 0 0 0 0

Quest’ultimo fatto dipende chiaramente dall’esperimento eseguito. Va notato comunqueche su 20 test a livello α = 0.99, la probabilita che nessuno dia come risultato di rigettarel’ipotesi nulla e circa 0.82.Eseguiamo ora un test χ2 di indipendenza con i dati dell’Esercizio 14.11.

>> A=[30 40; 58 67; 37 18]A =

30 4058 6737 18

Ricordando che sum somma i dati per colonne, calcoliamo SR (somme righe) e SC (sommacolonne) per mezzo dell’operatore ’ (trasposta di una matrice).

>> SR=sum(A’)’, SC=sum(A),SR =

7012555

SC =125 125

Calcoliamo quindi NT il numero totale di osservazioni e quindi la matrice del numero attesodi osservazioni stimato sotto l’ipotesi nulla. Il calcolo diviene in MATLAB algebricamentesemplice, moltiplicando il vettore colonna SR per il vettore riga SC e divendo il totale perNT.

NT =250

>> VA=SR*SC/NTVA =

84 Metodi statistici inferenziali

35.0000 35.000062.5000 62.500027.5000 27.5000

Il nostro test di indipendenza, a priori, rigettera l’ipotesi di adattamento se D > d1−α,con D ∼ χ2

(r−1)(c−1). Scegliamo α = 0. 95 Dalle tavole, d0.95 ∼ 5. 991.

>> D=sum(sum((A-VA).^2./VA))D =

8.6402

Rifiutiamo l’ipotesi di indipendenza, come si puo piu facilmente osservare con l’operatorelogico >

>> (D > 5.991)ans =

1

Indice Analitico

catene di Markov, 69

frequenze naturali, 56

istogramma di frequenza, 75

leggidi Keplero, 25

MATLABbackslash, 49comandi

abs, 5acos, 5asin, 5atan, 5axis, 12, 14, 63bar, 65, 75cd, 15clear, 7colormap, 47compose, 22conv, 20copyfile, 15cos, 4, 5deconv, 20delete, 15det, 51diag, 9diff, 27dir, 15dsolve, 42eig, 55erf, 73expand, 15exp, 5eye, 9ezplot, 11, 20figure, 62find, 65finverse, 22floor, 72fminbnd, 28fmin, 28for...end, 24format, 4fplot, 20

fzero, 25ginput, 61grid, 14histc, 75hist, 75hold off, 13hold on, 13image, 47inline, 20inv, 52legend, 13, 14length, 6limit, 23linspace, 7, 9log10, 5log2, 5logspace, 9log, 5mean, 75median, 75meshgrid, 68mesh, 63mkdir, 15mod, 79norm, 57null, 50num2str, 78ode45, 43ones, 9plot3, 38plot, 12polyval, 19pwd, 15quad8, 37quad, 37quiver3, 63quiver, 63randn, 72rand, 9, 71rank, 50real, 70root, 32sign, 5simplify, 32sin, 4, 5solve, 30, 32sqrt, 5std, 75

86 Indice Analitico

str2double, 78subplot, 27subs, 22, 51surf, 63symsum, 41syms, 15tan, 5taylor, 41text, 14title, 14, 63trapz, 39type, 15var, 66, 75whos, 6who, 6xlabel, 14, 63ylabel, 14, 63zeros, 9zlabel, 63

funzioni vettoriali, 9M-file, 14notazione due punti, 7numeri complessi, 3operatori\, 49*, 4+, 4-, 4.*, 4./, 4.ˆ, 4/, 4ˆ, 4

operazioni puntuali, 10Symbolic Mathematics Toolbox, 15variabili

Inf, 3NaN, 3ans, 2, 3computer, 3i,j, 3pi, 3realmax, 3realmin, 3version, 3

matricedi stiffness, 56di transizione, 69di Vandermonde, 52

probabilitadi transizione, 69

risonanza, 56

Laboratorio Matlab - €¦ · MATLAB pu`o essere utilizzato in modo diretto per calcolare semplici...

Documents

Transcript of Laboratorio Matlab - €¦ · MATLAB pu`o essere utilizzato in modo diretto per calcolare semplici...