MODELLI MATEMATICI PERL’AMBIENTE

dr. Franco Pasquarelli

2005/2006

Indice

1 Equazioni differenziali ordinarie 41.1 Definizioni sulle equazioni differenziali ordinarie. . . . . . . . . 4

1.1.1 Cos’e la soluzione di un’equazione differenziale. . . . . 41.1.2 Curve soluzione e campi di direzione. . . . . . . . . . . 5

1.2 Problema di Cauchy. . . . . . . . . . . . . . . . . . . . . . . . 61.3 Metodi numerici per il problema di Cauchy. . . . . . . . . . . 8

1.3.1 Metodo di Eulero esplicito. . . . . . . . . . . . . . . . . 81.3.2 Sviluppi in serie e metodi Runge-Kutta. . . . . . . . . 91.3.3 Metodi impliciti e metodi multipasso. . . . . . . . . . . 101.3.4 Metodi Predictor-Corrector. . . . . . . . . . . . . . . . 12

1.4 Problemi ai limiti. . . . . . . . . . . . . . . . . . . . . . . . . . 131.5 Metodi numerici per problemi ai limiti: metodo alle differenze

finite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141.6 Modelli matematici. . . . . . . . . . . . . . . . . . . . . . . . . 17

1.6.1 Oscillatore armonico. . . . . . . . . . . . . . . . . . . . 171.6.2 Circuito elettrico. . . . . . . . . . . . . . . . . . . . . . 181.6.3 Decadimento radioattivo. . . . . . . . . . . . . . . . . . 191.6.4 Equilibrio di una sbarra elastica vincolata. . . . . . . . 201.6.5 Equazione della diffusione. . . . . . . . . . . . . . . . . 21

2 Equazioni alle derivate parziali. 232.1 Propagazione delle onde. . . . . . . . . . . . . . . . . . . . . . 23

2.1.1 Approssimazione numerica. . . . . . . . . . . . . . . . 262.1.2 Caso non lineare. . . . . . . . . . . . . . . . . . . . . . 28

2.2 Equazione della diffusione. . . . . . . . . . . . . . . . . . . . . 312.2.1 Approssimazione numerica. . . . . . . . . . . . . . . . 31

2.3 Equazione di Laplace. . . . . . . . . . . . . . . . . . . . . . . 352.3.1 Approssimazione numerica nel caso bidimensionale. . . 35

2

Prefazione

Le scienze applicate, come, per esempio, alcuni settori della matematica,della fisica, della chimica e dell’ingegneria, si occupano della formulazionedi modelli matematici per sistemi cosı detti dinamici, cioe che evolvono neltempo. Questi modelli sono di particolare interesse anche per sistemi legati aproblematiche ambientali, basta pensare a problemi quali la diffusione di in-quinanti nell’aria o nell’acqua, o all’assorbimento di sostanze nocive da partedi piante o animali. Chi dovra studiare queste problematiche deve pertantoconoscere le basi teoriche delle equazioni differenziali che costituiscono unodegli strumenti fondamentali per la stesura dei modelli matematici.L’interesse applicativo non e pero rivolto solo a ricavare modelli che descri-vano bene il sistema studiato, ma anche nel saperne calcolare una soluzione.Solo in questo modo, infatti, si puo da un lato verificare la concordanza trala soluzione del modello e il comportamento reale del sistema, dall’altro pre-vedere il comportamento del sistema in particolari situazioni. La soluzionedi problemi differenziali puo , pero , essere calcolata in modo analitico solo inpochissimi casi, si deve allora ricorrere a tecniche di discretizzazione numericache permettono di calcolare una soluzione approssimata del modello.Scopo del corso e allora quello di illustrare le principali caratteristiche delleequazioni differenziali dal punto di vista matematico/teorico e di fornire alcu-ne tecniche di approssimazione numerica. L’argomento e pero molto ampioe non puo essere affrontato in modo dettagliato in un unico corso, si cer-chera pertanto di fornire la maggior parte dei risultati teorici solo in formaintuitiva e di far capire le problematiche legate all’approssimazione numeri-ca delle equazioni differenziali fornendo esempi di buon/mal comportamentodegli schemi descritti.In ultima analisi l’intento e quello di far acquisire le basi per poter utiliz-zare con cognizione di causa codici di approssimazione prefatti e per potergiudicare in modo critico i risultati ottenuti.

3

Capitolo 1

Equazioni differenziali ordinarie

1.1 Definizioni sulle equazioni differenziali or-

dinarie.

Per equazione differenziale si intende una equazione in cui viene espressoun legame tra una funzione incognita e le sue derivate. Si parla di equazionidifferenziali ordinarie o ODE se le derivate si riferiscono ad una sola incognitae di equazioni a derivate parziali o PDE in caso di piu variabili.

Per ordine di una equazione differenziale si intende l’ordine massimo del-la derivata della funzione incognita che compare nell’equazione, quindi unafunzione differenziale di ordine n sara del tipo:

G(x, y(x), y′(x), y′′(x), ..., y(n)(x)) = 0

dove G e una funzione di n+2 variabili.

Quando e possibile esplicitare la derivata n-esima in funzione delle altre inco-gnite si ottiene la cosı detta forma normale dell’equazione differenziale chesi presta meglio sia allo studio teorico che numerico.

1.1.1 Cos’e la soluzione di un’equazione differenziale.

Per soluzione di una equazione differenziale di ordine n si intende una fun-zione y(x) che sostituita nell’equazione con le sue derivate fino all’n-esima laverifica per ogni x in un certo intervallo (a, b). Questo ovviamente implica

4

che y(x) deve essere sufficientemente regolare per ammettere derivate finoalla n-esima in (a, b).

Si dice inoltre che una soluzione y(x) e massimale quando non esiste un’altrasoluzione su un intervallo (α, β) che contiene (a, b) e che coincide con y(x) in(a, b).

Notiamo ora che con la sostituzione zi(x) = y(i)(x) i = 0, 1, · · · , n − 1una equazione differenziale ordinaria di ordine n puo essere trasformata nelseguente sistema di n equazioni differenziali del primo ordine:

G(z0, z1, · · · , z′n−1) = 0z′0 = z1

z′1 = z2...z′n−2 = zn−1

(1.1)

Possiamo allora limitarci a studiare le equazioni del primo ordine, intese even-tualmente anche nel caso vettoriale. D’ora in avanti considereremo pertantouna equazione differenziale del tipo:

y′(x) = f(x, y(x)) (1.2)

1.1.2 Curve soluzione e campi di direzione.

Per curva soluzione si intende il grafico di una soluzione dell’equazione dif-ferenziale. La sua utilita e legata ad una piu facile ed immediata interpreta-zione del comportamento della soluzione stessa, ma puo essere ottenuta solose si conosce un’espressione analitica della soluzione. Quando questa non edisponibile, sono necessari dei metodi numerici che approssimano il valoredella soluzione in un certo numero di punti fissati, per poi tracciare il graficodi una curva che passa per i valori ottenuti.

Ricordiamo ora che l’equazione differenziale 1.2 permette di conoscere il va-lore della derivata prima di y nella regione R del piano x-y in cui e definitaf e quindi di stabilire la tangente ad una soluzione che passa per un certopunto di R. Puo allora essere utile considerare il cosı detto campo di direzio-ni ottenuto rappresentando dei vettori di lunghezza fissata, applicati in unareticolazione di R e aventi per direzione la pendenza f(x, y).

5

Osservando il campo di direzioni possiamo allora sapere se due diverse solu-zioni che passano per punti vicini tendono ad avvicinarsi o ad allontanarsi;infatti se ∂f

∂y> 0 le soluzioni tenderanno ad allontanarsi e si dira che le solu-

zioni sono instabili, viceversa se ∂f∂y< 0 allora le soluzioni si avvicinano e si

dicono stabili. Questa terminologia e giustificata dal fatto che se le soluzionisono stabili allora piccoli errori iniziali tendono ad annullarsi (o comunquerimangono limitati) mentre nell’altro caso si amplificano.

Figura 1.1: Campo di direzioni per f(x, y) = .8y(1− y).

1.2 Problema di Cauchy.

Spesso nei problemi che si vogliono modellizzare si deve trovare come evolveun sistema a partire da dei dati iniziali; cio si traduce imponendo che lasoluzione del problema differenziale passi per un determinato punto del pianox-y. Si ottiene allora il cosı detto problema ai valori iniziali o problema diCauchy:

y′(x) = f(x, y)y(x0) = y0

(1.3)

Il nostro problema ha quindi come dati la funzione f e il valore iniziale (x0, y0)e come risultato deve fornire una funzione y definita in un intervallo I taleche x0 ∈ I e y(x0) = y0.Dal punto di vista analitico per trovare la soluzione del problema di Cauchybisogna trovare tutte le soluzioni massimali dell’equazione differenziale (cosıdetto integrale generale) per poi trovare la soluzione che verifica la condizioneiniziale. Questa tecnica, pero, puo essere utilizzata solo in pochi casi perche

6

esistono solo alcune techiche per calcolare l’integrale generale e queste ultimenon sempre possono essere usate. Pertanto e indispensabile utilizzare deimetodi numerici che approssimano il valore della soluzione su un insiemefissato di valori dell’incognita x.Prima di parlare di questi metodi e pero necessario richiamare alcune pro-prieta che devono essere verificate dal problema di Cauchy perche la soluzionedel problema possa essere approssimata in modo corretto:

• esistenza: il problema ammette esistenza della soluzione?

• unicita : la soluzione e unica?

• dipendenza continua dai dati: come varia la soluzione al variaredei dati del problema?

I primi due punti trovano risposta nel seguente:

Teor 1 Sia f(x, y) una funzione continua in S = (x, y)|a ≤ x ≤ b, y ∈ Re lipschitziana rispetto a y, allora per ogni (x0, y0) ∈ S esiste ed e unica lasoluzione del problema di Cauchy 1.3.

La richiesta di lipschitzianita di f si traduce nell’esistenza di una costante L(detta costante di Lipshitz) tale che:

|f(x, y1)− f(x, y2)| ≤ L|y1 − y2| ∀x ∈ [a, b] e ∀y1, y2 ∈ R

Applicando il teorema del valor medio, e facile verificare che se la derivataparziale ∂f

∂yesiste ed e limitata (cioe ∂f

∂x≤ L ) allora f e lipschitziana in y.

Per quanto riguarda la dipendenza dai dati consideriamo il seguente problemaperturbato:

Y ′(x) = f(x, y) + δ(x)Y (x0) = y0 + ε0

(1.4)

Se f verifica le condizioni del teorema 1 e se δ(x) e ε0 sono limitate allora sipuo dimostrare che esistono due costanti ε e k tali che:

7

|y(x)− Y (x)| ≤ kε ∀x ∈ [a, b]

La condizione precedente esprime una dipendenza continua della soluzionedai dati che e molto importante per le approssimazioni numeriche percheinevitabilmente i dati numerici sono affetti da errori.Tale condizione, pero , non e sufficiente a garantire che gli errori numericinon tendano ad amplificarsi perche la costante k dipende dalla costante diLipschitz e dall’ampiezza (b − a) dell’intervallo su cui si vuole calcolare lasoluzione. Il valore di k puo quindi diventare grande.

1.3 Metodi numerici per il problema di Cau-

chy.

In questa sezione ci occuperemo dell’approssimazione numerica del problemadi Cauchy 1.3. Ovviamente daremo per scontato che il problema soddisfialle condizioni del teorema 1 nell’intervallo (x0, T ) dove vogliamo calcolarela soluzione.L’idea di base di tutti i metodi che vedremo e di introdurre dei punti xi+1 =xi + hi i = 1, 2, ..., n − 1 in (x0, T ) e di calcolare un valore approssimato ui

della soluzione esatta y(xi). Tali punti prendono il nome di nodi e i valori hi

sono detti passi della discretizzazione. Piu piccoli sono i passi e maggiore sarail numero di nodi in cui calcolo la soluzione. Questo significa sicuramente unmaggior costo computazionale che si spera venga compensato da una maggiorprecisione nell’approssimazione ottenuta.Il modo piu semplice per calcolare i nodi e che siano equispaziati, cioe taliche h = hi ∀i; studieremo pertanto principalmente questo caso. Ricordiamocomunque che approssimazioni con hi diversi possono essere indispensabili inalcuni problemi concreti.

1.3.1 Metodo di Eulero esplicito.

Iniziamo questo capitolo ricordando che l’equazione della retta tangente aduna curva in un punto (x0, y0) e la seguente:

y = y0 + y′(x0)(x− x0)

8

Tenendo conto che la nostra funzione y soddisfa il problema 1.3 otteniamoallora la seguente approssimazione u1

∼= y(x1):

u1 = y0 + f(x0, y0)h

Una volta ottenuto u1 posso utilizzare la stessa tecnica per calcolare il valoreapprossimato u2. In generale si avra lo schema

ui+1 = ui + f(xi, yi)h (1.5)

detto di Eulero esplicito.Sotto opportune ipotesi si puo dimostrare che vale la seguente stima dell’er-rore:

|ui − y(xi)| ≤ Ch (1.6)

dove C dipende dall’ampiezza dell’intervallo (x0, T ), dalla costante di lipschi-tzianita di f e dalla regolarita di y.Dalla formula 1.6 si ricava che l’errore tra la soluzione calcolata e quellaesatta e di ordine h, cioe se dimezzo h l’errore si dimezza. Quindi piu he piccolo e piu l’errore e piccolo. In realta cio e vero solo in parte, perchevalori troppo piccoli di h tendono ad amplificare gli errori di arrotondamentonumerico del computer e questi comportano una perdita di accuratezza dellasoluzione. Per avere una maggiore precisione conviene pertanto ricorrere ametodi che abbiano una convergenza dell’errore a zero rispetto a potenze piualte di h. Infatti, se per un certo metodo l’errore ha un’espressione del tipo|ui − y(xi)| ≤ Chp, allora, dimezzando h, l’errore si riduce di un fattore 2p equindi, in proporzione, si riduce in modo molto piu sensibile.

1.3.2 Sviluppi in serie e metodi Runge-Kutta.

Quando la funzione f(x, y) e sufficientemente regolare si possono otteneredei metodi di ordine superiore utilizzando lo sviluppo in serie di Taylor di y;ricordiamo, infatti, che una funzione y regolare in x0 puo essere approssimatain un intorno di x0 col seguente polinomio:

y(x) = y(x0) + y′(x0)(x− x0) + y′′(x0)2

(x− x0)2 + ...+ y(n)(x0)

n!(x− x0)

n + ...

+y(N)(ξ)N !

(x− x0)N

9

Possiamo allora definire una approssimazione del tipo:

un+1 = un + h(u′n + u′′n2h+ u′′′n

6h2 + ...)

e dunque l’approssimazione sara tanto migliore quanti piu termini della suc-cessione si considerano. Le derivate di y possono essere ottenute tenendoconto del problema 1.3:

y′(x) = f(x, y)y′′(x) = f ′(x, y) = fx + fyfy′′′(x) = f ′′(x, y) = fxx + 2fxyf + fyyf

2 + fxfy + f 2y f

(1.7)

Ovviamente per poter usare questi metodi f deve essere regolare e deve esserepossibile calcolarne tutte le derivate parziali.Si preferisce allora sostituire i termini in 1.7 con delle combinazioni di valoridi f calcolata in opportuni punti in modo da non perdere l’ordine in h delloschema; in questo modo si ottengono i cosı detti metodi di Runge-Kuttaespliciti.Esempi:

un+1 = un + h2(f(xn, un) + f(xn+1, un + hf(xn, un))) Metodo di Heun

un+1 = un + hf(xn + h2, un + h

2f(xn, un)) Metodo di Eulero Modificato

un+1 = un + h6(k1 + 2k2 + 2k3 + k4)

k1 = f(xn, un)k2 = f(xn + h

2, un + h

2k1)

k3 = f(xn + h2, un + h

2k2)

k4 = f(xn+1, un + hk3)

I primi due metodi sono del secondo ordine mentre il terzo e del quarto ordineed e il metodo Runge-Kutta piu usato.

1.3.3 Metodi impliciti e metodi multipasso.

I metodi che abbiamo visto fino ad ora sono tutti di tipo esplicito, ovvero adogni passo la soluzione nel punto xn+1 viene facilmente calcolata utilizzano

10

in modo opportuno i valori gia noti nei nodi precedenti. Questi schemi so-no in generale poco costosi e facili da implementare, ma presentano spessodei fenomeni di instabilita numerica soprattutto quando la soluzione esattadel problema che si sta approssimando presenta dei forti gradienti. Spessoper ottenere una soluzione non oscillatoria si dovrebbe utilizzare un passoh cosı piccolo da rendere proibitivo l’uso di metodi espliciti. In questi casie indispensabile ricorrere a dei metodi impliciti, cioe a dei metodi in cui lafunzione f deve essere valutata anche in un+1. Vediamo allora come possonoessere ricavati dei metodi di questo tipo.Dal teorema fondamentale del calcolo integrale e da 1.2 otteniamo:

yn+1 = yn +

∫ xn+1

xn

f(t, y(t))dt (1.8)

Se ora usiamo degli schemi numerici per approssimare l’integrale definito in1.8 otteniamo degli schemi per approssimare il problema di Cauchy.Se per esempio approssimiamo l’integrale con l’area del rettangolo di baseh ed altezza f(xn, un), riotteniamo il metodo di Eulero esplicito, mentre seusiamo l’area del rettangolo di base h ed altezza f(xn+1, un+1) otteniamo ilmetodo di Eulero implicito:

un+1 = un + hf(xn+1, un+1) (1.9)

Se invece l’integrale viene approssimato con l’area del trapezio di basi f(xn, un)e f(xn+1, un+1) e altezza h otteniamo il metodo dei trapezi:

un+1 = un +h

2(f(xn+1, un+1) + f(xn, un)) (1.10)

Utilizzando la stessa tecnica si possono ottenere anche metodi cosı detti apiu passi, cioe in cui il valore di un+1 dipende non solo da un, ma anchedagli altri valori un−1, un−2, ..., un−p calcolati precedentemente. Sostituendoall’integrale di f in 1.8 l’integrale del polinomio che interpola f nei nodixn−p, xn−p+1, ..., xn si ottengono i cosı detti metodi di Adams-Bashforth chesono di tipo esplicito, mentre sostituendo all’integrale di f in 1.8 l’integrale

11

del polinomio che interpola f nei nodi xn−p, xn−p+1, ..., xn, xn+1 si ottengonoi cosı detti metodi di Adams-Multon che sono impliciti.

Sottolineiamo il fatto che per poter innescare un metodo multipasso non epiu sufficiente il solo dato iniziale del problema di Cauchy, ma sono necessarip valori iniziali che andranno calcolati con un metodo, per esempio, ad unpasso. Inoltre questi dati iniziali devono essere calcolati con una precisioneche dipende dall’ordine del metodo multipasso che si vuole utilizzare per nonpregiudicarne l’approssimazione.

Riportiamo di seguito alcuni esempi di questi metodi con il loro ordine:

metodi di Adams-Bashforth:

un+1 = un + h2[3fn − fn−1] secondo ordine

un+1 = un + h12

[23fn − 16fn−1 + 5fn−2] terzo ordine

metodi di Adams-Multon:

un+1 = un + h2[fn+1 + fn] secondo ordine

un+1 = un + h12

[5fn+1 + 8fn − fn−1] terzo ordine

1.3.4 Metodi Predictor-Corrector.

Abbiamo gia detto che i metodi impliciti sono piu stabili rispetto a quelliespliciti, ma costano di piu dal punto di vista computazionale perche ad ognipasso si deve risolvere una equazione non lineare.

Cio comporta di dover utilizzare un metodo numerico opportuno per questotipo di problemi (bisezione, corde, Newton,...) oppure la forma stessa delloschema numerico suggerisce l’introduzione dei cosı detti metodi Predictor-Corrector con i quali si ottiene una successione di valori che convergono versoun+1 con un numero di iterate molto basso. Uno schema implicito, infatti, edella forma:

un+1 = φ(un+1)

che suggerisce di definire la seguente successione di numeri che, se converge,converge al valore di un+1:

12

vk+1 = φ(vk)v0 = dato assegnato

(1.11)

Scegliendo in modo opportuno il dato iniziale si puo allora pensare di poterfare poche valutazioni dello schema 1.11 (tipicamente una sola) per ottenereuna soluzione che ha lo stesso errore rispetto alla soluzione esatta yn+1 dellasoluzione esatta un+1 dello schema implicito.Se per calcolare il dato iniziale v0 di 1.11 si usa un metodo esplicito, si otten-gono i cosı detti metodi Predictor-Corrector. Una scelta tipica e di accoppiarei metodi di Adams-Bashforth e Adams-Multon con lo stesso numero di passi.Per esempio il metodo di Heun e di questo tipo, un altro esempio e il seguente

u0n+1 = un + h

2[3f(xn, un)− f(xn−1, un−1)]

un+1 = un + h12

[5f(xn+1, u0n+1) + 8f(xn, un)− f(xn−1, un−1)]

(1.12)

che complessivamente risulta essere un metodo del terzo ordine come loschema implicito che lo definisce (lo schema esplicito e invece di ordine 2).Per quanto riguarda la stabilita , in generale i metodi Predictor-Correctorhanno una regione di stabilita che si colloca tra quelle dei due metodi (espli-cito ed implicito) che lo definiscono.

1.4 Problemi ai limiti.

Questo tipo di problemi differenziali si differenzia dal problema di Cauchy peril fatto che vengono imposte delle condizioni sulla soluzione in punti distinti.Esempio tipico per un problema del secondo ordine e quello in cui si richiedeche la soluzione assuma valori assegnati in due punti:

y′′(x) = f(x, y(x), y′(x)) a < x < by(a) = α ; y(b) = β

(1.13)

Piu in generale le cosı dette condizioni al bordo o condizioni ai limiti possonoessere delle combinazioni lineari o non lineari dei valori della soluzione e dellasua derivata prima in punti assegnati e si vuole approssimare la soluzione delproblema cosı ottenuto nell’intervallo determinato dai punti dati.

13

Per problemi ai limiti non esiste un teorema che garantisca che un certoproblema ammette esistenza ed unicita della soluzione in funzione della re-golarita dei dati. In questo caso, infatti, non e sufficiente la regolarita deidati perche il problema sia ben posto, ma ogni problema deve essere studiatoed analizzato con tecniche opportune. Prima di procedere alla discretizza-zione di problemi di questo tipo e dunque necessario uno studio teorico delmodello che ne validi la correttezza e che garantisca il grado di regolaritadella soluzione del problema.

1.5 Metodi numerici per problemi ai limiti:

metodo alle differenze finite.

Esistono vari metodi numerici per problemi ai limiti, tra i principali possia-mo citare le differenze finite, volumi finiti, elementi finiti, metodi spettrali.Ognuno di questi metodi ha una sua formulazione ed uno studio teorico chene garantisce l’accuratezza e la precisione. Per semplicita descriviamo solo ilmetodo delle differenze finite, per gli altri metodi si possono consultare testispecifici di Analisi Numerica per problemi ai limiti.

Per costruire il metodo alle differenze finite si deve innanzitutto definire unadiscretizzazione x0 = a, x1, x2, ..., xn = b dell’intervallo su cui e defini-to il problema ai limiti. Come in precedenza supporremo che i nodi sianoequispaziati, quindi xi = a+ i ∗ h con h = b−a

n.

Ora non resta che approssimare le derivate che compaiono nel problema condelle opportune combinazioni lineari di valori di y che possono essere ottenutedallo sviluppo di Taylor. Per esempio, se scriviamo lo sviluppo di Taylor diy centrato in x fermato al quarto termine e lo valutiamo in x + h e x − h,otteniamo:

y(x+ h) = y(x) + hy′(x) + h2

2y′′(x) + h3

6y′′′(x) + h4

y

(IV )(ξ)

y(x− h) = y(x)− hy′(x) + h2

2y′′(x)− h3

6y′′′(x) + h4

y

(IV )(ξ)

Se ora sottraiamo le due equazioni otteniamo:

y′(x) =y(x+ h)− y(x− h)

2h+O(h2) (1.14)

14

da cui si ricava la cosı detta differenza finita centrata del primo ordine:y′(x) ∼= y(x+h)−y(x−h)

2hche approssima la derivata prima di y con un errore

dell’ordine di h2.Se invece sommiamo le due equazioni otteniamo:

y′′(x) =y(x+ h)− 2y(x) + y(x− h)

h2+O(h2) (1.15)

da cui si ricava la cosı detta differenza finita centrata del secondo ordine:y′′(x) ∼= y(x+h)−2y(x)+y(x−h)

h2 che approssima la derivata seconda con un erroredell’ordine di h2.Direttamente dallo sviluppo di Taylor si possono invece ricavare le differenzefinite in avanti e in indietro del primo ordine date rispettivamente da y′(x) ∼=y(x+h)−y(x)

he y′(x) ∼= y(x)−y(x−h)

hche approssimano la derivata prima di y con

un errore dell’ordine di h.Utilizzando queste espressioni si possono discretizzare tutti i problemi diffe-renziali del secondo ordine. Per problemi di ordine superiore si ricavano delleespressioni analoghe.Vediamo ora con un esempio come puo essere applicato il metodo delledifferenze finite. Consideriamo il seguente problema:

−y′′(x) + q(x)y(x) = f(x) a < x < by(a) = α ; y(b) = β

(1.16)

dove q(x) e f(x) sono funzioni note continue su [a, b]. Si puo dimostrareche se q(x) ≥ 0 in [a, b] il problema 1.16 ammette esistenza ed unicita dellasoluzione e quest’ultima e continua con le sue derivate di ordine 1 e 2. Inoltre,una maggiore regolarita nelle funzioni q(x) e f(x) porta ad una soluzione ycon maggiore regolarita .Se ora imponiamo che l’equazione differenziale venga verificata nei nodi in-terni della discretizzazione e se approssimiamo le derivate con le formule alledifferenze viste in precedenza, otteniamo:

−ui+1−2ui+ui−1

h2 + qiui = fi i = 1, 2, ..., n− 1.

dove, in generale, gi = g(xi) e ui e l’approssimazione di y(xi). Aggiungendoora allo schema precedente le due condizioni al bordo otteniamo il problemadiscreto:

15

u0 = α−ui−1 + (2 + h2qi)ui − ui+1 = h2fi i = 1, 2, ..., nun = β

(1.17)

Risolvendo il sistema lineare 1.17, si calcolano dunque i valori approssimatidi y(xi). In forma matriciale 1.17 puo essere riscritto come:

Au = f

dove:

A =

1 0 0 0−1 2 + q1h

2 −1 00 −1 2 + q2h

2 −1 0...

0 −1 2 + qn−1h2 −1

0 0 1

u =

u0

u1...

un−1

un

f =

αh2f1

...h2fn−1

β

Il sistema lineare puo essere risolto con uno dei metodi numerici proposti perquesto tipo di problemi.Per quanto riguarda la convergenza dello schema, si puo dimostrare il se-guente :

Teor 2 sia y(x) la soluzione del problema 1.16 con f(x) e q(x) ≥ 0 suffi-cientemente regolari da garantire che |y(IV )| ≤M in [a, b], allora:

|yi − ui| ≤Mh2

24(xi − a)(b− xi) i = 1, 2, ..., n− 1

Lo schema proposto fornisce quindi una soluzione approssimata a meno diun errore dell’ordine di h2.

16

1.6 Modelli matematici.

In questo paragrafo presentiamo alcuni esempi di modelli matematici che siesprimono mediante delle equazioni differenziali ordinarie. In particolare, iprimi tre casi riguardano modelli che si riconducono alla risoluzione di unproblema di Cauchy, mentre gli ultimi due portano ad equazioni ai limiti.

1.6.1 Oscillatore armonico.

Vogliamo modellizzare il moto di un punto materiale di massa m vincolato amuoversi lungo una retta e legato ad un filo elastico come indicato in figura1.2.

Figura 1.2: Oscillatore armonico

Se indichiamo con y(t) lo spostamento del punto in funzione del tempo, alloray′(t) e y′′(t) saranno rispettivamente la velocita e l’accelerazione del punto.Sul punto agiscono due forze, una di tipo elastico (e quindi proporzionaleallo spostamento) ed una dipendente dall’attrito (e quindi proporzionale allavelocita ). Quindi la forza totale F che viene applicata avra la forma F =−k1y

′(t)− k2y(t). Inoltre, dalla legge di Newton, sappiamo che F = my′′(t)e quindi uguagliando le due espressioni otteniamo:

my′′(t) + k1y′(t) + k2y(t) = 0 (1.18)

che esprime la legge del moto del punto. Per conoscere come si muoverail punto, e ora necessario specificare i dati iniziali che potranno essere lospostamento e la velocita iniziali.

L’equazione 1.18 puo essere risolta analiticamente calcolando le radici delpolinomio caratteristico dell’equazione differenziale, quindi le radici di x2 +k1

mx+ k2

m= 0

17

Se il determinante ∆ =k21

m2 − 4 ∗ k2

me positivo allora il polinomio ammette

due soluzioni reali e distinte λ1 e λ2 e la soluzione di 1.18 si esprime comey(t) = c1e

λ1t + c2eλ2t dove le due costanti vanno determinate in funzione dei

dati iniziali. Per esempio, se lo spostamento iniziale e y0 e la velocita inizialee nulla, allora la soluzione e y(t) = y0λ2

λ2−λ1eλ1t − y0λ1

λ2−λ1eλ2t ed essendo λ1 e λ2

negativi y(t) tende a zero all’aumentare del tempo. Quindi il punto si spostalentamente verso il punto di equilibrio.Comportamento analogo si ottiene quando ∆ = 0.

Se invece ∆ < 0, allora la soluzione generale e y(t) = e−k12m

t[c1 cos(√−∆2t) +

c2 sin(√−∆2t)] il punto quindi oscilla intorno al punto di equilibrio con oscil-

lazioni che diventano sempre piu piccole.Nel caso in cui l’attrito e nullo (o trascurabile), allora nell’equazione 1.18non c’e’ il termine di primo grado e il polinomio caratteristico ha sempreradici immaginarie. In questo caso, pero , il termine esponenziale legato allaparte reale delle soluzioni e costante ed uguale ad 1 e quindi non si ha piul’effetto di smorzamento delle oscillazioni. Si ottiene,allora, il cosı detto motoarmonico.

1.6.2 Circuito elettrico.

Riportiamo ora un altro esempio che, pur partendo da un problema fisi-co completamente diverso dal precedente, porta ad un modello matematicomolto simile a quello del moto armonico. Consideriamo il seguente circuitoelettrico:

C

R

~

L

E

Figura 1.3: Circuito elettrico

Supponiamo di voler determinare la carica Q(t) del condensatore. Dobbia-mo allora imporre la seconda legge di Kirchoff, che afferma che in un circuitochiuso la differenza di potenziale fornita da una sorgente E e uguale allasomma delle differenze di potenziale generate nel resto del circuito. Nel no-stro caso queste differenze di potenziale sono causate da una resistenza R

18

(dunque uguale a RI per la legge di Ohm dove I e l’intensita di corren-te), dall’induttanza L (dunque uguale a LdI

dt) e dal condensatore C (dunque

uguale a Q(t)C

).Tenendo conto che l’intensita di corrente e la variazione della carica nel tempoI = dQ

dtotteniamo la seguente equazione:

Ld2Qdt2

+ RdQdt

+ QC

= E(t)

che e dello stesso tipo dell’equazione del moto armonico. In questo caso,pero , l’equazione non e omogenea, ma presenta il termine forzante E(t). Perrisolvere analiticamente quest’ultima equazione, oltre all’integrale generaledell’equazione di secondo grado, occorre pertanto trovare anche una soluzioneparticolare del problema.

1.6.3 Decadimento radioattivo.

Come ultimo esempio, vediamo come puo essere modellizzato il fenomenodel decadimento radioattivo. In natura esistono degli elementi di naturainstabile, cioe che tendono a decadere in isotopi diversi rilasciando alcuneparticelle (alpha, beta, o fotoni). Il decadimento di un singolo nucleo disostanza radioattiva non puo essere previsto, ma da osservazioni sperimentalisi puo ricavare la seguente legge:in un campione con un numero elevato dinuclei radioattivi, la dimunuzione del numero dei nuclei durante un intervallodi tempo e direttamente proporzionale alla lunghezza dell’intervallo temporalee al numero di nuclei presenti all’inizio dell’intervallo. Se allora indichiamocon N(t) il numero di nuclei al tempo t, la legge puo essere scitta come:

N(t+ ∆t)−N(t) = −kN(t)∆t

dove k e una costante positiva di proporzionalita . Se ora supponiamo cheN(t) sia una funzione a valori reali (ovviamente questa e una richiesta chenon puo essere verificata in natura, ma se l’ultilizzo del modello e solo di tipoqualitativo puo essere accettata), dividiamo l’equazione per ∆t e ne consi-deriamo il limite per ∆t che tende a zero otteniamo l’equazione differenzialedel primo ordine:

19

N ′(t) + kN(t) = 0

che puo essere facilmente risolta in modo analitico. In particolare, se siconosce la concentrazione N0 al tempo T0, si ha N(t) = N0e

−k(t−t0) che puoessere utilizzata per determinare il cosı detto tempo di dimezzamento, cioe iltempo richiesto perche il numero di nuclei radioattivi si dimezzi.

1.6.4 Equilibrio di una sbarra elastica vincolata.

Studiamo ora l’equilibrio di una sbarra elastica fissata ad un estremo e sotto-posta a delle forze di trazione (esempio la forza peso). Indichiamo con u(x)lo spostamento di un punto inizialmente posto alla distanza x dal lato fissatodella sbarra.

(x+ x)

δ

u=o

δ ∆

(x)x

x+ x∆

Figura 1.4: Equilibrio di una sbarra elastica.

L’allungamento della sbarra alla distanza x e misurato dalla variazione diu(x) e quindi e calcolabile come ε = du

dx. Se l’allungamento non e nullo, nella

sbarra si genera una forza interna σ(x) detta stress che dipende da ε.Nella sua forma piu semplice, il legame tra σ e ε (detto equazione costitutivadel materiale) e della forma σ(x) = E(x)du

dx, cosı detta legge di Hooke, dove

il fattore di proporzionalita E(x) e detto modulo di Young.Se ora consideriamo un elemento infinitesimo ∆x, ed imponiamo che le forzeche vi agiscono siano in equilibrio, otteniamo l’equazione:

E(x+ ∆x)du

dx(x+ ∆x)− E(x)

du

dx(x) + f∆x = 0

dividendo ora per ∆x e passando al limite per ∆x→ 0, otteniamo:

20

− d

dx(E(x)

du

dx(x)) = f(x)

Per quanto riguarda le condizioni al bordo, evidentemente u(0) = 0 inquanto x = 0 corrisponde all’estremo bloccato, mentre in corrispondenzadell’estremo libero x = L si ha che lo stress deve essere nullo, e quindiE(L)du

dx(L) = 0.

1.6.5 Equazione della diffusione.

Il processo della diffusione e legato al trasporto di materia dovuto al movi-mento caotico delle molecole. Da un punto di vista macroscopico tale effettoha un comportamento analogo a quello del trasferimento di calore ed in effettii due fenomeni hanno delle modellizzazioni matematiche molto simili.Consideriamo innanzitutto un mezzo isotropo, cioe una sostanza in cui ladiffusione avviene nello stesso modo lungo tutte le direzioni. Si puo alloraipotizzare che la velocita di trasferimento di materia dovuta alla diffusioneattraverso un’area unitaria sia proporzionale alla variazione di concentra-zione misurata lungo la direzione normale all’area stessa. Il coefficiente diproporzionalita D viene detto coefficiente di diffusione e puo essere costanteo dipendere dalla concentrazione.Per ricavare l’equazione della diffusione in un mezzo isotropo consideriamoun parallelepipedo rettangolo con spigoli paralleli agli assi di riferimento edi lunghezza rispettivamente di AA′ = 2dx,AB = 2dy,AD = 2dz comeindicato in figura 1.5.

A

BB’

C’

D D’

C

A’

P

Figura 1.5: Schema per lo studio della diffusione.

La variazione della concentrazione di materia nel parallelepipedo sara datadal bilancio della materia entrante ed uscente dal parallelepipedo. Dobbiamo

21

allora cercare di esprimere la velocita di trasferimento di materia attraversole singole facce.Cominciamo dalle facce (ABCD) e (A’B’C’D’) che consideriamo ortogonaliall’asse x. Indicando con Fx la velocita di trasferimento per unita di area or-togonale all’asse x e passante per il punto centrale del parallelepipedo P, pos-siamo esprimere la velocita di ingresso in (ABCD) come 4dydz(Fx − ∂Fx

∂xdx).

e la velocita di uscita in (A’B’C’D’) come 4dydz(Fx + ∂Fx

∂xdx).

In modo analogo, si possono esprimere le velocita su tutte le facce del paral-lelepipedo. Se ora indichiamo con C(x, y, z, t) la concentrazione di materianel punto P = (x,y, z) all’istante t, otteniamo la seguente equazione che legala variazione di concentrazione al flusso di materia attraverso le facce:

∂C

∂t+∂Fx

∂x+∂Fy

∂y+∂Fz

∂z= 0

Tenendo ora conto della legge di proporzionalita tra la concentrazione e lavelocita di diffusione otteniamo la seguente equazione della diffusione:

∂C

∂t=

∂

∂x(D

∂C

∂x) +

∂

∂y(D

∂C

∂y) +

∂

∂z(D

∂C

∂z) (1.19)

L’equazione 1.19 nel caso D costante o solo dipendente dal tempo diventa:

∂C

∂t= D(

∂2C

∂x2+∂2C

∂y2+∂2C

∂z2) (1.20)

All’equazione 1.19 andra poi aggiunto un termine forzante f(x, y, z, t) nelcaso in cui delle forze esterne al sistema causino creazione o rimozione dimateria all’interno del sistema stesso.Solitamente questi problemi vengono poi completati assegnando la concen-trazione sul bordo di un certo dominio assegnato (per esempio se si studia ladiffusione di un gas in un tubo, sperimentalmente si puo misurare la concen-trazione sulla superficie del tubo) oppure si assegna la quantita di materia chetransita attraverso delle superfici (sempre nel caso precedente tale condizionepuo essere assegnata all’ingresso e all’uscita del tubo). La prima condizionee del tipo C(x, y, z, t) = dato la seconda e del tipo ∂C

∂n= dato dove n indica

la direzione normale alla superficie.

22

Capitolo 2

Equazioni alle derivate parziali.

In questo capitolo verranno esaminati alcuni tipi di equazioni alle derivateparziali e per ognuno verranno indicate delle opportune tecniche di approssi-mazione. Si tratta di equazioni la cui incognita e una funzione a piu variabilied essa compare nell’equazione assieme ad alcune sue derivate parziali. L’or-dine dell’equazione e dato dal grado massimo di derivazione che comparenell’equazione.

Per questo tipo di equazioni non esiste una teoria generale che garantiscal’esistenza e l’unicita della soluzione, pertanto, per ogni problema e necessariouno studio teorico specifico.

2.1 Propagazione delle onde.

Per onda si intende un segnale od una perturbazione in un mezzo che sipropaga nel tempo trasportando dell’energia. Non necessariamente, quindi,con un’onda vi e trasporto di materiale, semplicemente si ha la propagazionedi una perturbazione che, trasportando energia, puo eventualmente causaredel trasporto di materiale.

Un semplice modello matematico che rappresenta un’onda e il seguente:

u(x, t) = f(x− ct) (2.1)

dove u indica l’intensita della perturbazione, x indica la posizione, t il tempoe c > 0 e la velocita di propagazione dell’onda. In questo caso, dunque,

23

h

t=0;f(x) t>0; f(x−ct)

Figura 2.1: Propagazione di un’onda f(x) a velocita costante.

l’onda procede con velocita costante e dunque se inizialmente (t = 0) il suoprofilo e u = f(x) al tempo t il suo profilo sara traslato a destra di ct.

Calcoliamo ora le derivate parziali di u rispetto a t e x (che indicheremorispettivamente con ut e ux), otteniamo:

ut = −cf ′(x− ct) ux = f ′(x− ct)

e dunque u verifica l’equazione:

ut + cux = 0 x ∈ R; t > 0 (2.2)

che e un’equazione differenziale a derivate parziali del primo ordine.Un’equazione di questo tipo viene detta equazione del trasporto o advec-tion equation perche puo essere utilizzata per studiare la propagazione dellaconcentrazione di una sostanza in un fluido che si muove con velocita c.Se invece considerassimo l’equazione ut − cux = 0 allora la soluzione sarebbedel tipo u(x, t) = f(x+ ct), cioe si tratterebbe di un’onda che si sposta versosinistra con velocita costante c.Per studiare queste equazioni, e necessario introdurre la definizione di curvecaratteristiche.Consideriamo allora l’equazione 2.2 con la condizione iniziale:

u(x, 0) = ψ(x) (2.3)

La soluzione del problema e allora u(x, t) = ψ(x − ct). Le rette x − ct = dcon d dato costante, rappresentano le linee lungo le quali il dato inizialeviene propagato con valore costante. Si tratta pertanto delle linee del pianospazio-temporale lungo le quali viene trasportato il segnale. Lungo tali linee

24

la derivata temporale di u e nulla poiche e uguale all’equazione differenzialeche definisce il problema.

Infatti, se calcoliamo la derivata temporale di u lungo la retta C data dax(t) = ct+ d otteniamo:

du(x(t), t)

dt= ux(x(t), t)

dx

dt+ ut(x(t), t) = ux(x(t), t)c+ ut(x(t), t) = 0

La famiglia di rette x = ct+d con d costante prende il nome di famiglia dellecurve caratteristiche dell’equazione differenziale.

Vediamo ora come la definizione puo essere estesa al caso di equazioni ditrasporto in cui il termine c non e costante.

Indichiamo con C la famiglia di curve definite dall’equazione differenziale:

dx

dt= c(x, t) (2.4)

Con conti analoghi ai precedenti troviamo che lungo una curva di C laderivata di u rispetto a t e nulla e dunque u e costante.

Conoscendo le curve caratteristiche si puo facilmente trovare l’espressione del-la soluzione dell’equazione differenziale come vediamo nel prossimo esempio.Consideriamo il problema:

ut + 2tut = 0 x ∈ R, t > 0

u(x, 0) = e−x2x ∈ R

Le curve caratteristiche, definite dall’equazione dxdt

= 2t, sono le parabolex = t2 + d con d costante. Se ora vogliamo conoscere il valore della soluzioneu(x, t), sappiamo che il suo valore e uguale al valore di u(ξ, 0) dove ξ e taleche x = t2 + ξ. Ricaviamo cosı :

u(x, t) = e−ξ2

= e−(x−t2)2

si tratta pertanto di un’onda che si muove uguale a se stessa, ma con velocitache aumenta nel tempo.

25

Figura 2.2: Rappresentazione dell’onda u(x, t) = e−(x−t2)2 ai tempi t = 0,t = 1 e t = 2

2.1.1 Approssimazione numerica.

Vogliamo ora approssimare la soluzione del problema 2.2 con le condizioniiniziali 2.3, utilizzando schemi alle differenze finite. Introduciamo innanzitut-to dei nodi equispaziati di discretizzazione xj = jh |j| = 0, 1, 2, ... in spazioe tn = nk n = 0, 1, 2, ... in tempo. h e k saranno, rispettivamente, i passsi didiscretizzazione in spazio e tempo.Utilizzando la notazione uj,n

∼= u(xj, tn) e collocando l’equazione 2.2 nelpunto u(xj, tn), otteniamo le seguenti equazioni utilizzando lo schema alledifferenze finite in avanti del primo ordine sia per lo spazio che per il tempo:

uj,n+1 − uj,n

k+ c

uj+1,n − uj,n

h= 0 |j| = 0, 1, 2, ... n = 0, 1, 2, ...

Calcolando ora i valori iniziali uj,0 come approssimazione del dato ψ(xj),resta da risolvere il problema:

uj,n+1 = (1 + λ)uj,n + λuj+1,n |j| = 0, 1, 2, ... n = 0, 1, 2, ...

dove λ = ckh. Quindi per conoscere la soluzione approssimata nel nodo xj+1

al tempo (n + 1)k devo conoscere la soluzione approssimata nei nodi xj+1 exj al tempo precedente.Analogamente, usando per la derivata spaziale lo schema alle differenze finiteall’indietro del primo ordine otteniamo:

uj,n+1 = (1− λ)uj,n + λuj−1,n |j| = 0, 1, 2, ... n = 0, 1, 2, ...

26

dove ora per conoscere la soluzione approssimata nel nodo xj+1 al tempo(n + 1)k devo conoscere la soluzione approssimata nei nodi xj−1 e xj altempo precedente.Infine, se utilizziamo una differenza centrata per approssimare la derivataspaziale, otteniamo:

uj,n+1 = uj,n +λ

2uj+1,n −

λ

2uj−1,n |j| = 0, 1, 2, ... n = 0, 1, 2, ...

dove per conoscere la soluzione approssimata nel nodo xj+1 al tempo (n +1)k devo conoscere la soluzione approssimata nei nodi xj−1 e xj+1 al tempoprecedente.Schematicamente nella figura 2.3 riportiamo il dominio di dipendenza discre-ta (cioe i nodi in cui si deve conoscere la soluzione per poter calcolare il valoreuj,n+1) per i tre metodi proposti.

n+1

n

j j+1 j+2 j−1 j j+1 j−1 j j+1

Figura 2.3: Dominio di dipendenza discreta per gli schemi a differenze finite.

Se supponiamo che c > 0, allora sicuramente lo schema basato sulle differenzein avanti non puo essere utilizzato perche in questo caso l’onda si sposta versodestra mentre il dominio di dipendenza discreta si muove verso sinistra. Lasoluzione in uj,n+1 verrebbe dunque a dipendere da valori che fisicamente nonla influenzano. Lo schema basato sulla differenza all’indietro, invece, potraessere usato a patto che i suo dominio di dipendenza discreta contenga lalinea caratteristica che passa per (xj, tn+1) e quindi se k

h≤ 1

c. Dunque per

delle opportune scelte dei passi di discretizzazione spaziale e temporale.Analogamente, se c < 0 e lo schema basato sulla differenza all’indietro a nonpoter essere utilizzato.Lo schema basato sulla differenza centrata e apparentemente migliore deiprecedenti perche sicuramente per opportune scelte dei parametri di discre-tizzazione il dominio discreto contiene la linea caratteristica che passa per(xj, tn+1) indipendentemente dal segno di c. Purtroppo, pero , lo schema non

27

e convergente, nel senso che esistono dei dati iniziali ψ(x) tali che, facendotendere a zero h e k in modo che il loro rapporto sia costante e minore di 1

c,

la soluzione non converge a quella esatta.Per ottenere un metodo convergente occorre utilizzare lo schema:

uj,n+1 =1− λ

2uj+1,n +

1 + λ

2uj−1,n |j| = 0, 1, 2, ... n = 0, 1, 2, ...

che e uno schema con un errrore dell’ordine di k+h2 e che converge se λ ≤ 1.Quest’ultima condizione e fondamentale per lo studio della convergenza ditutti i metodi numerico per le equazione delle onde.Se si vogliono degli schemi con un errore di ordine 2 sia in spazio che intempo, si possono usare i seguenti:

uj,n+1 = (1− λ2)uj,n −λ(1− λ)

2uj+1,n +

λ(1 + λ)

2uj−1,n Schema LaxWenroff

uj,n+1 = uj,n−1 − λ(uj+1,n − uj−1,n) Schema leap − frog

L’ultimo schema richiede l’uso di tre passi temporali e dunque e piu costosoper quanto riguarda l’uso della memoria del calcolatore.

2.1.2 Caso non lineare.

Analizziamo ora brevemente il caso di equazioni non lineari ed in partico-lare il caso di equazioni di conservazione. Piu precisamente consideriamoun’equazione del tipo

ut + Fx(u) = 0 (2.5)

che esprime il fatto che la variazione di una quantita fisica u e uguale allavariazione del flusso −F (u) attraverso una superficie. Il significato di u,di flusso e dell’espressione di quest’ultimo in funzione di u dipendono dal

28

problema che si sta modellizzando. In particolare, l’espressione di F (u) espesso il rsultato di legami trovati sperimentalmente.Per esempio, in fluidodinamica u puo essere la densita di un fluido e il flussoe la quantita di fluido che passa in un punto x al tempo t. Sperimentalmente,allora, si possono ricavare equazioni del tipo:

F (u) = ku(1− u)F (u) = ku flusso lineareF (u) = ku2 flusso quadratico

dove k e una costante e 1 rappresenta la densita massima.Il flusso quadratico e un modello che si adatta alla rappresentazione di unflusso automobilistico, se supponiamo dunque di voler risolvere il seguenteproblema:

ut + 2uux = 0 −∞ < x < +∞, 0 < t < +∞u(x, 0) = ψ(x) −∞ < x < +∞ (2.6)

la soluzione u rappresentera come evolve la densita di automobili a partiredallo stato iniziale ψ(x).Cerchiamo ora la forma delle curve caratteristiche del nostro problema. Dob-biamo allora risolvere il problema dx

dt= 2u che puo essere interpretato come

studiare il moto di una particella che si muove con velocita 2u. Al tempo t laparticella sara dunque nella posizione x = x0 + 2u(x, t)t (in generale, se l’e-quazione non lineare da risolvere e ut + g(u)ux = 0 l’equazione caratteristicasara x = x0 + g(u)t). Ora, siccome sulle linee caratteristiche la quantita u siconserva, ovvero u(x, t) = u(x0, 0), otteniamo che

x = x0 + 2u(x0, 0)t

e l’equazione della linea caratteristica uscente da (x0, 0) (nel caso generalesara x = x0 + g(u(x0, 0))t).Supponiamo ora che la densita iniziale sia

u(x, 0) =

1 x ≤ 01− x 0 < x < 10 x ≥ 0

(2.7)

29

otteniamo allora le seguenti linee caratteristiche:

x = x0 + 2t x ≤ 0x = x0 + 2(1− x0)t 0 < x < 1x = x0 x ≥ 0

(2.8)

quindi le linee caratteristiche sono delle rette con inclinazione diversa in fun-zione del punto da cui partono. In particolare (si veda la figura 2.4) le lineecaratteristiche uscenti dai punti dell’intervallo 0 ≤ x ≤ 1 si intersecano altempo t = 1

2.

Cio e dovuto al fatto che i valori di u piu grandi viaggiano con una velocitasuperiore a quelli piu piccoli e quindi il profilo dell’onda si distorce. Se, comenel nostro caso, le linee caratteristiche si incrociano, si ha il cosı detto feno-meno dell’onda d’urto con una soluzione che diventa discontinua e presentaun cosı detto fronte d’onda.

Quando si forma la discontinuita , l’equazione differenziale non e piu veri-ficata ed il fenomeno deve essere studiato in modo diverso. In generale, sipuo osservare che il fronte d’onda si sposta con una velocita che e legataall’ampiezza del salto del fronte.

x

t=0.5

t

1

Figura 2.4: Linee caratteristiche del flusso di automobili.

t=0 t=0.5

x x x

u u u

Figura 2.5: Profilo di densita delle automobili.

30

2.2 Equazione della diffusione.

Abbiamo gia descritto come si ricava l’equazione per problemi di diffusione,vogliamo ora presentare alcuni risultati sulla risoluzione numerica di questimodelli utilizzando schemi alle differenze finite.

Consideriamo dunque il seguente problema in una dimensione spaziale:ut = Duxx a ≤ x ≤ b, 0 < t ≤ T

u(x, 0) = f(x) a ≤ x ≤ bu(a, t) = l(t) 0 < t ≤ Tu(b, t) = r(t) 0 < t ≤ T

(2.9)

La soluzione di problemi di diffusione ad un certo istante temporale, contra-riamente a quanto succedeva nel moto delle onde, e influenzata dal valoredella soluzione al tempo 0 in tutti i punti dell’intervallo (a, b). In questocaso, si dice che i segnali si propagano con velocita infinita (e quindi l’in-formazione che parte dai punti di (a, b) arriva contemporaneamente in unpunto dello spazio-tempo). Cio comporta che gli schemi numerici che pos-sono essere utilizzati devono essere sufficientemente veloci nel trasportare leinformazioni passando da un’istante temporale al successivo, altrimenti lasoluzione numerica risutera instabile.

2.2.1 Approssimazione numerica.

Introduciamo nel dominio computazionale (0, T )X(a, b) una griglia di nodiequispaziati. Precisamente, fissati due numeri N e M , si calcolano i passidi discretizzazione h = b−a

Ne k = T

Mrispettivamente in spazio e tempo e si

definiscono gli N + 1 nodi xj = a + j ∗ h j = 0, 1, ..., N e gli M + 1 noditn = nk n = 0, 1, ...,M .

Vogliamo ora calcolare il valore approssimato uj,n della soluzione del nostroproblema nei nodi (xj, tn) j = 0, 1, ..., N n = 0, 1, ...,M .

Dalle condizioni al bordo e dalle condizioni iniziali ricaviamo che:

uj,0 = f(xi) j = 0, 1, ..., N

u0,n = l(tn), uN,n = r(tn) n = 1, ...,M

31

Restano quindi da trovare delle equazioni che permettono di calcolare lasoluzione nei nodi interni al dominio di discretizzazione.La derivata seconda uxx in 2.9 puo essere approssimata con lo schema in dif-ferenze finite centrate del secondo ordine, cambiando lo schema di differenzefinite per la derivata temporale si possono allora ottenere schemi diversi comepresentiamo nelle prossime sezioni.

Metodo alle differenze esplicito.

Se scegliamo di approssimare la derivata temporale con una differenza finitadel primo ordine in avanti otteniamo le equazioni:

uj,n+1 − uj,n

k−D

uj+1,n − 2uj,n + uj−1,n

h2= 0

j = 1, 2, ..., N − 1n = 1, 2, ...,M − 1

ovvero

uj,n+1 = (1− 2r)uj,n + r(uj+1,n + uj−1,n)j = 1, 2, ..., N − 1n = 1, 2, ...,M − 1

dove r = Dkh2 .

x x

t n

xj−1 j j+1

n+1t

Figura 2.6: Dominio di dipendenza discreta per uj,n+1, schema esplicito.

Nella figura seguente, riportiamo il grafico della soluzione ottenuta con loschema proposto del problema della diffusione del calore in una sbarra portatainizialmente ad avere un salto di temperatura.Si nota che la soluzione numerica presenta delle oscillazioni nel primo grafico,mentre e corretta nel secondo grafico. In effetti lo schema proposto e stabilesolo se r ≤ 0.5; in questo caso il passo temporale e sufficientemente piccolo

32

rispetto al passo spaziale da garantire la propagazione del segnale discreto conuna velocita accurata. Lo schema e dunque condizionatamente convergente.L’errore tra la soluzione esatta e quella approssimata e dell’ordine di O(k +h2).

-5 -4 -3 -2 -1 0 1 2 3 4 5 6 0

0.2

0.4

0.6

0.8 1

Figura 2.7: Metodo esplicito, soluzioni con r = .6 e r = .4.

Metodo alle differenze implicito.

Approssimiamo ora la derivata temporale con lo schema a differenze finitedel primo ordine all’indietro. Otteniamo allora:

uj,n+1 − uj,n

k−D

uj+1,n+1 − 2uj,n+1 + uj−1,n+1

h2= 0

j = 1, 2, ..., N − 1n = 1, 2, ...,M − 1

ovvero

−ruj−1,n+1 + (1 + 2r)uj,n+1 − ruj+1,n+1 = uj,nj = 1, 2, ..., N − 1n = 1, 2, ...,M − 1

che, completato con le condizioni iniziali e al bordo, risulta essere un sistemalineare con matrice tridiagonale simmetrica e diagonalmente strettamentedominante. Il sistema e dunque sempre risolubile. In questo caso l’erroretra soluzione esatte ed approssimata si comporta come O(k) + O(h2) e loschema e incondizionatamente stabile, cioe la soluzione numerica non oscillaindipendentemente dall’ampiezza del passo temporale.

33

x x

t n

xj−1 j j+1

n+1t

Figura 2.8: Dominio di dipendenza discreta uj,n+1, schema implicito.

0

0.2

0.4

0.6

0.8 1

Figura 2.9: Metodo implicito, soluzione con r = 1.

Metodo di Crank-Nicolson e θ -metodi.

Possiamo ora ottenere degli altri schemi prendendo delle opportune combi-nazioni lineari dello schema esplicito ed implicito. Precisamente, sia θ unnumero reale compreso tra 0 e 1. Sommiamo allora l’equazione del metodoimplicito moltiplicata per θ con l’equazione del metodo esplicito moltiplicataper 1− θ, otteniamo cosı :

uj,n+1−uj,n

k−D [θ

uj+1,n+1−2uj,n+1+uj−1,n+1

h2

+(1− θ)uj+1,n−2uj,n+uj−1,n

h2 ] = 0j = 1, 2, ..., N − 1n = 1, 2, ...,M − 1

Al variare di θ si ottengono schemi diversi, in particolare per θ = 0 si ottieneil metodo esplicito, per θ = 1 si ottiene il metodo implicito presentato prima,per θ = 1

2si ottiene il cosı detto metodo di Crank-Nicolson che risulta essere

del secondo ordine sia in spazio che in tempo.

34

Per tutti gli altri valori del parametro si ottengono dei metodi impliciti edel primo ordine in tempo ma con un fattore moltiplicativo di h2 + k nellaformula dell’errore che diminuisce per θ → 1

2.

Per quanto riguarda la stabilita , si tratta di metodi incondizionatamentestabili per θ ≥ 1

2, mentre per θ < 1

2si ha la seguente condizione di stabilita :

r ≤ 1

2(1− 2θ)

2.3 Equazione di Laplace.

In quest’ultimo capitolo ci occupiamo dell’approssimazione dell’equazione diLaplace e precisamente dell’equazione:

−∆u = −∂2u

∂x2− ∂2u

∂y2− ∂2u

∂z2= 0 in Ω (2.10)

dove Ω e una regione dello spazio R3.Piu in generale l’equazione non omogenea −∆u = f prende il nome di equa-zione di Poisson; dove il termine noto f tiene conto di eventuali sorgenti inΩ. Si tratta di equazioni di tipo cosı detto ellittico e sono alla base dellamodellizzazione di molti problemi fisici.Per esempio sono equivalenti allo stato stazionario, cioe quando la soluzionenon dipende dal tempo, di problemi di diffusione, ma compaiono anche inproblemi di fluidodinamica e elettrostatica.L’equazione 2.10 puo essere completato o con delle condizioni di Dirichlet(u = g sul bordo di Ω ) o con delle condizioni di Newmann (∂u

∂νsul bordo di

Ω dove ν indica la normale uscente da Ω).

2.3.1 Approssimazione numerica nel caso bidimensio-nale.

Sia Ω = (x, y)|a ≤ x ≤ b ; c ≤ y ≤ d il dominio su cui vogliamo risolvereil problema:

−∆u(x, y) = −∂2u∂x2 (x, y)− ∂2u

∂y2 (x, y) = f(x, y) (x, y) ∈ Ω

u(x, y) = g(x, y) (x, y) ∈ ∂Ω

35

Introduciamo allora la discretizzazione di Ω con nodi equispaziati (xi, yj) =(a+ ih, c+ jk) i = 0, · · · , N − 1, j = 0, · · · ,M − 1 dove h = b−a

N−1, k = d−c

M−1,

ed approssimiamo l’operatore di Laplace con le differenze finite centrate delsecondo ordine; otteniamo cosı :

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

h2−ui,j+1 − 2ui,j + ui,j−1

k2= f(xi, yj)

i = 1, · · · , N − 2j = 1, · · · ,M − 2

(2.11)

dove ui,j∼= u(xi, yj).

Dalle condizioni al bordo, otteniamo inoltre le seguenti equazioni:

u0,j = g(a, yj) 0 < j < M − 1uN−1,j = g(b, yj) 0 < j < M − 1ui,0 = g(xi, c) 0 ≤ i ≤ N − 1ui,M−1 = g(xi, d) 0 ≤ i ≤ N − 1

Introduciamo ora un’unica numerazione dei nodi della discretizzazione. Peresempio, se supponiamo di numerare i nodi da sinistra a destra e dal bassoverso l’alto, otteniamo che il numero jj = j ∗N + i+ 1 corrisponde al nodo(xi, yj) e il nostro problema si riscrive come:

ujj = g(xi, c) jj = i+ 1, i = 0, · · · , N − 1ujj = g(a, yj) jj = j ∗N + 1, j = 1, · · · ,M − 2ujj = g(b, yj) jj = j ∗N, j = 2, · · · ,M − 1ujj = g(xi, d) jj = (M − 1) ∗N + i+ 1, i = 0, · · · , N − 1

−ujj−N

k2 − ujj−1

h2 + ujjh2+k2

h2k2 − ujj+1

h2 − ujj+N

k2 =

= f(xi, yj) jj = j ∗N + i+ 1,i = 1, · · · , N − 2j = 1, · · · ,M − 2

Si tratta pertanto di risolvere un sistema lineare di dimensione N ∗M conmatrice simmetrica, diagonalmente dominante e definita positiva. Si puoallora dimostrare che la matrice e non singolare e quindi il sistema e risolubile.In particolare la matrice ha solo 5 diagonali non nulle (la diagonale principale,la prima sottodiagonale, la prima sovradiagonale, la N-esima sottodiagonalee la N-esima sovradiagonale); si tratta quindi di un matrice sparsa (cioe

36

con un numero di elementi non nulli piccolo rispetto alla dimensione dellamatrice). Per la risoluzione del sistema e pertanto consigliato utilizzare deimetodi iterativi.

x xx

y

y

yj−1

j

j+1

i+1i−1 i

Figura 2.10: Dominio di dipendenza discreta per ui,j.

Condizioni al bordo di Newmann

Supponiamo ora che l’equazione di Laplace sia completata con delle condi-zioni al bordo di tipo Newmann. Innanzitutto, se questo tipo di condizioniviene imposto su tutto il bordo di Ω, il problema che ne deriva non ha unasola soluzione, ma infinite e tali che la differenza tra due soluzioni qualunquee la funzione costante.Se si vuole risolvere il problema numericamente, e allora necessario sostituireuna delle equazioni con una condizione che permetta di fissare una dellepossibili soluzioni; per esempio si puo assegnare il valore della soluzione inun punto qualunque della discretizzazione.Per descrivere come trattare le condizioni al bordo di Newmann, consideriamoil seguente problema:

−∆u(x, y) = f(x, y) (x, y) ∈ Ω

u(a, y) = g(a, y) c ≤ y ≤ du(x, c) = g(x, c) a ≤ x ≤ bu(x, d) = g(x, d) a ≤ x ≤ b∂u∂x

(b, y) = g1(y) c < y < d

Aggiungiamo allora alla discretizzazione proposta per il caso di Diricheletuna nuova colonna (xN , yj) j = 0, 1, ...,M − 1 e imponiamo le equazioni di

37

2.11 relative ai nodi interni anche per x = N − 1. In questo modo abbiamoaggiunto delle nuove incognite uN,j che ci servono per poter imporre la condi-zione di Newmann che, approssimata con lo schema alle differenze centrate,diventa:

uN,j − uN−2,j

2h= g1(yj) j = 1, 2, ...,M − 2

da cui ricaviamo uN,j = uN−2,j + 2hg1(yj) j = 1, 2, ...,M − 2. Questo valorepuo ora essere sostituito nella discretizzazione dell’operatore di Laplace neinodi del lato verticale di destra del dominio; ricaviamo cosı le equazioni:

−ujj−N

k2−2

ujj−1

h2+ujj

h2 + k2

h2k2−ujj+N

k2= f(xi, yj)−2

g1(yj)

h

jj = j ∗Nj = 2, · · · ,M − 1

Queste ultime equazioni andranno a sostituire le corrispondenti equazioniche nel problema di Dirichelet assegnavano il valore della soluzione nei nodi(xN−1, j), j = 1, · · · ,M − 2.

38

Bibliografia

[1] Valeriano Comincioli, “Metodi Numerici e Statistici per le ScienzeApplicate”, Casa Editrice Ambrosiana, Milano.

[2] Alfio Quarteroni, Riccardo Sacco, Fausto Saleri, “MatematicaNumerica”, Springer-Verlag.

[3] Alfio Quarteroni, “Elementi di Calcolo Numerico”, Editrice Esculapio.

39