[Mat Numerica] [RISOLTO] Risoluzione esercizio con Matlab

terminal1
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.

Risposte
Raptorista1
"terminal":

I coefficienti del polinomio sono memorizzati nella matrice C, però non ho capito il procedimento che mi suggerisci.
Grazie

Non hai cpaito perché non c'è nulla da capire :)
Hai un polinomio di terzo grado, devi trovare il massimo. Punto.
Lo puoi fare pure a mano, come si fa al liceo.

terminal1
E in Matlab come potrei esprimerlo?

Raptorista1
Ma non è che se non usi Matlab muori :)
Anzi, se usi matlab otterrai un risultato approssimato, quindi peggiore!

In ogni caso, se trovi il massimo di un polinomio di terzo grado lasciando indicate tutte le lettere, poi arrivi ad una qualche formula che puoi dare direttamente in pasto a matlab, forse XD

terminal1
Hai frainteso... il "problema" era come tradurlo in codice...

Raptorista1
Non penso di aver frainteso, ma forse non riesco a spiegarmi bene: non è necessario usare matlab per risolvere il problema.

Anche ammesso che tu voglia usare matlab, fargli risolvere numericamente l'equazione è sconveniente: tanto vale che tu faccia il conto simbolicamente a mano e trovi la formuletta generale che poi puoi dargli in pasto.

terminal1
Il probelma è appunto risolverlo in ambiente Matlab. Nel mio caso, sò che per calcolare la quota devo:
1 - o trovare il punto di flesso della curva temperatura-profondità
2 - o il punto in cui si annulla la derivata seconda della temperatura in funzione della profondità
3 - o il punto in cui è massima la derivata prima
Finq qui non ci piove.


Il mio problema è come assegnare i coefficienti di ogni singolo polinomio (memorizzati nella matrice) nella formula risolutiva...
$ (-b \pm sqrt(b^{2} - 4ac )) / (2a) $

Raptorista1
Quella non è la formula giusta!
Se hai una funzione polinomiale del tipo \(f(x) = ax^3 + bx^2 + cx + d\) e vuoi trovare il punto di flesso, dovresti azzerare la derivata seconda, non la derivata prima!

P.s. Nella formula quadratica che hai scritto manca comunque un segno meno.

terminal1
Io ho pensato di calcolare il punto di massimo della derivata I, i cui polinomi sono di 2° grado.
Ho escluso i punti dove si annulla della derivata II, perchè individuano la fascia del termocline e non la quota esatta.
Per la formula..si, mi è sfuggito un - ... sistemato ;)

Raptorista1
Guarda che nei punti di massimo della derivata prima [in questo caso in cui è continua] necessariamente si annulla la derivata seconda. Da questo non si scappa!

terminal1
:oops: hai ragione.. ci ho ragionato bene.
Mi sono lasciato "imbrogliare" dal grafico.

Faccio il punto della situazione:
nella mtrice C ci sono memorizzati i coefficienti dei polinomi della Spline; devo trovare il massimo della derivata della spline, i cui coefficienti dei polinomi sono memorizzati nella matrice Cder.

Le mie funzioni di partenza, quindi, si trovano nelle righe di Cder:
$ f ( x_1 ) = a_1x^2 + b_1x + c_1 $
$ f ( x_2 ) = a_2x^2 + b_2x + c_2 $
...
...
$ f ( x_7 ) = a_7x^2 + b_7x + c_7 $


Devo calcolare i punti di eventuale massimo (che mi interessa), minimo e felsso;

1. Calcolo la derivata della funzione $ f ( x_1 ) = a_1x^2 + b_1x + c_1 $ :
$ f ' ( x_1 ) = a_1x + b_1 $

2. Pongo la derivata $ f ' ( x_1 ) = 0 $

3. Calcolo le soluzioni dell'equazione;

4. Sostituisco il valore ottenuto in $ x_1 $ nella funzione di partenza $ f ( x_1 ) = a_1x^2 + b_1x + c_1 $ e il valore che ottengo è la componente y del punto;

5. La coppia ( x, y ) è un punto di eventuale massimo, minimo o flesso;

6. Trovo la derivata seconda:
$ f ' ' ( x_1 ) = a_1 $

7. Ora dovrei calcolare il valore della derivata II° sostituendo in $ f ' ' ( x_1 ) = a_1 $ la $ x_1 $ con il valore ottenuto al punto 3, e poi vado a vedere se $ f ' ' ( x_1 ) $ è < o > o = a 0

Nello schema che ho riportato ho considerato solo la funzione f relativa alla prima riga (o primo polinomio), quindi tutto ciò va applicato per tutte e 7 le $ f ( x_i )$.

Il ragionamento è corretto dal punto di vista matematico?

Grazie

Raptorista1
Secondo me stai mischiando un po' di cose... Vediamo di fare ordine!

Stai usando un'interpolazione di tipo spline con una spline cubica \(s_3\). Questo significa che in ciascuno degli intervalli la tua funzione approssimante è un polinomio di terzo grado, e come tale è individuato da 4 [quattro] coefficienti, che sono memorizzati [a quel che ho capito dal tuo primo post] nella matrice C, che quindi dovrebbe essere di dimensioni \(n \times 4\) dove \(n\) è il numero di intervalli di approssimazione.
"Scopo del gioco" è allora trovare i flessi di questi polinomi di terzo grado, ok?
Per fare ciò, presa la funzione polinomiale dell'i-esimo intervallo \([x_i, x_{i+1}]\)
\[
f_i(x) = C_{i1}x^3 + C_{i2}x^2 + C_{i3}x + C_{i4}, \qquad x \in [x_i, x_{i+1}]
\]
devi risolvere
\[
f_i''(x) = 0 : 6 C_{i1}x + 2C_{i2} = 0, \qquad x \in [x_i, x_{i+1}].
\]

Non è difficile, e non sono in grado di spiegartelo in maniera più semplice di così.

Una piccola nota sulla validità fisica e numerica di quello che stai facendo, indipendente da quali siano le tue conoscenze di matematica [a proposito, puoi dirmi a che corso di studi sei iscritto?].
Introduco la mia nota riassumendo quanto ho capito [correggimi dove sbaglio]: tu parti da una funzione \(T = T(z)\) [temperatura in funzione della profondità] di cui a priori non sai nulla. Allora la valuti in alcuni punti \(z_i\) per avere un'idea di come sia fatta; successivamente, ti chiedi se e dove ci siano punti di flesso \(\tilde z_j\) nella tua curva \(T = T(z)\).
Prosegui poi approssimando la funzione con una spline cubica [che è, sostanzialmente, un patchwork di polinomi] e trovi i flessi \(\bar z_j\) di questa spline cubica.
Va da sé, ma ci tengo a ribadirlo, che tu non hai alcuna garanzia, nella maniera più assoluta, che gli \(\tilde z_j\) siano legati in qualche modo agli \(\bar z_j\). È chiaro questo concetto?
Gli \(\tilde z_j\) potrebbero esserci o non esserci, essere in numero ed in valore diversi dagli \(\bar z_j\), ed a priori non hai alcuna informazione su questa cosa, non la puoi controllare.

Il mio monito è quindi questo: non so se questo esercizio ti serva solo per imparare ad usare le spline o se sia parte di qualcosa di più grosso, non andare però ciecamente a dire che il termoclino si trova dove ci sono i flessi della spline.

terminal1
"Raptorista":
Secondo me stai mischiando un po' di cose... Vediamo di fare ordine!

Stai usando un'interpolazione di tipo spline con una spline cubica \(s_3\). Questo significa che in ciascuno degli intervalli la tua funzione approssimante è un polinomio di terzo grado, e come tale è individuato da 4 [quattro] coefficienti, che sono memorizzati [a quel che ho capito dal tuo primo post] nella matrice C, che quindi dovrebbe essere di dimensioni \(n \times 4\) dove \(n\) è il numero di intervalli di approssimazione.

Perfetto, precisamente n è il numero delle mie misurazioni, quindi C = 7 x 4 .

"Scopo del gioco" è allora trovare i flessi di questi polinomi di terzo grado, ok?

Perfetto, trovare il flesso della curva temperatura-profondità (rappresntata dai polinomi di terzo grado) che equivale a trovare:
1 - o il punto di flesso della curva temperatura-profondità;
2 - o il punto in cui si annulla la derivata seconda della temperatura in funzione della profondità
3 - o il punto in cui è massima la derivata prima

Per fare ciò, presa la funzione polinomiale dell'i-esimo intervallo \([x_i, x_{i+1}]\)
\[
f_i(x) = C_{i1}x^3 + C_{i2}x^2 + C_{i3}x + C_{i4}, \qquad x \in [x_i, x_{i+1}]
\]
devi risolvere
\[
f_i''(x) = 0 : 6 C_{i1}x + 2C_{i2} = 0, \qquad x \in [x_i, x_{i+1}].
\]

Non è difficile, e non sono in grado di spiegartelo in maniera più semplice di così.

Qui forse non ci siamo capiti: tu vuoi calcolare il flesso della curva temperatura-profondità, mentre io volevo calcolare il punto in cui è massima la derivata prima; che sarebbe la stessa cosa, non credo di sbagliare.

Una piccola nota sulla validità fisica e numerica di quello che stai facendo, indipendente da quali siano le tue conoscenze di matematica [a proposito, puoi dirmi a che corso di studi sei iscritto?].
Introduco la mia nota riassumendo quanto ho capito [correggimi dove sbaglio]: tu parti da una funzione \(T = T(z)\) [temperatura in funzione della profondità] di cui a priori non sai nulla. Allora la valuti in alcuni punti \(z_i\) per avere un'idea di come sia fatta; successivamente, ti chiedi se e dove ci siano punti di flesso \(\tilde z_j\) nella tua curva \(T = T(z)\).

Ok

Prosegui poi approssimando la funzione con una spline cubica [che è, sostanzialmente, un patchwork di polinomi] e trovi i flessi \(\bar z_j\) di questa spline cubica.


Ok


Va da sé, ma ci tengo a ribadirlo, che tu non hai alcuna garanzia, nella maniera più assoluta, che gli \(\tilde z_j\) siano legati in qualche modo agli \(\bar z_j\). È chiaro questo concetto?
Gli \(\tilde z_j\) potrebbero esserci o non esserci, essere in numero ed in valore diversi dagli \(\bar z_j\), ed a priori non hai alcuna informazione su questa cosa, non la puoi controllare.


Qui non ti seguo, \(\bar z_j\) sono i flessi mentre \(\tilde z_j\)?
Puoi spiegarmi meglio?


Il mio monito è quindi questo: non so se questo esercizio ti serva solo per imparare ad usare le spline o se sia parte di qualcosa di più grosso, non andare però ciecamente a dire che il termoclino si trova dove ci sono i flessi della spline.

E' ovvio che le spline non sono una cosa a se stante, non possono essere trattate senza un background; è il testo dell'esercizio che indica come soluzione il calcolo del flesso della spline.
Comunque ti cito lo scopo dell'esercizio:
L’esercizio consiste nel costruire innanzitutto il grafico riportato di seguito e trovare grazie, eventualmente, a grafici più dettagliati e/o a tabelle la quota del termocline.


Ti ringrazio davvero per la disponibilità che dimostri.
:smt023

P.s. seguo un corso di laurea in informatica

Raptorista1
"terminal":

Qui forse non ci siamo capiti: tu vuoi calcolare il flesso della curva temperatura-profondità, mentre io volevo calcolare il punto in cui è massima la derivata prima; che sarebbe la stessa cosa, non credo di sbagliare.

Sono la stessa cosa, ma alla fine l'unico modo operativo di procedere è annullare la derivata seconda e verificare col segno che sia un flesso.

"terminal":

Qui non ti seguo, \(\bar z_j\) sono i flessi mentre \(\tilde z_j\)?
Puoi spiegarmi meglio?

Gli \(\tilde z_j\) sono i flessi della funzione \(T = T(z)\); gli \(\bar z_j\) sono i flessi della spline.

"terminal":

Comunque ti cito lo scopo dell'esercizio:
L’esercizio consiste nel costruire innanzitutto il grafico riportato di seguito e trovare grazie, eventualmente, a grafici più dettagliati e/o a tabelle la quota del termocline.


CVD: attento a fare deduzioni :)

terminal1
Aggiorno la situazione, con i passaggi che faccio:
- ho una curva interpolante, una spline cubica, che interpola una serie di polinomi di 3° grado, ossia $ f(x)=ax^3+bx^2+cx+d $
- devo calcolare il felsso della curva
- allora calcolo la derivata II dei polinomi della spline
- impongo la derivata II = 0 , e ricavo la x per ciascun polinomio
- sostituisco le x nel rispettivo polinomio iniziale di 3° grado, risolvo e trovo la y per ciascun polinomio

Di questi punti, nessuno appartiene alla spline... come mai?
Il ragionamento non sembra essere sbagliato...

:''''(

Raptorista1
"terminal":

Di questi punti, nessuno appartiene alla spline... come mai?
Il ragionamento non sembra essere sbagliato...
:''''(

Ecco, questo è abbastanza strano, perché la spline cubica è una curva globalmente \(C^2([a,b])\), quindi un flesso da qualche parte DEVE esserci [lo dico perché il grafico del primo post suggerisce una certa forma della spline].

Vediamo se ho capito: stai dicendo che se prendi la funzione \(i\)-esima \(f_i(x) = a_i x^3 + b_i x^2 + c_i x + d_i\) nell'intervallo \([x_i,x_{i+1}]\), questa ha sempre [i.e. indipendentemente da \(i\)] un flesso che non appartiene all'intervallo \([x_i,x_{i+1}]\)?

terminal1
Sicuramente commetto qualche errore nella valutazione della funzione nell'intervallo...perchè, come dici te, un flesso deve esserci, e (sembra) essere compreso nel 4 tratto / polinomio-interpolante.

Con:
x1rad, x2rad, ..., x7rad
indico le radici dei polinomi di primo grado (contenuti nella derivata II della spline, ax + b )

Poi sostituisco le x calcolate nelle rispettive funzioni $ f_i(x)=a_ix^3+b_ix^2+c_ix+di $
valutando però in $(x_i,x_(i+1))$ ossia nell'intervallo del tratto di polinomio.


Somo certo di sbagliare, e quindi , come si valuta un polinomio in un intervallo?

Raptorista1
I polinomi si valutano in quel modo, esattamente come ogni funzione XD
Fai così: disegna su un piano le derivate seconde \(f_i''(x)\) che ottieni negli intervalli \([x_i, x_{i+1}]\) e guarda se almeno graficamente uno zero c'è.
Indicativamente, se i conti sono giusti dovresti ottenere il terzo grafico del tuo secondo post [quello con una linea spezzata] però "orientato in orizzontale" [i.e. sulle ascisse metti la profondità, sulle ordinate metti la temperatura].

Se già questo non ti esce, allora c'è un errore più a monte.

terminal1
Questo è il grafico delle derivate seconde dei polinomi:


Devo valutare ogni polinomio nel suo intervallo:
ad es. il primo polinomio (o la prima riga della matrice dei coefficienti della derivata II)
Cder2(1,:)
nell'intervallo $(x_1, x_2)$
che corrisponde ai punti $x_1 = 0$ e $x_2=2.3$

$x=[0 , 2.3 , 4.9 , 9.1 , 13.7 , 18.3 , 22.9 , 27.2];$

Giusto?

Raptorista1
Bene, la curva è esattamente quella del primo post, ed uno zero si vede ad occhio.
A questo punto non serve valutare ogni polinomio: prendi solo quello in cui vedi che c'è uno zero e risolvi l'equazione di primo grado che ottieni.

terminal1
Ok, allora prendo il polinomio che "contiene" lo 0, lo impongo = 0 e lo risolvo;
La x che ottengo è la componente x del punto di flesso; corretto?
Sostituendo la x calcolata, nel polinomio di 3° grado della spline (la funzione di partenza) trovo la y;
Ora questo punto non ricade sulla retta:



In cosa sbaglio?

Grazie infinite per la pazienza :smt023

PS. ...
Ma la x che calcolo dalla (derivata seconda = 0) di un polinomio, non deve essere interna all'intervallo di interpolazione?
ad es. se x1 (la x del primo polinomio) non deve essere compresa tra 0 e 2.3 (considerando $x=[0 , 2.3 , 4.9 , 9.1 , 13.7 , 18.3 , 22.9 , 27.2]$ )?

Rispondi
Per rispondere a questa discussione devi prima effettuare il login.