Ciao, vorrei condividere un tutorial su un esempio di simulazioni fatto con Python. Il caso che propongo è una traduzione di un post  dei professori L.A. Barba, G.F. Forsyth, C.Cooper della George Washington University che trovate in originale a questo link.

L’articolo propone una simulazione dell’evoluzione nel tempo di un profilo termico, risolvendo numericamente l’equazione di diffusione (del calore) in una sola dimensione. Il risultato non sarà quindi bello come quello della simulazione in 3D che vedete all’inizio ma è come se di questa vedessimo l’intersezione con un piano “verticale”.

Visti gli argomenti trattati ho diviso il post in tre parti, la prima è un cenno, a volo d’uccello, sulle derivate ordinarie, parziali e sulle relative equazioni. La seconda e terza parte sono invece il post originale (una traduzione). La seconda parte riguarda il metodo di discretizzazione dell’equazione di diffusione, la terza il codice python per la simulazione. Nella seconda e terza parte mi sono permesso qualche commento personale che potrete identificare in quanto in corsivo e fra parentesi.

Una “svolazzata” su derivate ed equazioni differenziali.

Partiamo dalla derivata ordinaria. L’idea nasce da funzioni di una sola variabile, reale, f(x): se x varia, varia anche la funzione f(x); a pari Δx la variazione di f(x) dipende dalla funzione stessa; per avere un’idea di come varia la funzione al variare di x si introduce il rapporto incrementale Δf(x)/Δx. Se immaginiamo di ridurre il Δx a intervalli sempre più piccoli (al limite infinitesimi) il rapporto incrementale tende ad un particolare valore che è la derivata della funzione nel punto (Δx si chiude su un punto generico).

Solitamente la derivata si indica con $$df(x)/dx$$ oppure $$f'(x)$$ o $$f_{x}$$. Non tutte le funzioni sono derivabili in tutti i punti del proprio dominio e l’analisi fornisce precise ipotesi al riguardo. Per le funzioni derivabili esistono delle regole molto chiare e semplici che rendono l’operazione di derivazione quasi meccanica. I libri di analisi fornisco tutte le dimostrazioni ed ipotesi necessarie. Spesso sono disponibili negli stessi libri le tabelle di derivazioni di molte funzioni. Alcuni esempi di derivate sono $$ f(x)=x^n $$ allora $$f_{x}=nx^{n-1} $$ oppure $$ f(x)=e^x $$ allora $$f_{x}=e^x $$ e poi $$ f(x)=cos(x) $$ allora $$f_{x}=-sen(x) $$, e poi $$ f(x)=sen(x) $$ allora $$f_{x}=cos(x) $$.

Il concetto di derivata non è solo un’idea matematica ma ha un forte legame con la realtà fisica: infatti la derivata si può leggere come la variazione infinitesimale (piccolissima) di una variabile fisica che è legata ad un’altra variabile fisica. Ad esempio, se un corpo è in movimento, lo spazio percorso da questo corpo rispetto al punto iniziale cambia con il tempo. Possiamo allora pensare alla variazione dello spazio in un intervallo di tempo: questa variazione è la velocità del corpo.

Se consideriamo un intervallo di tempo misurabile (finito) ed il relativo spazio percorso, abbiamo una velocità media (se ad esempio in 1h abbiamo percorso 105Km si può dire che siamo andati in media a 105 Km/h) se pensiamo che la variazione di tempo è infinitesimale (e quindi anche la variazione di spazio è infinitesimale) otteniamo la velocità puntuale del corpo: matematicamente possiamo dire allora che la velocità di un corpo è la derivata dello spazio rispetto al tempo; analogamente la corrente in un filo di rame è la quantità di carica passata attraverso una sezione in un intervallo di tempo. Quando collegate un qualunque carico elettrico alle comuni prese di casa, la corrente nel filo non sarà costante ma varierà con il tempo seguendo la variazione di tensione: la corrente sarà esattamente la variazione di carica che attraversa la sezione del filo in un intervallo di tempo infinitesimo.

Se proviamo ad affrontare in generale un problema di fisica le variabili da considerare sono certamente più di una. Basti pensare alla temperatura dell’aria intorno ad una stufa: supponiamo di avere in una stanza una stufa collocata in uno dei vertici della stanza. Se potessimo misurare la temperatura dell’aria con elevata precisione (si pensi a qualche decimo di grado) la temperatura si ridurrà allontanandoci dalla stufa; invece fissato un punto sul pavimento della stanza la temperatura dell’aria tenderà ad aumentare andando dal pavimento verso il tetto. Se immaginiamo che nella stanza sia fissato un sistema di riferimento cartesiano tridimensionale con il vertice dove si trova la stufa potremmo dire che la temperatura nella stanza sarà una funzione dello spazio cioè delle tre coordinate x,y,z.

E’ allora naturale chiedersi come varia la grandezza in questione rispetto a tutte le altre grandezze che la influenzano. La derivata parziale è proprio questo: la derivata di una funzione a più variabili rispetto ad una sola di queste, supponendo costanti tutte le altre. La nomenclatura sarà $$\frac{\partial u}{\partial x},\partial _{x} u$$ oppure $$u_{x}$$.

Qualche esempio: $$u(x,y)=x^2+y^2$$ possiamo calcolare la derivata $$\frac{\partial u}{\partial x}$$ poi la $$\frac{\partial u}{\partial y}$$ ed ancora $$\frac{\partial^2 u}{\partial y\partial x}$$ nel caso specifico:

$$\frac{\partial u}{\partial x}=2x$$

$$\frac{\partial u}{\partial y}=2y$$

$$\frac{\partial^2 u}{\partial x\partial y}=\frac{\partial^2 u}{\partial y\partial x}=\frac{\partial }{\partial x}(\frac{\partial u}{\partial y})=\frac{\partial }{\partial y}(\frac{\partial u}{\partial x})=0$$

Come appena visto, nelle derivate parziali si presenta il caso delle derivate miste cioè la derivata della funzione u(x,y) prima rispetto ad x e poi rispetto ad y. In generale l’ordine di derivazione cambia il risultato: ($$\frac{\partial^2 u}{\partial x\partial y}\neq\frac{\partial^2 u}{\partial y\partial x} $$). Il teorema di Schwarz fornisce le condizioni perché l’ordine di derivazione non sia rilevante.

A livello operativo, conoscendo le regole di derivazione delle funzioni di una sola variabile si calcola facilmente la derivata parziale: si calcola la derivata rispetto una sola variabile considerando le altre costanti.

 Le derivate delle funzioni sono esse stesse delle funzioni e possono essere quindi derivate successivamente ottenendo la derivata seconda, terza e cosi via: il numero di volte con il quale una funzione viene derivata si chiama ordine della derivata.

Le equazioni alle derivate parziali (PDE, acronimo spesso usato nei testi anglosassoni per partial differential equation) sono equazioni in cui compaiono una o più derivate parziali di una funzione non conosciuta che, chiaramente, è funzione di più variabili. Formalmente:

$$ f(x_{1},x_{2}….x_{m}, \frac{\partial u}{\partial x_{1}},…,\frac{\partial^n u}{\partial x_{m}^n})=0 $$

dove la funzione u dipende dalle variabili indipendenti $$x_{1},x_{2}….x_{m}$$.

Esempio:

(1) $$-{\partial^2_{t} u}+(1+ cos u)\partial^3_{x}u=0 $$ equazione differenziale (PDE) non lineare del terzo ordine (massimo ordine della derivata presente), oppure

(2) $$-{\partial_{t} u}+2\partial_{x}u+u=t$$ PDE lineare del secondo ordine (massimo ordine della derivata presente)

L’ordine di una equazione differenziale è il massimo ordine delle derivate che in essa compaiono: nella (1) ad esempio abbiamo una derivata del terzo ordine per cui l’equazione è del terzo ordine. Una equazione differenziale è detta lineare se la funzione incognita u non è elevata a potenza, le sue derivate sono del primo ordine e non ci sono prodotti di u per le sue derivate: la (2) è ad esempio lineare.

Il problema in una equazione differenziale è trovare la funzione u che soddisfa l’equazione stessa. Trovare queste soluzioni non è semplice ed infatti non esistono metodi generali per risolvere le equazioni. In genere si prova a ricondurre l’equazione a casi noti di equazioni simili (ad esempio altre equazioni lineari o dello stesso grado o con le stesse derivate).

Inoltre nel risolvere un’equazione differenziale si introducono necessariamente dei gradi di libertà nelle soluzioni: questi gradi di libertà possono essere delle costanti arbitrarie o delle funzioni arbitrarie. Matematicamente le soluzioni trovate vanno bene ma restano del tutto generali. Faccio un esempio per essere più chiaro: data l’equazione $$ u_{xx}-u=0$$ con u funzione di x e y. La soluzione di questa equazione è la funzione: $$u(x,y)=A(y)e^x+B(y)e^{-x})$$. Qualunque siano le funzioni in y, A(y) e B(y) la u trovata è una soluzione dell’equazione iniziale: è la soluzione generale dell’equazione.

Se volessimo essere più precisi su u(x,y) ci serve qualche condizione in più, oltre l’equazione differenziale. Le condizioni aggiuntive sono dette condizioni al contorno ed infatti indicano dei vincoli che la funzione deve assumere in determinati punti del suo dominio. Ad esempio nella soluzione trovata possiamo imporre che u(0,y)=2y e u(L,y)=2 così facendo otteniamo:

$$ A(y)=\frac{2(ye^{-L}-1)}{e^{-L}-e^L}$$

$$ B(y)=\frac{2(1-ye^{-L})}{e^{-L}-e^L}$$

Le condizioni al contorno fissate ci permettono quindi di trovare una specifica funzione che risolve l’equazione differenziale iniziale e soddisfa le condizioni date: in questo caso si parla di soluzione particolare.  Le condizioni al contorno sono molto importanti nella soluzione delle PDE, infatti, al mutare delle condizioni al contorno una equazione può essere risolvibile o meno. Così come le equazioni differenziali descrivono fenomeni fisici è utile sottolineare che anche le condizioni al contorno hanno un significato fisico indicando solitamente lo stato del sistema nella fase iniziale del processo o nella fase finale o altri vincoli del sistema studiato.

Dovendo parlare nei prossimi post dell’equazione del calore, proviamo ad analizzare questa equazione nel caso monodimensionale: $$u_t=c^2u_{xx} $$ dove c indica la diffusività termica e dipende dal particolare materiale. Questo tipo di equazioni si possono risolvere tramite separazione delle variabili, assumendo cioè che la soluzione abbia una forma del tipo u(x,y)=F(y)G(x). questo metodo è spesso usato nella soluzioni di equazioni differenziali. Come visto prima, per trovare una soluzione all’equazione di diffusione occorre fissare le condizioni al contorno. Se ne possono considerare diverse, barra finita o infinita isolata o meno, alcune condizioni al contorno che permettono di trovare una soluzione analitica dell’equazione di diffusione sono date da u(0,t)=0 u(L,t)=0 e u(x,0)=f(x). Fisicamente le condizioni al contorno considerate impongono che data una barra di lunghezza finita L, ai suoi estremi la temperatura resti fissa nel tempo (ed uguale a zero) e che all’inizio del nostro esame la distribuzione di temperatura lungo il corpo abbia un certo andamento f(x). Quest’ultima funzione (per ragioni in cui non entrerò) non può essere una funzione arbitraria ma deve avere certe caratteristiche (continuità) affinché l’equazione di diffusione sia risolvibile. A queste condizioni la soluzione è:

$$ u(x,t)=\sum_{n=1}^{+\infty}B_{n}sin\frac{n{\pi}x}{L}e^{\frac{cn{\pi}t}{L}}$$

$$ B_{n}=\int_0^L f(x)sin\frac{n{\pi}x}{L}dx $$

Il risultato analitico, va ben oltre lo scopo del post; aggiungo solo due considerazioni su questa soluzione: 1) la funzione u(x,t) è data dalla somma di infinite sinusoidi di ampiezza $$B_{n}$$ e frequenza $$\frac{n{\pi}x}{L}$$. I coefficienti $$B_{n}$$ sono determinati dalla funzione iniziale f(x). Se f(x) fosse una sinusoide u(x,t) sarebbe formata da una sola sinusoide (esisterebbe un solo $$B_{1}$$ mentre $$B_{2}=B_{3}=….=0$$). 2) ogni sinusoide è moltiplicata per un esponenziale negativo del tempo. Ciò significa che tutte le componenti di u(x,t) avranno un’ampiezza decrescente nel tempo. Visto che u(x,0)=f(x), la soluzione trovata ci indica che il grafico della f(x) si spalmerà su x, appiattendosi al passare del tempo.

Conoscendo la soluzione analitica dell’equazione del calore la si può utilizzare per sviluppare delle simulazioni. Si tratterebbe di sviluppare i coefficienti $$B_n$$ a partire da f(x) ed utilizzarli per riscrivere u(x,t). Questa strada condurrebbe ad un risultato molto accurato ma impegnativo come elaborazione. Per sviluppare una simulazione non è indispensabile conoscere la soluzione analitica dell’equazione ma si può procedere con funzioni approssimate. Questo e’ quello che descriverò nella seconda parte del post per poi scrivere il codice in Python per la definitiva simulazione del profilo di temperatura.

Ciao alla prossima!

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