Universita del SalentoFacolta di Scienze MM.FF.NN.

Corso di Laurea in Fisica

Tesi di Laurea triennale

Aspetti computazionalidi equazioni differenziali

di interesse fisico.

RelatoreClaudio Coriano

CandidatoLuigi Manni

Lecce, 20 dicembre 2007

Anno accademico 2006/2007

Prefazione

Lo studio numerico di equazioni differenziali di interesse fisico e, piu in ge-nerale, scientifico, e uno tra i settori piu interessanti sia nell’ambito teorico-formale che nelle applicazioni di carattere tecnologico. Basti pensare che lamodellizzazione della dinamica dei fluidi o lo studio di processi dissipativinell’ambito di complessi sistemi dinamici o ancora in fisica statistica richiedenon solo una corretta formulazione matematica nel continuo delle corrispon-denti equazioni differenziali, ma anche un’opportuna discretizzazione su reti-colo delle medesime. Questa complementarieta o convergenza di intenti perla risoluzione di problemi che hanno una forte rilevanza scientifica e praticaimpone lo sviluppo di strumenti sofisticati di analisi ma anche l’elaborazionedi specifiche tecniche numeriche.Mentre lo studio (analitico/formale) del problema di Cauchy di equazionidi propagazione (cioe di equazioni non ellittiche) permette di distinguere traproblemi ben posti e problemi mal posti nella ricerca dell’esistenza ed unicitadelle corrispondenti soluzioni, lo studio numerico da la possibilita di studiarecasi molto complessi per i quali una soluzione analitica e molto difficile da ot-tenere. Per queste ragioni una vera e propria branca della fisica/matematicaapplicata che ha un forte connotato informatico si e venuta affermando inambito internazionale, caratterizzandosi come un’attivita scientifica multi-disciplinare che riceve nelle migliori Universita straniere una collocazioneaccademica specifica: la Scienza Computazionale o ComputationalScience. Sotto questo nome si inquadrano gli studi che confinano con lamatematica applicata, la fisica e l’informatica (Computer Science) e chesi avvalgono dell’uso massiccio delle capacita di calcolo offerte dalle modernereti di calcolatori e dal calcolo distribuito.L’obiettivo di questa tesi e quello di passare brevemente in rassegna alcunidegli aspetti fondamentali di questo processo e di descrivere le basi salientie piu elementari di questo tipo di studi, cosı da fornire un punto di partenzaper analisi piu estese che l’autore spera di poter mandare avanti nella secondaparte del corso di studi.Nell’ambito di questo lavoro che, come abbiamo appena accennato, e es-

1

senzialmente di rassegna, in cui introduciamo la metodologia essenziale perpassare dalla descrizione di equazioni al continuo al caso discreto, abbiamoanche avuto modo di affrontare una tematica originale e molto complessa cheillustra in modo molto diretto le difficolta che si incontrano a livello nume-rico. Questa consiste nello studio di un’equazione molto specifica introdottada Eduardo Pascali in alcuni lavori recenti [14, 15, 16] che forse porta aduno tra i casi piu complessi di ricerca di soluzioni numeriche di equazionidifferenziali di tipo diffusivo e che motiveremo in seguito. Tali equazioni, chepossono essere formulate sia in forma scalare che matriciale, sono caratte-rizzate da un lato destro di tipo non-locale e, come chiariremo nel corso dellavoro, auto-referenziale. Equazioni di questo tipo sono particolarmentedifficili da analizzare numericamente perche, come vedremo nel corso dellanostra analisi, sono caratterizzate da una soluzione marciante che presentauna regione di influenza che tende a diventare troppo ampia nella direzionespaziale a mano a mano che si avanza nella direzione temporale. Equazionidi questo tipo, che possono essere utili in vari contesti applicativi, sono, sottoun profilo fisico, classificabili come processi deterministici non-markoviani incui la auto-referenzialita dell’equazione svolge, in un certo senso, il ruolo dimemoria del processo stesso.Prima di passare ad uno studio dell’equazione auto-referenziale di Pascaliintrodurremo alcuni dei concetti di base che riguardano la classificazione deitipi piu semplici di equazioni differenziali. Successivamente ci occuperemodella loro discretizzazione, definendo l’ordine di accuratezza ed i parametriche sono essenziali per una caratterizzazione della stabilita (o instabilita) diun certo schema discreto. Nel terzo capitolo introdurremo e giustificheremol’equazione di Pascali nel caso scalare e cercheremo di definire le proprietadella corrispondente soluzione numerica marciante. Introdurremo un algo-ritmo che permette di risolvere questa equazione nell’ambito di uno schemamarciante e ne descriveremo l’implementazione.Abbiamo optato per una descrizione literate dei programmi numerici, se-condo la definizione data da D. Knuth, che altro non sono se non programmiscritti in un linguaggio di programmazione specifico (nel nostro caso in lin-guaggio C) interleaved, cioe intermezzati con il testo usuale. Questo stiledi presentazione (detto anche di dissezione, o by dissection) rende i li-terate programs piu facili da leggere e permette di comunicare al lettoreche conosce il linguaggio la struttura dell’implementazione.

2

Indice

1 Introduzione alle equazioni alle derivate parziali 51.1 Classificazione di equazioni alle derivate parziali quasi lineari. 5

1.1.1 Metodo di Cramer . . . . . . . . . . . . . . . . . . . . 61.1.2 Metodo agli autovalori . . . . . . . . . . . . . . . . . . 10

1.2 Tipi di equazioni alle derivate parziali quasi lineari. . . . . . . 111.2.1 Equazioni iperboliche . . . . . . . . . . . . . . . . . . . 111.2.2 Equazioni paraboliche . . . . . . . . . . . . . . . . . . 131.2.3 Equazioni ellittiche . . . . . . . . . . . . . . . . . . . . 14

2 Introduzione alla discretizzazione 162.1 Il processo di discretizzazione . . . . . . . . . . . . . . . . . . 162.2 Metodo delle differenze finite . . . . . . . . . . . . . . . . . . . 172.3 Equazioni alle differenze . . . . . . . . . . . . . . . . . . . . . 242.4 Approcci alle soluzioni . . . . . . . . . . . . . . . . . . . . . . 262.5 Analisi degli errori e studio della stabilita . . . . . . . . . . . . 312.6 Risoluzione numerica dell’equazione di

propagazione del calore in una dimensione . . . . . . . . . . . 40

3 Equazioni autoreferenziali 493.1 Equazioni evolutive in fisica fondamentale ed in statistica . . . 503.2 Un’equazione con memoria . . . . . . . . . . . . . . . . . . . . 523.3 Interpretazione e discretizzazione . . . . . . . . . . . . . . . . 533.4 Codice sorgente e grafici . . . . . . . . . . . . . . . . . . . . . 56

4 Verso la CFD 704.1 Cenni preliminari . . . . . . . . . . . . . . . . . . . . . . . . . 704.2 Le equazioni della CFD . . . . . . . . . . . . . . . . . . . . . . 71

4.2.1 L’equazione di continuita . . . . . . . . . . . . . . . . . 714.2.2 L’equazione del momento . . . . . . . . . . . . . . . . . 734.2.3 L’equazione dell’energia . . . . . . . . . . . . . . . . . 754.2.4 Equazioni di Navier-Stokes . . . . . . . . . . . . . . . . 78

3

4.2.5 Equazioni di Eulero . . . . . . . . . . . . . . . . . . . . 794.3 Il metodo di Lax - Wendroff . . . . . . . . . . . . . . . . . . . 81

5 Conclusioni 86

A Dimostrazioni relative alle equazioni differenziali alle deriva-te parziali 87A.1 Dimostrazione: la (2.28) e parabolica . . . . . . . . . . . . . . 87A.2 Dimostrazione: la (2.28) e iperbolica . . . . . . . . . . . . . . 88

B Il metodo di Gauss-Legendre 89

Bibliografia 91

4

Capitolo 1

Introduzione alle equazioni allederivate parziali

In questo primo capitolo vogliamo avvicinarci allo studio delle equazioni dif-ferenziali alle derivate parziali. Introduciamo quindi alcuni metodi che con-sentono di classificare tali equazioni e procediamo quindi ad una loro analisipiu particolareggiata mettendo in luce le loro caratteristiche salienti.

1.1 Classificazione di equazioni alle derivate

parziali quasi lineari.

Definiamo sistema quasi lineare un sistema di equazioni differenziali allederivate parziali in cui le derivate all’ordine piu alto si presentano ne molti-plicate tra loro ne tanto meno sottoforma di esponenziali, bensı linearmentemoltiplicate per coefficienti che risultano essere funzioni delle variabili dipen-denti.Dallo studio delle equazioni alle derivate parziali, notiamo che esse presen-tano dei comportamenti matematici molto diversi tra loro. E’ dunque utileclassificare tali tipi di equazioni. Queste possono essere di tipo iperbolico,parabolico, ellittico o misto (se non presentano uno dei primi tre andamen-ti). Di seguito esponiamo due diversi metodi che permettono di analizzare unsistema di equazioni quasi lineari allo scopo di determinarne la classificazione.

5

1.1.1 Metodo di Cramer

Consideriamo il seguente sistema di equazioni alle derivate parziali quasilineari

a1∂u

∂x+ b1

∂u

∂y+ c1

∂v

∂x+ d1

∂v

∂y= f1 (1.1)

a2∂u

∂x+ b2

∂u

∂y+ c2

∂v

∂x+ d2

∂v

∂y= f2 (1.2)

dove u e v sono le variabili dipendenti (funzioni di x e y), e i coefficientia1, a2, b1, b2, c1, c2, d1, d2, f1 ed f2 possono essere funzioni di x, y, u e v. u e vsono funzioni continue di x e y. Quindi, il loro differenziale totale puo esserescritto come

du =∂u

∂xdx +

∂u

∂ydy (1.3)

dv =∂v

∂xdx +

∂v

∂ydy. (1.4)

Notiamo anche che ad ogni punto nel piano xy e associato un unico valoredi u e un unico valore di v, mentre le loro derivate risultano essere ivi finite.Con riferimento alla figura 1.1, fissato un punto P nel piano xy, cerchiamole linee o direzioni passanti per P lungo le quali le derivate di u e v sonoindeterminate e attraverso le quali possono presentarsi delle discontinuita.Tali direzioni sono dette linee caratteristiche.Consideriamo ora il sistema di quattro equazioni in quattro incognite (∂u/∂x,∂u/∂y,∂v/∂x,∂v/∂y) costituito dalle equazioni lineari (1.1), (1.2), (1.3) e(1.4). Possiamo scrivere il tutto in forma matriciale

a1 b1 c1 d1

a2 b2 c2 d2

dx dy 0 00 0 dx dy

∂u/∂x∂u/∂y∂v/∂x∂v/∂y

=

f1

f2

dudv

(1.5)

Denotiamo con [A] la matrice dei coefficienti, cioe

[A] ≡

a1 b1 c1 d1

a2 b2 c2 d2

dx dy 0 00 0 dx dy

(1.6)

Risolviamo ora, con la regola di Cramer, l’equazione (1.5) per l’incognita∂u/∂x. Introduciamo quindi la matrice [B], costruita a partire dalla matrice

6

[A] sostituendo la prima colonna col vettore colonna dei termini noti della(1.5)

[B] =

f1 b1 c1 d1

f2 b2 c2 d2

du dy 0 0dv 0 dx dy

(1.7)

Utilizzando la regola di Cramer scriviamo

∂u

∂x=|A||B| (1.8)

in cui con |A| e |B| abbiamo indicato rispettivamente i determinanti dellematrici [A] e [B]. Per ottenere un valore dall’equazione precedente dobbia-mo pero assegnare alle quantita dx, dy, du e dv dei valori ben definiti.

Figura 1.1: Esempio di curva caratteristica.

Con riferimento alla figura 1.1 immaginiamo di muoverci lungo una linea abpassante per il punto P nel piano xy. Supponiamo di effettuare uno sposta-mento infinitesimo ds verso il punto Q. In questo caso sia x che y subirannouna variazione infinitesima data da dx = xQ − xP e dy = yQ − yP . Ana-logo discorso puo esser fatto per le quantita u e v, quindi du = uQ − uP edv = vQ − vP . Una volta note queste quantita, possiamo porle all’internodelle matrici [A] e [B], e attraverso la (1.8), calcolare una delle quattro inco-gnite del sistema (1.5).Il valore cosı trovato risulta essere indipendente dalla direzione lungo la qualeci allontaniamo dal punto P, sempre nell’ipotesi di un cammino infinitesimo.Se ad esempio ci muoviamo lungo la linea cd, i valori di dx, dy, du e dv sononaturalmente diversi dai precedenti, tuttavia svolgendo i calcoli si ottieneun identico valore della variabile ∂u/∂x. Questo valore risulta essere fissato,

7

potremmo intenderlo come valore puntuale, mentre la scelta della curva pas-sante per P sulla quale ci muoviamo e del tutto arbitraria.Definiamo come curva caratteristica quella curva, ef nella figura 1.1,muovendosi lungo la quale le variazioni infinitesime sono tali da rispettare laseguente condizione

|A| = 0. (1.9)

Proprio ponendo questa condizione siamo in grado di trovare, se esistono, lecurve caratteristiche passanti per il generico punto P nel piano xy. Natural-mente la determinazione di tali linee e indipendente da quale delle quattroequazioni andiamo a risolvere dal momento che il denominatore |A| e comu-ne a tutte le equazioni analoghe alla (1.8). Quando queste particolari lineeesistono, possiamo calcolare le loro equazioni ed in particolare la pendenzanel punto P, infine tracciarle nel piano xy, come nel caso della linea ef nellafigura di cui sopra. Riscriviamo allora la (1.9) e risolviamo il determinante.

∣∣∣∣∣∣∣∣

a1 b1 c1 d1

a2 b2 c2 d2

dx dy 0 00 0 dx dy

∣∣∣∣∣∣∣∣= 0 (1.10)

(a1c2 − a2c1)dy2 − (a1d2 − a2d1)dxdy + (b1d2 − b2d1)dx2 = 0. (1.11)

Dividiamo a questo punto ambo i membri per dx2, ottenendo

(a1c2 − a2c1)

(dy

dx

)2

− (a1d2 − a2d1)

(dy

dx

)+ (b1d2 − b2d1) = 0. (1.12)

Notiamo che quest’ultima e un’equazione quadratica in (∂y/∂x). La soluzio-ne dell’equazione (1.12) ci da le pendenze delle curve caratteristiche passantiper il generico punto P. Poniamo ora

a = (a1c2 − a2c1) (1.13)

b = −(a1d2 − a2d1 + b1c2 − b2c1) (1.14)

c = (b1d2 − b2d1) (1.15)

di modo che la (1.12) assuma la forma seguente

a

(∂y

∂x

)2

+ b

(∂y

∂x

)+ c = 0. (1.16)

8

Quest’ultima equazione puo essere integrata per ottenere l’equazione y =y(x) della curva caratteristica, oppure puo, attraverso la determinazione dellapendenza

∂y

∂x=−b±√b2 − 4ac

2a, (1.17)

indicarci la direzione della stessa curva caratteristica nel piano xy.Inoltre, la (1.17) ci permette, attraverso la determinazione del valore del suodiscriminante, di classificare il comportamento delle varie linee caratteristi-che.Sia D il discriminante della (1.17), cioe

D = b2 − 4ac. (1.18)

Se D > 0 esistono due caratteristiche reali e distinte passanti per ogni puntodel piano xy. Il sistema costituito dalle equazioni (1.1) e (1.2) e detto iper-bolico.Se D = 0 il sistema in questione e detto parabolico.Se D < 0 le linee caratteristiche sono immaginarie ed il sistema e definitoellittico.L’origine degli appellativi iperbolico, parabolico ed ellittico, con i quali ven-gono definiti i sistemi di equazioni alle derivate parziali quasi lineari, traeorigine da un’analogia con le sezioni coniche, la cui equazione generale e datada

ax2 + bxy + cy2 + dx + ey + f = 0. (1.19)

In questo caso, se

• b2 − 4ac > 0 si ha un’iperbole;

• b2 − 4ac = 0 si ha una parabola;

• b2 − 4ac < 0 si ha un’ellisse.

Torniamo a considerare l’equazione (1.8) e notiamo che qualora fosse rispet-tata la condizione (1.9), la derivata (∂u/∂x) risulterebbe infinita, contrav-venendo all’ipotesi stessa che riguarda le curve caratteristiche. La derivatadeve essere indeterminata ma non infinita. Per ovviare a tutto cio occorreche anche il determinante di [B] sia nullo in modo tale che la (1.8) si presentinella forma

∂u

∂x=|A||B| =

0

0, (1.20)

9

mentre

|B| =

∣∣∣∣∣∣∣∣

f1 b1 c1 d1

f2 b2 c2 d2

du dy 0 0dv 0 dx dy

∣∣∣∣∣∣∣∣= 0 (1.21)

Lo sviluppo di tale determinante genera un’equazione differenziale ordinarianelle variabili du e dv, mentre dx e dy sono limitati a restare lungo la lineacaratteristica. Procedendo con la sua risoluzione si ottiene un’equazione perle variabili dipendenti u e v che e detta equazione di compatibilita. Taleequazione risulta avere una dimensione in meno di quelle alle derivate par-ziali di partenza, ed essere, in generale, piu facile da risolvere. Cio conducead una tecnica risolutiva del problema originario nota come metodo dellecaratteristiche, secondo cui si risolve la semplice equazione di compatibi-lita solo sulle linee caratteristiche individuate nel piano xy. Tale metodo,richiedendo almeno due curve caratteristiche, e applicabile al solo caso delleequazioni iperboliche. Rimandiamo il lettore al [12] per l’applicazione di taletecnica al caso di fluidi supersonici non viscosi.

1.1.2 Metodo agli autovalori

Questo metodo, basato sugli autovalori del sistema, e piu generale e legger-mente piu sofisticato del precedente.Consideriamo il sistema formato dalle equazioni (1.1) e (1.2) del paragrafoprecedente, in cui per semplicita supponiamo che sia f1 che f2 siano nulli,quindi

a1∂u

∂x+ b1

∂u

∂y+ c1

∂v

∂x+ d1

∂v

∂y= 0 (1.22)

a2∂u

∂x+ b2

∂u

∂y+ c2

∂v

∂x+ d2

∂v

∂y= 0. (1.23)

Utilizzando il formalismo matriciale, ed introducendo il vettore colonna W ,dato da

W =

[uv

](1.24)

possiamo scrivere il tutto come

[a1 c1

a2 c2

]∂W

∂x+

[b1 d1

b2 d2

]∂W

∂y= 0. (1.25)

10

Ponendo quindi [K] e [M ] uguali alle matrici che moltiplicano rispettivamente(∂W/∂x) e (∂W/∂y) si ha

[K]∂W

∂x+ [M ]

∂W

∂y= 0. (1.26)

Moltiplicando la precedente equazione per [K]−1, inversa1 della matrice [K],otteniamo

∂W

∂x+ [K]−1[M ]

∂W

∂y= 0 (1.27)

o equivalentemente∂W

∂x+ [N ]

∂W

∂y= 0 (1.28)

in cui [N ] = [K]−1 · [M ]. Proprio gli autovalori di quest’ultima matrice deter-minano la classificazione del sistema in quanto non sono altro che le pendenzedelle linee caratteristiche. Infatti, se gli autovalori sono tutti reali le equazio-ni sono iperboliche, mentre se sono immaginari le equazioni sono ellittiche2.Talvolta pero gli autovalori di [N ] possono presentarsi come un misto travalori reali e complessi. In questi casi non e possibile classificare in una dellecategorie viste prima le equazioni date poiche queste presentano un compor-tamento eterogeneo.

1.2 Tipi di equazioni alle derivate parziali qua-

si lineari.

Attraverso l’utilizzo dei metodi presentati nel paragrafo precedente possia-mo classificare le varie equazioni, oggetto dei nostri studi, in tre categorieprincipali che rispecchiano altrettanti comportamenti a livello matematico.

1.2.1 Equazioni iperboliche

Consideriamo un’equazione iperbolica nelle variabili indipendenti x e y.Come gia osservato un’equazione iperbolica ammette due curve caratteristi-che reali e distinte per ogni punto P nel piano xy. Denotiamo tali curve conil nome di direzione destra e sinistra. Per fare cio immaginiamo di porci nel

1Per ottenere la matrice inversa di una data se ne calcola prima la matrice dei cofattori,quindi la trasposta, ed infine si moltiplica quest’ultima per il determinante della matricedi partenza.

2[13] per maggiori dettagli ed esempi.

11

punto P del piano xy con lo sguardo volto verso la direzione positiva dell’as-se x. La curva che si trova alla nostra destra viene quindi chiamata curvacaratteristica con direzione destra, e viceversa per quella alla nostra sinistra.Queste due linee che passano per P ne determinano la relativa area d’in-fluenza. Una lieve perturbazione che coinvolge tale punto produce una certamodificazione nei punti che si trovano nella regione I della figura 1.2 che epercio detta regione d’influenza di P. Al contrario la parte di piano indi-viduata dall’asse y e dalle due caratteristiche passanti per P ne determinanoil cosiddetto dominio di dipendenza di P (regione III) poiche il compor-tamento nel punto P dipende solo da cio che avviene in tale regione, in cuivengono fissate le condizioni al contorno.

Figura 1.2: Dominio e contorni della soluzione di un’equazione iperbolica indue dimensioni.

Un generico punto c che non ricade in nessuna delle due regioni appenaesposte non e legato in alcun modo al punto P, ma avra una propria areadi influenza (come mostrato in 1.2 dalla regione II), ed anche un dominio didipendenza se non e caratterizzato da x = 0, cioe se non si trova sul contornocostituito dall’asse y.Queste equazioni vengono risolte con degli algoritmi basati su metodi iterati-vi che prendono le mosse dalle condizioni iniziali fissate, ad esempio, sull’assey.

12

1.2.2 Equazioni paraboliche

Consideriamo un’equazione parabolica nelle variabili indipendenti x e y.Ricordiamo che un’equazione parabolica presenta un’unica curva caratteri-stica passante per il generico punto P del piano xy. Analizziamo questasoluzione schematizzata nella figura 1.3.

Figura 1.3: Dominio e contorni della soluzione di un’equazione parabolica indue dimensioni.

Supponiamo che le condizioni iniziali del problema siano date sull’asse y nel-l’intervallo ac e che siano note le condizioni al contorno lungo le curve ab ecd. La linea caratteristica, denotata con un tratteggio nella figura, e datada una linea verticale passante per P e congiungente le linee di contorno. Laregione di influenza di P si estende allora tra tale linea e le curve sulle qualisono definite le condizioni al contorno nella direzione delle x crescenti.Questo tipo di equazioni, cosı come le precedenti, presenta delle soluzionimarcianti nel senso che partendo dalle condizioni iniziali, assegnate sulla li-nea ac, la soluzione, tra i contorni cd e ab, e ottenuta procedendo nel versopositivo dell’asse x.La figura 1.4 mostra un esempio di una possibile estensione in 3D, dove levariabili indipendenti risultano essere tre (ad esempio x,y, z). Le condizioniiniziali non sono piu date lungo una linea (1D) bensı sull’area con contornoabcd (2D) nel piano yz. Risulta invece essere ancora valida la propagazionedell’informazione lungo la direzione positiva delle x. In questo caso la regio-ne di influenza del generico punto P non e un’area (2D) ma il volume (3D)tratteggiato.

13

Figura 1.4: Dominio e contorni della soluzione di un’equazione parabolica intre dimensioni.

1.2.3 Equazioni ellittiche

Consideriamo un’equazione ellittica nelle variabili indipendenti x e y.Dal momento che in questo caso le due curve caratteristiche sono immagi-narie non possiamo ricorrere al metodo delle caratteristiche per risolvere ilproblema. Per le equazioni ellittiche le regioni di influenza e i domini di di-pendenza risultano essere non limitati. Inoltre, l’informazione si propaga daogni punto in tutte le direzioni.Con riferimento alla figura 1.5 immaginiamo che il dominio sia definito dalquadrato abcd e che il generico punto P sia al suo interno. La peculiarita delsistema considerato e che il punto P influenza tutti gli altri punti del domi-nio e, a sua volta, viene influenzato dall’intero contorno abcd. La risoluzionedell’equazione nel punto P dovrebbe avvenire contemporaneamente con tuttigli altri punti del dominio.Cio e in netto contrasto con i processi iterativi visti fin ora per le equazio-ni iperboliche e paraboliche. I problemi che coinvolgono equazioni ellittichesono detti problemi di giuria (jury problems) in quanto la soluzionedipende dalle condizioni applicate all’intero contorno del dominio. Questecondizioni al contorno possono essere ottenute

1. assegnando il valore delle variabili dipendenti u e v lungo il contor-no. Questo tipo di condizione al contorno e detto condizione allaDirichlet;

14

2. assegnando il valore delle derivate delle variabili dipendenti lungo ilcontorno. Questo secondo tipo di condizione al contorno e detto con-dizione alla Neumann;

3. con una sovrapposizione di condizioni alla Dirichlet e alla Neumann.

Figura 1.5: Dominio e contorni della soluzione di un’equazione ellittica indue dimensioni.

In questo primo capitolo ci siamo occupati della classificazione delle equazio-ni differenziali alle derivate parziali attraverso due metodi: quello di Cramere quello degli autovalori. Siamo poi passati a esplicitare le caratteristiche re-lative alle soluzioni di ciascun tipo di equazioni, illustrando con degli esempigrafici (figure 1.2, 1.3, 1.4, 1.5) i diversi domini.Ulteriori esempi e approfondimenti in merito alle equazioni differenziali pos-sono essere trovati in [3], [10] e [13].

15

Capitolo 2

Introduzione alladiscretizzazione

Passiamo ora allo studio del processo di discretizzazione che permette dideterminare il valore di una funzione continua a partire da quello che essaassume su un insieme discreto di punti. Nel fare cio introduciamo diversitipi di stencil, cioe di configurazioni che possono essere utilizzate, nonche leproblematiche relative ai diversi tipi di approccio alle soluzioni e allo studiodella stabilita delle stesse.

2.1 Il processo di discretizzazione

La discretizzazione delle equazioni alle derivate parziali e detta alle dif-ferenze finite, quella delle forme integrali invece ai volumi finiti.La discretizzazione consiste in un processo che permette di risolvere un’equa-zione differenziale o integrale, che coinvolge delle funzioni definite su dominicon la potenza del continuo, attraverso delle approssimazioni, queste ultimedefinite su insiemi discreti di punti o volumi all’interno del dominio origina-rio.Ricordiamo che le soluzioni analitiche costituiscono soluzioni continue de-finite su tutto il dominio, mentre le soluzioni numeriche possono esserecalcolate solo su un insieme discreto di punti che formano un reticolo, comeindicato in figura 2.1.La spaziatura del reticolo non deve essere per forza uniforme. ∆x puo avereun passo diverso da ∆y. Inoltre, anche la distanza tra le stesse coppie di pun-ti in una direzione puo cambiare. In quest’ultimo caso i calcoli non vengonoeseguiti sul piano di partenza bensı su uno ottenuto dal primo attraverso unatrasformazione dei punti in modo tale da ottenere un reticolo regolare, molto

16

Figura 2.1: Esempio di reticolo bidimensionale.

piu facile da trattare dal punto di vista computazionale.Se consideriamo delle equazioni alle derivate parziali possiamo cercare diottenere due tipi di soluzioni. La soluzione analitica, come gia detto, ci per-mette di calcolare i valori delle incognite su ogni punto del dominio, tuttavianon sempre si puo pervenire ad un tale tipo di soluzione, o comunque nonfacilmente. D’altro canto si possono usare delle approssimazioni delle fun-zioni originarie, soppiantando cosı l’equazione alle derivate parziali con unnuovo sistema di equazioni algebriche che possono essere risolte solo sui puntidel reticolo considerato. Il metodo di discretizzazione appena introdotto edetto metodo delle differenze finite in quanto le derivate parziali sonoapprossimate con dei quozienti che fanno riferimento a un numero finito dipunti del reticolo.

2.2 Metodo delle differenze finite

Vediamo ora come e possibile approssimare una derivata parziale con il quo-ziente di una differenza algebrica, cioe usando il metodo alle differenze finite.Generalmente per fare cio si usa lo sviluppo in serie di Taylor della funzioneconsiderata. Con riferimento alla figura 2.1 se supponiamo che u(i, j) sia ilvalore della funzione u nel punto P di coordinate (i, j), allora il valore di unel punto (i + 1, j) puo essere scritto come

ui+1,j = ui,j +

(∂u

∂x

)

i,j

∆x+

(∂2u

∂x2

)

i,j

(∆x)2

2!+

(∂3u

∂x3

)

i,j

(∆x)3

3!+ . . . . (2.1)

Naturalmente nel caso limite in cui i termini della sommatoria sono infinitie l’intervallo ∆x tende a zero, l’equazione (2.1) fornisce il valore esatto della

17

funzione nel punto di coordinate (i+1, j). Risolvendo la precedente equazionein termini della pendenza della funzione u si ottiene

(∂u

∂x

)

i,j

=ui+1,j − ui,j

∆x−

(∂2u

∂x2

)

i,j

∆x

2−

(∂3u

∂x3

)

i,j

(∆x)2

6+ . . . . (2.2)

Consideriamo ora i termini scritti sulla destra. Il primo di essi altro non eche la rappresentazione alle differenze finite della derivata parziale scritta asinistra nella (2.2), mentre i restanti termini corrispondono all’errore di tron-camento che si commette nell’approssimazione. Dal momento che il primotermine ivi presente e dell’ordine di ∆x, allora possiamo dire che un’equa-zione alle differenze finite, come la seguente (2.3), e approssimata al primoordine. (

∂u

∂x

)

i,j

≈ ui+1,j − ui,j

∆x. (2.3)

Piu propriamente possiamo scrivere il tutto come

(∂u

∂x

)

i,j

=ui+1,j − ui,j

∆x+ O(∆x) (2.4)

in cui il termine O(∆x) sta ad indicare l’accuratezza della rappresentazione.

Figura 2.2: Differenza al primo ordine avanti.

Inoltre, la (2.4) e detta differenza al primo ordine avanti o first-orderforward difference in quanto, con riferimento alla figura 2.2, per la suarisoluzione si fa riferimento solo ai valori della funzione a destra del punto dicoordinate (i, j). Il disegno schematizza l’equazione di cui sopra illustrandoquali punti sono coinvolti nella rappresentazione alle differenze finite appenaillustrata, nonche il segno + o −, e l’eventuale peso, dei vari termini.Se, al contrario, sviluppiamo in serie di Taylor la funzione u a partire dal

18

punto di coordinate (i, j) per ricavarne il valore immediatamente alla suasinistra, otteniamo

ui−1,j = ui,j−(

∂u

∂x

)

i,j

∆x+

(∂2u

∂x2

)

i,j

(∆x)2

2−

(∂3u

∂x3

)

i,j

(∆x)3

6+ . . . . (2.5)

Quindi, risolvendo rispetto a (∂u/∂x)i,j

(∂u

∂x

)

i,j

=ui,j − ui−1,j

∆x+ O(∆x). (2.6)

Figura 2.3: Differenza al primo ordine retro.

Con analogo ragionamento visto per la (2.4), facendo attenzione al fatto checi spostiamo verso sinistra, possiamo dire che la precedente equazione costi-tuisce una differenza al primo ordine retro o first-order rearwarddifference.Dalle rappresentazioni di cui sopra e possibile ottenere facilmente un’appros-simazione con un’accuratezza del secondo ordine. Sottraendo la (2.5) dalla(2.1) si ha

ui+1,j − ui−1,j = 2

(∂u

∂x

)

i,j

∆x + 2

(∂3u

∂x3

)

i,j

(∆x)3

6+ . . . , (2.7)

da cui (∂u

∂x

)

i,j

=ui+1,j − ui−1,j

2∆x+ O(∆x)2. (2.8)

Notiamo come in questo caso, a differenza della (2.4) e della (2.6), per calco-lare il valore approssimato della pendenza di u nel punto P si fa riferimentoai due punti ad esso adiacenti mentre l’errore di troncamento risulta essere

19

Figura 2.4: Differenza centrale al secondo ordine.

dell’ordine di (∆x)2. La (2.8) e appunto detta differenza centrale al se-condo ordine o second-order central difference.Espressioni analoghe alla (2.4), (2.6) e (2.8) si ottengono se invece di spo-starci lungo l’asse x nella figura 2.1, lo facciamo lungo l’asse y.Talvolta pero si ha a che fare con delle derivate parziali di ordini superiorial primo. Di seguito vediamo come possono essere rappresentate quelle delsecondo ordine col metodo delle differenze finite.Sommando la (2.1) e la (2.5) si ha

ui+1,j + ui−1,j = 2ui,j +

(∂2u

∂x2

)

i,j

(∆x)2 +

(∂4u

∂x4

)

i,j

(∆x)4

12+ . . . (2.9)

e risolvendo rispetto alla derivata seconda di u rispetto a x

(∂2u

∂x2

)

i,j

=ui+1,j − 2ui,j + ui−1,j

(∆x)2+ O(∆x)2. (2.10)

Figura 2.5: Differenza al secondo ordine centrale per la derivata seconda.

20

Quest’ultima e detta differenza centrale al secondo ordine per la de-rivata seconda, o second-order central second difference, di u rispet-to a x. Il problema e ora quello di trovare le rappresentazioni per le derivatemiste del secondo ordine. Per far cio deriviamo rispetto alla variabile y la(2.1) e la (2.5) ottenendo

(∂u

∂y

)

i+1,j

=

(∂u

∂y

)

i,j

+

(∂2u

∂x∂y

)

i,j

∆x

+

(∂3u

∂2x∂y

)

i,j

(∆x)2

2+

(∂4u

∂3x∂y

)

i,j

(∆x)3

6+ . . . (2.11)

(∂u

∂y

)

i−1,j

=

(∂u

∂y

)

i,j

−(

∂2u

∂x∂y

)

i,j

∆x

+

(∂3u

∂2x∂y

)

i,j

(∆x)2

2−

(∂4u

∂3x∂y

)

i,j

(∆x)3

6+ . . . . (2.12)

Sottraendo quest’ultima alla (2.11) si ha

(∂u

∂y

)

i+1,j

−(

∂u

∂y

)

i−1,j

= 2

[(∂2u

∂x∂y

)

i,j

∆x +

(∂4u

∂3x∂y

)

i,j

(∆x)3

6+ . . .

]

(2.13)ed infine, risolvendo rispetto alla derivata mista di secondo ordine

(∂2u

∂x∂y

)

i,j

=(∂u/∂y)i+1,j − (∂u/∂y)i−1,j

2∆x−

(∂4u

∂3∂y

)(∆x)2

6+ . . . . (2.14)

Quest’ultima puo essere esplicitata tenendo conto delle formule di rappre-sentazione delle derivate parziali del primo ordine viste sopra. Utilizzandoformule analoghe alla (2.10), scritte nella variabile y per i punti di coordinate(i + 1, j) e (i− 1, j), e sostituendole nella (2.14) risulta

(∂2u

∂x∂y

)

i,j

=ui+1,j+1 − ui+1,j−1 − ui−1,j+1 + ui−1,j−1

4∆x∆y+ O[(∆x)2, (∆y)2].

(2.15)

L’errore di troncamento e del secondo ordine tanto rispetto a ∆x quantoa ∆y. Si parla allora di differenza centrale al secondo ordine perla derivata mista (∂2u/∂x∂y)i,j, o anche second-order central mixeddifference.E’ ovvio che le rappresentazioni viste sin ora non sono ne le uniche ne tan-to meno le piu accurate. Come detto in precedenza, possiamo migliorare

21

Figura 2.6: Differenza centrale al secondo ordine per la derivata mista.

la precisione dell’approssimazione andando a considerare termini di ordinevia via maggiore nello sviluppo di Taylor della funzione u. Ad esempio, ladifferenza finita centrale al quarto ordine, o fourth-order centraldifference, per ∂2u/∂x2 e

(∂2u

∂x2

)

i,j

= −−ui+2,j + 16ui+1,j − 30ui,j + 16ui−1,j − ui−2,j

12(∆x)2+ O(∆x)4.

(2.16)Vediamo come nella precedente, a differenza della (2.10), si fa riferimento acinque punti del reticolo, e non a tre.Un’approssimazione del comportamento della funzione sui punti che costi-tuiscono il contorno nel caso di un reticolo unidimensionale puo presentaredei problemi nel caso in cui ci sia un’unica direzione possibile e cioe quellad’allontanamento dal contorno.Con riferimento alla figura 2.7, l’approssimazione al primo ordine e datasemplicemente da (

∂u

∂y

)

1

=u2 − u1

∆y+ O(∆y). (2.17)

Per ottenere una rappresentazione della derivata del secondo ordine col me-todo del punto centrale dovremmo conoscere il valore della funzione in 2’,simmetrico rispetto a 2. Tuttavia nel caso considerato il punto 2’ e al di fuo-ri del dominio della nostra funzione. Per bypassare tale problema si puo farricorso alla condizione al contorno di riflessione secondo la quale la funzioneassume gli stessi valori nei punti simmetrici rispetto al contorno considerato.Ovviamente tale posizione puo risultare anche troppo azzardata.Consideriamo quindi un nuovo approccio, detto approccio polinomiale.Supponiamo quindi che la funzione u sul contorno possa essere espressa informa polinomiale, cioe

u = a + by + cy2. (2.18)

22

Figura 2.7: Esempio di reticolo unidimensionale in allontanamento dal bordo.

E’ facile dimostrare che al punto 1 y = 0. Allora

u1 = a . (2.19)

Al punto 2 y = ∆y quindi

u2 = a + b∆y + c(∆y)2 . (2.20)

Al punto 3 y = 2∆y e

u3 = a + b(2∆y) + c(2∆y)2. (2.21)

Da queste ultime tre equazioni possiamo calcolare b come

b =−3u1 + 4u2 − u3

2∆y. (2.22)

Deriviamo ora la (2.18) rispetto a y ottenendo

∂u

∂y= b + 2cy. (2.23)

Quindi nel punto 1, in cui y = 0, si otterra(

∂u

∂y

)

1

= b (2.24)

e ricordando la (2.22) possiamo scrivere(

∂u

∂y

)

1

=−3u1 + 4u2 − u3

2∆y. (2.25)

23

Quest’espressione e detta differenza finita lungo una direzione poicheper la sua risoluzione si fa riferimento solo ai punti che si trovano da un latodel contorno. Concludiamo determinandone il grado di approssimazione.Scriviamo lo sviluppo in serie di Taylor della (2.18)

u(y) = u1 +

(∂u

∂y

)

1

y +

(∂2u

∂y2

)

1

y2

2+

(∂3u

∂y3

)

1

y3

6+ . . . (2.26)

i cui primi tre termini possono essere posti in corrispondenza proprio conquelli della (2.18). Di conseguenza quest’ultima (cosı come u1, u2, u3) eun’approssimazione al terzo ordine. Tenendo pero conto del ∆y presente aldenominatore della (2.25) si conclude che questa rappresenta un’approssi-mazione del secondo ordine che puo essere messa, piu correttamente, nellaseguente forma

(∂u

∂y

)

1

=−3u1 + 4u2 − u3

2∆y+ O(∆y)2. (2.27)

2.3 Equazioni alle differenze

Nel precedente paragrafo abbiamo visto com’e possibile approssimare unaderivata parziale con un quoziente alle differenze finite. Nel momento in cuiin una data equazione differenziale sostituiamo tutte le derivate con dellerappresentazioni alle differenze finite otteniamo un’equazione algebrica cheprende il nome di equazione alle differenze. Cerchiamo di tradurre inpratica quanto esposto finora studiando un’equazione alle derivate parzialipropria della fluidodinamica.Consideriamo allora l’equazione della conduzione del calore messa nella forma

∂T

∂t= α

∂2T

∂x2. (2.28)

Da quanto espresso nel precedente capitolo e facile dimostrare che si trattadi un equazione parabolica 1. Essa ammette dunque una soluzione marcianterispetto al parametro t.Consideriamo una griglia costruita come in figura 2.8.Siano cioe i ed ∆x l’indice e l’intervallo sulle ascisse, n e ∆t quelli sulle ordina-te. Evidenziamo il fatto che n caratterizza la variabile marciante ponendolo,per convenzione, come apice nell’espressione del quoziente alle differenze fi-nite. Rappresentiamo la derivata parziale di T fatta rispetto al tempo con

1Si rimanda all’appendice A.1 per tale dimostrazione

24

Figura 2.8: Griglia bidimensionale.

un’espressione di tipo forward analoga alla (2.4)

(∂T

∂t

)n

i

=T n+1

i − T ni

∆t−

(∂2T

∂t2

)n

i

∆t

2+ . . . , (2.29)

mentre sostituiamo la derivata rispetto alla x con una differenza centrale deltipo della (2.10), cioe

(∂2T

∂x2

)n

i

=T n

i+1 − 2T ni + T n

i−1

(∆x)2−

(∂4T

∂x4

)n

i

(∆x)2

12+ . . . . (2.30)

A partire dalla (2.28), portando ambo i termini al primo membro si ottiene

∂T

∂t− α

∂2T

∂x2= 0 (2.31)

da cui, sfruttando la (2.29) e la (2.30), si ha

∂T

∂t− α

∂2T

∂x2= 0 =

T n+1i − T n

i

∆t− α

T ni+1 − 2T n

i + T ni−1

(∆x)2

+

[(∂2T

∂t2

)n

i

∆t

2+ α

(∂4T

∂x4

)n

i

(∆x)2

12+ · · ·

].(2.32)

In ultima analisi, troncando la serie, la (2.32) puo essere scritta come

T n+1i − T n

i

∆t= α

T ni+1 − 2T n

i + T ni−1

(∆x)2. (2.33)

Quest’ultima e un’equazione alle differenze che approssima la (2.28).E’ questo un concetto fondamentale poiche l’equazione alle differenze e solo

25

un’equazione algebrica che e risolvibile soltanto su un certo insieme di punticostituenti un reticolo. I valori ivi calcolati possono non essere identici aquelli ricavati per mezzo di una soluzione analitica, tuttavia tendono a que-sti ultimi al diminuire del passo di discretizzazione dal momento che l’erroredi troncamento, come si nota dalla (2.32), dipende da ∆x. A tal proposi-to se il passo e opportunamente piccolo si dice che la rappresentazione alledifferenze finite e consistente con l’equazione differenziale data entro l’erroredi troncamento. Tuttavia non basta sostituire le derivate con dei quozientialle differenze finite ed eventualmente far tendere a zero il passo per otteneredei buoni risultati. Bisogna infatti tener conto della stabilita dell’algoritmousato per calcolare il valore della funzione sulla griglia, nonche del correttotrattamento delle condizioni al contorno del problema. Naturalmente qualo-ra ∆x fosse nullo otterremmo un risultato esatto. Cio non puo mai verificarsiin pratica dal momento che il calcolatore ha comunque una precisione finita,legata al numero finito di bit in cui vengono salvati i numeri e agli errori diarrotondamento. Potremmo allora pensare di mantenere costante l’interval-lo e contemporaneamente far diminuire rapidamente il passo. Questo perocomporterebbe un aumento altrettanto repentino del numero di mesh (magliedella griglia), quindi dei tempi macchina ed anche della memoria necessariaper salvare i valori.Capiamo allora come sia opportuno trovare, di volta in volta, un punto d’e-quilibrio tra la precisione che vogliamo ottenere e le problematiche relativeall’utilizzo dei calcolatori.

2.4 Approcci alle soluzioni

Un ulteriore passo in avanti nello studio dei metodi di risoluzione delle equa-zioni differenziali attraverso processi di discretizzazione consiste nel conside-rare i vari possibili approcci al problema. I due metodi piu usati consistonoin un approccio di tipo esplicito o implicito. Per capire in cosa consistonoconsideriamo nuovamente la (2.28), che esprime la conduzione del calore nelcaso monodimensionale.Nel precedente paragrafo abbiamo visto come la (2.28) possa essere messasottoforma di equazione alle differenze come

T n+1i − T n

i

∆t= α

T ni+1 − 2T n

i + T ni−1

(∆x)2(2.34)

da cui, banalmente, si ottiene

T n+1i = T n

i + α∆t

(∆x)2(T n

i−1 − 2T ni + T n

i−1). (2.35)

26

A questo punto ricordiamo che la (2.28) e un’equazione di tipo parabolicoche ammette una soluzione di tipo marciante. In questo caso la variabile cheevolve e costituita dal tempo t. Per l’appunto, come evidenziato dall’equa-zione (2.35), i valori della variabile T sulla mesh con indice temporale pari an + 1 sono determinati a partire dai valori della stessa funzione sulla meshprecedente. Schematizziamo questa situazione nella figura 2.9.

Figura 2.9: Esempio di approccio esplicito.

Nell’area tratteggiata sono evidenziati i punti interessati dall’equazione dicui sopra.Considerando ancora la (2.35), notiamo che tutti i valori (noti) della meshcon indice n compaiono a destra dell’uguaglianza, mentre l’unica incognitapresente e rappresentata dal valore di T nel punto di coordinate (i, n + 1).Siamo di fronte quindi ad un’equazione con un’unica incognita che puo essererisolta conoscendo i valori di T negli istanti di tempo precedenti.Per definizione, in un approccio esplicito ogni equazione alle differenzecontiene una sola incognita ed inoltre puo essere risolta esplicitamente inmaniera abbastanza semplice.Tuttavia la (2.34) non e l’unica equazione alle differenze che si puo otteneredalla (2.28). Ad esempio, possiamo svilupparla in termini di proprieta medietra i livelli con indice temporale n ed n + 1, ottenendo

T n+1i − T n

i

∆t= α

12(T n+1

i+1 + T ni+1) + 1

2(−2T n+1

i − 2T ni ) + 1

2(T n+1

i−1 + T ni−1)

(∆x)2.

(2.36)Questo tipo di differenziazione, usato spesso per risolvere problemi che coin-volgono equazioni paraboliche, e detto Crank-Nicolson [6]. In tal casosono presenti tre incognite legate ai valori della funzione sui punti del livel-

27

lo n + 1. Capiamo allora perche, in questo secondo caso, occorre che tuttele equazioni che interessano i punti del reticolo a n + 1 devono essere risol-te simultaneamente. La (2.36) altro non e che un esempio di un approccioimplicito al problema considerato. Per definizione infatti, in un approccioimplicito le incognite devono essere ottenute dalla risoluzione simultaneadelle equazioni alle differenze applicate a tutti i punti della griglia dispostisu una data mesh. E’ cosı spiegato il motivo per cui nell’approccio implicitosi fa ricorso al calcolo matriciale.Nella figura 2.10 e schematizzata l’area di influenza di un approccio implicitosu una griglia avente sette punti su ogni livello.

Figura 2.10: Esempio di approccio implicito.

Esplicitiamo ora la (2.36) portando a sinistra tutti i termini che si trovanosul livello temporale n + 1 e a destra quelli che appartengono alla mesh pre-cedente.

α∆t

2(∆x)2T n+1

i−1 −[1 +

α∆t

(∆x)2

]T n+1

i +α∆t

2(∆x)2T n+1

i+1 =

− T ni −

α∆t

2(∆x)2(T n

i+1 − 2T ni + T n

i−1). (2.37)

Ponendo dunque

A =α∆t

2(∆x)2, (2.38)

B = 1 +α∆t

(∆x)2, (2.39)

28

Ki = −T ni −

α∆t

2(∆x)2(T n

i+1 − 2T ni + T n

i−1), (2.40)

si ottieneAT n+1

i−1 −BT n+1i + AT n+1

i+1 = Ki (2.41)

in cui Ki e una quantita nota in quanto fa riferimento ai valori di T al passoprecedente della griglia. Con riferimento alla figura 2.10, applichiamo la(2.41) ai punti della griglia da 2 a 6 sul livello temporale n + 1.Per il punto 2 si ottiene

AT1 −BT2 + AT3 = K2. (2.42)

Dal momento che l’equazione parabolica e definita sul contorno per x = 1 ex = 7, allora sia K2 che T1 sono quantita note, che possiamo portare perciotutte al secondo membro

−BT2 + AT3 = K2 − AT1 (2.43)

e sostituire col valore K ′2

−BT2 + AT3 = K ′2. (2.44)

Analogamente, procedendo con gli altri punti nel punto 3 si ha

AT2 −BT3 + AT4 = K3; (2.45)

nel punto 4AT3 −BT4 + AT5 = K4; (2.46)

nel punto 5AT4 −BT5 + AT6 = K5; (2.47)

nel punto 6AT5 −BT6 + AT7 = K6. (2.48)

Per l’ultima equazione vale il discorso fatto per la (2.42) in quanto T7 ecalcolato sul contorno, quindi e noto. Possiamo allora riscrivere la (2.48)come

AT5 −BT6 = K6 − AT7 = K ′6. (2.49)

Le equazioni dalla (2.44) alla (2.47), con la (2.49) formano un sistema di cin-que equazioni in altrettante incognite; in forma matriciale possiamo scrivere

−B A 0 0 0A −B A 0 00 A −B A 00 0 A −B A0 0 0 A −B

T2

T3

T4

T5

T6

=

K ′2

K3

K4

K5

K ′6

(2.50)

29

La matrice dei coefficienti e una matrice tridiagonale avente, cioe, ele-menti non nulli solo lungo le tre diagonali centrali. Questo tipo di matricepuo essere risolta facendo ricorso all’algoritmo di Thomas 2.Notiamo che l’equazione di cui ci siamo occupati finora, cioe la (2.28), e un’e-quazione differenziale alle derivate parziali di tipo lineare. Anche le equazionialle differenze dedotte, sia con l’approccio esplicito (2.35) che con quello im-plicito (2.36), sono di tipo lineare.Assumiamo ora che il fattore di diffusione termica α dipenda dalla tempera-tura T , di modo che l’equazione (2.28) non risulti essere lineare. Scriviamoallora

∂T

∂t= α(T )

∂2T

∂x2. (2.51)

Analogamente alla (2.35), per l’approccio esplicito otteniamo

T n+1i = T n

i + α(T ni )

∆t

(∆x)2(T n

i+1 − 2T ni + T n

i−1). (2.52)

Quest’equazione risulta essere ancora lineare nell’unica incognita T n+1i poi-

che α e stimata a partire dal valore di T sul livello n, che e ben noto.Se applichiamo alla (2.51) la Crank-Nicolson otteniamo un’equazione analogaalla (2.36), tenendo conto che α e ora rappresentato da 1/2[α(T n+1

i +α(T ni )].

In questo caso, risolvendo l’equazione, otteniamo dei termini dati dal pro-dotto delle variabili dipendenti α e T calcolate al passo n + 1. Ne consegueun’equazione alle differenze non lineare. Per risolvere allora un tale problemadovremmo risolvere un sistema formato da equazioni non lineari, il che non ecertamente semplice. Per superare tale ostacolo si puo approssimare il valoredi α sulla mesh n+1 a quello ch’essa assume nella precedente.Dopo aver abbondantemente esposto ambedue gli approcci facciamo delleconsiderazioni finali. Notiamo allora come gli incrementi ∆x e ∆t appaionosia nella (2.35) che nella (2.36). Per quanto riguarda l’approccio esplicito,una volta fissato ∆x il valore che ∆t puo assumere non e arbitrario, ma deveessere minore o uguale di un certo valore prescritto dal criterio di stabilita.Se ∆t risulta maggiore di questo limite il processo diviene subito instabile.Per ovviare a cio si preferisce un passo molto piccolo con un conseguenteaumento del numero di processi che la macchina deve eseguire per ottenere

2Il metodo di Thomas consiste nel risolvere un sistema di equazioni lineari in formatridiagonale attraverso un processo iterativo. In particolare si trasforma il sistema origi-nario in uno bidiagonale eliminando con opportune moltiplicazioni e sottrazioni i terminidella diagonale inferiore e facendo delle sostituzioni di variabile. A questo punto, partendodall’ultima equazione, avente una sola incognita, si procede a ritroso risolvendo il sistema.Per maggiori informazioni si consulti [18].

30

la soluzione. Al contrario i processi impliciti non prevedono la presenza dialcuna restrizione sul valore da assegnare a ∆t, anzi talvolta si incontrano deisistemi che risultano essere incondizionatamente stabili, nel senso che sonostabili per qualsiasi valore di ∆t. In questo caso il tempo impiegato dai cal-colatori sara dunque minore. Ci sono pero anche dei risvolti negativi legatialla scelta di un ∆t molto ampio. Innanzi tutto con l’aumentare di ∆t crescel’errore di troncamento legato alla relativa equazione alle differenze per laderivata temporale; inoltre un aumento del passo va a discapito della preci-sione della descrizione del sistema. Capiamo bene che non c’e un approcciomigliore dell’altro in quanto entrambe, come esposto, hanno dei pregi e deidifetti; percio bisogna far ricorso all’uno o all’altro in base alle situazioni incui ci si trova.

2.5 Analisi degli errori e studio della stabilita

Come gia osservato, utilizzando un approccio esplicito si puo pervenire aduna soluzione che non risulta essere stabile qualora il valore del passo dellavariabile marciante supera un determinato valore soglia. Quest’ultimo e bendeterminato dallo studio della stabilita effettuato sull’equazione alle dif-ferenze sotto esame. Il comportamento di una data equazione e strettamentelegato agli errori concernenti il processo tramite il quale si ottiene una suasoluzione. E’ facilmente intuibile allora che se tale errore durante il processocresce in continuazione da un passo al successivo, la soluzione non puo esse-re stabile. Viceversa se esso resta costante, o ancora meglio diminuisce, lasoluzione puo avere un comportamento stabile.Gli errori che riguardano il processo di discretizzazione sono di due diffe-renti tipi. Il primo e l’errore di discretizzazione propriamente detto inquanto riguarda l’errore che si commette approssimando la soluzione esatta,ottenuta analiticamente, con quella derivante dall’equazione alle differenzefinite. Questo e dato dall’ errore di troncamento e da altri errori che possonoriguardare le condizioni al contorno. Il secondo tipo di errore e detto erroredi arrotondamento ed e relativo agli arrotondamenti alle cifre significativecompiuti dal calcolatore.Per capire meglio l’argomento della discussione, studiamo la stabilita dellagia nota equazione differenziale lineare (2.28), riportata qui di seguito

∂T

∂t= α

∂2T

∂x2.

31

La relativa equazione alle differenze ottenuta col metodo esplicito e data dalla(2.34)

T n+1i − T n

i

∆t= α

T ni+1 − 2T n

i + T ni−1

(∆x)2.

Se poniamo

A = soluzione analitica dell’equazione alle derivate parzialiD = soluzione esatta dell’equazione alle differenze

N = soluzione numerica data da un computer reale avente precisione finita

allora

Errore di discretizzazione = A - DErrore di arrotondamento = ε = N - D.

Dall’ultima si ricava cheN = D + ε. (2.53)

Questa, essendo una soluzione dell’equazione alle differenze, rispetta la (2.34).Quindi, da tale equazione, con N al posto di T , si ottiene

Dn+1i + εn+1

i −Dni − εn

i

α∆t=

Dni+1 + εn

i+1 − 2Dni − 2εn

i + Dni−1 + εn

i−1

(∆x)2. (2.54)

A questo punto, ricordando che D e per definizione la soluzione esatta del-l’equazione alle differenze, possiamo scrivere

Dn+1i −Dn

i

α∆t=

Dni+1 − 2Dn

i + Dni−1

(∆x)2. (2.55)

Essendo l’equazione differenziale di partenza lineare, possiamo sottrarre in-fine la (2.55) dalla (2.54)

εn+1i − εn

i

α∆t=

εni+1 − 2εn

i + εni−1

(∆x)2(2.56)

notando come anche l’errore (di arrotondamento) ε rispetta l’equazione alledifferenze.Solitamente gli errori εi sono presenti in ogni fase del processo risolutivodell’equazione alle differenze per cui la relativa soluzione risulta essere stabilese tali εi non aumentano, o meglio diminuiscono, procedendo da un passo al

32

successivo. Possiamo allora dire che una soluzione risulta essere stabile seviene rispettata la seguente condizione

∣∣∣∣εn+1i

εni

∣∣∣∣ ≤ 1. (2.57)

Dobbiamo allora cercare di vedere quando questa condizione viene rispetta-ta. Per fare cio consideriamo il fatto che il valore dell’errore εi puo variarein maniera del tutto casuale lungo il dominio sul quale si vuole risolvere l’e-quazione.Denotando tale lunghezza con L possiamo immaginare che l’andamento di εin funzione di x venga schematizzato come in figura 2.11.

Figura 2.11: Esempio di andamento di ε(x)

Come si vede, l’errore relativo ai punti caratterizzati da x = ±L/2 e nullo dalmomento che ci troviamo agli estremi dove il valore della funzione e fissatodalle condizioni al contorno.Sfruttando poi la teoria di Fourier possiamo pensare di poter esprimereanaliticamente il comportamento di ε(x) attraverso una serie, cioe

ε(x) =∑m

Ameikmx (2.58)

in cui Am rappresenta la fase e km il vettore d’onda dato da

km =2π

λ(2.59)

o equivalentemente, ricordando che il dominio L deve essere pari a un numerointero di λ,

km =2π

(L/m)=

(2π

L

)m. (2.60)

Capiamo quindi come l’indice m, presente nel numero d’onda k, indichi ilnumero delle onde che ricadono all’interno del dominio L.

33

Passiamo ora a considerare la massima e la minima lunghezza d’onda permes-sa. La prima, come ci si aspetta, e data dalla lunghezza stessa del dominiodella funzione, quindi λMAX = L (con m = 1). La minima lunghezza d’ondae quella che si annulla in tre punti adiacenti del reticolo. Con riferimentoalla figura 2.12 possiamo allora dire che λmin = 2∆x.

Figura 2.12: Lunghezza d’onda massima e minima

Immaginiamo di avere un reticolo di N + 1 punti distanziati con spaziatura∆x = L/N , avendo N intervalli sulla lunghezza L del dominio. Allora nelcaso in cui λ = λmin otteniamo

km =2π

(2L/N)=

(2π

L

)N

2. (2.61)

Ne consegue che nella sommatoria della (2.58) corre da 1 a N/2

ε(x) =

N/2∑m=1

Ameikmx. (2.62)

La precedente equazione fornisce una stima dell’errore in funzione della solavariabile spaziale x, tuttavia siamo interessati al valore di ε in due istanticonsecutivi percio immaginiamo che l’ampiezza Am sia funzione anche di t.Supponendo che tale ampiezza abbia un andamento esponenziale

ε(x, t) =

N/2∑m=1

eateikmx (2.63)

in cui a e una costante.Notiamo che, essendo l’equazione alle differenze (2.34) di tipo lineare e l’er-rore ε una soluzione della stessa, sostituendo la (2.63) nella (2.56), il com-portamento di ogni singolo termine e equivalente a quello della stessa serie,percio

εm(x, t) = eateikmx. (2.64)

34

Senza perdere in generalita, sostituiamo l’espressione di ε data dalla prece-dente (2.64) nell’equazione (2.56)

ea(t+∆t)eikmx − eateikmx

α∆t=

eateikm(x+∆x) − 2eateikmx + eateikm(x−∆x)

(∆x)2. (2.65)

Dividendo ambo i membri per eateikmx si ha

ea∆t − 1

α∆t=

eikm∆x − 2 + e−ikm∆x

(∆x)2, (2.66)

quindi

ea∆t = 1 +α∆t

(∆x)2(eikm∆x − 2 + e−ikm∆x). (2.67)

Ricordando poi l’espressione del coseno

cos δ =eδ + e−δ

2

scriviamo

ea∆t = 1 +2α∆t

(∆x)2[cos(km∆x)− 1]. (2.68)

Facendo ricorso a questo punto ad un’altra equazione trigonometrica

sin2 δ

2=

1− cos δ

2,

dalla (2.68) otteniamo

ea∆t = 1− 4α∆t

(∆x)2sin2

(km∆x

2

). (2.69)

Inoltre ea∆t puo anche essere visto come

ea∆t =ea(t+∆t)eikmx

eateikmx=

εn+1i

εni

. (2.70)

Sfruttando queste due ultime equazioni possiamo mettere la (2.57) nellaforma seguente

∣∣∣∣εn+1i

εni

∣∣∣∣ = |ea∆t| =∣∣∣∣1−

4α∆t

(∆x)2sin2

(km∆x

2

)∣∣∣∣ ≤ 1 (2.71)

in cui indichiamo col termine di fattore di amplificazione G la quantita

G ≡∣∣∣∣1−

4α∆t

(∆x)2sin2

(km∆x

2

)∣∣∣∣ .

Affinche la soluzione numerica della (2.28), ottenuta con un approccio espli-cito al metodo delle differenze finite, sia stabile occorre allora che G ≤ 1.Possono presentarsi due casi:

35

1. 1− 4α∆t(∆x)2

sin2(

km∆x2

) ≤ 1 ⇒ 4α∆t(∆x)2

sin2(

km∆x2

) ≥ 0.

La condizione espressa dall’ultima equazione risulta essere sempre ri-spettata dal momento che il fattore 4α∆t/(∆x)2 e sempre positivo;

2. 1− 4α∆t(∆x)2

sin2(

km∆x2

) ≥ −1 ⇒ 4α∆t(∆x)2

sin2(

km∆x2

)− 1 ≤ 1

da cui si ricava cheα∆t

(∆x)2≤ 1

2. (2.72)

La precedente (2.72) esprime in termini di ∆x e ∆t il criterio di stabilitaper la soluzione della (2.34) nel senso che finche tale condizione vieneverificata l’errore ε non cresce e la soluzione ′′non esplode′′; viceversanel caso in cui il ∆t scelto non la rispetti.

Il procedimento appena presentato prende il nome di metodo di stabilitaalla von Neumann ed e spesso utilizzato per lo studio della stabilita delleequazioni alle differenze lineari.

Presentiamo di seguito, attraverso un altro esempio, un metodo che per-mette di studiare la stabilita delle soluzioni di equazioni iperboliche comequella delle onde al primo ordine data da

∂u

∂t+ c

∂u

∂x= 0. (2.73)

Utilizzando per la derivata spaziale una differenza centrale del tipo

∂u

∂x=

uni+1 − un

i−1

2∆x(2.74)

e per la derivata temporale una differenza al primo ordine avanti, alloral’equazione delle differenze per la (2.73) sara data dalla seguente

un+1i − un

i

∆t= −c

uni+1 − un

i−1

2∆x. (2.75)

L’equazione trovata e anche detta forma esplicita di Eulero. Applicandoa tale equazione il metodo di von Neumann per lo studio della stabilita si notache essa e incondizionatamente instabile, cioe risulta essere instabile perogni valore di ∆t.Immaginiamo allora di poter esprimere u(t) nel punto di coordinate (i, n)come media dei valori di u sui due punti ad esso adiacenti sulla stessa mesh,cioe

u(t) =1

2(un

i+1 + uni−1).

36

Quindi riscriviamo la derivata temporale come differenza al primo ordine conil valore di u(t) = un

i ottenuto precedentemente

∂u

∂t=

un+1i − 1

2(un

i+1 + uni−1)

∆t. (2.76)

Sostituendo la (2.74) e la (2.76) nella (2.73), e risolvendo si perviene a

un+1i =

uni+1 + un

i−1

2− c

∆t

∆x

uni+1 − un

i−1

2. (2.77)

Tale metodo di approssimazione e detto metodo di Lax. Se poi supponiamoche l’errore abbia la gia nota forma εm(x, t) = eateikmt ed andiamo a sostituirlonella (2.77) troviamo un’espressione del fattore di amplificazione del tipo

ea∆t = cos(km∆x)− iC sin(km∆x) (2.78)

in cui C = c∆t/∆x e il numero di Courant.Il requisito di stabilita |ea∆t| ≤ 1 diviene allora

C = c∆t

∆x≤ 1. (2.79)

Capiamo bene che, secondo quanto affermato dalla precedente (2.79), affin-che una soluzione della (2.77) risulti essere stabile occorre che ∆t ≤ ∆x/c.La (2.79) esprime allora la condizione di stabilita ed e definita come condi-zione di Courant-Friedrichs-Lewy (CFL).La CFL risulta essere per giunta anche la condizione di stabilita per lesoluzioni dell’equazione delle onde al secondo ordine

∂2u

∂t2= c2∂2u

∂x2. (2.80)

Proprio a partire da quest’ultima equazione possiamo ricavare un’interpreta-zione fisica della (2.79). Consideriamo a tal proposito le linee caratteristichedella (2.80). Esse sono date da

x =

{ct cammino destro

−ct cammino sinistro(2.81)

Naturalmente si possono presentare due situazioni, schematizzate nelle im-magini della figura 2.13, che dipendono dal valore del numero di Courant C.Nella figura 2.13 si nota come il punto b sia dato, in entrambe i grafici, dal-l’intersezione tra la linea caratteristica destra, passante per il punto i− 1, equella sinistra, passante per i + 1.

37

Figura 2.13: Esempi di situazione stabile (a) ed instabile (b).

38

Nel disegno (a) della figura 2.13 supponiamo C < 1 ed in quello (b) C > 1.E’ ovvio che se C = 1 allora dalla (2.79) segue che

∆tC=1 =∆x

c. (2.82)

Dalle equazioni (2.81) si ricava invece, considerando un incremento finito dix e t

∆t = ±∆x

c. (2.83)

A questo punto notiamo che la (2.82) ha a che fare con la condizione distabilita espressa dalla (2.79), mentre la (2.83) prende le mosse dalle lineecaratteristiche indicate nelle (2.81).Dal momento che ∆t e il medesimo allora capiamo come ∆tC=1 sia propriola distanza del punto b dall’asse x.Se C < 1, allora dalla (2.79) segue che ∆tC<1 < ∆tC=1 ed il punto d, dicoordinate (i, t + ∆tC<1) si trova al di sotto di b. In questo caso il dominionumerico, dato dal triangolo abc, contiene quello analitico, dato dall’interse-zione della retta con t = n con le parallele alle linee caratteristiche passantiper il punto d.Se C > 1 invece si ha che ∆tC>1 > ∆tC=1 ed il punto d, di coordinate (i,t + ∆tC>1) si trova al di sopra di b. In tal caso il dominio numerico noncontiene piu quello analitico.Perveniamo cosı all’interpretazione della CFL cercata, secondo la quale af-finche una funzione sia stabile il dominio numerico deve contenere quelloanalitico.In ultima analisi consideriamo l’accuratezza del processo. Il comportamentodella funzione nel punto b dipende solo dai punti compresi tra le due lineecaratteristiche tuttavia la soluzione numerica attinge alle informazioni deipunti i− 1 e i+1 che si trovano al di fuori del dominio analitico. Anche se ilprocesso risulta essere stabile la soluzione numerica puo essere poco accuratarispetto al valore analitico cosicche si preferisce considerare un C ≤ 1, mamolto prossimo all’unita.

39

2.6 Risoluzione numerica dell’equazione di

propagazione del calore in una dimensio-

ne

001. #include <stdio.h>002. #include <math.h>003. #include <stdlib.h>

Per prima cosa immettiamo nel codice le librerie che contengono alcunefunzioni standard utilizzate.

004. #define PI 3.141592653858979005. #define P 9006. #define K 50007. #define N 1000

Definiamo il valore della costante π nonche il numero di mesh che vo-gliamo campionare e salvare in singoli file. Introduciamo poi come variabiliglobali i punti della griglia da prendere lungo l’asse delle ascisse e delle ordi-nate. Naturalmente questi, cosı come le dimensioni dei rispettivi intervalli,devono essere scelti opportunamente in modo tale da verificare la condizionedi stabilita dell’algoritmo. Nel nostro caso deve risultare

∆t

(∆x)2≤ 1

2.

Se ad esempio prendessimo dei valori troppo grandi per N, mantenendo fissitutti gli altri parametri, noteremmo come i risultati, fornitici dal program-ma, scarterebbero in maniera quanto mai brusca rispetto quelli vicini poichel’algoritmo risulterebbe ivi instabile. Per ovviare a cio facciamo in modoche l’algoritmo si blocchi nel caso in cui la condizione di stabilita non vienerispettata.

008. double H0(double x);009. double HSx(double t);010. double HDx(double t);011. double ** allocate(int N_, int K_);

Vengono dichiarati i prototipi delle funzioni H0, HSx, HDx e **allocateusate nel codice. Le espressioni delle stesse funzioni sono riportate per co-modita subito sotto il main.

40

012. int main(void)013. {014. double **H, *X, *T;015. double deltax, deltat;016. double x0, x1, t0, t1;017. int k,n,i;018. int M[P] = {0, 99, 199, 299, 399, 499, 699, 799, 999};019. double alpha;020. double parametro;021. char mesh[50];

Si inizia il main definendo tutte le varie variabile presenti nel corpo delprogramma:**H, *X e *T sono puntatori a double;deltax e deltat sono double e rappresentano il passo della griglia rispettiva-mente sulle ascisse e sulle ordinate;x0, x1, t0 e t1 sono gli estremi della griglia, quindi essendo dei reali assegna-mo loro la dimensione dei double;k, n ed i sono degli indici di tipo intero;M[ ] e un vettore di interi di dimensione P i cui elementi sono riportati traparentesi graffe;alpha e un double che rappresenta la costante di diffusione termica;parametro e dato dal rapporto di deltat su deltax2 e serve per controllare lastabilita dell’algoritmo;mesh e un vettore che puo contenere fino a 50 elementi di tipo char, costi-tuisce dunque una stringa di caratteri.

022. X = (double*)calloc(K, sizeof(double));023. T = (double*)calloc(N, sizeof(double));

Con queste due linee di codice si alloca memoria per i vettori in cui sal-vare i punti della griglia che si trovano sull’asse spaziale e temporale.

024. x0 = 0.;025. x1 = 6.283185307;026. t0 = 0.;027. t1 = 1.;028. alpha=1.0;

029. deltax=(x1 -x0)/(K-1);030. deltat=(t1 -t0)/(N-1);

Si definiscono gli estremi degli intervalli lungo l’asse x e t in modo taleche venga rispettata, alla luce di K ed N fissati precedentemente, la condi-zione di stabilita data da (deltat/deltax2) ≤ 1/2. Il parametro alpha, cherappresenta la costante di diffusione termica, e fissato inizialmente a uno.

41

Vengono quindi calcolati i passi deltax e deltat.

031. parametro = deltat/((deltax)*(deltax));032. if (parametro > 0.5)033. {034. printf("L’algoritmo non risulta essere stabile poiche’ il"035. "parametro deltat/deltax^2 e’ %lf>1/2.\n\n",parametro);036. system("PAUSE");037. return 0;038. }

Con questo pezzo di codice si effettua un controllo sulla stabilita in quan-to si verifica che il parametro sia minore o uguale a 1/2. Se cio e verificatoallora il compilatore passa alle linee che seguono, altrimenti stampa in videoun messaggio di errore e termina l’applicazione.

039. for(k=0; k<K; ++k) X[k]= x0 + k*deltax;040. for(n=0; n<N; ++n) T[n]= t0 + n*deltat;

041. H = allocate(N, K);

Le prime due linee di codice di cui sopra servono a riempire i vettori X[k]e T[n] con i rispettivi punti del reticolo lungo i rispettivi assi.Nell’ultima linea si richiama la funzione allocate in modo tale da creare unopportuno spazio di memoria per i valori di H definiti sulla griglia di N*Kpunti.

042. for(k=0;k<K;++k) H[0][k] = H0(X[k]);

043. for(n=0;n<N;++n) H[n][0] = HSx(T[n]);044. for(n=0;n<N;++n) H[n][K-1] = HDx(T[n]);

Lo scopo del primo ciclo for e quello di riempire la prima mesh con ivalori fornitici dalla condizione iniziale, mentre gli altri due servono a fissarele condizioni al contorno lungo i bordi verticali della griglia.Da quanto riportato nel codice delle funzioni, si e scelto un andamento sinu-soidale per la variabile H(x,t) al tempo t=0, mentre supponiamo che essa siannulli ai bordi per qualsiasi istante di tempo.

045. for(n=0; n<N-1; ++n)046. {047. for(k=1;k<K-1;++k)048. {

42

049. H[n+1][k] = H[n][k]+alpha*deltat*pow(deltax,-2)*(H[n][k+1]-2*H[n][k]+H[n][k-1]);

050. }051. }

I due cicli for innestati permettono di muoverci su tutta la griglia primalungo una mesh e poi passando alla successiva assegnando ad ogni punto ilvalore della funzione H.Si noti come questo e legato ai valori che H assume sui tre punti ad essosottostanti nella mesh precedente, come mostrato dall’area tratteggiata infigura 2.9.Bisogna fare particolare attenzione nel riempire in maniera opportuna la gri-glia, cioe bisogna farlo rispettando la condizione iniziale e i bordi. Eccospiegato il motivo per cui l’indice k corre da 1 a K-1 (riga 047).

052. FILE *pfile;053. pfile = fopen("dati.txt","w");054. for(i=0;i<n+1;++i)055. {056. for(k=0;k<K;++k)057. {058. fprintf(pfile,"%f %f %f\n",X[k],T[i],H[i][k]);059. }060. }061. fclose(pfile);

Tramite un puntatore a file si apre, o crea nel caso in cui non sia statogia generato, il file ′′dati tot.txt′′ in cui vengono riportati i vari punti dellagriglia ed i relativi valori della funzione H. Il file di testo cosı generato puopoi essere usato per graficare la funzione H(x, t) in 3D con Mathematica 5.

062. for(i=0;i<P;++i)063. {064. sprintf(mesh, "dati_%d_meshL.txt", M[i]+1);065. pfile = fopen(mesh,"w");066. printf("Mesh %d al tempo t = %lf\n", M[i]+1, T[M[i]]);067. for(k=0;k<K;++k)068. {069. fprintf(pfile,"%f %f\n",X[k], H[M[i]][k]);070. }071. fclose(pfile);072. }073. printf("\n");074. system("PAUSE");075. return 0;076. }

43

Abbiamo fatto uso all’istruzione sprintf per scrivere nella stringa deno-tata con l’appellativo ′′mesh′′ il nome del file che andiamo di volta in volta agenerare. Tale stringa viene usata quindi dall’istruzione successiva per apri-re tale file. Dopo di che il compilatore passa ad eseguire un’operazione discrittura sul video e un ciclo for. Viene dunque stampato in video il tempot corrispondente alla mesh definita a partire dal vettore M[ ]. Lo scopo delciclo for e quello di salvare nel file appena aperto i punti della griglia lungol’asse x ed i relativi valori di H(x). I file cosı generati possono essere usatiper generare dei grafici bidimensionali con GNUplot. Si chiude infine il file.txt ed il ciclo for piu esterno, quindi il main.

Il codice segue con la definizione operativa delle funzioni usate nel mainche non si trovano nelle librerie standard.

077. double H0(double x)078. {079. return sin(x);080. }

Tale funzione definisce la condizione iniziale. Alla funzione H0 si passaun numero reale salvato come double ed essa restituisce un altro double, chenel nostro caso e dato dal corrispondente valore della funzione seno.

081. double HSx(double t)082. {083. double temp;084. temp = 0;085. return temp;086. }

087. double HDx(double t)088. {089. double temp;090. temp = 0;091. return temp;092. }

HSx e HDx definiscono le condizioni sui bordi sinistro e destro della gri-glia. Nel nostro caso stiamo supponendo che la funzione H si annulli in tuttii punti con x=x0 e x=x1 per qualsiasi valore di t.

093. double ** allocate(int N_, int K_)094. {095. double **A;

44

096. A = (double **)calloc(N_,sizeof(double *));097. for(int i=0;i<N_;++i) *(A+i)=(double *)calloc(K_,sizeof(double));098. return A;099. }

Con questa funzione si alloca memoria dove salvare i dati relativi ai puntidella griglia. Si genera cosı un vettore di dimensione N che fa riferimentoai diversi istanti di tempo e i cui elementi rimandano ad un altro vettore,di dimensione K , che interessa la mesh puntata da n. Si crea dunque unagriglia N *K .

Riportiamo di seguito i grafici ottenuti fittando i dati relativi alla soluzionenumerica dell’equazione del calore per due diversi valori della costante didiffusione α. Si e usato il programma Mathematica 5 per fare i grafici in 3De GNUplot 4.2 per quelli in 2D per mettere in risalto le differenze tra i dueandamenti. Per entrambe abbiamo considerato le condizioni iniziali ed i variparametri come riportato nel codice precedente.

Figura 2.14: Grafico 3D per α=0.5.

45

Figura 2.15: Sovrapposizione del grafico di mesh temporali differenti perα=0.5.

46

Figura 2.16: Grafico 3D per α=1.

Figura 2.17: Sovrapposizione del grafico di mesh temporali differenti perα=1.

47

Dai grafici 3D, che riportano l’andamento delle soluzioni sull’intero inter-vallo di tempo considerato, si nota come queste abbiano una forma sinusoi-dale nei primi istanti di tempo tuttavia essa tende man mano ad appiattirsi.Nei grafici ottenuti con GNUplot si osserva come la condizione iniziale siaidentica per entrambe i casi ma ben presto la soluzione relativa ad α = 1 equasi del tutto appiattita sull’asse delle ascisse. Non e cosı per il grafico incui α = 0.5 in cui si considera cioe un fattore di diffusione minore.

Dopo aver introdotto i vari tipi di equazioni alle derivate parziali e le loropeculiarita siamo passati in questo secondo capitolo a presentare il processodi discretizzazione, descrivendo vari tipi di stencil e i problemi relativi allastabilita delle relative soluzioni.Il processo di discretizzazione risulta essere molto utile nello studio di equa-zioni differenziali per cui in molti testi, come [9], [11] e [8], si affronta ampia-mente tale tematica ricavando diversi tipi di griglia.Notiamo infine che nel capitolo, si e fatto spesso ricorso all’equazione delcalore (2.28) per tradurre in pratica i vari temi affrontati. Inoltre e statoscritto un codice, si veda 2.6, per risolvere l’equazione suddetta ricorrendoad un approccio esplicito alla discretizzazione.

48

Capitolo 3

Equazioni autoreferenziali

Prima di procedere oltre, va precisato che una delle caratteristiche princi-pali delle soluzioni marcianti di una certa equazione al continuo consiste nelfatto che la stabilita di un certo schema evolutivo e legata alla localita del-l’informazione nei tempi immediatamente precedenti a quello in esame. Lastruttura della stencil - la sua regione di influenza - deve essere connessa alpunto di aggiornamento. Questo, chiaramente, avviene per equazioni locali,cioe equazioni che sul lato destro fanno comparire un numero finito di deri-vate. E’ naturale chiedersi come si puo estendere questa procedura al caso incui la stencil prende informazione da punti molto distanti nel reticolo. Se ilnumero di derivate cresce la stencil tende ad espandersi, raggiungendo puntimolto distanti nel passato per un dato punto di aggiornamento. E’ evidenteche il limite non locale di una certa implementazione numerica viene ad es-sere raggiunto per un numero infinito di derivate. In taluni casi l’equazionepuo essere descritta da un nucleo evolutivo non locale, come e tipico in certeequazioni integro-differenziali, in altri casi la serie delle derivate rimane indi-cata in quanto tale, cioe non puo essere risommata. Nella prossima sezionecercheremo di effettuare una distinzione tra localita e non localita nell’am-bito di certe equazioni di diffusione che appaiono in fisica statistica e chesono essenziali per una corretta descrizione di processi di natura stocastica.Chiariremo che la presenza di memoria in un processo e una caratteristicadel processo stesso che non solo porta a stencil non locali, ma addirittura astencil che non sono connesse. Questo ultimo punto verra discusso nell’am-bito della equazione autoreferenziale di Pascali, che descrive un processo conmemoria.

49

0 1x y

+−

Figura 3.1: Rappresentazione grafica della equazione master nel caso uni-dimensionale. Abbiamo scelto l’intervallo 0 < x < 1 nella variabileindipendente.

3.1 Equazioni evolutive in fisica fondamenta-

le ed in statistica

Non c’e forse modo piu semplice per introdurre il concetto di assenza dimemoria in un processo diffusivo (caratteristica detta Markoviana [19], daMarkov, che per primo analizzo questo tipo di processi) che mediante uncosiddetto random walk o cammino aleatorio. In un cammino aleatorioil camminante (walker) si muove da un sito all’altro con una probabilita localedi transizione. In altri termini il salto aleatorio puo essere assunto esserelocalizzato. Questo significa che la probabilita di transizione decade in modosufficientemente rapido, ad esempio in modo esponenziale, a partire dal puntodi partenza. Chiaramente questa non e l’unica possibilita. Esistono tipi dicammini aleatori in cui si ammettono anche salti non locali ed ad intervallidi tempo non prefissati, tra i quali i continuous time random walksche portano alle distribuzioni di Levy. E’ noto, invece, che cammini aleatoricon salti localizzati ed ad intervalli di tempo fissato portano a distribuzionidi probabilita di tipo gaussiano. Queste distribuzioni, ad esempio, nel casounidimensionale funzioni della posizione x e del tempo t, q(x, t), soddisfanoad equazioni di diffusione del tipo dell’equazione del calore nella forma

∂p(x, t)

∂t= D

∂2p(x, t)

∂x2. (3.1)

Ad esempio, si puo assumere, nel nostro caso, che il camminatore, ad ogniistante di tempo, abbia una probabilita p di spostarsi a destra dal puntopresente e 1− p di andare a sinistra. In tal caso avremo

q ≡ 1− p D = 2pq(∆x)2

∆t. (3.2)

Nel limite continuo la densita di probabilita che si trovi un camminatore in(x, x + ∆x) e data da

p(x, t) =

(1

2π2Dt

)1/2

exp

[−1

2

(x− vt)2

2Dt

]. (3.3)

50

Va osservato che nella caratterizzazione del cammino aleatorio abbiamo bi-sogno sia di una distribuzione iniziale, che ci dica quale e la densita di pro-babilita di trovare il camminatore in (x, x + ∆x), che di una probabilita ditransizione che ci dica quale e la probabilita che il camminatore si muova daun certo sito ad un altro. Il lettore interessato puo trovare una discussionedettagliata della sua derivazione in [20]. Qui ci limitiamo a dire che nel casotipico di moto browniano la probabilita di transizione medesima da (x0, t0)ad (x, t), cioe w(x, t|x0, t0), soddisfa una equazione del tipo (3.1), ossia

∂w(x, t|x0, t0)

∂t= D

∂2w(x, t|x0, t0)

∂x2(3.4)

con la condizione di localita

w(x, t, |x0, t0) → δ(x− x0), t → t0. (3.5)

Non e difficile mostrare che e possibile riscrivere l’equazione di diffusionenella forma di rate equation o master equation

∂

∂τf(x, τ) =

∫ +∞

−∞[w(x|x′)f(x′, τ)− w(x′|x)f(x, τ)] dx′ (3.6)

∂

∂τq(x, τ) =

∫ +∞

−∞[w(x|x′)q(x′, τ)− w(x′|x)q(x, τ)] dx′. (3.7)

Inoltre sotto condizioni opportune di decrescita del nucleo w, questa equa-zione puo essere riscritta come equazione differenziale di ordine arbitrario

∂

∂τq(x, τ) =

∫ +∞

−∞dy

∞∑n=1

(−y)n

n!

∂n

∂xn(w(x + y|x)q(x, τ)) . (3.8)

Il risultato puo essere riscritto nella forma

∂

∂τq(x, τ) =

∞∑n=1

(−y)n

n!

∂n

∂xn(an(x)q(x, τ)) (3.9)

usando la cosiddetta espansione di Kramers-Moyal dove

an(x) =

∫ +∞

−∞dy(y − x)nw(y|x). (3.10)

Per noi e importante osservare che, in generale, un’equazione con memoria,non ha una master equation o rate equation ad essa associata e quindi, ingenerale, non e equivalente ad una equazione di diffusione con un numero

51

infinito di derivate sul lato destro dell’equazione di diffusione. Per questaragione, lo studio della stencil per queste equazioni e molto complesso ed erichiesta una approfondita indagine teorica/numerica per arrivare alla descri-zione di queste equazioni. Chiaramente, sarebbe interessante vedere se, comenel caso del moto browniano, esiste una descrizione microscopica di questeequazioni. Nel caso dell’equazione che analizzeremo in seguito questo fattonon e noto.

3.2 Un’equazione con memoria

Passiamo adesso allo studio di un’equazione differenziale che ha espressamen-te un termine di memoria (implicito) nella sua definizione. Questa equazionee stata introdotta da Eduardo Pascali in alcuni lavori recenti [14, 15, 16].In essi l’autore discute e prova l’unicita della soluzione di questa equazioneper condizioni iniziali Lipschiziane. Il nostro studio rappresenta un primo ap-proccio nella direzione di determinare numericamente una possibile soluzionedi questa equazione. L’equazione in questione e definita da

u(x, t) = u0(x) +

∫ t

0

u

(∫ τ

0

u(x, s)ds, τ

)dτ, t ≥ 0, x ∈ R (3.11)

che corrisponde al sistema{

∂∂t

u(x, t) = u(∫ t

0u(x, s)ds, t

), t ≥ 0, x ∈ R

u(x, 0) = u0(x), x ∈ R(3.12)

dove naturalmente u0 esprime la condizione iniziale.Cerchiamo di dare una giustificazione statistica di questa equazione.L’equazione descrive l’evoluzione di una variabile deterministica, la u(x, t), ilcui aggiornamento temporale dipende da tutta la sua storia, cioe dai valoridella funzione negli istanti precedenti, misurati a partire da un tempo inizialet0, qui scelto essere nullo. Tale dipendenza e di tipo integrale. Una realiz-zazione pratica di una evoluzione di questo tipo puo avvenire nell’ambitodei prodotti finanziari derivati (derivative products) e quindi puo giocare unruolo nell’econofisica. In genere, in un accordo finanziario, il valore moneta-rio di un contratto, la u(x, t) ad un certo tempo t dipende da una variabile,ad esempio un’azione, in questo caso rappresentata dalla x che puo variarealeatoriamente. Nel nostro caso assumeremo che l’azione stessa abbia unaimplicita dipendenza dal rendimento nei tempi precedenti al tempo in esame,cioe si puo formalmente ridefinire

〈x〉 ≡∫ t

0

d t′u(x, t′) (3.13)

52

ed assumere che l’aggiornamento del valore del prodotto derivato al tempo ttenga conto del suo rendimento precedente. Una variante possibile di questaequazione consisterebbe nel sostituire alla (3.13) la sua media temporale

〈x〉 ≡ 1

t

∫ t

0

d t′u(x, t′). (3.14)

In tal caso, l’aggiornamento del contratto al tempo t sarebbe dipendente dal-la media temporale del valore monetario dell’azione sottostante negli istantiprecedenti. Va osservato come contratti di questo tipo sono molto piu dif-ficili da prezzare rispetto a contratti basati su prodotti derivati ordinari,quali le calls1 e le puts2 (rimandiamo il lettore a [20] per una introduzio-ne alla econofisica dei mercati finanziari). E’ chiaro che questa equazione,nonostante la sua complessita, descrive ancora un processo deterministico.Una sua estensione stocastica richiederebbe un’analisi specifica che va al dila degli obiettivi del nostro studio. In tal caso, comunque, la variabile x,l’azione o share, dovrebbe essere interpretata come una variabile aleatoria erichiederebbe l’introduzione del calcolo stocastico [20].

3.3 Interpretazione e discretizzazione

Cerchiamo di capire meglio, con riferimento alla figura 3.2, il funzionamentodella (3.12) per passare poi a discretizzarla. Innanzi tutto questa indica comel’evoluzione temporale della funzione u(x, t) dipenda dal valore che la stessaassume, nello stesso istante di tempo, in un punto la cui ordinata e datadall’integrale all’interno delle parentesi tonde.Discretizziamo la derivata fatta rispetto al tempo con una differenza finitaal primo ordine avanti, cioe

∂

∂tu(x, t) =

un+1i − un

i

∆t(3.15)

in cui i ed n sono gli indici che individuano i punti della griglia, rispettiva-mente sull’asse x e t, e ∆t il passo lungo quest’ultimo asse.Denotiamo poi il membro di destra nella (3.12) come un

j . Notiamo che l’a-pice n fa riferimento alla mesh temporale mentre il pedice j e relativo allacoordinata spaziale, il cui valore e dato dall’integrale della funzione u(x, t)

1Nel gergo finanziario con l’appellativo di call si indica l’opzione di voto o opzione diacquisto (diritto di prelazione degli azionisti ad acquistare nuove azioni entro un tempodeterminato).

2Una put o put option indica la scelta di vendita cedibile ad un prezzo piu alto di quellodeciso in precedenza allo scopo di evitare la caduta nel tasso di cambio.

53

Figura 3.2: Schema relativo all’equazione (3.12).

54

fatto sull’ordinata per τ compreso tra 0 e t.Discretizzando la (3.12) otteniamo dunque la seguente equazione alle diffe-renze

un+1i − un

i

∆t= un

j (3.16)

da cui e facile ricavareun+1

i = uni + ∆tun

j . (3.17)

Dall’ultima equazione notiamo come il valore che la variabile u assume nelpunto di coordinate (i, n+1) e costruito a partire dal valore della stessa fun-zione nel punto temporalmente precedente, di coordinate (i, n), e in quellocaratterizzato da (j, n). Inoltre si puo osservare come l’incognita un+1

i siapresente solo nel lato sinistro dell’equazione il che ci indica che l’approccioda noi seguito e di tipo esplicito per cui dovremmo fare delle considerazionisulla stabilita della soluzione.Per risolvere il sistema riportato nella (3.12) abbiamo fatto ricorso a dei meto-di numerici di integrazione e differenziazione. Abbiamo inoltre supposto chela condizione iniziale fosse data dalla semplice equazione u0(x) = sinxcosxin quanto tale funzione risulta essere Lipschitziana3. Tale condizione e allabase del teorema 4 che assicura l’esistenza di una ed una sola soluzione perl’equazione considerata.Particolare attenzione deve essere posta sul fatto che non si e riusciti a deter-minare la stabilita dell’algoritmo alla base del codice riportato in appendicein quanto essendo l’equazione sotto esame non di tipo lineare ad essa nonpossono essere applicati i metodi di indagine visti precedentemente. Nonpossiamo dunque provare la stabilita della soluzione numerica.Dal momento che non possiamo sapere se l’eventuale carattere divergentedella soluzione e ad essa intrinseco o dovuto agli errori numerici abbiamopensato di ovviare a tale impedimento considerando un intervallo temporalemolto ridotto. Si e dunque studiato l’andamento della soluzione fornitacidall’algoritmo ipotizzando che i relativi errori fossero cosı piccoli da poteressere trascurati. In altre parole abbiamo assunto l’andamento in questo pic-colo intervallo come proprio della funzione stessa e indipendente dagli erroridovuti all’algoritmo e agli arrotondamenti attuati dal calcolatore.Si nota allora come in questo caso la soluzione sembra avere un caratteredivergente gia per piccoli intervalli di tempo, il che ci lascia pensare che lasoluzione dell’equazione (3.12) segua un andamento di tipo esponenziale etenda subito a ∞.

3Dato un insieme F ⊆ Rn, una funzione f : F → R si dice Lipschitziana su F , seesiste una costante L > 0 (chiamata costante di Lipschitz) tale che per ogni x, y ∈ F si ha|f(y)− f(x)| ≤ L‖y − x‖ [17].

4Si consulti [14, pags 257-262].

55

3.4 Codice sorgente e grafici

001. #include <stdio.h>002. #include <math.h>003. #include <stdlib.h>

Vengono passate al pre-compilatore le librerie (di input/output, matema-tica e standard) che contengono i prototipi delle funzioni usate nel corpo delprogramma.

004. #define P 8005. #define K 500006. #define N 500007. #define NGP 50

Si definiscono le variabili a livello globale. La prima indica la dimensionedella stringa mesh[ ]. Le due seguenti costituiscono il numero di punti daprendere per il campionamento della funzione nell’intervallo rispettivamentesull’asse x e y, mentre con NGP definiamo il numero di punti gaussiani.

008. double U0(double x); double ** allocate(int N_, int K_);009. doubleIntegra1(double *T, double *XG, double *WG, double **AA, int n,

int k, double h);010. void gauleg(double x1, double x2, double x[ ],double w[ ], int n);011. void map(double *Y, double *WY, double *YY,double *WWY, int NG,

double RMIN, double RMAX, double RMED);012. double phit(double *T, int k, double **AA, double yy, double h);013. double phix(double *X, int n, double **AA, double xp, double h);

Prima del main si riportano infine i prototipi delle funzioni che vengonorichiamate nel resto del codice. Ricordiamo che la prima parola di ogni rigadi testo fa riferimento al tipo di variabile restituita dalla funzione, la secondaal nome stesso della funzione, mentre tra le parentesi tonde sono riportati iparametri da cui essa dipende.

014. int main(void)015. {016. double **U, *X, *T;017. double deltax, deltat;018. double x0, x1, t0, t1;019. double xp;020. double PI = 3.141592653858979;021. int k,n,i;022. int M[P] = {0, 2, 4, 9, 19, 29, 49, 54};

56

023. char mesh[50];

024. X = (double*)calloc(K, sizeof(double));025. T = (double*)calloc(N, sizeof(double));

Si apre il main definendo la tipologia (double, int o char) delle diversevariabili. Si passa quindi ad allocare memoria in maniera continua, tramitela funzione calloc [1], creando una matrice K*N in cui salvare i dati di tipodouble.

026. x0 = -3*PI/2;027. x1 = 3*PI/2;028. t0 = 0;029. t1 = 10;

030. deltax=(x1 -x0)/(K-1);031. deltat=(t1 -t0)/(N-1);

032. for(k=0; k<K; ++k) X[k]= x0 + k*deltax;033. for(n=0; n<N; ++n) T[n]= t0 + n*deltat;

Vengono definiti a questo punto gli estremi degli intervalli su ambo gliassi. Si calcolano dunque i rispettivi passi e si riempiono i vettori X[k] e T[n]con i valori dei punti della griglia cosı ottenuti.

034. U = allocate(N, K);

035. for(k=0;k<K;++k) U[0][k] = U0(X[k]);

Con queste due linee di codice si crea lo spazio necessario per salvare ivalori dalla variabile U e si riempie la prima mesh con i dati relativi allacondizione iniziale.

036. double XG[NGP], WG[NGP];

037. bool exit = false;038. double LL = x1-x0;039. double TT = t1-t0;

Si definiscono i vettori XG e WG con dimensione NGP e le variabili LLe TT rappresentano delle costanti di normalizzazione.La variabile exit e di tipo booleano e viene settata su false.

57

040. for(n=0; n<N-1; ++n)041. {042. for(k=0;k<K;++k)043. {044. xp = Integra1(T, XG, WG, U, n, k, deltat);045. if (xp>x1 || xp<x0)046. {047. printf("DEBUG:\n Siamo fuori dall’intervallo:

x’=%f\n",xp);048. exit = true;049. n--;050. break;051. }

052. U[n+1][k] = U[n][k] + deltat*phix(X, n, U, xp, deltax);053. }054. if (exit) break;055. }056. printf("\n");

I due cicli for innestati permettono di muoverci su tutta la griglia primalungo una mesh e poi passando alla successiva assegnando ad ogni punto ilvalore della funzione U. Quest’ultimo e calcolato a partire dal valore ad essosottostante e facendo ricorso alla funzione phix.Per calcolare il valore di xp si usa invece la funzione Integra1 a cui vengonopassati i valori della variabili T, XG, WG, U, n, k e deltat. Per un’analisidettagliata della funzione si rimanda al codice seguente il main.Tramite il comando if si controlla che xp ricada all’interno del dominio. Incaso negativo viene stampato in video un messaggio d’errore, si ritorna allamesh precedente e si esce dal ciclo.

057. FILE *pfile;058. pfile = fopen("dati.txt","w");059. for(i=0;i<n+1;++i)060. {061. for(k=0;k<K;++k)062. {063. fprintf(pfile,"%f %f %f\n",X[k],T[i],U[i][k]);064. }065. }066. fclose(pfile);

Nell’ultima parte del main si fa ricorso ai puntatori a file in modo taleda creare un file txt di nome ′′dati′′ nel quale vengono salvati i valori dellagriglia, ricorrendo ai vettori X[ ] e T[ ], e quelli della funzione U.

58

067. for(i=0;i<P;++i)068. {069. sprintf(mesh, "dati_%d_mesh.txt", M[i]+1);070. pfile = fopen(mesh,"w");071. printf("Mesh %d al tempo t = %lf\n", M[i]+1, T[M[i]]);072. for(k=0;k<K;++k)073. {074. fprintf(pfile,"%f %f\n",X[k], U[M[i]][k]);075. }076. fclose(pfile);077. }078. printf("\n");079. system("PAUSE");080. return 0;081. }

Usiamo dei cicli for innestati e l’istruzione sprintf per creare tanti filequanti sono gli elementi del vettore M[ ]. Ciascun file ha nel nome il numerodella mesh considerata e all’interno i valori di U(x, t) al variare di x. Fac-ciamo inoltre stampare in video il valore della variabile t a cui ogni meshconsiderata si riferisce.

Riportiamo ora le funzioni usate. Solitamente queste vengono scritte su-bito dopo il main, oppure inserite in un’apposita libreria che deve esserepero richiamata all’inizio del file. Abbiamo optato per la prima possibilitache rende l’intero codice piu semplice da leggere e modificare.

082. double U0(double x)083. {084. return sin(x)*cos(x);085. }

Questa funzione definisce le condizioni iniziali della funzione U sulla primamesh, caratterizzata dall’indice i pari a 0. In questo caso stiamo supponendoche U0 sia del tipo sin(x)cos(x).

086. double phit(double *T,int k, double **AA, double yy, double h)087. {088. double zp,zz;089. int np=0;090. while((np*h)<yy)091. ++np;092. if (T[np]==yy)093. {094. return 0;095. }

59

096. zz = AA[np][k]*( 1. - (yy -T[np])/h) + (yy - T[np])/h*AA[np+1][k];097. return zz;098. }

La funzione precedente ha lo scopo di interpolare i valori lungo l’asse delleordinate restituendo il valore della funzione cui punta **A, cioe U, nel puntoindicato con XXG[j] in Integra1.

099. double phix(double *X,int n, double **AA, double xp, double h)100. {101. double zp,zz;102. int np=0;103. while((np*h)<xp)104. ++np;105. zz= AA[n][np]*( 1. - (xp -X[np])/h) + (xp - X[np])/h*AA[n][np+1];106. return zz;107. }

Anche questa funzione ha il compito di interpolare i valori della funzionecui punta **A. A differenza della precedente viene restituito il valore dellafunzione in questione, quindi U, nel punto d’ascissa xp che non coincide apriori con uno dei punti del reticolo.

108. double Integra1(double *T, double *XG, double *WG, double **AA,int n, int k, double h)

109. {110. int j;111. double zed,rmed;112. double XXG[NGP], WWG[NGP];

113. rmed = T[n]/2;114. map(XG,WG,XXG,WWG,NGP,0.,T[n],rmed);

115. gauleg(0., T[n], XXG, WWG, NGP);

116. zed=0.;117. for(j=0;j<NGP;++j)118. {119. double prova = phit(T,k,AA,XXG[j], h);120. zed += prova *WWG[j];121. }122. return zed;123. }

La funzione Integra1 calcola l’integrale che costituisce poi il valore di xp.Per far cio si ricorre al metodo di Gauss-Legendre 5 calcolando il valore della

5Si consulti B per maggiori dettagli sul processo di integrazione gaussiana.

60

U solo su un insieme discreto di punti, nel nostro caso NGP=50. Si noti ilricorso alle funzioni gauleg e map, esposte di seguito, per generare e map-pare sull’intervallo [0,t] i punti gaussiani.

124. double ** allocate(int N_, int K_)125. {126. double **A;127. A = (double **)calloc(N_,sizeof(double *));

128. for(int i=0;i<N_;++i) *(A+i)=(double *)calloc(K_,sizeof(double));129. return A;130. }

Con questa funzione si alloca memoria in cui salvare i dati relativi ai puntidella griglia. Si genera per prima cosa un vettore di dimensione N che fariferimento ai diversi istanti di tempo e i cui elementi rimandano ad un altrovettore, di dimensione K , che interessa la parte spaziale. Si crea cosı unagriglia N *K .

131. void gauleg(double x1,double x2,double x[], double w[], int n)132. {133. double PI=3.141592653858979;134. double EPS = 3.0e-11;

135. int m,j,i;136. double z1,z,xm,xl,pp,p3,p2,p1;137. m = (n+1)/2;138. xm=0.5 *(x2 + x1);139. xl= 0.5 *(x2 - x1);

140. for(i=1; i<m+1; ++i)141. {142. z = cos(PI*(i-0.25)/(n+0.5));143. do144. {145. p1=1.0;146. p2=0.0;147. for(j=1;j<=n;j++)148. {149. p3=p2;150. p2=p1;151. p1=((2.0*j - 1.0)*z*p2 -(j - 1.0)*p3)/j;152. }

153. pp=n*(z*p1 - p2)/(z*z - 1.0);154. z1=z;

61

155. z=z1 - p1/pp;156. }

157. while( fabs(z - z1)>EPS);158. x[i-1]=xm - xl*z;159. x[n-i]=xm + xl*z;160. w[i-1]=2.0*xl/((1.0 - z*z)*pp*pp);161. w[n-i]=w[i-1];162. }163. }

Come detto in precedenza tale funzione serve a generare i punti ed i pesigaussiani. I valori in ingresso sono costituiti dagli estremi dell’intervallo inquestione, da due vettori in cui salvare la posizione dei singoli punti gaussianied il relativo peso, e infine da un intero che determina il numero complessivodei punti gaussiani da determinare. Si noti inoltre come il parametro d’usci-ta di gauleg [4] sia impostato su void in quanto la funzione non restituiscealcun valore bensı calcola i punti gaussiani ed i rispettivi pesi salvandoli neivettori x[ ] e w[ ].

164. void map(double *Y, double *WY, double *YY, double *WWY, int NG,double RMIN, double RMAX, double RMED)

165. {166. double PI=3.141592653858979;167. double DELTA = RMAX-RMIN;168. double SN,CS,RDIST=(RMED-RMIN)*(RMAX-RMIN)/(RMAX-RMED);169. int i;

170. for(i=0;i<NG;++i)171. {172. SN=sin(PI/4*(1+Y[i]));173. CS=cos(PI/4*(1+Y[i]));174. WWY[i]=WY[i]*RDIST*PI/4/(CS+RDIST/DELTA*SN)/

(CS+RDIST/DELTA*SN);175. YY[i]=RMIN+RDIST*SN/(CS+RDIST/DELTA*SN);176. }177. }

La funzione map ridistribuisce opportunamente i punti gaussiani, trovatiper un intervallo unitario, sull’intero dominio d’integrazione.

62

Figura 3.3: Grafico 3D della soluzione numerica per U0(x) = sin(x)cos(x).

Figura 3.4: Grafico 3D della soluzione numerica per U0(x) = sin(x)cos(x).

63

Figura 3.5: Sezione bidimensionale della soluzione a diversi istanti di tempo.

Da questi grafici si nota come la funzione sia caratterizzata da un an-damento cosinusoidale lungo l’asse x, dovuto alla condizione iniziale, anchesulle mesh successive alla prima. Inoltre la funzione U(x, t) sembra decaderein maniera abbastanza rapida, quasi in modo esponenziale sebbene questoandamento potrebbe essere osservato meglio per un intervallo temporale mag-giore. Tuttavia la funzione U va a pescare ben presto dei punti che si trovanoal di fuori dell’intervallo considerato. Ipotizziamo allora che essa sia caratte-rizzata da un’instabilita che risulta essere intrinseca al problema (3.12) cheabbiamo tentato di risolvere numericamente.Supponiamo di fissare i parametri t0, t1, x0, x1 e K=300, come riportatonel codice di cui al paragrafo 3.4. Vogliamo studiare il comportamento dellasoluzione al variare del parametro α dato da

α =deltat

deltax=

(t1− t0

x1− x0

)(K − 1

N − 1

). (3.18)

Si nota allora che α e funzione della sola variabile N che rappresenta il nume-ro di mesh considerate sull’intervallo, costante, t1-t0. Nella seguente tabella

64

riportiamo, al variare di N, il valore di α(N) e di ts, tempo relativo alla meshin cui la soluzione considera punti esterni al dominio.

N α(N) ts400 0.795110 1.203008500 0.635769 1.182365600 0.529631 1.185309700 0.453861 1.187411800 0.397057 1.1889861000 0.317566 1.1911915000 0.063462 1.186237

Grafichiamo ora il comportamento delle varie soluzioni ottenute, per diversivalori di α, sull’ultima mesh.

Figura 3.6: Grafico relativo all’andamento della soluzione numericasull’ultima mesh al variare di N.

65

Naturalmente si ottengono dei risultati del tutto confrontabili con i pre-cedenti se si considera α come α(K). In questo caso si tiene dunque fisso ilnumero N e si fa variare K. Ovviamente cio che importa e il comportamentodella funzione al variare del rapporto (deltat/deltax), qualsiasi sia la varia-bile N o K che andiamo a cambiare. Ecco spiegato il motivo del perche irisultati, che omettiamo per ovvi motivi, risultano essere equivalenti.Proviamo poi a cambiare la condizione iniziale. Supponiamo, ad esempio,che la U0 sia una funzione lineare con coefficiente angolare pari ad un ven-tesimo dell’unita. Modificando opportunamente il codice e ponendo N=500,K=300, x1=200 otteniamo l’andamento mostrato nella figura 3.7.

Figura 3.7: Grafico 3D della soluzione numerica per U0(x) = 0.05 · x.

66

Figura 3.8: Andamento della soluzione per U(x) = 0.05 · x a diversi istantidi tempo.

Abbiamo notato che la soluzione, in questo caso, usciva dall’intervallo ats = 8.016032. Naturalmente aumentando il coefficiente angolare della con-dizione iniziale si otteneva una soluzione per un intervallo temporale minore,mentre facendo tendere a zero la costante che moltiplica x si notava un al-lungamento dei tempi in cui la funzione si aggiorna considerando punti neldominio. Inoltre per la condizione limite in cui U0 sia identicamente nulla siottiene una soluzione, sebbene banale, definita su ogni intervallo temporaleconsiderato.Consideriamo infine un’altra condizione iniziale, U0(x) = 1/(x+3),che risul-ta essere ancora Lipschitziana.Fissando i parametri N=600, K=300, e facendo variare x e t rispettivamentein [0,50] e [0,5] otteniamo

67

Figura 3.9: Grafico 3D della soluzione numerica per U0(x) = 1/(x + 3).

Figura 3.10: Andamento della soluzione per U0(x) = 1/(x + 3) a diversiistanti di tempo.

68

Si noti come in questo caso la funzione viene rappresentata, nel grafico3D, su tutto il dominio temporale in quanto per la particolare configurazioneconsiderata la funzione non esce mai dall’intervallo lungo le ascisse. Sebbenesiamo costretti ad osservare la funzione solo su un intervallo molto ridottodella variabile t, possiamo comunque notare come la funzione assuma uncomportamento quasi esponenziale. Purtroppo non e possibile osservare ilcomportamento della soluzione per t maggiori, in quanto aumentando l’in-tervallo temporale la funzione va a considerare punti esterni al dominio dellex.

In questo capitolo abbiamo introdotto il concetto di processi markoviani enon markoviani, nonche l’equazione di Pascali. Abbiamo tentato di ottener-ne una soluzione numerica attraverso un processo di discretizzazione tuttaviadata la sua complessita non si puo condurre uno studio sulla stabilita. Inoltrenon si puo neanche tentare, come indicato in [19], di rendere non markovia-no il processo cercando di ottenere poi una soluzione approssimata. Tuttocio che abbiamo fatto non e che un possibile approccio al problema, il qualesicuramente deve essere studiato piu a fondo non solo dal punto di vista ma-tematico, ma anche dal punto di vista fisico per capire se equazioni come la(3.12) possono essere usate per lo studio di fenomeni fisici o finanziari. Anchel’algoritmo per lo studio della (3.12) puo essere migliorato per fare in modoche la funzione resti nel dominio il piu a lungo possibile. Questo probabil-mente puo essere fatto normalizzando l’integrale di u(x, s), moltiplicandoloper (1/t), oppure inserendo opportuni fattori di scala.

69

Capitolo 4

Verso la CFD

Continuiamo il nostro studio considerando alcune equazioni differenziali pro-prie della Fluidodinamica computazionale (CFD). Ripercorriamo quin-di la strada che porta alla loro derivazione prima di passare a consideraredei metodi di discretizzazione che permettono di risolvere alcune di esse inparticolari casi.

4.1 Cenni preliminari

Un primo approccio e quello di considerare il fluido in ciascuna sua parte,cioe facendo riferimento alle caratteristiche di tutte le particelle che lo com-pongono. Dato pero il numero molto elevato di queste possiamo far ricorsoai metodi statistici per determinare i vari parametri in termini di grandezzemedie. Alternativamente si puo spostare l’attenzione passando a considera-re il fluido non nel suo insieme bensı restringendo il campo d’indagine aduna regione finita o addirittura infinitesima. Nel primo caso si ottengonodirettamente delle equazioni di tipo integrale, nel secondo invece delleequazioni differenziali. Sia le prime che le seconde possono poi essereconservative o non conservative a seconda del fatto che l’elemento con-siderato si pensi fisso all’interno del fluido oppure che si possa muovere conle linee di flusso.Iniziamo con l’introduzione di alcuni concetti chiave che verranno spesso usa-ti nel proseguo.Sia

df

dt=

∂f

∂t+

∂f

∂x

∂x

∂t+

∂f

∂y

∂y

∂t+

∂f

∂z

∂z

∂t(4.1)

la derivata totale di una generica funzione f(t, x, y, z) fatta rispetto al tempo.Notiamo che ∂x/∂t, ∂y/∂t, ∂z/∂t possono essere interpretate come le com-

70

ponenti, lungo i tre assi cartesiani, di un vettore V che esprime la velocita.Ricordando poi l’espressione dell’operatore divergenza, data da

∇ ≡ i∂

∂x+ j

∂

∂y+ k

∂

∂z, (4.2)

possiamo scrivere la (4.1) nella forma seguente

df

dt≡ ∂f

∂t+ V · ∇f. (4.3)

Definiamo allora l’operatore derivata sostanziale come

D

Dt≡ ∂

∂t+ (V · ∇). (4.4)

Possiamo interpretare la derivata sostanziale come la velocita di variazionedella funzione f(t, x, y, z), mentre ∂/∂t altro non e che la derivata localefatta rispetto al tempo e V · ∇ la derivata convettiva. Si puo pensareche la prima sia dovuta ad una variazione istantanea, mentre la seconda alcambiamento della posizione.Possiamo considerare, in ultima analisi, la derivata sostanziale come la deri-vata totale fatta rispetto al tempo di una funzione che dipende direttamenteanche da quest’ultimo parametro e interpretare fisicamente l’espressione di∇ ·V come la velocita della variazione del volume di un elemento del fluido,di volume unitario, che si muove con esso.

4.2 Le equazioni della CFD

4.2.1 L’equazione di continuita

L’equazione che si ricava dal principio fisico della conservazione della massae detta equazione di continuita. Questa puo assumere diverse forme aseconda del particolare sistema fisico al quale si applica il principio fisicoconsiderato. Inoltre, le equazioni ottenute sono tra loro equivalenti e si puopassare dall’una all’altra in maniera abbastanza semplice.Procediamo allora a considerare il sistema fisico costituito da un elementoinfinitamente piccolo del fluido che si muove con esso. Capiamo bene cheil nostro sistema puo intendersi caratterizzato da una massa infinitesimacostante

δm = ρδV (4.5)

in cui denotiamo con ρ la densita, che date le dimensioni supponiamo co-stante, e δV l’elemento infinitesimo di volume.

71

Ricordando che stiamo assumendo la conservazione della massa, deduciamoche

D(δm)

Dt= 0, (4.6)

cioe che la velocita con cui varia la massa dell’elemento del fluido consideratoe nulla.Sostituendo la (4.5) nella (4.6) e dividendo tutto per δV otteniamo

Dρ

Dt+ ρ

[1

δV

D(δV )

Dt

]= 0 (4.7)

e, ricordando che∇·V esprime proprio la velocita della variazione del volumedell’elementino per unita di volume, possiamo scrivere

Dρ

Dt+ ρ∇ ·V = 0. (4.8)

Siamo cosı pervenuti ad una forma differenziale alle derivate parziali dell’e-quazione di continuita.Facciamo subito due osservazioni. L’ipotesi di aver a che fare con un sistemainfinitamente piccolo si riflette nella forma differenziale della (4.8), mentre ilfatto che l’elemento del fluido e supposto mobile si ripercuote sul caratterenon conservativo dell’equazione di cui sopra.Considerando differenti tipi di sistemi di partenza, finiti o infinitesimi, fissio in movimento, si ottengono quattro differenti equazioni che si presentanorispettivamente in forma integrale o differenziale, conservativa o non conser-vativa. Possiamo pero passare dall’una all’altra attraverso semplici passaggiche sfruttano la definizione di derivata sostanziale, di divergenza di un pro-dotto di uno scalare per un vettore ed il teorema della divergenza.Si ottengono cosı le restanti equazioni

∂

∂t

∫∫∫

V

ρdV +

∫∫

S

ρV · dS = 0 (4.9)

partendo da un elemento finito di fluido fisso nello spazio;

D

Dt

∫∫∫

V

ρdV = 0 (4.10)

nel caso in cui l’elemento finito di fluido si muova con esso;

∂ρ

∂t+∇ · (ρV) = 0 (4.11)

considerando un elemento infinitesimo fisso nello spazio con il fluido che gliscorre attraverso.

72

4.2.2 L’equazione del momento

Ci occupiamo ora dell’equazione del momento, cioe dell’equazione che siottiene considerando la seconda legge di Newton (F=ma). Come nel casoprecedente consideriamo il sistema formato da un elemento infinitesimo difluido che si muove unitamente ad esso.Notiamo a questo punto che l’equazione che esprime la seconda legge diNewton, riportata tra parentesi, e di tipo vettoriale. Consideriamone laproiezione lungo l’asse x

Fx = max (4.12)

in cui Fx e dovuta sia a forze a distanza che a forze che agiscono direttamentesulla superficie del sistema considerato.Se indichiamo con f la forza a distanza per unita di massa che agisce sulsistema, allora

ρfx(dxdydz) (4.13)

e la forza che agisce sull’intero volumetto (dxdydz) nella direzione x.Le forze che agiscono sulla superficie sono invece dovute sia alla pressione delfluido circostante che a deformazioni e stiramenti della superficie stessa.

Figura 4.1: Schematizzazione delle forze agenti, lungo l’asse x, sull’elementoinfinitesimo in movimento.

Con riferimento alla figura 4.1 segue che la forza di superficie che agiscesull’elemento di fluido considerato nella direzione positiva dell’asse x e datada [

p−(

p +∂p

∂xdx

)]dydz +

[(τxx +

∂τxx

∂xdx

)− τxx

]dydz

73

+

[(τyx +

∂τyx

∂ydy

)− τyx

]dxdz +

[(τzx +

∂τzx

∂zdz

)− τzx

]dxdy. (4.14)

Quindi, la forza totale agente nella direzione delle x positive e data dallasomma della (4.13) e della (4.14), cioe da

Fx =

[−∂p

∂x+

∂τxx

∂x+

∂τyx

∂y+

∂τzx

∂z

]dxdydz + ρfxdxdydz. (4.15)

La precedente (4.15) altro non e che il membro sinistro della (4.12).Ricordiamo poi che la massa e data dal prodotto tra la densita (ρ) ed ilvolume (dxdydz) e che l’accelerazione lungo l’asse x e la derivata sostanzialedi u. Alla luce di cio e della precedente (4.15), riscriviamo la (4.12) dividendotutto per dxdydz, ottenendo

ρDu

Dt= −∂ρ

∂x+

∂τxx

∂x+

∂τyx

∂y+

∂τzx

∂z+ ρfx. (4.16)

Analoghe espressioni si ottengono per le componenti lungo y e z

ρDv

Dt= −∂ρ

∂y+

∂τxy

∂x+

∂τyy

∂y+

∂τzy

∂z+ ρfy, (4.17)

ρDw

Dt= −∂ρ

∂z+

∂τxz

∂x+

∂τyz

∂y+

∂τzz

∂z+ ρfz. (4.18)

Le equazioni (4.16),(4.17) e (4.18) sono dette equazioni di Navier-Stokesnella forma non conservativa. Notiamo che quest’ultima caratteristicaderiva dal fatto che l’elementino di partenza era supposto in movimento conil fluido; tuttavia si puo passare alla forma conservativa osservando che

∂(ρu)

∂t= ρ

∂u

∂t+ u

∂ρ

∂t⇒ ρ

∂u

∂t=

∂(ρu)

∂t− u

∂ρ

∂t(4.19)

e

∇· (ρuV) = u∇· (ρV)+(ρV) ·∇u ⇒ (ρV) ·∇u = ∇· (ρuV)−u∇· (ρV).(4.20)

Possiamo allora riscrivere il membro di sinistra della (4.16) esplicitando laderivata sostanziale

ρDu

Dt= ρ

∂u

∂t+ ρV · ∇u (4.21)

e, sfruttando la (4.19) e la (4.20), si ottiene

ρDu

Dt=

∂(ρu)

∂t− u

∂ρ

∂t− u∇ · (ρV) +∇ · (ρuV)

=∂(ρu)

∂t− u

[∂ρ

∂t+∇(ρV)

]+∇ · (ρuV). (4.22)

74

A questo punto si osserva che la somma presente all’interno delle parentesigraffe e il membro sinistro dell’equazione di continuita (4.11) per cui la (4.22)diventa

ρDu

Dt=

∂(ρu)

∂t+∇ · (ρuV). (4.23)

Sostituendo infine la (4.23) nell’equazione (4.16) si ha

∂(ρu)

∂t+∇ · (ρuV) = −∂ρ

∂x+

∂τxx

∂x+

∂τyx

∂y+

∂τzx

∂z+ ρfx. (4.24)

Analogamente per gli assi y e z si ottiene

∂(ρv)

∂t+∇ · (ρvV) = −∂ρ

∂y+

∂τxy

∂x+

∂τyy

∂y+

∂τzy

∂z+ ρfy, (4.25)

∂(ρw)

∂t+∇ · (ρwV) = −∂ρ

∂z+

∂τxz

∂x+

∂τyz

∂y+

∂τzz

∂z+ ρfz. (4.26)

Le equazioni (4.24), (4.25) e (4.26) costituiscono le equazioni di Navier-Stokes nella forma conservativa.

4.2.3 L’equazione dell’energia

Cerchiamo ora di ricavare l’equazione dell’energia che si ottiene impo-nendo la legge di conservazione dell’energia.Come nei casi precedenti consideriamo un elemento infinitesimo di fluido chesi muove con esso. Ci aspettiamo di ottenere dunque l’equazione dell’energianella sua forma differenziale e non conservativa.Dalla legge fisica suddetta si ricava che la velocita di variazione dell’energiaall’interno del sistema (che indichiamo con A) e data dalla somma del flussodi calore (B) che interessa il medesimo elemento del fluido e della variazionedel lavoro (C) compiuto dalle forze a distanza e superficiali, che agisconosullo stesso. Quindi, possiamo scrivere

A = B + C. (4.27)

Il membro di sinistra e dato dal prodotto della derivata sostanziale dell’ener-gia totale del sistema, calcolata per unita di massa, per la massa complessivadell’elemento considerato, quindi

A = ρD

Dt

(e +

V 2

2

)dxdydz (4.28)

75

in cui con e e V 2/2 indichiamo rispettivamente l’energia interna e cineticadel sistema per unita di massa. 1

Passiamo ora a determinare un’espressione per il termine indicato nella (4.27)con la lettera B.Il flusso del calore e dovuto all’eccitazione che interessa l’intero volume (comel’emissione o l’assorbimento di radiazione) ed alla conduzione termica cheavviene attraverso la superficie.Se indichiamo con q la variazione del calore volumetrico per unita di massa,allora il primo dei due contributi a B e dato semplicemente da

ρq(dxdydz) (4.29)

mentre a causa della conduzione termica, considerando le due facce dell’ele-mentino della figura 4.2 perpendicolari all’asse x, si ha

[qx −

(qx +

∂qx

∂xdx

)]dydz = −∂qx

∂xdxdydz. (4.30)

Ricordando poi che q e proporzionale al gradiente della temperatura possiamoriscrivere la (4.30) e sommarla alla (4.29) ottenendo

B =

[ρq +

∂

∂x

(k∂T

∂x

)+

∂

∂y

(k∂T

∂y

)+

∂

∂z

(k∂T

∂z

)]dxdydz. (4.31)

Determiniamo infine il valore di C, cioe della variazione del lavoro esercitatosull’elemento considerato. Come visto in precedenza bisogna tener conto siadelle forze che agiscono a distanza sull’intero sistema che di quelle che siesercitano solo sulla superficie. Le prime contribuiscono al valore di C con laquantita

ρf ·V(dxdydz) (4.32)

in cui indichiamo con ρ la densita dell’elemento, V la sua velocita e f la forzaa distanza agente su di esso.Le forze superficiali sono dovute alla pressione ed alle sollecitazioni esterne.Con riferimento alla figura 4.2 e con la convenzione che le forze agenti lungola direzione delle x crescenti producono un lavoro positivo, e viceversa perquelle che vanno nel verso contrario, possiamo scrivere

[up−

(up +

∂(up)

∂xdx

)]dydz = −∂(up)

∂xdxdydz (4.33)

1Per una discussioni piu articolata dell’energia di sistemi molecolari e atomici si consulti[11].

76

Figura 4.2: Schematizzazione delle forze agenti, lungo l’asse x, sull’elementoinfinitesimo in moto col fluido.

che rappresenta la variazione del lavoro dovuta alla pressione lungo la di-rezione x, mentre aggiungendovi le tensioni che agiscono sulle varie facce,sempre lungo l’asse x si ottiene

[−∂(up)

∂x+

∂(uτxx)

∂x+

∂(uτyx)

∂y+

∂(uτzx)

∂z

]dxdydz. (4.34)

Il termine indicato con C e allora dato dalla somma delle forze agenti sul-l’intero volume con quelle superficiali, agenti lungo le direzioni x, y e z. Perqueste ultime si ottengono delle espressioni analoghe alla (4.34).Possiamo dunque scrivere

C =

[−

(∂(up)

∂x+

∂(vp)

∂y+

∂(wp)

∂z

)+

∂(uτxx)

∂x+

∂(uτyx)

∂y

+∂(uτzx)

∂z+

∂(vτxy)

∂x+

∂(vτyy)

∂y+

∂(vτzy)

∂z+

∂(wτxz)

∂x

+∂(wτyz)

∂y+

∂(wτzz)

∂z

]dxdydz + ρf ·Vdxdydz. (4.35)

Sostituendo le equazioni (4.28), (4.31) e (4.35) nella (4.27) e dividendo amboi membri per il volume dxdydz si ottiene

ρD

Dt

(e +

V 2

2

)= ρq +

∂

∂x

(k∂T

∂x

)+

∂

∂y

(k∂T

∂y

)+

∂

∂z

(k∂T

∂z

)

− ∂(up)

∂x− ∂(vp)

∂y− ∂(wp)

∂z+

∂(uτxx)

∂x

77

+∂(uτyx)

∂y+

∂(uτzx)

∂z+

∂(vτxy)

∂x+

∂(vτyy)

∂y+

∂(vτzy)

∂z

+∂(wτxz)

∂x+

∂(wτyz)

∂y+

∂(wτzz)

∂z+ ρf ·V. (4.36)

Quest’ultima rappresenta l’equazione dell’energia nella forma non conserva-tiva che stavamo cercando.Naturalmente questa non e che una delle possibili forme che si possono ot-tenere. In particolare essa e relativa all’energia totale data dalla sommadi quella interna e di quella cinetica. Altre forme possono essere ottenuteper la sola energia interna o passando dalla forma non conservativa a quel-la conservativa, con un procedimento analogo a quello visto nel paragrafoprecedente.

4.2.4 Equazioni di Navier-Stokes

Vogliamo riassumere ora le leggi che governano un fluido viscoso per ilquale si presentano dei fenomeni dissipativi, quali l’attrito o la conduzionetermica, che ne aumentano l’entropia. Trascurando i fenomeni legati alladiffusione della massa, otteniamo

Equazioni di continuita

Forma non conservativaDρ

Dt+ ρ∇ ·V = 0 (4.37)

Forma conservativa∂ρ

∂t+∇ · (ρV) = 0 (4.38)

Equazioni del momento

Forma non conservativa

ρDu

Dt= −∂p

∂x+

∂τxx

∂x+

∂τyx

∂y+

∂τzx

∂z+ ρfx (4.39)

ρDv

Dt= −∂p

∂y+

∂τxy

∂x+

∂τyy

∂y+

∂τzy

∂z+ ρfy (4.40)

ρDw

Dt= −∂p

∂z+

∂τxz

∂x+

∂τyz

∂y+

∂τzz

∂z+ ρfz (4.41)

78

Forma conservativa

∂(ρu)

∂t+∇ · (ρuV) = −∂p

∂x+

∂τxx

∂x+

∂τyx

∂y+

∂τzx

∂z+ ρfx (4.42)

∂(ρv)

∂t+∇ · (ρvV) = −∂p

∂y+

∂τxy

∂x+

∂τyy

∂y+

∂τzy

∂z+ ρfy (4.43)

∂(ρw)

∂t+∇ · (ρwV) = −∂p

∂z+

∂τxz

∂x+

∂τyz

∂y+

∂τzz

∂z+ ρfz (4.44)

Equazioni dell’energia

Forma non conservativa

ρD

Dt

(e +

V 2

2

)= ρq +

∂

∂x

(k∂T

∂x

)+

∂

∂y

(k∂T

∂y

)+

∂

∂z

(k∂T

∂z

)

− ∂(up)

∂x− ∂(vp)

∂y− ∂(wp)

∂z+

∂(uτxx)

∂x+

∂(uτyx)

∂y

+∂(uτzx)

∂z+

∂(vτzy)

∂z+

∂(vτxy)

∂x+

∂(vτyy)

∂y

+∂(wτxz)

∂x+

∂(wτyz)

∂y+

∂(wτzz)

∂z+ ρf ·V (4.45)

Forma conservativa

∂

∂t

[ρ

(e +

V 2

2

)]+ ∇ ·

[ρ

(e +

V 2

2

)V

]= ρq +

∂

∂x

(k∂T

∂x

)+

∂

∂y

(k∂T

∂y

)

+∂

∂z

(k∂T

∂z

)− ∂(up)

∂x− ∂(vp)

∂y− ∂(wp)

∂z+

∂(uτxx)

∂x

+∂(uτyx)

∂y+

∂(uτzx)

∂z+

∂(vτxy)

∂x+

∂(vτyy)

∂y

+∂(vτzy)

∂z+

∂(wτxz)

∂x+

∂(wτyz)

∂y+

∂(wτzz)

∂z+ ρf ·V (4.46)

4.2.5 Equazioni di Eulero

Se nello studio di un fluido trascuriamo i termini dissipativi otteniamo leseguenti equazioni, proprie quindi dei fluidi non viscosi.

79

Equazioni di continuita

Forma non conservativaDρ

Dt+ ρ∇ ·V = 0 (4.47)

Forma conservativa∂ρ

∂t+∇ · (ρV) = 0 (4.48)

Equazioni del momento

Forma non conservativa

ρDu

Dt= −∂p

∂x+ ρfx (4.49)

ρDv

Dt= −∂p

∂y+ ρfy (4.50)

ρDw

Dt= −∂p

∂z+ ρfz (4.51)

Forma conservativa

∂(ρu)

∂t+∇ · (ρuV) = −∂p

∂x+ ρfx (4.52)

∂(ρv)

∂t+∇ · (ρvV) = −∂p

∂y+ ρfy (4.53)

∂(ρw)

∂t+∇ · (ρwV) = −∂p

∂z+ ρfz (4.54)

Equazioni dell’energia

Forma non conservativa

ρD

Dt

(e +

V 2

2

)= ρq − ∂(up)

∂x− ∂(vp)

∂y− ∂(wp)

∂z+∇f ·V (4.55)

Forma conservativa

∂

∂t

[ρ

(e +

V 2

2

)]+∇·

[ρ

(e +

V 2

2

)V

]= ρq−∂(up)

∂x−∂(vp)

∂y−∂(wp)

∂z+∇f·V

(4.56)

80

4.3 Il metodo di Lax - Wendroff

Alla luce di quanto esposto nei primi capitoli, e dopo aver introdotto le equa-zioni proprie della fluidodinamica passiamo a dare un esempio di discretiz-zazione delle suddette equazioni. Vi sono diversi metodi per fare cio, ognunodei quali e legato al particolare fluido con cui si ha a che fare, ad esempioviscoso, stazionario, ecc. Riportiamo allora un semplice esempio con cui siimplementa la discretizzazione delle equazioni della fluidodinamica, che pren-de il nome di metodo di Lax - Wendroff.Il metodo considerato e applicabile a soluzioni marcianti sia nello spazio chenel tempo ed e particolarmente usato per la risoluzione di equazioni iperbo-liche.Consideriamo dunque le equazioni di Eulero in forma non conservativa perun fluido non stazionario, bidimensionale e non viscoso

∂ρ

∂t= −

(ρ∂u

∂x+ u

∂ρ

∂x+ ρ

∂v

∂y+ v

∂ρ

∂y

)Eq. di continuita (4.57)

∂u

∂t= −

(u∂u

∂x+ v

∂u

∂y+

1

ρ

∂p

∂x

)Eq. del momento lungo x(4.58)

∂v

∂t= −

(u

∂v

∂x+ v

∂v

∂y+

1

ρ

∂p

∂y

)Eq. del momento lungo y (4.59)

∂e

∂t= −

(u

∂e

∂x+ v

∂e

∂y+

p

ρ

∂u

∂x+

p

ρ

∂v

∂y

)Eq. dell’energia (4.60)

Tali equazioni risultano essere per l’appunto iperboliche rispetto al tempo.Inoltre stiamo assumendo che siano nulle sia le forze a distanza che il cosid-detto calore volumetrico, cioe f = 0 e q = 0.Cosı come esposto nel capitolo 2 consideriamo una griglia nel piano xye con riferimento alla figura 4.3 sviluppiamo in serie di Taylor una dellevariabili dipendenti, ad esempio la densita ρ

ρt+∆ti,j = ρt

i,j +

(∂ρ

∂t

)t

i,j

∆t +

(∂2ρ

∂t2

)t

i,j

(∆t)2

2+ . . . . (4.61)

Analoghe equazioni possono essere ottenute per le altre variabili, quindi

ut+∆ti,j = ut

i,j +

(∂u

∂t

)t

i,j

∆t +

(∂2u

∂t2

)t

i,j

(∆t)2

2+ . . . , (4.62)

vt+∆ti,j = vt

i,j +

(∂v

∂t

)t

i,j

∆t +

(∂2v

∂t2

)t

i,j

(∆t)2

2+ . . . , (4.63)

81

Figura 4.3: Esempio di griglia bidimensionale nel piano xy

et+∆ti,j = et

i,j +

(∂e

∂t

)t

i,j

∆t +

(∂2e

∂t2

)t

i,j

(∆t)2

2+ . . . . (4.64)

Da tali equazioni capiamo di essere in presenza di un metodo esplicito checi consente di ottenere il valore delle variabili considerate all’istante t + ∆ta partire dai valori delle stesse sulla mesh precedente. Il passo seguente eallora quello di discretizzare le derivate presenti nelle equazioni dalla (4.61)alla (4.64) in modo tale da poter studiare le funzioni su un reticolo.Occupiamoci ad esempio della (4.61). Notiamo subito che il valore delladerivata prima di ρ fatta rispetto al tempo e data proprio dall’equazione dicontinuita (4.57). Approssimando le derivate con delle differenze centrali alsecondo ordine, riscriviamo la (4.57) come

(∂ρ

∂t

)t

i,j

= −(

ρti,j

uti+1,j − ut

i−1,j

2∆x+ ut

i,j

ρti+1,j − ρt

i−1,j

2∆x

+ρti,j

vti,j+1 − vt

i,j−1

2∆y+ vt

i,j

ρti,j+1 − ρt

i,j−1

2∆y

). (4.65)

Ovviamente tutte le quantita al secondo membro della (4.65) sono note poichesi trovano sulla mesh corrispondente al tempo t.Per ottenere ora ∂2ρ/∂t2 deriviamo la (4.57) rispetto al tempo ottenendo

∂2ρ

∂t2= −

(ρ

∂2u

∂x∂t+

∂ρ

∂t

∂u

∂x+ u

∂2ρ

∂x∂t+

∂u

∂t

∂ρ

∂x

82

+ρ∂2v

∂y∂t+

∂ρ

∂t

∂v

∂y+ v

∂2ρ

∂y∂t+

∂v

∂t

∂ρ

∂y

). (4.66)

Si rimpiazzano poi le derivate prime con le relative differenze centrali alsecondo ordine, mentre le derivate miste presenti nella (4.66) sono ottenutedifferenziando in maniera opportuna le equazioni da (4.57) a (4.60). Adesempio, per ottenere ∂2u/∂x∂t deriviamo la (4.58) rispetto a x

∂2u

∂x∂t= −u

∂2u

∂x2−

(∂u

∂x

)2

− v∂2u

∂x∂y− ∂u

∂y

∂v

∂x− 1

ρ

∂2p

∂x2+

1

ρ2

∂p

∂x

∂ρ

∂x(4.67)

e discretizzando con differenze centrali al secondo ordine si ottiene

(∂2u

∂x∂t

)t

i,j

= −uti,j

uti+1,j − 2ut

i,j + uti−1,j

(∆x)2−

(ut

i+1,j − uti−1,j

2∆x

)2

− vti,j

uti+1,j+1 + ut

i−1,j−1 − uti−1,j+1 − ut

i+1,j−1

4(∆x)(∆y)

− uti,j+1 − ut

i,j−1

2∆y

vti+1,j − vt

i−1,j

2∆x− 1

ρti,j

pti+1,j − 2pt

i,j + pti−1,j

(∆x)2

+1

(ρti,j)

2

pti+1,j − pt

i−1,j

2∆x

ρti+1,j − ρt

i−1,j

2∆x. (4.68)

In questo modo si discretizza l’equazione (4.61) ed e possibile ricavare il va-lore della densita nel punto di coordinate (i, j) al tempo t+∆t. Per ottenerecio si fa ricorso all’informazione relativa ai punti della griglia segnati con uncerchietto nella figura 4.3. L’accuratezza di tale metodo e del secondo ordine,sia nello spazio che nel tempo. Inoltre esso puo essere velocizzato attraversoaltre tecniche. Vediamone una di queste, detta tecnica di MacCormack.Secondo questo secondo metodo il valore della variabile ρ nel punto di coor-dinate (i, j) al tempo t + ∆t e dato da

ρt+∆ti,j = ρt

i,j +

(∂ρ

∂t

)

MED

∆t (4.69)

in cui con (∂ρ/∂t)MED si indica una stima del valore di tale variabile tra gliistanti di tempo t e t+∆t. Inoltre confrontando la (4.69) con la (4.61) si notache in quest’ultimo caso stiamo trascurando il termine relativo alla derivataseconda.Logicamente espressioni analoghe alla (4.69) si ottengono per le variabile u,v ed e.Procediamo dunque alla determinazione di (∂ρ/∂t)MED. Cio viene fatto in

83

due passi.Consideriamo nuovamente la (4.57) e sostituiamo le derivate prime con delledifferenze al primo ordine di tipo forward, ottenendo

(∂ρ

∂t

)t

i,j

= −(

ρti,j

uti+1,j − ut

i,j

∆x+ ut

i,j

ρti+1,j − ρt

i,j

∆x

+ ρti,j

vti,j+1 − vt

i,j

∆y+ vt

i,j

ρti,j+1 − ρt

i,j

∆y

). (4.70)

Possiamo allora ottenere una stima della densita (ρ)t+∆ti,j considerando lo

sviluppo in serie di Taylor di questa funzione e fermandoci al primo ordine

(ρ)t+∆ti,j = ρt

i,j +

(∂ρ

∂t

)t

i,j

∆t. (4.71)

Quest’equazione fornisce un possibile valore della variabile ρ; dunque, e daconsiderarsi soltanto come una stima accurata al primo ordine. Come sem-pre le espressioni relative alle variabili u, v ed e sono del tutto analoghe alla(4.71).Le espressioni cosı ottenute vengono poi impiegate nel passo successivo perottenere una stima del valore della derivata di ρ calcolata nel punto di coor-dinate (i, j) al tempo t + ∆t. In particolare si ricorre ancora una voltaalla (4.57), la si discretizza facendo uso delle differenze questa volta di tiporearward, e si ottiene

(∂ρ

∂t

)t+∆t

i,j

= −[(ρ)t+∆t

i,j

(u)t+∆ti,j − (u)t+∆t

i−1,j

∆x+ (u)t+∆t

i,j

(ρ)t+∆ti,j − (ρ)t+∆t

i−1,j

∆x

+ (ρ)t+∆ti,j

(v)t+∆ti,j − (v)t+∆t

i,j−1

∆y+ (v)t+∆t

i,j

(ρ)t+∆ti,j − (ρ)t+∆t

i,j−1

∆y

].(4.72)

Infine il valore della derivata presente nella (4.69) e dato dalla media aritme-tica delle espressioni che si ricavano dalla (4.70) e dalla (4.72), quindi

(∂ρ

∂t

)

MED

=1

2

[(∂ρ

∂t

)t

i,j

+

(∂ρ

∂t

)t+∆t

i,j

]. (4.73)

La (4.73) sostituita nella (4.69) ci consente di ottenere un valore correttodella densita al tempo t + ∆t con molti meno passaggi rispetto al metodo diLax-Wendroff nonostante abbia lo stesso grado di accuratezza.

84

In quest’ultimo capitolo abbiamo passato in rassegna le principali equazioniche governano la fluidodinamica. Come abbiamo avuto modo di dire esseassumono forme diverse a seconda del sistema sotto esame. Possono dunqueessere discretizzate e usate per la risoluzione di vari problemi fisici. Per ilmomento il lavoro termina qui, con il proposito di continuare lo studio nelladirezione segnata da questa tesi, e cioe riguardante lo studio di funzioni ela loro implementazione numerica. Ai lettori piu interessati alle applicazionifisico-numeriche delle leggi esposte in questo capitolo segnaliamo comunque[7], [8], [10].

85

Capitolo 5

Conclusioni

In questo lavoro di tesi abbiamo cercato di mettere in evidenza alcuni degliaspetti piu salienti dell’approccio alla discretizzazione di equazioni differen-ziali alle derivate parziali, di vario tipo, che sono di rilevanza in fisica classica.In particolare, abbiamo analizzato il problema della stabilita delle discretiz-zazioni, introdotto le soluzioni marcianti per equazioni paraboliche e discussole varie tecniche di aggiornamento su reticolo caratterizzandone la regione diinfluenza. Siamo poi passati a discutere uno dei casi piu complessi di questaprocedura, prendendo in esame un’equazione con memoria della quale abbia-mo cercato di derivare delle soluzioni numeriche. Infine abbiamo introdottole equazioni di Navier-Stokes e discusso alcune loro discretizzazioni. Sembraquasi superfluo sottolineare che gli aspetti computazionali che abbiamo ana-lizzato si limitano solo a scalfire la superficie di un campo di studio che siarricchisce ogni anno di nuovi importanti contributi essendo l’interfaccia travarie branche scientifiche fondamentali.

86

Appendice A

Dimostrazioni relative alleequazioni differenziali allederivate parziali

A.1 Dimostrazione: la (2.28) e parabolica

Consideriamo l’equazione di propagazione del calore lungo una direzionemessa nella forma

∂T

∂t= α

∂2T

∂x2. (A.1)

Portando tutto al secondo membro otteniamo

α∂2T

∂x2− ∂T

∂t= 0. (A.2)

Notiamo una certa analogia tra la precedente e la (1.19)

ax2 + bxy + cy2 + dx + ey + f = 0 (A.3)

che esprime l’equazione generale delle sezioni coniche.Nel nostro caso possiamo porre:

a = α, b = 0, c = 0, d = 0, e = −1, f = 0. (A.4)

QuindiD = b2 − 4ac = 0 ∀a ⇔ ∀α (A.5)

dalla definizione data nel paragrafo 1.1.1 segue che la (A.1) e un’equazioneparabolica.

87

A.2 Dimostrazione: la (2.28) e iperbolica

Consideriamo l’equazione delle onde (2.73)

∂u

∂t+ c

∂u

∂x= 0.

E’ ovvio che il differenziale totale della funzione u(x, t) e dato da

du =∂u

∂tdt +

∂u

∂xdx. (A.6)

Scriviamo quindi il sistema dato dalla (2.73) e dalla (A.6). Questo puo esserescritto come

{∂u∂t

+ c∂u∂x

= 0∂u∂t

dt + ∂u∂x

dx = du(A.7)

o equivalentemente in forma matriciale[

1 cdt dx

] [∂u/∂t∂u/∂x

]=

[0du

](A.8)

In questo caso la matrice dei coefficienti e data da

[A] ≡[

1 cdt dx

](A.9)

La condizione da imporre per trovare le curve caratteristiche e

|A| = 0

che corrisponde a diredx− cdt = 0. (A.10)

Dal momento che le variazioni sono infinitesime ma non nulle possiamodividere tutto per dt ottenendo

dx

dt− c = 0. (A.11)

Confrontando infine la (B.4) con una tipica equazione di secondo grado echiamando con a, b e c rispettivamente i coefficienti del termine quadratico,lineare e noto otteniamo

D = b2 − 4ac = 1 > 0 ∀c, (A.12)

ne consegue che l’equazione parabolica e di tipo iperbolico.

88

Appendice B

Il metodo di Gauss-Legendre

La tecnica d’integrazione gaussiana permette di approssimare l’integrale diuna funzione, definito su un intervallo finito, attraverso la somma dei valoriche la stessa assume su un insieme discreto di punti.In generale questo metodo viene applicato, e fornisce degli ottimi risultati,nel caso in cui la funzione integranda non varia rapidamente, quando cioepuo essere ben approssimata da un polinomio. In caso contrario si puo sud-dividere l’intervallo di partenza in piccoli intervalli in cui la funzione variameno rapidamente, rendendo possibile l’approssimazione considerata.Per quanto riguarda la scelta dei punti sui quali calcolare la funzione pos-siamo fare l’ipotesi che essi costituiscano le radici di una particolare classedi polinomi, detti polinomi di Legendre, da cui il nome della relativa tecnicad’integrazione1. Proprio tale posizione viene sfruttata nella funzione gauleg.Inoltre, dal momento che le radici risultano essere simmetriche nell’intervalloconsiderato, il calcolatore ne puo calcolare direttamente solo una meta conun conseguente risparmio dei tempi macchina.Supponiamo di voler calcolare l’integrale

I =

∫ 1

−1

f(x)dx. (B.1)

In generale, se una formula d’integrazione basata sulle serie di Taylor usa Npunti, allora essa integra esattamente un polinomio di grado N − 1. Suppo-niamo quindi di avere a che fare con l’integrale di un polinomio di grado p,definito sull’intervallo [-1,1], e di voler usare a tale scopo proprio N punti.Oltre a questi ultimi e importante conoscere i relativi pesi, i quali rispettano

1Questa scelta non e l’unica possibile. I punti gaussiani possono essere scelti anche inaltro modo, come si puo vedere su [2] e [5].

89

la seguente condizione

∫ 1

−1

xpdx =N∑

n=1

wnxpn; p = 0, 1, . . . , N − 1. (B.2)

Abbiamo quindi a che fare con 2N parametri, dovuti alla scelta dei punti edei pesi Gaussiani e possiamo dunque considerare polinomi di grado 2N − 1.Come gia accennato supponiamo che i punti gaussiani coincidano con le radicidei polinomi di Legendre. Naturalmente Pi e un polinomio di grado i coni radici nell’intervallo [−1, 1]. Inoltre ogni polinomio di grado 2N − 1 puoessere scritto come

f(x) = Q(x)PN(x) + R(x) (B.3)

in cui Q(x) e R(x) sono polinomi al piu di grado N − 1.Possiamo allora riscrivere la (B.1) come

I =

∫ 1

−1

(QPN + R)dx =

∫ 1

−1

Rdx (B.4)

avendo fatto uso del fatto che PN e ortogonale a tutti i polinomi di gradoN − 1 o minori.Possiamo quindi ottenere una stima dell’integrale definito nella (B.1) consi-derando

I =N∑

n=1

wn[Q(xn)PN(xn) + R(xn)] =N∑

n=1

wnR(xn) (B.5)

i cui pesi gaussiani vengono dedotti dalla (B.2) nel caso in cui consideriamocome xn gli N zeri del polinomio di Legendre PN . In particolare si nota comewn dipende sia dai punti xn che dalla derivata di PN secondo la relazione chesegue

wn =2

(1− x2n)[P ′

N(xn)]2. (B.6)

Si tenga presente, infine, che il metodo di Gauss-Legendre puo essere ap-plicato a tutti gli integrali definiti su un generico dominio finito [a,b]. Perapplicare le formule di cui sopra bastera quindi fare un cambiamento linearedi variabile, dato da

t = −1 + 2(x− a)

(b− a). (B.7)

90

Bibliografia

[1] Al Kelley and Ira Pohl A Book on C: Programming in C (4th Edition).Addison-Wesley, Pearson Education, 1997.

[2] Steven E. Koonin Computational Physics. The Benjamin/CummingsPublishing Company, Inc., 1986.

[3] Hildebrand, Francis B. Advanced Calculus for Applications Prentice-Hall,1976

[4] S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery Numerical Recipesin C: The Art of Scientific Computing. Cambridge University Press, SerieVII, Vol. 25, 2005.

[5] Abramowitz M., Stegun I. A. Handbook of Mathematical Functions. DoverPublications, Applied Mathematics Series, Vol. 55, 1968.

[6] Morten Hjorth-Jensen Computational Physics. University of Oslo, Fall,2006.

[7] Dale A., John C. Tannehill, Richard H. Pletcher Computational FluidMechanics and Heat Transfer. McGraw-Hill, Inc., 1984.

[8] Fletcher,C. A. Computational Techniques for Fluid Dynamics, vol.I:Fundamental and General Techniques. Springer-Verlag, 1988.

[9] Thompson, Joe F., Z. V. A. Warsi, C. Wayne Mastin Numerical GridGeneration: Foundations and Applications. North-Holland, 1985.

[10] Moretti, G and M. Abbett A Time Dependent Computational Methodfor Blunt Body Flows.vol. 4, , 1966.

[11] Anderson, John D., Jr. Hypersonic and High Temperature GasDynamics. McGraw-Hill, Inc., 1989.

91

[12] Anderson, John D., Jr. Modern Compressible Flow: With HistoricalPerspective (2th edition). McGraw-Hill, Inc., 1990.

[13] Whitham, G. B. Linear and Nonlinear Waves. Wiley, 1974.

[14] Michele Miranda Jr. and Eduardo Pascali On a type of evolution of self-referred and hereditary phenomena. Aequationes mathematicae, Vol. 71,No. 3, pags. 253-268, 2006.

[15] Eduardo Pascali Existence of solutions to a self-referred and heredi-tary system of differential equations. Electronic Journal of DifferentialEquations, Vol. 2006, No. 7, pags. 1-7, 2006.

[16] Michele Miranda Jr. and Eduardo Pascali On a class of differentialequations with self-reference. Rendiconti di matematica.

[17] Choquet, Dewitt, Dillard Analysis, Manifolds and Physics. North-Holland, 1977.

[18] A. Quarteroni, R. Sacco, F. Saleri Matematica Numerica (2◦ edizione)Springer, 2000.

[19] Gunter Bolch, Stefan Greiner, Hermann de Meer, Kishor S. TrivediQueueing networks and Markov chains: modeling and performance eva-luation with computer science applications (2th Edition). John Wiley &Sons, Inc., 2006.

[20] W. Paul and J. Bachnagel Stochastic Processes from Physics to Finance,Springer, 1999.

92