Questo post prosegue sul tema del precedente: la determinazione delle linee di flusso da campi vettoriali in una modalità GIS-friendly. Software come Paraview consentono di effettuare questa determinazione, ma sia l'input sia l'output sono in formati non usuali per i GIS.
Queste tecniche numeriche rientrano nel tema dell'integrazione di ODE (Ordinary Differential Equations).
Le fonti bibliografiche per la creazione degli algoritmi sono Computational Physics di P.O.J. Scherer, Advanced engineering Mathematics di E. Kreyszig e A Primer on Scientific Programming with Python di H.P. Langtangen. Esistono inoltre molte fonti di informazione accessibili online (vedi p.e. Wikipedia).Nel precedente post ho presentato uno script Python che calcola le pathlines per campi vettoriali costanti nel tempo. Quello script preliminare si basava sul metodo più semplice fra quelli proposti in letteratura, cioè il metodo di Eulero esplicito (forward Euler), la cui formula è:
u k+1 = u k + Δt f( u k , t k) (eq. 9.23 in Langtangen, 2009)
u k+1 = u k + Δt f( u k , t k) (eq. 9.23 in Langtangen, 2009)
dove:
u k : valore della variabile allo step k,
u k+1 : il suo valore allo step successivo
Δt : incremento temporale considerato
f( u k , t k) : la funzione che ci fornisce la derivata di u k al tempo t k, nel caso specifico la velocità.
Essendo un campo vettoriale con componenti x e y, la formula soprastante verrà applicata sia all'una sia all'altra componente. La derivata viene stimata calcolando l'incremento della variabile durante l'intervallo di tempo considerato. Nell'algoritmo implementato, i valori delle componenti in un particolare punto vengono stimati per interpolazione bilineare.
I risultati ottenibili con questo metodo non sono stabili ed è preferibile utilizzare altri metodi, il più noto dei quali è il metodo di Runge-Kutta, con le sue varianti.
Prima di considerare le altre metodologie, è necessario stabilire come si possa confrontare la qualità dei risultati fra i diversi metodi. In Computational Physics è utilizzato come esempio quello di un campo vettoriale centrale (esempi sono il campo gravitazionale e quello elettrostatico) con orbite circolari (Fig. 1a). Un metodo di interpolazione di linee di flusso in un tale campo dovrebbe generare orbite strettamente circolari. Si può dimostrare sia teoricamente sia sperimentalmente che il metodo di Eulero esplicito (forward) non è stabile e determina delle orbite che man mano hanno un raggio sempre maggiore (Fig. 1b). Metodi con una maggiore stabilità sono l'improved Euler method (midpoint) e soprattutto i metodi di Runge-Kutta del terzo e quarto ordine ed il metodo di Runge-Kutta-Fehlberg (Fig. 1c). Questi producono risultati molto più stabili, in cui le orbite si degradano molto più lentamente con i vari cicli di iterazione.
 |
| Fig. 1. A) campo centrale con orbite circolari. B) L'orbita interpolata con il metodo di Eulero esplicito non mantiene un raggio costante ma cresce rapidamente con il tempo. C) Metodi come Runge-Kutta 4 e Runge-Kutta-Fehlberg sono stabili e le orbite determinate rimangono circolari al passare del tempo. Ridisegnato da Schroeder & Martin, fig. 1.17, in 'The visualization handbook', Academic Press. |
Dei metodi inseriti nella nuova versione dello script Python (mid-point, Runge-Kutta del terzo e quarto ordine, Runge-Kutta-Fehlberg) inseriamo le formule del più complesso e preciso fra questi, cioé il metodo di Runge-Kutta-Fehlberg (RKF) (riprese da Kreyszig, 2006):
yn+1 = yn + h ( γ1k1 + ... + γ6k6)
con vettore di coefficienti:
γ = [ 16/135 0 6656/12825 28561/56430 -9/50 2/55]
e stimatori:
k1 = f(xn, yn)
k2 = f(xn + (1/4) h, yn + h (1/4) k1)
k3 = f(xn + (3/8) h, yn + h ( (3/32) k1 + (9/32) k2) )
k4 = f(xn + (12/13) h, yn + h ( (1932/2197) k1 - (7200/2197) k2 + (7296/2197) k3) )
k5 = f(xn + h, yn + h ( (439/216) k1 - (8) k2 + (3680/513) k3 - (845/4104) k4) )
k6 = f(xn + (1/2) h, yn + h ( - (8/27) k1 + (2) k2 - (3544/2565) k3 + (1859/4104) k4 – (11/40) k5) )
Confrontando i risultati del metodo di Eulero esplicito con il Runge-Kutta-Fehlberg ottenuti per un campo con orbite circolari con la nuova versione dello script Python, si nota appunto l'instabilità del primo e la stabilità del secondo.
yn+1 = yn + h ( γ1k1 + ... + γ6k6)
con vettore di coefficienti:
γ = [ 16/135 0 6656/12825 28561/56430 -9/50 2/55]
e stimatori:
k1 = f(xn, yn)
k2 = f(xn + (1/4) h, yn + h (1/4) k1)
k3 = f(xn + (3/8) h, yn + h ( (3/32) k1 + (9/32) k2) )
k4 = f(xn + (12/13) h, yn + h ( (1932/2197) k1 - (7200/2197) k2 + (7296/2197) k3) )
k5 = f(xn + h, yn + h ( (439/216) k1 - (8) k2 + (3680/513) k3 - (845/4104) k4) )
k6 = f(xn + (1/2) h, yn + h ( - (8/27) k1 + (2) k2 - (3544/2565) k3 + (1859/4104) k4 – (11/40) k5) )
Confrontando i risultati del metodo di Eulero esplicito con il Runge-Kutta-Fehlberg ottenuti per un campo con orbite circolari con la nuova versione dello script Python, si nota appunto l'instabilità del primo e la stabilità del secondo.
 |
| Fig. 2. Confronto tra risultati del metodo di Runge-Kutta-Fehlberg e di Eulero esplicito su un punto iniziale (quadrato) in un campo centrale con orbite circolari: il circolo interno rappresenta vari cicli calcolati dal primo metodo, la spirale crescente i corrispondenti cicli prodotti dal secondo. |
Consideriamo ora l'applicazione di questi due metodi, agli estremi come precisione, su un dataset naturale, già analizzato in post precedenti: i flussi del ghiacciaio David in corrispondenza della sua zona di disancoraggio dal substrato per formare la lingua Drygalski (Antartide). Questo dataset deriva dalla tesi di dottorato in Scienze Polari di Debbie Biscaro: [http://www.mna.it/italiano/Didattica/dott_polare/Tesi_abstract/Biscaro_XXI_2010.pdf]
Il campo di velocità è stato dedotto confrontando due immagini Landsat fra loro coregistrate del 2001 e del 2003, con una separazione temporale di 400 giorni. Le velocità sono state determinate tramite il software open-source Imcorr: http://nsidc.org/data/velmap/imcorr.html .
Applichiamo lo script usando il metodo di Eulero esplicito e di Runge-Kutta-Fehlberg, confrontando i risultati fra di loro e rispetto alle lineazioni glaciali visibili dalle immagini satellitari Landsat.
Il prodotto tra i due metodi non differiscono in maniera sensibile. Alla scala della mappa sottostante sono indistinguibili (Fig. 3). Fra i punti terminali della traccia più settentrionale in Fig. 3, la differenza è di soli 9.3 m su un percorso totale di circa 15.3 km.
Queste basse differenze sono presumibilmente dovute alle traiettorie fondamentalmente rettilinee seguite dai flussi.
Il prodotto tra i due metodi non differiscono in maniera sensibile. Alla scala della mappa sottostante sono indistinguibili (Fig. 3). Fra i punti terminali della traccia più settentrionale in Fig. 3, la differenza è di soli 9.3 m su un percorso totale di circa 15.3 km.
Queste basse differenze sono presumibilmente dovute alle traiettorie fondamentalmente rettilinee seguite dai flussi.
 |
| Fig. 3. Mappa delle linee di flusso determinate con il metodo di Runge-Kutta-Fehlberg. A questa scala il risultato è indistinguibile da quello del metodo di Eulero esplicito. Sfondi: mosaico satellitare LANDSAT e grid della magnitudine delle velocità. Mappa creata con QGIS Copiaco. |
Quali sono le relazioni con le lineazioni di flusso visibili nel mosaico satellitare LANDSAT?
Confrontando le traiettorie da RKF con le lineazioni, si vedono chiaramente delle discordanze dove le traiettorie sono particolarmente variabili (Fig. 4). Questo potrebbe essere dovuto a varie cause: determinazione non ottimale del campo di velocità, inadeguatezza del metodo RKF, variazioni nel tempo del campo di velocità glaciale, tali per cui le streamlines determinate da RKF differiscono dalle pathlines che sono evidenti nelle immagini glaciali.
Non abbiamo indizi per discernere fra queste e altre ipotesi possibili. L'ultima ipotesi però ci sembra verosimile, in quanto è noto che i flussi glaciali possono variare su tempi anche abbastanza brevi ed inevitabilmente le lineazioni glaciali subiscono processi di deformazione che possono farle divergere dalle streamlines "istantanee". Metodi di analisi della deformazione glaciale potrebbero aiutare a valutare meglio questa possibilità.
Confrontando le traiettorie da RKF con le lineazioni, si vedono chiaramente delle discordanze dove le traiettorie sono particolarmente variabili (Fig. 4). Questo potrebbe essere dovuto a varie cause: determinazione non ottimale del campo di velocità, inadeguatezza del metodo RKF, variazioni nel tempo del campo di velocità glaciale, tali per cui le streamlines determinate da RKF differiscono dalle pathlines che sono evidenti nelle immagini glaciali.
Non abbiamo indizi per discernere fra queste e altre ipotesi possibili. L'ultima ipotesi però ci sembra verosimile, in quanto è noto che i flussi glaciali possono variare su tempi anche abbastanza brevi ed inevitabilmente le lineazioni glaciali subiscono processi di deformazione che possono farle divergere dalle streamlines "istantanee". Metodi di analisi della deformazione glaciale potrebbero aiutare a valutare meglio questa possibilità.
 |
| Fig. 4. Mappa delle traiettorie determinate con il metodo di Runge-Kutta-Fehlberg confrontabili con le lineazioni di flusso glaciale visibili nel mosaico stellitare LANDSAT. Mappa creata con QGIS Copiaco. |
Lo script Python è ancora in fase di test e inserimento di ulteriori funzionalità. Chi fosse interessato, può richiedere la versione attuale a alberti.m65@gmail.com.
No comments:
Post a Comment