Ciao, come scritto nel precedente post, condivido un tutorial scritto dai prof. L.A. Barba, G.F. Forsyth, C. Cooper che descrive come sviluppare una simulazione dell’equazione di diffusione del calore tramite python. Al seguente link potete trovare il tutorial in originale. Il tutorial si suddivide in due parti, la prima descrive come passare dall’equazione di diffusione in forma differenziale ad una equazione discreta ovvero in cui le variazioni delle variabili in gioco (lo spazio x ed il tempo t) sono finite: una equazione in forma discreta ci permette di capire come cambia il profilo di temperatura u(x,t) considerando variazioni del tempo o dello spazio finite ad esempio di un 1 secondo oppure variazioni dello spazio di 1 cm.

La seconda parte del tutorial (descritta nel prossimo post)  utilizza la relazione discreta ottenuta prima per scrivere il relativo codice python che partendo da una condizione iniziale opportuna (la stabilità del sistema è sempre un problema), calcola come cambia il profilo di temperatura nel tempo. La simulazione fornisce in fine un’animazione che mostra dinamicamente il profilo di temperatura nel tempo.

Partiamo quindi con la prima parte. Per colmare alcuni passaggi nel tutorial originale ho inserito dei commenti personali che potete identificare rispetto al testo originale in quanto  sono fra parentesi ed in corsivo. Saltare questi commenti non compromette la comprensione del tutorial.

Buona lettura e programmazione.

Nota Introduttiva

Inizio con una personale nota per approfondire il polinomio di Taylor, che verrà utilizzato poco dopo nel tutorial, per passare dalle derivate parziali a differenze finite. 

Data una funzione f(x) continua insieme alle proprie derivate fino a quella di ordine n in un intorno di $$x_0$$ si può scrivere che per qualunque x nell’intorno di $$x_0$$ vale l’equazione:

$$f(x)=f(x_0)+(x-x_0)f^{‘}(x_0)+(x-x_0)^{2}\frac{f^{”}(x_0)}{2!}+(x-x_0)^{3}\frac{f^{”’}(x_0)}{3!}+….+(x-x_0)^{n}\frac{f^{n}(x_0)}{n!}+R_n(x-x_0))$$

dove l’ultimo termine indica una funzione (resto) per cui vale la proprietà:

$$\lim_{x \to x_0}\frac {R_n (x-x_0)}{(x-x_0)^n} =0$$

La funzione f(x) si può quindi scrivere come somma di un polinomio di grado n più la funzione $$R_n(x-x_0)$$. Quest’ultima funzione va a zero più rapidamente di $$(x-x_0)^n$$ avvicinandosi a $$x_0$$. Se si scrivesse la funzione f(x) con il polinomio di grado n di Taylor si commetterebbe un errore che espresso dalla funzione $$R_n(x-x_0)$$. Questo errore sarà tanto più trascurabile quanto più stringo l’intervallo intorno $$x_0$$. Ma vale anche al contrario che fissata una certa x in un intorno di $$x_0$$ l’errore $$R_n(x-x_0)$$ sarà via via inferiore se aumento il grado n del polinomio. 

Questi risultati hanno due interessanti conseguenze:

1) una qualunque funzione f(x) di classe $$C^n$$ si può scrivere come polinomio di grado n commettendo un errore. L’errore commesso è tanto più piccolo quanto più alto è il grado del polinomio. La figura seguente mostra graficamente quanto detto: la funzione coseno di x (curva verde) viene approssimata con un polinomio di grado crescente, rispettivamente le curve blu rossa e nera.

taylor-expansion-003

2) Data una certa funzione di classe $$C^n$$ è possibile calcolarne il valore numero in qualunque punto grazie al polinomio di grado n. Ad esempio si potrebbe calcolare il valore del numero e oppure di $$e^2$$ tramite lo sviluppo del polinomio o ancora cos(0.5) etc..

Ai fini del seguente tutorial il punto 1 sarà il metodo che ci consente di passare dall’equazione differenziale del calore ad una forma alla differenze finite.

Diffusione in una dimensione

Bentornati! Questo è un IPython Notebook della serie Spazio e Tempo-Introduzione delle soluzioni alle differenze finite delle PDE, il secondo modulo di “Practical Numerical Methods with Python”.

Nel precedente IPython notebooks di questa serie abbiamo studiato la soluzione numerica delle equazioni di convezione lineare e non lineare usando il metodo delle differenze finite ed abbiamo imparato ad usare le condizioni CFL (“Courant-Friedrichs-Lewy“). Ora passiamo a studiare l’equazione mono dimensionale della diffusione.

$$\partial_{t}u-\nu\partial^2_{x}u=0$$

dove $$\nu$$ è una costante nota come coefficiente di diffusione. La prima cosa da notare è che questa equazione è una equazione differenziale del secondo ordine. Per prima cosa abbiamo da imparare cosa fare con questa equazione!

Discretizzazione della derivata seconda

(“una nota sulla nomenclatura: data una funzione u(x,t) nelle righe seguenti si indicherà con l’espressione $$u^n_{i}=u(i,n)$$ ovvero la funzione u calcolata nel punto x=i e nel tempo t=n“)

La derivata seconda può essere rappresentata geometricamente come la linea tangente alla curva data dalla derivata prima. Si può rendere discreta la derivata seconda con lo schema della differenza centrale: una combinazione di differenze in avanti ed indietro della derivata prima. Consideriamo lo sviluppo di Taylor di $$u_{i+1}$$ e di $$u_{i-1}$$

$$u_{i+1}=u_{i}+\Delta x\frac{\partial u}{\partial x}\mid_{i}+\Delta x^2\frac{\partial^2u}{\partial^2x}\mid_{i}+\Delta x^3\frac{\partial^3u}{\partial^3x}\mid_{i}+O(x^4)$$

 $$u_{i-1}=u_{i}-\Delta x\frac{\partial u}{\partial x}+\Delta x^2\frac{\partial^2 u}{\partial^2 x}-\Delta x^3\frac{\partial^3 u}{\partial^3 x}+O(x^4)$$

Se le due espressioni vengono sommate, le derivate di ordine dispari si eliminano. Trascurando ogni termine di orine $$O(x^4)$$ o superiore (e sono veramente piccoli) possiamo riadattare la somma di queste due espansioni per trovare un’espressione per la derivata seconda:

$$u_{i+1}+u_{i-1}=2u_{i}+\Delta x^2\frac{\partial^2 u}{\partial^2 x}+O(x^4)$$

e quindi in fine

$$\frac{\partial^2 u}{\partial^2 x}=\nu \frac{u_{i+1}+u_{i-1}-2u_{i}}{\Delta x^2} +O(x^2)$$

L’approssimazione della derivata parziale è accurata al secondo ordine.

Tornando all’equazione di Diffusione

Possiamo quindi riscrivere l’equazione di diffusione del calore in una direzione nella versione discreta :

$$ \frac{u^{n+1}_{i}-u^n_{i}}{\Delta t}=\frac{u^n_{i+1}+u^n_{i-1}-2 u^n_{i}}{\Delta x^2}$$

Come osservato prima, una volta ottenuta una condizione iniziale il solo termine non noto è $$u^{n+1}_{i}$$ (ricordiamo che la notazione $$u^{n+1}_{i}$$ indica il valore della funzione u(x,t) nel punto i e nell’istante n+1 ) l’equazione quindi si può riscrivere:

$$ u^{n+1}_{i}=u^{n}_{i}+\frac{\nu \Delta t}{\Delta x^2}{}$$

Questa equazione discreta ci permette di scrivere un programma che mostra l’avanzamento nel tempo della funzione u(x,t), partendo però da una condizione iniziale. Continuiamo usando la condizione iniziale favorita: la funzione rettangolare. E quindi per t=0, u=2 nell’intervallo $$0.5\leq x\leq 1$$ ed u=1 altrove.

Stabilità dell’equazione di diffusione

L’equazione di diffusione è limitata da vincoli di stabilità. Come per l’equazione della convezione lineare e non lineare esistono insiemi di parametri di discretizzazione $$ \Delta x$$ e $$\Delta t$$ che rendono l’equazione divergente (Ovvero non siamo liberi di scegliere il passo in x ed in t che vogliamo, affinchè l’equazione precedente porti ad una soluzione). Per l’equazione di diffusione ed il metodo di discretizzazione utilizzato qui, la stabilità per l’equazione di diffusione è data da:

$$\nu\frac{\Delta x}{\Delta t}\leq \frac{1}{2} $$

(A questo punto abbiamo un’equazione che ci indica per ogni punto x quanto sarà la temperatura nell’istante successivo (n+1), noto l’istante precedente. Quindi fissate le condizioni iniziali ed i passi di integrazione opportuni possiamo scrivere un programma che mostri questa evoluzione. Ciò è quanto vedremo nel prossimo post.)

CC BY-NC-SA 4.0
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.