Nuovo Blog

L’ultimo articolo pubblicato su questo blog, risale oltre ad un anno fa.

Gia da un po’ di tempo pensavo di riprendere in mano il blog, per motivazioni varie però l’hosting su wordpress mi andava stretto, molte funzionalità di wordpress stesso sono bloccate (se non ne sapevate nulla leggete qua)), ma nonostante i diversi vantaggi (non doversi preoccupare della gestione del sito, gli aggiornamenti e quant’altro), ho deciso di provare a ricominciare da zero.

Al momento ho deciso di lasciare anche online questo blog, quotidianamente ci accedono ancora molti visitatori, piano a piano aggiornerò gli articoli scritti e li posterò sul nuovo blog, finchè non deciderò di cancellare il vecchio sito.

Il sito è appena stato installato, passerà un po’ di tempo prima che sarà completamente a posto, ma tanto non ho fretta .

Questo è il link per il nuovo blog

Fekir

Interpolazione polinomiale

Cominciamo a spiegare cosa è l’interpolazione polinomiale. Immaginate di avere una curva analitica da studiare y=”qualcosa di complicato”, oppure una serie di dati in coppie (x1,y1), (x2,y2)…. con tutti gli xi diversi.
Come potete ben immaginare non sarà sempre un oggetto semplice da studiare.

Un tipo di curve particolarmente semplici da studiare sono i polinomi, e quindi ci piacerebbe trovare dei polinomi che assomiglino alle funzioni che vogliamo studiare, in modo da avere un lavoro più semplice da svolgere.
Per chi conosce relativamente bene la matematica avrà magari pensato ai polinomi di Taylor o al caso particolare McLaurin, ma noi vogliamo prendere polinomi che siano anche molto semplici da costruire. Fare le derivate non è sempre una operazione semplice e poi questi due polinomi sono adatti per fare una analisi locale in un punto, non sempre approssimano una funzione su un intervallo.

Il metodo è utilizziamo è il seguente (supponendo di dover sudiare una funzione in una variabile o qualcosa di analogo)
Scegliamo sull’intervallo che vogliamo studiare n+1 punti. Fissati questi punti sull’ascissa delle x li valutiamo e otteniamo la coppie (x1,y1),(x2,y2),…

A questo punto si può dimostrare che esiste un solo polinomio di grado n (ad esempio per due punti distinti passa una sola retta, per tre punti distinti una sola parabola…). Il problema (ma neanche tanto complicato) è come costruirlo.

Il primo metodo è quello di Vandermonde. Si considera un generico polinomio di grado n e si impongono tutte le condizioni. Dovremo quindi poi risolvere un sistema lineare, che coincide con una matrica, guarda a caso, chiamata matrice di Vandermonde…

Il secondo metodo, forse più “diretto” e che uso più spesso (ma siete liberi di usare il metodo che volete) è quello di Lagrange.
Si costruiscono n polinomi di grado n tali che $latexL_i(x_j)=0$ per ogni i diverso da j, mentre . Il polinomio interpolante sarà quindi semplicemente . Provare per credere se inserite nel polinomio il valore otterrete il valore .

Facciamo un esempio, dati i nostri punti da interpolare:
(1,2), (2,3), (3,4)
Costruiamo

E quindi abbiamo il nostro polinomio interpolante!

Il terzo metodo, concettualmente più complicato ma da un punto di vista operativo molto più conveniente.
Infatti se adesso volessimo il punto (4,5) ai nostri punti precedenti, dovremmo rifare tutti i conti per ottenere il polinomio interpolante. Con il metodo delle differenze divise i conti si riducono drasticamente se vogliamo aggiungere dei punti.

Sono un poco più complicate da spiegare e mi porterebbero via tanto spazio, rischiando di rendere poco leggibile il tutto, pertanto vi rimando a Wikipedia, purtroppo l’articolo in italiano non esiste, e purtroppo la spiegazione in inglese non è (a mio parere) molto buona, trovo quella in tedesco più comprensibile, qui l’algoritmo viene anche chiamato di Neville-Aitken ma temo che non tutti conosciate abbastanza questa lingua.

In ogni caso, il metodo di Lagrange sarà quello che useremo per interpolare un polinomio, vediamo un esempio!

Ho scritto le function lagr.m per il metodo di lagrange, e per chi volesse anche diffdiv_coeff.m, che servono per creare i polinomi interpolanti.
Vogliamo, ad esempio, interpolare con un polinomio di grado 5 la funzione sin(2*pi*x) nell’intervallo [-1,1].
Creiamo quindi uno script-file (considerate che da ora in poi qualsiasi esercizio lo dovete fare su uno script-file, per poterlo correggere, rileggere, salvare…) e cominciamo a scrivere:

f=inline('sin(2*pi*x)');
xnodi=linspace(-1,1,5+1); %ci servono n+1 punti per avere un polinomio di grado n!
ynodi=f(xnodi); %qui stiamo valutando la nostra funzione nei nodi d'interpolazione
x=linspace(-1,1);
p_lag = lagr(xnodi,ynodi,x);
y=f(x);
plot(x,y,x,p_lag)

Dopo aver eseguito lo script-file verrà visualizzato un grafico contenente le due funzioni. È auspicabile pensare che all’aumentare del grado n del polinomio aumenti la somiglianza tra le due funzioni. Provate infatti ad aumentare il numero di nodi a 11.

La funzione lagr.m si occupa di valutare il polinomio interpolanti sull’asse delle x, in sostanza quello che ci restituisce è un vettore con la valutazioni del polinomio interpolante.

Purtroppo però non è vero che prendendo nodi equispaziati (come li abbiamo scelti nell’esempio precedente) la funzione interpolante aumenti il grado di precisione all’aumentare del grado del polinomio stesso. Esiste infatti un famoso controesempio, il controesempio di Runge. Per visualizzarlo potete scaricare il file runge.m e provare ad eseguirlo con Octave.

Possiamo notare che al crescere di n le due funzioni non sono poi tanto simili e si può dimostrare che queste differenze diventano sempre più evidenti.

Risolto quindi il problema di come studiare una funzione generica attraverso un polinomio ci dobbiamo porre il problema di come scegliere i nodi d’interpolazione, in modo che al crescere del loro numero il polinomio interpolante assomigli alla funzione da studiare. Dal controesempio di Runge appare ovvio che scegliendo più nodi ai bordi si costringerebbe il polinomio interpolante a stare più vicino alla funzione da approssimare.

Si può dimostrare, che scegliendo i nodi di Tschebyscheff (valori che sono soluzione di un particolare polinomio), la funzione interpolante tenderà sempre al polinomio da studiare. Per fortuna questi nodi non sono difficili da calcolare, se stiamo lavorando nell’intervallo [-1,1] questi sono tutti nella forma: , dove n è il grado di interpolazione e j assume i valori da 1 a n.

Se lavoriamo su un generico intervallo [a,b] basterà effettuare la trasposizione e valutare il problema con i nodi .
Proviamo quindi a eseguire runge2.m. Avremo però bisogno di mettere tra le nostre function nodicheb.m che si occupa di creare i giusti nodi.
Possiamo osservare che all’aumentare dei nodi il polinomio interpolante diventa sempre più simile alla nostra funzione.

Ovviamente non è strettamente necessario usare le mie function per interpolare un polinomio, le ho messe a disposizione soltanto per scopo didattico. Octave integra già questa possibilità attraverso il comando polyfit(xnodi, ynodi, n), dove n è il grado del polinomio che vogliamo ottenere e xnodi, ynodi sono due vettori di lunghezza n+1.

A questo punto potete passare agli esercizi. Ovviamente potete interpolare anche altri tipi di funzioni!

Funzioni matematiche e grafici

Ecco un nuovo capitolo a questa guida su Octave

Grafici (2D)
In Octave è possibile creare grafici in modo molto semplice utilizzando il comando plot. Un esempio è sicuramente più che esplicativo:

x=linspace(0,10); %genera valori su cui valutare la funzione, vettore delle ascisse
y=x.*cos(x)/2; %valutiamo la funzione nei punti definiti prima, vettore delle ordinate
plot(x,y);

Ovviamente è possibile passare diversi paramatri con il comando plot, se volevamo una linea rossa è tratteggiata avremmo dovuto dare il comando

plot(x,y,'r--') %per altre opzioni date il comando
doc plot

(così imparate un poco a usare il doc, che può davvero salvare in molte situazioni).
Vale la pena menzionare però un altro paio di comandi, che si possono dare subito dopo aver disegnato la funzione.
Il comando “grid”, dato dopo il comando plot(…) aggiunge le linee di griglia, mentre il comando “hold on” permette di disegnare più grafici sullo stesso foglio. Per smettere di disegnare grafici sullo stesso foglio basta dare “hold off”. Oppure, sempre per disegnare più grafici sullo stesso foglio, li si può passare entrambi sotto lo stesso plot(…), ecco un paio di esempi:

x=linspace(0,10);
y1=x.*cos(x);
y2=sin(x);
plot(x,y1,'r' , x,y2,'b');

oppure:

x=linspace(0,10);
y1=x.*cos(x);
plot(x,y1,'r');
hold on
y2=sin(x)
plot(x,y2,'b')
hold off

Funzioni
Abbiamo già visto per il disegnare i grafici, come possiamo creare delle funzioni da disegnare. In sostanza abbiamo creato dei vettori con i valori valutati nei punti sul quale era definita la funzione.
Se però, oltre a disegnarla con una funzione dobbiamo anche lavorarci, sapere per quali punti passa potrebbe non bastare.
Attraverso il comando inline è possibile definire una generica funzione, ad esempio:

f=inline('sin(x).*x', 'x')

Se diamo il comando “whos” potremo infatti vedere che a f non è associato un numero, bensì una funzione.
Se adesso, come prima vogliamo disegnare la nostra funzione, possiamo operare nel seguente modo:

x=linspace(-5,5);
y=f(x);
plot(x,y)

Grafici (3D)
Quando ci ritroviamo a disegnare un grafico in 3 dimensioni dobbiamo distinguere se stiamo per andare a disegnare una curva, oppure una superficie, due oggetti che vengono trattati in modo distinto.

Se vogliamo disegnare una curva, possiamo adottare la seguente sintassi

t=linspace(0,20*pi,1000);
plot3(cos(t).*t, sin(t).*t,t);

Se invece vogliamo disegnare una superficie abbiamo bisogno di due variabili una x e una y. I grafici possiamo disegnarli come nel se

guente esempio con la funzione mesh

x=y=linspace(-2*pi,2*pi,40);
[X,Y]=meshgrid(x,y);
Z=sin(X).*sin(Y);
mesh(X,Y,Z)

oppure, al posto di mesh(X,Y,Z) possiamo dare il comando surf(X,Y,Z)

Polinomi
Abbiamo visto come definire una funzione generica attraverso il comando inline, nel caso di polinomi è però possibile considerare un’altra sintassi.
Supponiamo di avere un polinomio di grado n (cioè con il primo coefficente diverso da 0)

allora possiamo rappresentarlo attraverso il vettore dei coefficienti di dimensione n+1.
Se vogliamo valutare il polinomio in un punto, possiamo fare così:

p=[2 3 4 0]
x=2;
w=polyval(p,x)

Se vogliamo invece disegnare un grafico, allora possiamo valutarlo in un vettore di punti, ad esempio:

p=[-3 0 -1 2];
x=linspace(-5,5);
plot(x, polyval(p,x));

Mentre per la somma di due polinomi basta sommare due vettori (osservando che devono avere la stessa lunghezza, si può usare ad esempio la function polyadd), per il prodotto esiste una function apposita:

p=[-3 0 -1 2];
q=[1 -3];
conv(p,q)

Per calcolare la funzione derivata e la funzione integrale esistono i comandi appositi polyint e polyder.

Come al solito, ecco alcuni esercizi e le relative soluzioni.

Script file e function in Octave

Comincia qui la seconda lezione di Octave
Script-file
Normalmente è possibile navigare nella cronologia dei comandi effettuati premendo il tasto su o giu, risulta però molto comodo salvare su dei file tutte le istruzioni date, con tanto di commenti, in modo da poter ritirare fuori, a distanza di giorni, mesi o anni, il proprio lavoro.

I file che introduciamo ora si chiamano script-file, per saranno file con estensione .m, e il loro contenuto è una sequenza di istruzioni, scritte come se fossero digitate direttamente in Octave.
Per eseguirli, da Octave, se ci si trova nella stessa cartella dello script-file, basta digitare il nome del file (senza l’estensione).

Notate che normalmente si deve ad uno script-file un nome diverso dalle funzioni già presenti in Octave, altrimenti Octave utilizzerà la funzione al posto del vostro script-file.
Per controllare se un nome esiste già potete dare il comando

exist('nome_file')

Se il comando esiste vi verrà riportato un numero diverso da 0.
Il seguente script-file calcola, ad esempio, la somma dei primi n numeri:
somma.m

%Questo script-file (somma.m) calcola la somma da 1 a n.

n=input('inserisci un numero intero positivo')
v=[1:n];
s=sum(v);
disp('somma= ')
disp(s)

I commenti sono normalmente preceduti da %, oppure, per chi conosce già il C o un minimo di bash, #. Tuttavia in Matlab non è possibile, quindi in questi tutorial, userò la sintassi compatibile con Matlab, appunto per mostrare quanto siano intercambiabili (almeno fino ad un certo livello) i due programmi. Come potete infatti notare, nella documentazione di Octave i commenti vengono dati antecedendo il simbolo ‘#’, mentre nella documentazione di Matlab si usa il simbolo ‘%’ che può appunto venir anche usato in Octave.

Function
Le function sono, normalmente, degli “script-file” che accettano dei parametri e danno un risultato. Normalmente vengono richiamate dagli script-file, inoltre la loro sintassi è leggermente diversa, in modo da essere facilmente interpretati da altri file.
L’intestazione deve avere la struttura

Il file deve avere il nome “nomefunction.m” e un file può contenere una e una sola function.
Per interrompere la function è consigliato usare la parola return, e i valori output1,output2,…outputn devono essere stati tutti assegnati.
Una function viene richiamata semplicemente nel seguente modo:

[value1, value2, ...valuen]=nomefunction(input1, input2, ...inputn)

A titolo di esempio scriviamo un function che, dato un vettore, riporta il minimo e massimo e elemento

function [min, max]=minmax(v)
%Questa function riporta il minimo e massimo elemento di un vettore
min=min(v);
max=max(v);
return

Come già scritto più volte in precedenza, uso questa sintassi per il Matlab, mentre per Octave sarebbero possibile anche altre soluzioni più specifiche che non presento, appunto per mantenermi compatibili con ambedue i programmi. Qui trovate ad esempio la documentazione per Octave, non sono però sicuro che Matlab riconosca il comando endfunction (purtroppo non possiedo Matlab e quindi non ho potuto controllare), per evitare questo problema vi do’ il link alla documentazione di Matlab.

In uno dei primi articoli ho scritto che avremmo dovuto creare una cartella dove avremmo posizionato le nostre function. Se non lo avete ancora fatte fatelo, e salvate la function minmax.m (assicuratevi di salvarla con il nome giusto!)

Cicli in Octave
Come in tutti i linguaggi di programmazione anche in Octave è possibile definire dei cicli. In sostanza si dice a Octave di ripetere un certo calcolo (in cui ovviamente i numeri dati varieranno) finché non accade una certa condizione preimpostata.

Cicli while
La sintassi è la seguente
while (condizione==true)
istruzione
end

un esempio sciocco

>n=0;
>while (n<5) >n=n+1
>end

Questo ciclo aumenta di 1 il numero n (che abbiamo impostato a 0 prima) finché esso è minore di 5, per n=0,1,2,3,4 questa condizione è vera, quando n=5 non lo è più e il ciclo si ferma.

Cicli for
La sintassi è la seguente:
for contatore=inizio:passo:fine
istruzioni
end
oppure al contatore possiamo assegnare i valori di un vettore.

In ambedue i casi nel ciclo for un istruzione viene eseguito un numero predefinito di volte per i valori che il contato assume.
ad esempio

>for i=1:10 %(oppure i=1:1:10)
>disp('ciao')
>end

Qui il ciclo viene eseguito 10 volte, e in ogni ciclo verrà visualizzata la scritta “ciao”

Normalmente in Octave è sconsigliato, quando si può, scrivere cicli.
Ovviamente gli esempi di prima sono poco significativi, supponiamo però di voler calcolare la media di un vettore v=[1 2 3 4 5 6]
con un ciclo for scriveremmo

>n=0
>for i=v %v è il vettore, i assume i valori presenti nel vettore
>n=n+i
>end
>media=n/lenght(v)

Un modo molto più elegante è quello di scrivere

>media=sum(v)/length(v)

oppure

>mean(v)

visto che in questo caso esiste la funzione apposita.
Evitare cicli, ma effettuare operazioni vettoriali o matriciali permette di risparmiare tempo, su vettori con pochi componenti, come tutti quelli tratti in esempio, la differenza sarà praticamente nulla, ma se si deve lavorare con centinaia di centinaia di valori, dati che vengono da un esperimento, allora anche risparmiare pochi decimi di secondo ad ogni operazione può significare di risparmiare ore e ore di elaborazione dei dati per un codice scritto in maniera più pulita e efficiente.

Operatori logici
In informatica vengono spesso usato operatori logici (oppure booleani), che riportano il valore 0 se una cosa è falsa, un numero diverso da 0 (per convenzione 1) se una cosa è vera.
Gli operatori logici in Octave sono:

>,<,>=,<=,==,~=

potete provare a definire un vettore (o un numero) a e b provare i comandi

>a==b
>a~=b

e via dicendo e osservare i risultato. Le operazioni vengono fatte per ogni componente.

Vi sono poi i seguenti operatori

& , ||, ~, xor

a & b vi da una condizione vera se gli vengono passate due condizioni vere (cioè due numeri diversi da 0, ad esempio
>6 & 3)
a||b da una condizione vera se almeno una delle due condizione che gli vengono passate è vera
xor(a,b) da vero se solo una delle due condizioni è vera/falsa
~a vi da la condizione opposta; da falso se gli passata un qualcosa di vero, da vero se gli passate qualcosa di falso.

Istruzioni if
la sintassi generica è
if (condizione1==true)
istruzione1
elseif (condizione2==true)
istruzione2
elseif …
…
else
istruzione3
end

In pratica Octave controlla se una condizione è vera (cioè se nella parentesi ha un numero diverso da 0), se lo è esegue una operazione, se non lo è va a controllare se è vera un’altra condizione per eseguire un’altra istruzione.
Esempio

if (n==5) %si vedano gli operatori logici per il significato
disp('n è uguale a cinque')
elseif (n>5)
disp('n è maggiore di cinque')
else
disp('n è minore di cinque')
end

Finita anche questa lezione, eccovi gli esercizi e le soluzioni. Ovviamente è un’ottima idea salvare anche tutti gli esercizi del capitolo precedente sotto forma di script-file, in modo da averli tutti salvati insieme

Scalari, vettori e matrici in Octave

Comincia ora la prima, di una serie, guida di Octave.
In questo capitolo parlerò dei valori scalari, vettori e matrici, come anche delle operazioni principali che possiamo effettuare con esse. Per chi avesse ravanato un poco nei miei articolu più vecchi avrebbe trovato questo articolo su gnuplot, dove avevo messo una tabella riassuntiva di molte funzioni analitiche, ovviamente anche presenti in Octave.

Valori scalari
Con valori scalari si intendono le cifre, reali o complesse è indifferente. In questo modo si assegnano ad una variabile i valori scalari:

>a=5
a=
 	5
>b=5-1
b=
 	1

Se non viene specificata una variabile, il valore viene assegnato in automatico alla variabile ans

>3
ans=
 	3
>0*2
ans=
 	0

Le operazioni elementari +,-,*,/,^ sono definite nel modo naturale, attraverso le parentesi “(” e “)” è possibile dare la precedenza alle operazioni. (In Octave è possibile scrivere anche ** al posto di ^, ma in Matlab questa sintassi non è accettata.) Se scriviamo “;” alla fine di un comando, viene omesso l’output:

>z=5;

Se abbiamo bisogno di sapere quali variabili abbiamo già usato, possiamo usare il comando “who”

>who

Con il seguente comando cancelleremo tutte le variabili in memoria

>clear all

Sono definite inoltre molte altre funzioni, come ad esempio: abs(x), sqrt(x), exp(x), log(x), sin(x),cos(x), tan(x), asin(x), …
E anche alcune variabili: pi, i,j, exp(1) che ovviamente possono venir sovrascritte con una assegnazione

I vettori
I vettori vengono assegnati nel seguente modo: vettore riga

>a=[1 2 3]
a=
 	1	2	3

vettore colonna

>b=[1;2;3]
b=
 	1
 	2
	3

Si possono inoltre assegnare vettori nel seguente modo:

>c=[1:4]
c=
 	1	2	3	4
>d=[1:0.5:4]
d=
 	1.0000	1.5000	2.0000	2.5000	3.0000	3.5000	4.0000

qui 0.5 segna il “passo”, cioè la distanza tra due numeri adiacenti

>e=linspace(0,1,5) e= 	0.00000	0.25000	0.50000	0.75000	1.00000

mentre qui abbiamo l’intervallo [0,1] diviso in 5 sottointervalli larghi uguali.
Esistono svariate operazioni (non solitamente usate in matematica) che si possono fare con i vettori: per conoscere la lunghezza di un vettore:

>length(c)
ans= 	3

per accedere ad un singolo componente:

>d(2)
ans= 	1.5000

l’ultima posizione del vettore è sempre end:

>d(end)
ans= 	4.0000

Per chi sa un poco programmare: In Octave, al contrario che, ad esempio in C, l’indicizzazione comincia da 1, non da 0, infatti proviamo a vedere cosa conterrebbe il componente 0 di un vettore:

>d(0)
error: subscript indices must be either positive integers or logicals.

Ovviamente un indice non naturale non è di alcun significato.
Si può inoltre trasporre un vettore:

>c=c'
c=
 	1
 	2
 	3
 	4

Ovviamente le operazioni matematiche sono tutte definite, in modo particolare la norma euclidea:

>norm(e)
ans= 	1.3693

oppure altri tipi di norme con la seguente sintassi: “norm(v,p)” Dove p è un valore numerico positivo oppure inf. In Octave è definita anche la norma -inf (anche se non si tratta propriamente di una norma, come tutte le norme negative).
La somma e la differenza vengono trattati in modo analogo ai scalari, si deve solo fare attenzioni che i vettori abbiano le stesse dimensioni e siano tutti e due orizzontali o verticali.
Per il prodotto scalare invece ci si deve preoccupare di avere due vettori di ugual dimensioni, il primo orizzontale, il secondo verticale. In alternativa si può usare il comando “dot(v,w)” con v,w due vettori di ugual lunghezza. Octave provvederà da solo ad assicurarsi che uno sia orizzontale e l’altro verticale.
È definito anche il prodotto vettoriale attraverso il comando “cross(v,w)”, con v e w due vettori di ugual lunghezza.
È possibile sommare tutti gli elementi di un vettore, attraverso il comando “sum(v)”, con v vettore, mentre con “prod(v)” verrà eseguito il prodotto delle componenti. Le funzioni “max(v)” e “min(v)” vi daranno come output l’elemento massimo e minimo del vettore. La funzione “sort(v)” ordinerà gli elementi del vettore in ordine di grandezza mentre il comando “diff(v)” effettuerà la differenza a due a due di tutti i componenti del vettore (secondo componente meno il primo, terzo componente meno il secondo…) e come output vi darà un vettore di con n-1 entrate (supponendo che la dimensione di v fosse n) .
Inoltre è possibile fare anche alcuni operazioni “componente per componente”. Queste non sono operazioni che vengono di solito definite/usate in matematica, ma sono di largo uso in ambito numerico e permettono di risparmiare (per chi conoscesse un poco di programmazione) l’uso di cicli for.
Abbiamo ad esempio il prodotto componente per componente:

>a=[1 2 3];
>b=[4 5 6];
>c=a.*b
c=
 	4	10	8

Allo stesso modo si definisce la divisione, e l’elevamento a potenza tra due vettori. Operazioni come la somma, la differenza, il prodotto/divisione con uno scalare, sono già definite componente per componente come normalmente lo sono in matematica, mentre funzioni come modulo, seno, coseno, logaritmo, exp… oppure la somma/sottrazione con uno scalare sono anch’esse, in Octave, definite componente per componente (si può pensare che lo scalare viene prima moltiplicato con un vettore di 1 prima di effettuare la somma/sottrazione).
In sostanza le uniche funzioni per la quale serviva una definizione a parte per le operazioni “componente per componente” sono il prodotto, la divisione e l’elevamento a potenza tra due vettori.

Matrici
Come gli scalari sono un caso particolare dei vettori (si trattano, a tutti gli effetti, di vettori con una riga (o una colonna, dipende come si pone il vettore)), i vettori sono matrici particolari, nuovamente con una sola riga, oppure colonna. In Octave è possibile definire matrici n*m, e normalmente si lavora con esse. Ovviamente sono anche definite le operazioni “det(A)” (il determinante) e inv(A) (per calcolare l’inversa di una matrice), supponendo che A sia una matrice quadrata. Tutti i discorsi sulle operazioni “componente per componente” introdotte sui vettori, si applicano anche alle matrici, mentre per le funzioni che sono state introdotte ai vettori si estendono alle matrici, ricordandosi che le operazioni vengono effettuate sulle righe, e quindi l’output non sarà uno scalare, bensì un vettore.
Ad esempio

>A=[1 2 3;7 8 9; 4 5 6];
>sum(A)
ans=
 	12	15	18

La somma viene quindi fatta sulle colonne (si può pensare che le matrici siano vettori di vettori…)

>prod(A)
ans=
 	28	80	162
>sort(A)
ans=
	1	2	3
	4	5	6
 	7	8	9

Anche la funzione “max” e “min” lavorano allo stesso modo, su colonne. se invece diamo il comando

>a=diag(A)
a=
 	1 	8 	6

otteniamo in output la diagonale di A, se diamo invece

>A2=diag(a)
A2=
 Diagonal Matrix
 	1	0	0
	0	8	0
 	0	0	6

Dando “diag(a,n)” o “diag(A, n)” si sposta di n righe (verso l’alto o il basso dipende dal segno di n) la diagonale interessata. Al posto di lavorare direttamente sulle matrici, si può lavorare anche su sottomatrici, ad esempio:

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

Se vogliamo, ad esempio eliminare una colonna dobbiamo scrivere

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

Un’altra operazione importante che si può fare con le matrici (e quindi anche i vettori), è la concatenazione. Scrivete ad esempio

>B=[A A]

oppure

>B=[A;A]

Esistono poi matrici particolari che vengono spesso usate in calcolo numerico:

eye(n) %crea la matrice identità di dimensione n*n
ones(m,n) %crea una matrice di 1 di dimensione m*n
zeros(m,n) %crea una matrice di 0 di dimensione m*n
rand(n,m) %crea una matrice di dimensione n*m di valori pseudocasuali

Dopo questo lungo capitolo, eccovi un po’ di esercizi, giusto per vedere se avete capito tutto e per eventualmente chiarire dubbi!
I risultati sono salvati sui file con l’estensione .m, per ora vi basta sapere che dovete copiare-incollare le singole istruzioni sulle righe, la prossima volta vedremo invece meglio cosa sono.

Octave, installazione, personalizzazione e i numeri sul calcolatore

Voglio fare una piccola premessa.
Io il programma lo ho testato su un sistema Debian ma tutte le indicazioni fatte, rimangono invariate per tutti i sistemi (escluse le considerazioni fatte sui percorsi, su Windows saranno un po’ diversi).
Per l’installazione dovete usare i package manager presenti nella vostra distro, come aptitude oppure yum. È consigliato inoltre installare oltre al pacchetto octave anche octave-doc e octave-info per avere la documentazione. Se usate Windows o se Octave non è presente nei vostri repository recatevi al sito di Octave, i pacchetti che ho installato sul mio sistema terminale sono octave octave-doc octave-info
Come editor di testo (vedremo più avanti perché ce ne servirà uno) uso personalmente gedit, in quanto è in grado di riconoscere la sintassi di Octave, ma vi sono tantissimi editor in grado su farlo. Su Windows di default insieme ad Octave verrà installato noteppad++, ovviamente in grado di riconoscere la sintassi.

Octave, durante questo laboratorio, verrà usato dal terminale (non abbiamo infatti installato nessuna interfaccia grafica, ma niente vi vieta di cercarne una).
Ogni riga comincia con

>

I comandi scritti verranno interpretati solo dopo essere andati a capo (premere tasto enter).
Nel caso di errori di sintassi Octave produrrà un output segnando dove si trova l’errore oppure il tipo di errore commesso
Per chiudere il programma Octave basta dare il comando

>quit

oppure

>exit

Alcuni comandi che vengono normalmente usati in un sistema Gnu che si ritrovano all’interno di Octave sono:
pwd (mostra dove ci troviamo all’interno del sistema)
mkdir (serve per creare cartelle
ls (serve per visualizzare il contenuto della cartella in cui ci troviamo)
cd (serve per muoversi da una cartella all’altra)
Esempio:

>pwd
/home/fekir
>ls
Documenti	Multimedia
>mkdir Laboratorio
>ls
Documenti	Multimedia	Laboratorio
>cd Laboratorio
>pwd
/home/fekir/Laboratorio

Attenzione perché i comandi sono case sensitive, se avessimo scritto

>cd laboratorio
error: laboratorio: No such file or directory

Avremmo ottenuto un errore, che ci informa, che tale cartella (o file, ma non in questo caso) non esiste.

Spesso si inoltre a che fare con diverse funzioni predefinite e può darsi di non sapere bene come funzionano.
I comandi help e doc servono a capire come funzionano, ad esempio:

>help sqrt

>doc sin

Inoltre in Octave è possibile impostare una cartella dove posizionare le function che scriveremo. Vedremo più avanti cosa sono e come funzionano, ma possiamo già impostare una cartella.
Se abbiamo i permessi di amministratore andiamo a modificare il file

/usr/share/octave/3.2.4/m/startup/octaverc

(notate che se la vostra versione è diversa dalla 3.2.4 dovrete cambiare il percorso), dovreste trovare una linea del tipo

addpath (genpath ("/usr/local/share/octave/site-m"), "-begin");

sotto di questa aggiungiamo semplicemente

addpath("~/.octave/function")

. Dobbiamo poi creare nella nostra home la cartella .octave/function, e qui andremo, nei laboratori successivi, a inserire le function che vi darò o scriverete.
Altrimenti (sopratutto se non avete i permessi di amministratore) potete farlo direttamente da Octave dando il comando

addpath("~/.octave/function")

Il vantaggio di fare l’operazione da amministratore era nel caso aveste un computer con più utenti, in tal caso la modifica avrebbe interessato tutti gli utenti, in questo modo invece soltanto il vostro profilo.
Se volete rimuovere tale percorso basta rimuovere la riga scritta oppure dare da Octave il comando

rmpath("~/.octave/function")

Questa nota vale agli utenti che usano Windows, in quanto sui Sistemi Microsoft il terminale appare molto povero. Mentre sui sistemi Gnu/Linux il programma octave parte con le impostazioni predefinite del vostro terminale, su Windows (almeno nel mio caso quando ho fatto una prova), partiva con dei caratteri un po’ bruttini e troppo piccoli.
È per fortuna possibile personalizzare l’ambiente, facendo click destro sul collegamento a Octave (o direttamente con l’eseguibile) si può impostare i font da usare, la dimensione dei caratteri, i colori, e molto altro.

I numeri a virgola mobile
Vediamo adesso di dare una veloce (molto veloce) rappresentazione dei numeri sul calcolatore.
Sui calcolatori, si sa, la memoria disponibile è finita, quindi è ad esempio impossibile scrivere numeri irrazionali, che dopo la virgola hanno infinite cifre. Dal momento che quindi è impossibile rappresentare la maggior parte dei numeri reali, sui calcolatori si lavora con i numeri floating point, altrimenti chiamati numeri a virgola mobile.
Questi numeri rappresentano i numeri reali fino ad una certa cifra, dopo di essa vi è un troncamento. Normalmente gli errori sono molto piccoli, “invisibili” all’occhio umano, però, un errore, per quanto piccolo sia, si può propagare e aumentare, portando a dei risultati “sbagliati”.
Vediamo un semplice esempio, scriviamo in Octave:

>n=2;
>while (n>0)
>n=n/2;
>end
>n

Vediamo cosa abbiamo pressapoco scritto. Dapprima abbiamo assegnato ad n il valore due, poi diciamo ad Octave che finché n non è pari a 0, n verrà dimezzato.
Ovviamente questo ciclo dovrebbe risultare infinito. Nel senso che se dividete un numero in due, la sua metà sarà ancora diversa da zero. Ovviamente se pensate al significato analitico di quello che abbiamo scritto, stiamo facendo tendere n a 0, ma vi ricordo che questo è un limite per quando facciamo infinite operazioni. Eppure Octave si ferma dopo poco e valuta n come se fosse nullo proprio perché l’errore di cui vi ho parlato prima rende n e il numero 0 uguali.
Se riscrivete il ciclo di prima, senza il “;” dopo n=n/2, vi verrà mostrata ogni iterazione, e potrete notare come viene dimezzato n di volta in volta, fino a diventare un numero piccolissimo: 4.9407e-324
Vuol dire che ci sono 324 zeri prima che appaia una cifra diversa da zero. L’errore è senz’altro minuscolo, però rischia appunto di aumentare se il numero n venisse poi usato per altri calcoli, poiché (in questo caso specifico) moltiplicare per 0 o un numero diverso costituisce una grossa differenza. E sono state necessario “solo” 1076 iterazioni per annullare n.
Questi numeri possono apparire assurdi, ma durante la simulazioni di modelli matematici poter gestire numeri minuscoli o giganteschi mette a dura prova la precisione dei calcolatori.
Esistono poi dei numeri un poco particolari: Inf, -Inf, NaN, cioè infinito, meno infinito e Not a Number.
I primi due numeri si ottengono quando ad esempio state considerano numeri (in valore assoluto) troppo grandi. Visto che la memoria del calcolatore è finita, da un certo numero in poi (un numero a 308 cifre) i numeri verranno trattati come infinito. NaN è una espressione che si ottiene moltiplicando un numero infinito con un numero nullo, oppure facendo la differenza tra due infiniti, casi in cui, a priori, è impossibile dire quale sarà il risultato.
Ci sarebbe ancora tanto da dire sulla rappresentazione numerica sui calcolatori, ad esempio non valgono sempre le proprietà commutative o associative, ma penso che questa piccola parte sia sufficiente per capire un poco cosa succede dietro lo schermo e perché esistono gli errori di arrotondamento.

Octave al posto di Matlab

In questo Blog vi avevo già parlato una volta, di Octave (qui e qui), e ora torno a farlo nuovamente, per un motivo diverso.

In ambiente universitario può capitare che gli studenti debbano fare dei laboratori, ad esempio, di calcolo numerico. E il software che viene dato in dotazione è, purtroppo, troppo spesso Matlab. Ultimamente si parla sempre di crisi, che vogliono togliere i soldi alle università e via dicendo. Ora, non voglio fare polemica (ma parlare di software open) e analizziamo un attimo i fatti.

Personalmente non ho assolutamente nulla contro Matlab, probabilmente è un ottimo prodotto e di qualità, ma a mio parere non vale la pena spendere centinaia e centinaia di euro per tenere su dei laboratori, quando alternative libere come Octave, Scilab o altri ancora possono essere più che sufficienti. Anche perché, tra i vari problemi che incorrono nel software proprietario, escludendo il fattore puramente economico, si ha che (nel caso di Matlab, ma spesso si può generalizzare)
1) Le richieste minime dell’hardware sono assai più superiori delle controparti libere, e spesso i computer delle aule informatiche non sono recentissimi. Matlab richiede ad esempio 2 Gb di ram (minimo assoluto 1) e 10Gb di spazio sul disco. Ovviamente i computer di oggi sono molto più potenti, ma su computer di anche solo 5 anni fa i due Gb di ram cominciano a diventare pochi…
2) Vi è un incremento sostanziale alla pirateria, lo studente medio vorrebbe poter lavorare a casa con gli stessi programmi che si usano in laboratorio (giustamente)
3) Ovviamente mancano tutti i vantaggi del software libero (poter migliorare il programma alle proprie esigenze, offrire collaborazione per avere la documentazione italiana, poter aggiungere funzionalità intrinseche che possono servire anche per la ricerca all’interno di un dipartimento…)

Notate che non sto dicendo che le università dovrebbero smettere di usare Matlab dall’oggi al domani. Sto dicendo che è inutile usare Matlab se sono disponibili alternative libere, sopratutto per non mettere in difficoltà lo studente. Se poi il ricercatore universitario deve assolutamente usare Matlab perché gli servono delle funzioni che li sono presenti e in Octave/Scilab/altri software no e non si riescono a implementare facilmente, allora mi va bene che sul suo computer (e quindi una licenza, non centinaia…) ci sia Matlab. Ma dovrebbe essere appunto un poco l’eccezione e non la regola.

L’alternativa che voglio proporre è quindi Octave. Ovviamente la scelta poteva ricadere anche su altri programmi, ho scelto Octave perché ho già scritto su questo programma qualcosa in merito, e mi sono accorto che la sintassi tra Octave e Matlab, per quanto ho studiato, è pressapoco identica. Ho trovato soltanto un paio di cose che con Matlab avrei potuto fare in un certo modo e Octave no (almeno non al primo colpo), ma niente che potesse incidere sulla qualità dei due programmi. In sostanza se al posto di aver usato Matlab, gli studenti avessero usato Octave, non ci sarebbe stata nessuna differenza. I comandi, gli script-file, le funzioni…tutto avrebbe funzionato alla perfezione senza cambiare una virgola!
Ovviamente dire che tutto funzioni alla perfezione è quasi assurdo, alcune funzioni di Matlab in Octave non sono ancora presenti e la sintassi di Octave è molto più “ricca”  (poichè riprende in alcuni contesti anche la sintassi che si usa nel terminale) in confronto a Matlab, però la compatibilità è talmente alta da poter permettere di tenere la maggior parte dei laboratori inalterati dal punto di vista degli esercizi e dei contenuti.
Incuriosito dal fatto che si usasse un programma che difficilmente uno studente potesse avere a casa mi sono un poco informato sull’argomento. Il motivo ufficiale (o meglio: quello che mi è stato dato) è che Matlab sa gestire meglio i grafici.
La cosa, almeno per quello che ne so io, è abbastanza vera. Di default i grafici di Octave sono un po’ poveri e si deve prestare quindi un attimo di attenzione in più, attenzione che viene però a costare fino a 500 euro a licenza! Secondo me non ne vale quindi la pena, visto che non sembra portare altri vantaggi.

Per mostrare la veridicità di quello che sto dicendo ho deciso di “tenere un corso di Octave” (con calma, quando il tempo mi permette di scrivere qualcosa). Gli argomenti trattati saranno puramente matematici, spiegherò alcuni metodi di Calcolo Numerico su come risolvere determinati problemi e introdurrò la programmazione su Octave. Se avete la possibilità di riprovare tutto quello che ho scritto su Octave su Matlab, fate pure, non dovreste cambiare nemmeno una virgola e viceversa. Il corso non sarà incentrato su Octave stesso (altrimenti rischierei di perdere la compatibilità con Matlab visto che alcune differenze esistono!), quanto più su un suo metodo di utilizzo e alcuni metodi numerici per l’approssimazione di soluzioni, implementati poi ovviamente per i software come Octave e Matlab.
Il “corso” che intendo tenere vorrei dividerlo nelle seguenti lezioni:
1) Questa introduzione
2) Come installare Octave, personalizzazioni e come vengono rappresentati i numeri sul calcolatore
3) Scalari, Vettori, Matrici e funzioni elementari
4) Programmare in Octave, script-file, function, operatori logici, cicli
5) Disegnare grafici, rappresentare funzioni, polinomi
Dopo questa prima parte interamente dedicata a Octave, comincia quella più legata alla Matematica.
6) Interpolazione polinomiale
7) Spline interpolanti e Minimi quadrati discreti
8) Integrazione numerica
9) Risoluzione di equazioni non lineari
10) Localizzazione di Autovalori
11) Risoluzione di sistemi lineari
Vedrò di mettere online inoltre esercizi e soluzioni, in modo da poter appurare di aver capito le lezioni e di rendere il tutto un poco più interessante.