ElectroYou - la comunità dei professionisti del mondo elettrico

da **onire** » 5 mar 2019, 13:11

PietroBaima ha scritto:In pratica bisogna usare questo algoritmo: ......

Ho provato ad usare il tuo algoritmo, ed in effetti per gamma<0.5 funziona molto bene. Per gamma>0.5 mi restituisce un valore maggiore del diametro.

Non avendo capito come hai ottenuto la formula, ma piacerebbe capirlo, non sono stato in grado di correggerla.

da **onire** » 5 mar 2019, 13:32

Ho pensato, che potrei sottrarre al volume del serbatoio, il volume del liquido e quindi ricondurlo al caso 1.
Calcolata l'altezza, la sottraggo al diametro.

Ho verificato e mi sembra funzionare bene, ma come scritto prima mi piacerebbe capire meglio quanto scritto da

PietroBaima.

da **EcoTan** » 5 mar 2019, 13:43

LucaCassioli ha scritto:Riempilo con X litri e misura l'altezza del liquido

In effetti la taratura di un serbatoio fiscale l'abbiamo fatta così ma non è questo il caso.

da **PietroBaima** » 5 mar 2019, 16:54

onire ha scritto:
PietroBaima ha scritto:In pratica bisogna usare questo algoritmo: ......

Ho provato ad usare il tuo algoritmo, ed in effetti per gamma<0.5 funziona molto bene. Per gamma>0.5 mi restituisce un valore maggiore del diametro.

Non avendo capito come hai ottenuto la formula, ma piacerebbe capirlo, non sono stato in grado di correggerla.

hai sbagliato un segno in alpha2, quello del termine centrale.
Stasera scrivo due righe sul metodo.
Ciao
Pietro

da **onire** » 5 mar 2019, 17:09

Hai ragione sul segno :ok:

Script modificato.

Aspetto chiarimenti sul metodo.

da **PietroBaima** » 6 mar 2019, 0:07

Dato un serbatoio cilindrico di lunghezza L, raggio di base R, avendolo riposto a terra in modo che sia appoggiato solidalmente al terreno lungo la sua lunghezza e avendo inserito un volume di liquido al suo interno V, trovare una formula computazionalmente semplice che esprima l'altezza h raggiunta dal liquido nel serbatoio.

Come prima cosa conviene osservare che è molto più semplice trovare una formula che esprima il volume del liquido rispetto all'altezza, piuttosto che la funzione inversa, per cui troviamo il volume data l'altezza e poi invertiamo la formula.
Scrivo così perché per trovare il volume è sufficiente integrare la regione di area occupata dal liquido sulla base circolare e poi moltiplicarla per la lunghezza del cilindro, invece per trovare l'altezza dato il volume bisognerebbe calcolare l'integrale curvilineo della regione della calotta. Meglio essere pigri.

Se piazzo il cerchio con centro O tale che O=(0,0) e esprimo il versore della gravità concorde a x ottengo che il liquido occuperà questa sezione:

Sorgente FidoCadJ | Ingrandisci | Cos'è Fidocad

Se la formula del cerchio è x^2+y^2=R^2

il centro del cerchio è lo zero, quindi il punto di contatto con il terreno sarà x=R mentre la sommità della cisterna sarà x=-R.
Quindi se la cisterna è vuota (h=0) dovrò partire ad integrare da x=R e finire a terra, cioè a x=R. Ovviamente il punto di inizio e fine coincidono e l'area sarà nulla.
Se la cisterna è colma (h=2R) dovrò partire a integrare da x=-R (ricordo che lo zero coincide con il centro del cerchio della base della cisterna) e finire a x=R.
Quindi il punto generico, rispetto a questa convenzione, dell'altezza del liquido è R-h.
Devo quindi impostare questo integrale per calcolare l'area:

$A=2\ \int_{R-h}^{R}\sqrt{R^2-x^2}\text{d}x$

L'area è due volte quell'integrale, perché ricavando y dalla formula x^2+y^2=R^2

e prendendo solo la radice positiva calcolo solo metà area.

Sorgente FidoCadJ | Ingrandisci | Cos'è Fidocad

Ho calcolato l'integrale con Mathematica (sì, sono pigro, ma se volete calcolarlo a mano potete sostituire $z=R \cos x$ e ricondurlo all'integrale del seno al quadrato).

Codice: Seleziona tutto: Assuming[h > 0 && h <2* R && R > 0, Integrate[2*Sqrt[R^2 - x^2], {x, R - h, R}]]

Il comando qui sopra serve per far integrare la funzione facendo assumere a Mathematica che h sia positivo e minore del diametro della cisterna.
Avendo inoltre definito delle relazioni fra numeri Mathematica assume automaticamente che h ed R non siano complessi.
Mathematica risponde:

Codice: Seleziona tutto: ConditionalExpression[ R (h Sqrt[-((h (h - 2 R))/R^2)] + \[Pi] R - Sqrt[h (-h + 2 R)] - R ArcCos[-1 + h/R]), h < R]

In pratica ci risponde che ha trovato una formula per h<R, cioè che sappiamo come calcolare riempimenti fino a metà cisterna. Yuppie!

Non dobbiamo dimenticare di moltiplicare la formula per L, la lunghezza della cisterna, per ottenere il volume di liquido in essa contenuto.

La formula adesso deve essere invertita.
Questa formula ci permette di trovare il volume in funzione dell'altezza, noi dobbiamo trovare invece la funzione inversa.

Prima però facciamo una cosa importante.
Tutte le volte che vogliamo invertire (in questo caso numericamente, la formula è troppo complessa per cercare una forma chiusa) e approssimare DOBBIAMO avere a che fare con numeri puri.
Questa formula restituisce un volume, avendo in ingresso delle lunghezze (h ed R).
La normalizzazione più logica è quella che esprime l'altezza h in funzione del raggio R e il volume del liquido nel serbatoio in funzione del volume del serbatoio.
Abbiamo:

$V=R\ L \left(h \sqrt{-\frac{h (h-2 R)}{R^2}}-\sqrt{h (2 R-h)}-R \cos ^{-1}\left(\frac{h}{R}-1\right)+\pi R\right)$

Sapendo che il volume della cisterna è $V_{cist}=\pi R^2 L$ e chiamando $\alpha=h/R$ , con un po' di algebra, possiamo scrivere il rapporto fra il volume della cisterna e il volume V (che chiamo gamma) come:

$\gamma(\alpha)=\frac{2 (2 \alpha-1) \sqrt{\alpha (1-\alpha)}}{\pi }+\frac{\cos ^{-1}(1-2 \alpha)}{\pi }$

Quindi posso definire su Mathematica questa funzione:

Codice: Seleziona tutto: G[a_] := 2/Pi*Sqrt[a*(1 - a)]*(2*a - 1) + 1/Pi*ArcCos[1 - 2*a]

Che non posso invertire in forma chiusa, ma numericamente sì, col comando:

Codice: Seleziona tutto: a[G_] := InverseFunction[G][a]

Adesso la disegno col comando Plot e la guardo:

Codice: Seleziona tutto: Plot[a[x], {x, 0, 1}]

: alphadigamma.jpg (7.12 KiB) Osservato 11365 volte

Mi rendo conto che assomiglia ad una radice quadrata, che diventa lineare crescendo, che mi fa capire che a quella radice quadrata è stata sommata una retta.
Compilo quindi una tabella di valori con la funzione inversa (ricordo che la funzione inversa è stata ottenuta in via numerica, quindi non ho una formula) con:

Codice: Seleziona tutto: data = Table[{x, G[x]}, {x, 0, 0.5, 0.001}]

che produce una matrice di 500 punti x,G(x).
Ovviamente faccio calcolare i punti fino a gamma=0.5, perché successivamente la funzione dovrà essere simmetrica e posso poi scegliere di ricavarla per simmetria, senza sprecare risorse di calcolo.

Cerco quindi il best fitting della matrice con una radice, una retta e una costante (conviene sempre aggiungere una piccola costante, o meglio, se si usa anche un polinomio, fare il fitting sul polinomio a cui non sono stati tolti dei termini. In genere si ottengono risultati migliori).
Uso il comando:

Codice: Seleziona tutto: Fit[data, {1, x, Sqrt[x]}, x]

che mi restituisce:

Codice: Seleziona tutto: -0.0072414 + 0.355769 Sqrt[x] + 0.509338 x

Ho trovato quindi una interpolante. La disegno sopra la funzione invertita numericamente:

In rosso ho la funzione calcolata numericamente e in blu la mia interpolante.
Direi che fino a gamma=0.5 posso essere soddisfatto.
Dopo non ho dato istruzioni al modello, quindi non posso pretendere che sia fedele alla funzione.

Correggiamo quindi questo problema.
Se osservo nuovamente la funzione invertita numericamente:

: alphadigamma.jpg (7.12 KiB) Osservato 11365 volte

Posso osservare che da gamma=0.5 fino a gamma=1 la funzione (ovviamente) gode di simmetria rispetto all'intervallo in cui gamma va da 0 a 0.5.

Quindi è vera la relazione di simmetria $\gamma(\alpha)=1-\gamma(1-\alpha)$

La applico su:
$\gamma(\alpha)=-0.0072414 + 0.355769 \sqrt{\alpha} + 0.509338 \alpha$

ottenendo

$1-\gamma(1-\alpha)=1+0.0072414 - 0.355769\sqrt{1-\alpha} - 0.509338 (1-\alpha)$

cioè

$1-\gamma(1-\alpha)=0.4979034 - 0.355769 \sqrt{1-\alpha} + 0.509338 \alpha$

Ho quindi trovato un modello composto da due sottomodelli.
Valutiamo l'errore del modello, utilizzando un modello campione dato dalla inversione numerica. Faccio disegnare l'errore relativo fra modelli:

: errore.jpg (6.7 KiB) Osservato 11365 volte

Ottenendo un errore massimo di circa sei parti su mille. Direi sia accettabile.

Ecco come ho fatto questo semplice modello.
Spero possa tornarti utile

La ragione per la quale reputo sbagliato l'uso di un metodo iterativo è che si ucciderebbe una mosca con un cannone. Meglio sbagliare in modo accettabile, trovando una soluzione semplice, piuttosto che cercare una soluzione precisissima ma complicata che sprema le risorse di calcolo all'osso.

Ciao,
Pietro.

da **xyz** » 6 mar 2019, 19:42

Ho preparato una piccola guida per integrare all'ottimo lavoro fatto da

PietroBaima su come ha trovato una soluzione al problema in questione. Eviterò di ripetere cose già dette e mi concentrerò sulle differenze e estensioni delle possibili soluzioni.

Il software utilizzata questa volta è R:

http://www.r-project.org

sempre open source multi piattaforma. Nato come software per elaborare dati statistici, ha una serie di funzioni e molte librerie per eseguire le più complesse elaborazioni su molti tipi di dati.

Lo script per R è questo, basta leggere la documentazione su come eseguirlo (esistono più metodi):

Codice: Seleziona tutto: tank_volume <- function(height, diameter, volume_total) { radius <- diameter / 2.0 length <- volume_total / (pi * radius^2) length * (radius^2 * acos((radius - height) / radius) - (radius - height) * sqrt(2 * radius * height - height^2)) } DIAMETER <- 2 VOLUME_TOTAL <- 5 N <- 25 tank <- data.frame( height = seq(from = 0, to = DIAMETER, by = DIAMETER / (N - 1)) ) tank$volume <- tank_volume(tank$height, DIAMETER, VOLUME_TOTAL) # RSE: 0.04521 #model <- as.formula(height ~ a * volume + b) #init <- list(a = 1, b = 1) # RSE: 0.01697 model <- as.formula(height ~ a * volume^3 + b * volume^2 + c * volume + d) init <- list(a = 1, b = 1, c = 1, d = 1) # RSE: 0.009457 #model <- as.formula(height ~ a * volume^5 + b * volume^4 + c * volume^3 + d * volume^2 + e * volume + f) #init <- list(a = 1, b = 1, c = 1, d = 1, e = 1, f = 1) # RSE: 0.04171 #model <- as.formula(height ~ a * sqrt(volume) + b * volume + c) #init <- list(a = 1, b = 1, c = 1) # RSE: 0.04457 #model <- as.formula(height ~ a * volume^b + c ) #init <- list(a = 1, b = 1, c = 1) # RSE: 0.02036 #model <- as.formula(height ~ a * cos(volume) + b * sin(volume) + c * volume + d) #init <- list(a = 1, b = 1, c = 1, d = 1) fit <- nls(model, data = tank, start = init, control = nls.control(maxiter = 100, minFactor = 1E-6)) print(summary(fit)) print(coef(fit)) #svg(filename="plot.svg") #png(filename="plot.png") plot(tank$volume, tank$height, xlab = "volume", ylab = "height", main = toString(fit$m$formula()), cex.main=0.9, pch = 16, col = "blue", panel.first = grid()) lines(tank$volume, fitted(fit), col = "red", lwd = 2) #dev.off()

La funzione "tank_volume" calcola il volume di un liquido in un serbatoio cilindrico rovesciato di diametro e volume totale noto. E' come quella di

PietroBaima ho solo usato il diametro al posto del raggio. Se si vuole si può inserire una funzione normalizzata in base a qualche dato, in generale non ci sono problemi dell'uso di costanti con unità di misura, io ho fissato nel codice il diametro e il volume totale.

Il programma campiona la funzione "tank_volume" per un numero prefissato di punti "N". Io ho scelto 25, è un valore sufficiente per questo problema, se si vuole una maggiore precisione basta aumentarlo.

Per trovare la funzione inversa di "tank_volume" si è sfruttato una regrezione non lineare basata su questo comando "Nonlinear Least Squares":

https://stat.ethz.ch/R-manual/R-devel/l ... l/nls.html

è incluso direttamente in R. Esistono delle librerie per R più sofisticate per dati più complessi ma per questo caso va più che bene. Bisogna fornire il modello della funzione da trovare. In questo caso si tratta di trovare l'altezza in base al volume (funzione inversa di "tank_volume"), vengono usati dei coefficienti non noti ma bisogna fornire un loro valore iniziale, se si riesce a trovare una soluzione il programma riporta quali sono i migliori coefficienti che approssimano i punti della funzione inversa (invertendo i campioni) col minimo errore.

Come primo nodello uso quello di una retta:

Codice: Seleziona tutto: model <- as.formula(height ~ a * volume + b)

Nel linguaggio di R questo è equivalente a trovare i migliori coefficienti "a" e "b" di questa funzione:

$\mbox{height}(\mbox{volume} ) = a \cdot \mbox{volume} + b$

Questo è il risulto:

Codice: Seleziona tutto: Formula: height ~ a * volume + b Parameters: Estimate Std. Error t value Pr(>|t|) a 0.358678 0.005411 66.283 < 2e-16 *** b 0.103305 0.016272 6.349 1.77e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.04521 on 23 degrees of freedom Number of iterations to convergence: 1 Achieved convergence tolerance: 5.393e-08 a b 0.3586779 0.1033054

Sono stati stampati i valori di "a" e "b" e viene fornito un un valore Residual standard error (radice quadrata dell'errore quadratico medio) per valutare la qualità dei parametri trovati.

Il grafico è il seguente:

: a * volume + b; plot-1.png (17.79 KiB) Osservato 11273 volte

E' come avere risolto il problema dei minimi quadrati, la rette può essere una buona approssimazione ma non per i valori iniziali e finali.

Usiamo il modello usato da

PietroBaima, una funzione non lineare con la radice quadrata:

Codice: Seleziona tutto: Formula: height ~ a * sqrt(volume) + b * volume + c Parameters: Estimate Std. Error t value Pr(>|t|) a 0.12048 0.05378 2.24 0.0355 * b 0.31314 0.02093 14.96 5.18e-13 *** c 0.04349 0.03063 1.42 0.1697 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.04171 on 22 degrees of freedom Number of iterations to convergence: 1 Achieved convergence tolerance: 1.71e-07 a b c 0.12048449 0.31314431 0.04349396

: a * sqrt(volume) + b * volume + c; plot-2.png (18.46 KiB) Osservato 11273 volte

Vediamo che questa volta i valori iniziali vengono approssimati meglio ma non quelli finali. Si potrebbe specchiare e dividerla in due come ha fatto

PietroBaima ma vedendo l'andamento dei punti una funziona polinomiale del terzo grado sicuramente riesce a coprire meglio tutti i punti:

Codice: Seleziona tutto: Formula: height ~ a * volume^3 + b * volume^2 + c * volume + d Parameters: Estimate Std. Error t value Pr(>|t|) a 0.014681 0.001231 11.928 8.14e-11 *** b -0.110107 0.009363 -11.760 1.06e-10 *** c 0.568164 0.019331 29.392 < 2e-16 *** d 0.038369 0.009806 3.913 8e-04 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.01697 on 21 degrees of freedom Number of iterations to convergence: 2 Achieved convergence tolerance: 1.483e-07 a b c d 0.01468089 -0.11010670 0.56816354 0.03836909

: a * volume^3 + b * volume^2 + c * volume + d; plot-3.png (18.59 KiB) Osservato 11273 volte

Questa volta i valori finali vengono mappati meglio.

Possiamo usare anche funzioni non polinomiali come usare quelle trigonometriche, questo è un altro modello:

Codice: Seleziona tutto: Formula: height ~ a * cos(volume) + b * sin(volume) + c * volume + d Parameters: Estimate Std. Error t value Pr(>|t|) a 0.075725 0.009388 8.066 7.22e-08 *** b 0.101369 0.011215 9.039 1.10e-08 *** c 0.411895 0.006049 68.089 < 2e-16 *** d -0.029738 0.015665 -1.898 0.0715 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.02036 on 21 degrees of freedom Number of iterations to convergence: 1 Achieved convergence tolerance: 4.381e-07 a b c d 0.07572475 0.10136880 0.41189512 -0.02973778

: a * cos(volume) + b * sin(volume) + c * volume + d; plot-4.png (18.68 KiB) Osservato 11273 volte

Non ho tanto tempo per provare altre funzioni, se qualcuno ha in mente modelli migliori basta indicarli e vedo se si ottiene una soluzione.

Credo che per ora il modello con la funzione polinomiale cubica riesca a coprire meglio tutti i valori.

Questa tecnica può essere usata in altri ambiti basta modificare la funzione da campionare e scrivere il modello della funzione con i coefficienti da trovare.