[Mat Numerica] [RISOLTO] Risoluzione esercizio con Matlab
Ciao a tutti,
ho da poco iniziato un corso di matematica numerica e mi è stato assegnato un esercizio, e vorrei dei chiarimenti da parte di qualcuno che ha conoscenze consolidate in materia.
Faccio una piccola premessa per inquadrarvi lo scenario.
In pratica l'esercizio è incentrato sulla diffusine dell'inquinamento nei bacini d'acqua, in particolare nei laghi stratificati (ovvero in estate i laghi in zone temperate possono, da un punto di vista termico, diventare stratificati. Questa stratificazione divide il lago in due parti: l’epilimnio e l’ipolimnio separati da un piano detto Termoclino). È interessante studiare l’inquinamento in tale ambienti in quanto il termocline diminuisce molto lo scambio tra i due ambienti.
L'obiettivo dell'esercizio è conoscere la quota di profondità del Termoclino (lo stato intermedio).
Dalle basi teoriche sono giunto alla conclusione che per calcolare la posizione del termoclino occorre trovare il punto di flesso della curva temperatura-profondità, cioè è il punto in cui si annulla la derivata seconda della temperatura in funzione della profondità, che, poi, è anche il punto in cui è massima la derivata prima.
Ora parto dal principio così potete seguire:
Ho creato 2 vettori che contengono i valori di temperatura e profondità:
profondita= $ ( ( 20.6, 20.6, 20.6, 18.4, 12.7, 9.5, 8.9, 8.9 ) ) $;
temperatura=$ ( ( 0, 0.1, 2.7, 6.9, 11.5, 16.1, 20.7, 25.0 ) ) $;
Poi con il comando:
pp=spline(profondita,temperatura)
ho creato una variabile strutturata da cui è possibile estrarre i coefficienti su ogni intervallino della spline cubica interpolante i punti (profondità,temperatura).
Poi con il comando:
[x,C,l,k,d] = unmkpp(pp)
nella matrice C ho memorizzato di ogni intervallino i 4 coefficienti che individuano il polinomio di 3° grado:
Per ottenere i coefficienti della derivata prima della spline ho calcolato la derivata prima del polinomio di 3° grado :
Cder=[3*C(:,1) 2*C(:,2) C(:,3)];
Cder è quindi la matrice che contiene i coefficienti della derivata della spline.
Poi con il comando:
ppder=mkpp(x,Cder);
ottengo una nuova funzione ppder che rappresenta la derivata della spline nei nodi x.
Il procedimento fino ad ora è corretto oppure ho commesso qualche errore?
Se è corretto, vado avanti... altrimenti faccio mente locale e ricontrollo i passaggi.
Grazie a tutti coloro che vorranno intervenire.
ho da poco iniziato un corso di matematica numerica e mi è stato assegnato un esercizio, e vorrei dei chiarimenti da parte di qualcuno che ha conoscenze consolidate in materia.
Faccio una piccola premessa per inquadrarvi lo scenario.
In pratica l'esercizio è incentrato sulla diffusine dell'inquinamento nei bacini d'acqua, in particolare nei laghi stratificati (ovvero in estate i laghi in zone temperate possono, da un punto di vista termico, diventare stratificati. Questa stratificazione divide il lago in due parti: l’epilimnio e l’ipolimnio separati da un piano detto Termoclino). È interessante studiare l’inquinamento in tale ambienti in quanto il termocline diminuisce molto lo scambio tra i due ambienti.
L'obiettivo dell'esercizio è conoscere la quota di profondità del Termoclino (lo stato intermedio).
Dalle basi teoriche sono giunto alla conclusione che per calcolare la posizione del termoclino occorre trovare il punto di flesso della curva temperatura-profondità, cioè è il punto in cui si annulla la derivata seconda della temperatura in funzione della profondità, che, poi, è anche il punto in cui è massima la derivata prima.
Ora parto dal principio così potete seguire:
Ho creato 2 vettori che contengono i valori di temperatura e profondità:
profondita= $ ( ( 20.6, 20.6, 20.6, 18.4, 12.7, 9.5, 8.9, 8.9 ) ) $;
temperatura=$ ( ( 0, 0.1, 2.7, 6.9, 11.5, 16.1, 20.7, 25.0 ) ) $;
Poi con il comando:
pp=spline(profondita,temperatura)
ho creato una variabile strutturata da cui è possibile estrarre i coefficienti su ogni intervallino della spline cubica interpolante i punti (profondità,temperatura).
Poi con il comando:
[x,C,l,k,d] = unmkpp(pp)
nella matrice C ho memorizzato di ogni intervallino i 4 coefficienti che individuano il polinomio di 3° grado:
Per ottenere i coefficienti della derivata prima della spline ho calcolato la derivata prima del polinomio di 3° grado :
Cder=[3*C(:,1) 2*C(:,2) C(:,3)];
Cder è quindi la matrice che contiene i coefficienti della derivata della spline.
Poi con il comando:
ppder=mkpp(x,Cder);
ottengo una nuova funzione ppder che rappresenta la derivata della spline nei nodi x.
Il procedimento fino ad ora è corretto oppure ho commesso qualche errore?
Se è corretto, vado avanti... altrimenti faccio mente locale e ricontrollo i passaggi.
Grazie a tutti coloro che vorranno intervenire.
Risposte
Dovrei ottenere questi grafici:
Il primo, lo ottengo facilmente con:
plot(temperatura,-profondita)
Mentre per il 2° e il 3° incontro qualche difficoltà...
Il primo, lo ottengo facilmente con:
plot(temperatura,-profondita)
Mentre per il 2° e il 3° incontro qualche difficoltà...
"terminal":
Ora parto dal principio così potete seguire:
Ho creato 2 vettori che contengono i valori di temperatura e profondità:
profondita= $ ( ( 20.6, 20.6, 20.6, 18.4, 12.7, 9.5, 8.9, 8.9 ) ) $;
temperatura=$ ( ( 0, 0.1, 2.7, 6.9, 11.5, 16.1, 20.7, 25.0 ) ) $;
Poi con il comando:
pp=spline(profondita,temperatura)
ho creato una variabile strutturata da cui è possibile estrarre i coefficienti su ogni intervallino della spline cubica interpolante i punti (profondità,temperatura).
Non so se sia un errore di scrittura sul forum o direttamente nel tuo programma, ma:
pp=spline(profondita,temperatura)
dovrebbe dare errore, dato che la variabile profondità ha più elementi che assumono lo stesso valore. Dovresti invertire
pp=spline(temperatura,profondita).
Ciao Lory314,
intanto grazie per il tuo interesse; riguardo il tuo post:
- perchè dovrebbe dare errore se ci sono più elementi con lo stesso valore?
Il comando viene eseguito regolarmente.
intanto grazie per il tuo interesse; riguardo il tuo post:
- perchè dovrebbe dare errore se ci sono più elementi con lo stesso valore?
Il comando viene eseguito regolarmente.
La cosa mi sembra strana. Per due motivi:
1) A me provando, risluta un errore. Mi sembra strano che non lo dia a te. Che versione di MatLab stai usando?
2) Prova a pensare a quello che stai facendo. Se metti sull'asse orizzontale la profondità e sull'asse vertiale la temperatura, rappresentando i tuoi dati avrai 3 punti che hanno la stessa ascissa, ma diversa ordianata e quindi saranno allineati su una retta parallela all'asse della temperatura. Già di per sè, sarebbe impossibile costruire una funzione passante per quei tre punti; tantomeno è impossibile costruire una funzione interpolante. Nota infatti che i grafici che hai riportato hanno la temperatura come variabile indipendente e la profondità come variabile dipendente; usando il comando spline come lo stai usando te stai facendo il contrario. Spero di essere stato chiaro e, più che altro di non sbagliarmi.
1) A me provando, risluta un errore. Mi sembra strano che non lo dia a te. Che versione di MatLab stai usando?
2) Prova a pensare a quello che stai facendo. Se metti sull'asse orizzontale la profondità e sull'asse vertiale la temperatura, rappresentando i tuoi dati avrai 3 punti che hanno la stessa ascissa, ma diversa ordianata e quindi saranno allineati su una retta parallela all'asse della temperatura. Già di per sè, sarebbe impossibile costruire una funzione passante per quei tre punti; tantomeno è impossibile costruire una funzione interpolante. Nota infatti che i grafici che hai riportato hanno la temperatura come variabile indipendente e la profondità come variabile dipendente; usando il comando spline come lo stai usando te stai facendo il contrario. Spero di essere stato chiaro e, più che altro di non sbagliarmi.
1) Utilizzo la r2012a (7.14.0.739) - 64 bit - verisione per Mac;
2) A pensarci bene (e riguardando bene appunti e libro) non posso far altro che darti ragione... Gli algoritmi di interpolazione hanno come condizione iniziale che:
dati n+1 punti sperimentali
(x0,y0), (x1,y1), (x2,y2), ... , (xn,yn)
con x0 < x1 < x2 < ... < xn
si vuole determinare una funzione fn(x) tale che
fn (xi) = yi , y = 0 ,1 , 2, ..., n
Ovvero che le ascisse siano distinte.
Sottolineo che il MatLab non dà alcun errore e se eseguo il comando:
MatLab da errore:
2) A pensarci bene (e riguardando bene appunti e libro) non posso far altro che darti ragione... Gli algoritmi di interpolazione hanno come condizione iniziale che:
dati n+1 punti sperimentali
(x0,y0), (x1,y1), (x2,y2), ... , (xn,yn)
con x0 < x1 < x2 < ... < xn
si vuole determinare una funzione fn(x) tale che
fn (xi) = yi , y = 0 ,1 , 2, ..., n
Ovvero che le ascisse siano distinte.
Sottolineo che il MatLab non dà alcun errore e se eseguo il comando:
pp=spline(profondita,temperatura);
MatLab da errore:
[color=#FF0000]Error using chckxy (line 51) The data sites should be distinct. Error in spline (line 54) [x,y,sizey,endslopes] = chckxy(x,y);[/color]
"terminal":
Sottolineo che il MatLab non dà alcun errore e se eseguo il comando:
pp=spline(profondita,temperatura);
MatLab da errore:
[color=#FF0000]Error using chckxy (line 51) The data sites should be distinct. Error in spline (line 54) [x,y,sizey,endslopes] = chckxy(x,y);[/color]
Non ho capito questa parte; comunque l'errore che hai riportato è quello che ti dovrebbe dare se utilizzi un vettore con valori uguali.
Allora:
se eseguo:
Matlab non mi da errore, e crea una struct 1x1;
se eseguo:
MatLab da errore:
pp=spline(profondita,temperatura)
Matlab non mi da errore, e crea una struct 1x1;
pp=spline(temperatura,profondita);
MatLab da errore:
Error using chckxy (line 51)
The data sites should be distinct.
Error in spline (line 54)
[x,y,sizey,endslopes] = chckxy(x,y);
Ok. Quindi si comporta all'opposto di come si comporta il mio. Io ho una versione più vecchia e uso windows.
Prova a guardare l'help o la doc del comando spline e vedi se ne vieni a capo. Altrimenti postala e vediamo di capire se magari nelle versione successive alla mia hanno modificato qualcosa.
Prova a guardare l'help o la doc del comando spline e vedi se ne vieni a capo. Altrimenti postala e vediamo di capire se magari nelle versione successive alla mia hanno modificato qualcosa.
Ora sono a lezione, nel pomeriggio-sera lo studio meglio; per ora lo posto così (se ti colleghi prima di me) puoi già dargli spline un'occhiata
Cubic spline data interpolation.
PP = spline(X,Y) provides the piecewise polynomial form of the
cubic spline interpolant to the data values Y at the data sites X,
for use with the evaluator PPVAL and the spline utility UNMKPP.
X must be a vector.
If Y is a vector, then Y(j) is taken as the value to be matched at X(j),
hence Y must be of the same length as X -- see below for an exception
to this.
If Y is a matrix or ND array, then Y(:,...,:,j) is taken as the value to
be matched at X(j), hence the last dimension of Y must equal length(X) --
see below for an exception to this.
YY = spline(X,Y,XX) is the same as YY = PPVAL(spline(X,Y),XX), thus
providing, in YY, the values of the interpolant at XX. For information
regarding the size of YY see PPVAL.
Ordinarily, the not-a-knot end conditions are used. However, if Y contains
two more values than X has entries, then the first and last value in Y are
used as the endslopes for the cubic spline. If Y is a vector, this
means:
f(X) = Y(2:end-1), Df(min(X))=Y(1), Df(max(X))=Y(end).
If Y is a matrix or N-D array with SIZE(Y,N) equal to LENGTH(X)+2, then
f(X(j)) matches the value Y(:,...,:,j+1) for j=1:LENGTH(X), then
Df(min(X)) matches Y(:,:,...:,1) and Df(max(X)) matches Y(:,:,...:,end).
Example:
This generates a sine-like spline curve and samples it over a finer mesh:
x = 0:10; y = sin(x);
xx = 0:.25:10;
yy = spline(x,y,xx);
plot(x,y,'o',xx,yy)
Example:
This illustrates the use of clamped or complete spline interpolation where
end slopes are prescribed. In this example, zero slopes at the ends of an
interpolant to the values of a certain distribution are enforced:
x = -4:4; y = [0 .15 1.12 2.36 2.36 1.46 .49 .06 0];
cs = spline(x,[0 y 0]);
xx = linspace(-4,4,101);
plot(x,y,'o',xx,ppval(cs,xx),'-');
Class support for inputs x, y, xx:
float: double, single
See also interp1, pchip, ppval, mkpp, unmkpp.
Reference page in Help browser
doc spline
Cubic spline data interpolation.
PP = spline(X,Y) provides the piecewise polynomial form of the
cubic spline interpolant to the data values Y at the data sites X,
for use with the evaluator PPVAL and the spline utility UNMKPP.
X must be a vector.
If Y is a vector, then Y(j) is taken as the value to be matched at X(j),
hence Y must be of the same length as X -- see below for an exception
to this.
If Y is a matrix or ND array, then Y(:,...,:,j) is taken as the value to
be matched at X(j), hence the last dimension of Y must equal length(X) --
see below for an exception to this.
YY = spline(X,Y,XX) is the same as YY = PPVAL(spline(X,Y),XX), thus
providing, in YY, the values of the interpolant at XX. For information
regarding the size of YY see PPVAL.
Ordinarily, the not-a-knot end conditions are used. However, if Y contains
two more values than X has entries, then the first and last value in Y are
used as the endslopes for the cubic spline. If Y is a vector, this
means:
f(X) = Y(2:end-1), Df(min(X))=Y(1), Df(max(X))=Y(end).
If Y is a matrix or N-D array with SIZE(Y,N) equal to LENGTH(X)+2, then
f(X(j)) matches the value Y(:,...,:,j+1) for j=1:LENGTH(X), then
Df(min(X)) matches Y(:,:,...:,1) and Df(max(X)) matches Y(:,:,...:,end).
Example:
This generates a sine-like spline curve and samples it over a finer mesh:
x = 0:10; y = sin(x);
xx = 0:.25:10;
yy = spline(x,y,xx);
plot(x,y,'o',xx,yy)
Example:
This illustrates the use of clamped or complete spline interpolation where
end slopes are prescribed. In this example, zero slopes at the ends of an
interpolant to the values of a certain distribution are enforced:
x = -4:4; y = [0 .15 1.12 2.36 2.36 1.46 .49 .06 0];
cs = spline(x,[0 y 0]);
xx = linspace(-4,4,101);
plot(x,y,'o',xx,ppval(cs,xx),'-');
Class support for inputs x, y, xx:
float: double, single
See also interp1, pchip, ppval, mkpp, unmkpp.
Reference page in Help browser
doc spline
Allora, il comando fa le stesse identiche cose della mia versione. Il dubbio che mi era venuto è che magari con una versione più aggiornata (la tua) avessero cambiato qualcosa (ad esempio l'ordine degli input).
A questo punto direi che, se non hai ancora risolto, potresti provare a postare tutto il codice del tuo programma e vediamo cosa salta fuori.
A questo punto direi che, se non hai ancora risolto, potresti provare a postare tutto il codice del tuo programma e vediamo cosa salta fuori.
In questi giorni, sono andato comunque avanti col codice... anche se rimane sempre il dubbio delle ascisse dei punti...
Sono riuscito a plottare sia la spline, sia la funzione che rappresenta la derivata prima della spline nei nodi x, sia la funzione che rappresenta la derivata seconda della spline nei nodi x.
Quindi la costruzione dei grafici è riuscita
L'ultimo punto dell'esercizio chiede di individuare la quota del Termoclino; graficamente sarebbe il nodo "estremo sinistro" del grafico della derivata della spline; come lo potrei calcolare?
Sono riuscito a plottare sia la spline, sia la funzione che rappresenta la derivata prima della spline nei nodi x, sia la funzione che rappresenta la derivata seconda della spline nei nodi x.
Quindi la costruzione dei grafici è riuscita
% Vettore contenente i valori della temperatura
temp=[22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1];
% Vettore contenente i valori della profondità
prof=[0 2.3 4.9 9.1 13.7 18.3 22.9 27.2];
% Costruzione della spline cubica naturale interpolante memorizzando le informazioni in una variabile strutturata detta pp-form
pp=spline(prof,temp);
% Estrae le componenti della pp-form, essendo x il vettore dei nodi, C la matrice dei coefficienti di dimensione lxk la cui j-ma riga coincide con i coefficienti del j-esimo polinomio che rappresenta la spline
[x,C,l,k,d]=unmkpp(pp);
% Coefficienti della derivata prima della spline
Cder=[3*C(:,1) 2*C(:,2) C(:,3)];
ppder=mkpp(x,Cder);
zz=0:0.1:30;
spder=ppval(ppder,zz);
% Coefficienti della derivata seconda della spline
Cder2=[6*C(:,1) 2*C(:,2)];
ppder2=mkpp(x,Cder2);
zz2=0:0.1:30;
spder2=ppval(ppder2,zz2);
%Plot
subplot(1,3,1),plot(temp,-prof,'-o'), title('Z (T)'), xlabel('T (∞C)'), ylabel('z (m)'), axis([0 30 -30 0]);
subplot(1,3,2),plot(spder,-zz), title('dz/dT');
subplot(1,3,3),plot(spder,-zz), title('d^2/dT^2');
L'ultimo punto dell'esercizio chiede di individuare la quota del Termoclino; graficamente sarebbe il nodo "estremo sinistro" del grafico della derivata della spline; come lo potrei calcolare?
Scoperto il mistero!!!
Questo è quello che hai scritto nel primo post.
Questo è quello che hai scritto nel codice:
A parte i valori diversi, in questo caso i valori uguali sono sulle temperature, quindi sulle ordinate e questo non è un problema.
Questo è quello che hai scritto nel primo post.
"terminal":
profondita= $ ( ( 20.6, 20.6, 20.6, 18.4, 12.7, 9.5, 8.9, 8.9 ) ) $;
temperatura=$ ( ( 0, 0.1, 2.7, 6.9, 11.5, 16.1, 20.7, 25.0 ) ) $;
Questo è quello che hai scritto nel codice:
"terminal":
% Vettore contenente i valori della temperatura temp=[22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1]; % Vettore contenente i valori della profondità prof=[0 2.3 4.9 9.1 13.7 18.3 22.9 27.2];
A parte i valori diversi, in questo caso i valori uguali sono sulle temperature, quindi sulle ordinate e questo non è un problema.
"terminal":
L'ultimo punto dell'esercizio chiede di individuare la quota del Termoclino; graficamente sarebbe il nodo "estremo sinistro" del grafico della derivata della spline; come lo potrei calcolare?
Un'idea potrebbe essere quella di considerare il grafico con gli assi invertiti e prendere il massimo di quella funzione.
"Lory314":
[quote="terminal"]
L'ultimo punto dell'esercizio chiede di individuare la quota del Termoclino; graficamente sarebbe il nodo "estremo sinistro" del grafico della derivata della spline; come lo potrei calcolare?
Un'idea potrebbe essere quella di considerare il grafico con gli assi invertiti e prendere il massimo di quella funzione.[/quote]
Cavolo! Avevo sbagliato... colpa mia! Scusami
Si, per quanto riguarda quel punto avevo pensato anche io di invertire gli assi (avevo, però, pensato di individuare il minimo ruotando il sistema di +90)...
C'è qualche comando che effettua questo tipo di operazione?
Si be massimo o minimo dipende da come prendi gli assi. Ti butto lì un'idea...prova a usare la matrice di rotazione....
La matrice di rotazione +90 dovrebbe essere :
R=[cos(pi/4) sin(pi/4) 0; -sin(pi/4) cos(pi/4) 0; 0 0 1];
Giusto?
R=[cos(pi/4) sin(pi/4) 0; -sin(pi/4) cos(pi/4) 0; 0 0 1];
Giusto?
Quella sicuramente non può essere la rotazione che ti seve dato che è una matrice $3 x 3$ e tu sei in dimensione $2$....
La rotazione (antioraria) rispetto all'origini di un angolo $\theta$ è data da:
$[[\cos(\theta), -\sin(\theta)],[\sin(\theta),\cos(\theta)]]$
La rotazione (antioraria) rispetto all'origini di un angolo $\theta$ è data da:
$[[\cos(\theta), -\sin(\theta)],[\sin(\theta),\cos(\theta)]]$
Ok, grazie...
Ora, per ruotare la matrice dei coefficienti della spline (chimamola A) la devo moltiplicare per la amtrice di rotazione (R); Per calcolare A90° dovrei moltiplicare A*R; siccome le dimensioni sono divrese devo riempire la matrice R di rotazione con tanti 0, fino a quando le dimensioni sono uguali?
Ora, per ruotare la matrice dei coefficienti della spline (chimamola A) la devo moltiplicare per la amtrice di rotazione (R); Per calcolare A90° dovrei moltiplicare A*R; siccome le dimensioni sono divrese devo riempire la matrice R di rotazione con tanti 0, fino a quando le dimensioni sono uguali?
Secondo me con la matrice di rotazione vai un po' fuori strada..
Per iniziare a rispondere alla tua domanda, vorrei che mi spiegassi come è memorizzata la spline in matlab.
A quanto ne so, la spline è una funzione polinomiale a pezzi, quindi dovresti GIÀ avere l'espressione di quella curva in forma \(T = T(z)\) [cioè la temperatura in funzione della profondità]. Probabilmente in una delle variabili ci sono i coefficienti del polinomio di secondo grado che interpola la funzione in un intorno del punto [che ad occhio puoi stimare], quindi ti riduci a minimizzare un polinomio.
E questo è FACILE!
Per iniziare a rispondere alla tua domanda, vorrei che mi spiegassi come è memorizzata la spline in matlab.
A quanto ne so, la spline è una funzione polinomiale a pezzi, quindi dovresti GIÀ avere l'espressione di quella curva in forma \(T = T(z)\) [cioè la temperatura in funzione della profondità]. Probabilmente in una delle variabili ci sono i coefficienti del polinomio di secondo grado che interpola la funzione in un intorno del punto [che ad occhio puoi stimare], quindi ti riduci a minimizzare un polinomio.
E questo è FACILE!
"Raptorista":
Secondo me con la matrice di rotazione vai un po' fuori strada..
Per iniziare a rispondere alla tua domanda, vorrei che mi spiegassi come è memorizzata la spline in matlab.
Allora, con il comando:
pp=spline(prof,temp)
ho costruito la spline cubica naturale interpolante memorizzando le informazioni in una variabile strutturata detta pp-form.
In particolare le informazioni contenute sono:
Poi col comando:
[x,C,l,k]=unmkpp(pp);
ho estratto le componenti della pp-form, essendo x il vettore dei nodi, C la matrice dei coeffcienti di dimensione lxk la cui j-ma riga coincide con i coeffcienti del j-esimo polinomio che rappresenta la spline, ossia:
Per valutare una variabile nella pp-form bisogna utilizzare una opportuna routine, a cui passare un numero di valutazioni:
spder=ppval(ppder,zz);
A quanto ne so, la spline è una funzione polinomiale a pezzi, quindi dovresti GIÀ avere l'espressione di quella curva in forma \(T = T(z)\) [cioè la temperatura in funzione della profondità]. Probabilmente in una delle variabili ci sono i coefficienti del polinomio di secondo grado che interpola la funzione in un intorno del punto [che ad occhio puoi stimare], quindi ti riduci a minimizzare un polinomio.
E questo è FACILE!
I coefficienti del polinomio sono memorizzati nella matrice C, però non ho capito il procedimento che mi suggerisci.
Grazie
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo