Informatica per Scienze Geologiche LTa.a.2017-2018

Esercitazione Analisi dati PozzoTrebiciano

Parte B2Parte B2

Slides preparate da Dott. Pastorutti

Docente: Prof. Carla Braitenberg, Dipartimento Matematica e Geoscienze, Via Weiss 1, Università di Trieste

E-mail: berg@units.it Tel. 040 5582258

Esercizio livello acqua nell’Abisso di Trebiciano

• Anni di osservazione: 2006-2008, un campione ogni ora. Unita’ di misura: m sopra il livello del mare.

• Caricare il file. Nome dell’array: liv_Trebiciano

load('livello_AbissoTrebiciano_Tdec.mat');

• Determinare valore medio, valore massimo e data del valore massimo valore massimo

• Comandi necessari:

• max(livello_tdec): cerca il massimo sulla quinta colonna, comprendente il livello.

• mean(livello_tdec)

• std(livello_tdec)

%Max_TrebicianoWater

%read the data file and find extreme value.

load('livello_AbissoTrebiciano_Tdec.mat');

liv=livello_tdec;

% contains the workspace liv_Trebiciano with year, month, day hour, value

n=length(liv(:,1));

% extreme value

[xmax, kmax] = max(liv(:,6));[xmax, kmax] = max(liv(:,6));

tmax=liv(kmax,2:4);

disp([' max value at ' ,num2str(tmax),' value= ', num2str(liv(kmax,6)),'m']);

media=mean(liv(:,6)); scarto= std(liv(:,6));

disp(['average ' ,num2str(media),'m',' std= ', num2str(scarto),'m' ] );

Continua il programma: qui la parte che crea il grafico

figure;

plot(liv(:,1),liv(:,6))

title('Livello acqua nel pozzo Abisso Trebiciano');

xlabel('tempo in anni');

ylabel('livello (m)');ylabel('livello (m)');

Ricerca dei massimi localinella serie livello_AbissoTrebiciano

%% valori estremiload('livello_AbissoTrebiciano_Tdec_CORRETTO.mat');liv=livello_tdec;n=length(liv(:,1));[xmax, kmax] = max(liv(:,6));tmax=liv(kmax,2:4);disp(['max value at ' ,num2str(tmax),...

• Lo script fatto finora: • Risultato:

Parte B/2

' value= ', num2str(liv(kmax,6)),'m']);media = mean(liv(:,6));scarto = std(liv(:,6));disp(['average ' ,num2str(media),'m',...

' std= ', num2str(scarto),'m' ] );

%% graficofigura1 = figure;assi1 = axes('Parent',figura1);plot(assi1,liv(:,1),liv(:,6))hold ontitle('Livello acqua nel pozzo Abisso Trebiciano');xlabel('tempo in anni');ylabel('livello (m)');

Ricerca dei massimi localinella serie livello_AbissoTrebiciano

...[xmax, kmax] = max(liv(:,6));...

• Massimo assoluto: possiamo usare xmax, kmax già calcolati prima

scatter(assi1,liv(kmax,1),xmax)

x y

NB: la figura ottenuta prima non va chiusa,

altrimenti:

x y

Ricerca dei massimi localinella serie livello_AbissoTrebiciano

function [M,I] = maxsort(in,n)M = zeros(1,n); I = M;for i=1:n

[M(i),I(i)] = max(in);in(I(i)) = 0;

endend

• cercare i primi n valori più grandi: definiamo una funzione

1234567 end

• Salvate la function in un file con lo stesso nome della chiamata, quindi maxsort.m

«Ottieni il valore M e l’indice I degli n valori più grandi di in.»

Nella riga 2 vengono creati M ed I, vettori lunghi n, composti di soli zeri.

Le righe da 3 a 6 appartengono a un ciclo for: ad ogni iterazione viene calcolato l’n-esimo massimo e il suo indice, quindi quell’elemento di in viene sostituito con zero.

7

Ricerca dei massimi localinella serie livello_AbissoTrebiciano

[massimi, indici] = maxsort(liv(:,6),100);scatter(assi1,liv(indici,1),massimi);

• Usiamo la nuova funzione «maxsort» nel nostro script

i primi 100 valori

Ricerca dei massimi localinella serie livello_AbissoTrebiciano

[massimi, indici] = maxsort(liv(:,6),1000);scatter(assi1,liv(indici,1),massimi);

• Usiamo la nuova funzione «maxsort» nel nostro script

• Non è una strategia efficace.• Non è una strategia efficace.

• Potremmo calcolare i massimi in finestre di tempo, ma e’ una strategia soggettiva: devo definire inizio e lunghezza della finestra.

Ricerca dei massimi localinella serie livello_AbissoTrebiciano

• Un altro metodo: calcoliamo le derivate prima e seconda

• Cerchiamo i punti in cui

– la derivata prima è zero (entro una tolleranza)

– la derivata seconda è negativa (concavità verso l’alto)

• Produciamo tre grafici sovrapposti, usando subplot

Calcoliamo la derivata (numerica) con diff:Calcoliamo la derivata (numerica) con diff:

in questa maniera usiamo il passo corretto tra due campioni, che non è costante

step_livelli = liv(2:end,1) - liv(1:end-1,1);der1 = diff(liv(:,6)) ./ step_livelli;der2 = diff(der1) ./ step_livelli(1:end-1);der1 = der1(1:end-1);

dopo aver calcolato diff, il risultato è più corto di un elemento. Rendiamo la derivata 1° lunga quanto la derivata 2°, togliendo l’ultimo elemento

(fermandoci a end-1)

Dal help di matlab: diff(X), for a vector X, is [X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)].

Ricerca dei massimi localinella serie livello_AbissoTrebiciano

• Approssimiamo la derivata prima «uguale a zero» entro una tolleranza (più o meno tol_der1)

• Approssimiamo la derivata seconda «negativa» entro una tolleranza (minore di tol_der2)

• Teniamo solo i massimi locali trovati che sono superiori a tol_liv

tol_der1 = 1.25e3;tol_der2 = -0.001;tol_liv = 20;massimi_locali = find(and(and(and(der1<=tol_der1,der1>=-tol_der1),...

der2<=tol_der2),...der2<=tol_der2),...liv(1:end-2,6)>=tol_liv));

Ricerca dei massimi localinella serie livello_AbissoTrebiciano

• Tre grafici, con subplot (livello, prima derivata, seconda derivata)

figura2 = figure;AssiLivelli = subplot(4,1,[1 2]);plot(AssiLivelli,liv(:,1),liv(:,6))hold onscatter(AssiLivelli,liv(massimi_locali,1),liv(massimi_locali,6));

AssiDer1 = subplot(4,1,3);AssiDer1 = subplot(4,1,3);plot(AssiDer1,liv(1:end-2,1),der1)

AssiDer2 = subplot(4,1,4);plot(AssiDer2,liv(1:end-2,1),der2)

linkaxes([AssiLivelli,AssiDer1,AssiDer2],'x');