Informatica 5 Liceo Scientifico Scienze Applicate/Risolvere Numericamente una equazione differenziale

Equazioni differenziali risolte numericamente

Impariamo ora a risolvere numericamente delle equazioni differenziali. Per capire cosa si vuole calcolare procediamo in questo modo:

• ci ricordiamo della pendenza
• calcoliamo la pendenza in un punto di una curva (concetto di derivata)
• Molto spesso il legame fra 2 grandezze fisiche è dato da una derivata scopriamo di aver già scritto delle equazioni differenziali
• usiamo il metodo di Eulero per risolvere una equazione differenziale
• usiamo metodi più sofisticati (HEUN e Runge Kutta) per risolvere numericamente le equazioni differenziali

Pendenza

Basta andare in montagna e fare quattro passi per rispolverare il concetto di pendenza,

MonteNero

le strade di montagna a differenza di quelle di pianura sono in pendenza , la pendenza esprime il rapporto fra la variazione di altitudine di un certo percorso e la lunghezza del percorso stesso. In matematica una strada con una certa pendenza la posso rappresentare con una retta inclinata, tutti i punti della retta hanno la stessa pendenza.

Pendenza di una retta

Derivata di una funzione

Quando la strada non è approssimabile con una retta possiamo approssimarla con una spezzata, e per ogni tratto possiamo calcolare la pendenza , quando però la strada si fa decisamente curva ogni "tratto" incomincia ad assumere una pendenza diversa. Nel caso di una curva generica (descritta matematicamente dalla funzione f(x) )possiamo calcolare la pendenza in ciascun punto con questo metodo:

Concetto di derivata in un punto

pensiamo di dover calcolare la pendenza nel punto P di coordinate x e y=f(x), ingrandiamo la curva in un intorno del punto P e prendiamo un secondo punto P1 nelle vicinanze di coordinate x+Δx e f(x+Δx) , ora definiamo la pendenza nel punto P come la pendenza della retta passante per i punti P e P1 , naturalmente se prendiamo un Δx molto piccolo le cose vanno meglio perché la retta P P1 meglio approssima la curva attorno al punto P , se Δx tende a zero la retta considerata è la retta tangente alla curva f(x) nel punto P. I matematici definiscono questa pendenza nel punto P(x,y) come la derivata della funzione f(x) calcolata nel punto x , in matematica la definizione di derivata è $f'(x)=\lim _{\Delta x\rightarrow 0}{\frac {f(x+\Delta x)-f(x)}{\Delta x}}$ i simboli per indicare la derivata possono essere $f'(x)$ oppure ${\frac {\mathrm {d} }{\mathrm {d} x}}f(x)$

Quando si parla di derivata bisogna distinguere la derivata di una funzione calcolata in un particolare punto xo ( cioe' un numero che esprime la pendenza della retta tangente alla curva passante per quel punto) e si scrive

$f'(x_{0})$ oppure $\left.{\begin{matrix}{\frac {\mathrm {df(x)} }{\mathrm {d} x}}\end{matrix}}\right|_{x_{0}}$

e la derivata di una funzione f(x) ( cioe' una funzione indicata con f'(x) che esprime la pendenza di ciascun punto della curva al variare della x)che si scrive con $f'(x)$
oppure con ${\frac {\mathrm {d} }{\mathrm {d} x}}f(x)$

Spazio,Velocita' e Accelerazione

Ritornando alla fisica del primo anno si diceva che

l'accelerazione è la variazione di velocità nell'unita di tempo, cioè la funzione che descrive l'accelerazione è la derivata della funzione che ne esprime la velocità (in questo caso l'asse delle ascisse esprime invece di un generico x il tempo , e si ha f(t) , f'(t) etc)

acc(t) = v'(t)

la velocità è la variazione dello spazio nell'unità di tempo , cioè la funzione che descrive la velocità è la derivata della funzione che ne esprime lo spostamento (spazio al variare del tempo)

v(t)= s'(t)

Sono tantissime le grandezze fisiche legate ad altre mediante il concetto di derivata Ora queste leggi sono dei particolari casi di equazioni differenziali , cioe' delle leggi che legano fra loro delle funzioni (non singoli valori) mediante l'operazione di derivata ( e generalizzando mediante la derivata della derivata $f^{''}(t)$ , o derivata della derivata della derivata $f^{'''}(t)$ ).

noi vogliamo studiare delle funzioni differenziali particolari dove y'(t) = f( y(t),t).

Dalla y'(t) si vuole ricavare la y(t) , se la y' dipende solo dalla t come succede nel moto uniforme o nel moto uniformemente accelerato , allora è possibile mediante l'integrazione calcolare data la funzione acc(t) la funzione della velocità, oppure data la funzione velocità calcolare la funzione che descrive lo spazio percorso .

Quando si deve risolvere la y'(t) =f(t,y(t)) le cose si complicano, se non si è in grado di risolvere il problema possiamo ricavare la y(t) per via numerica discretizzando l'asse temporale (cioè solo per particolari valori del tempo).

In generale possiamo calcolare y(to+Δt) = y(to) + y'(to)*Δt cioè conoscendo la y(to) in un punto particolare to e la y'(to) è possibile calcolare la y(to+Δt) .

Premessa Notazione Matematica Funzioni

Una funzione è una legge che associa in modo univoco a ciascun elemento di un dominio un elemento di un secondo insieme detto codominio. Per esprimere questa legge possiamo farlo in modi diversi :

• tramite un grafico sull'asse cartesiano (che associa a ciascun punto dell'asse delle ascisse un punto dell'asse delle ordinate)
• tramite una funzione matematica y= 3*x+2;
• tramite una serie di regole (ad ogni alunno di questo istituto associo il computer nel laboratorio che ha come numero di postazione la posizione occupata dallo studente nell'elenco alfabetico degli studenti di ciascuna classe)

Quando si utilizza una funzione in matematica di solito la si scrive simbolicamente come y=f(x) , dove la x rappresenta un generico elemento del dominio, la f la legge che associa a un elemento del dominio uno del codominio , e la y il generico elemento del codominio associato al valore x del dominio.

Naturalmente possiamo anche cambiare i simboli y=f(t) oppure temperatura= riscaldamento(termosifoni) oppure d=g(m) Con la notazione y=f(x) si indica come variabile indipendente la x (dominio) e come variabile dipendente la y

In concreto una funzione del tipo y=f(x) potrebbe essere f(x)=3*sin(x)+x^2 che stabilisce che dato un particolare valore di x ad esempio x0 per trovare il corrispondente valore associato y0 , basta calcolare 3*sin(x0)+x0^2

Bisogna stare attenti che quando usiamo i simboli x,y,f(x) stiamo indicando la regola generale f(x) applicata a generici valori y e x , mentre quando usiamo i simboli y0,x0,f(x0) stiamo parlando di particolari valori numerici associati fra loro dalla medesima legge f.

Naturalmente si possono costruire funzioni che invece di una sola variabile indipendente lavorano su piu' variabili indipendenti e che associano a questo gruppo di variabili indipendenti un elemento del codominio, possiamo allora scrivere z=f(x,y) oppure temperatura=f(x,y,z) ad esempio pensiamo a una stanza nella quale la temperatura dell'aria non sia uniforme allora possiamo scrivere temperatura=f(x,y,z) intendendo che il valore della temperatura dipende dalla coordinate x,y,z dei diversi punti della stanza secondo una certa legge rappresentata simbolicamente da f.

Anche in questo caso possiamo calcolare il valore per particolari valori e allora la formula diventa z0=f(x0,y0) oppure z3=f(x3,y3) etc.

Quindi x è un valore generico x0 x1 x2 x3 sono valori particolari (anche se non ancora specificati numericamente).

• se scopriamo che la grandezza fisica y è legata alla grandezza fisica x nel seguente modo y=3*x+45 allora simbolicamente la possiamo indicare con y=f(x)
• se scopriamo che la grandezza fisica y dipende dalla grandezza fisica tempo nel seguente modo y=t^2-5 allora abbiamo y=f(t)
• se scopriamo che la grandezza fisica r dipende dalla grandezza fisica b e abbiamo già usato la lettera f per esprimere un'altra legge possiamo scrivere r=g(b)
• se scopriamo che una certa grandezza fisica z dipende da altre 2 grandezze fisiche x e y nel seguente modo z=3*x+y possiamo indicare che la nostra legge è del tipo z=f(x,y)

• ora se la legge la indichiamo con y=f(t) allora la legge f potrebbe essere f(t) = t+5* t^2
• se indichiamo la legge con y=f(a,b) allora la legge f potrebbe essere f(a,b)= a*b+4-a

certe volte si dà alla funzione lo stesso nome della variabile indipendente cioè si scrive y=y(t) creando un po' di ambiguità risolta facilmente perché se scrivo y intendo il valore generico del codominio e quando scrivo y(t) intendo la legge generale.

pensiamo che esista una funzione y(t) , possiamo scrivere z=f( y(t) , t) cosa intendiamo con questa scrittura? che la variabile generica del codominio è indicata con z, che la legge è la f e che questa dipende da due variabili una che è t e l'altra è una certa y=y(t) cioè invece di essere una semplice variabile y questa è a sua volta calcolata tramite una seconda legge y(t) è come aver scritto

y=y(t)= ad esempio 3*t^2
z=f(y,t)=f(y(t),t) = ad esempio  t+7*y(t) 

non si sostituisce nel calcolo della f il valore y(t) = 3*t^2 perché magari si è dato un particolare significato logico all'espressione 3*t^2 e si vuole nel calcolo della z fare riferimento a questo significato logico senza ridurlo alla sola dipendenza temporale, naturalmente se conosco la y(t) (non sempre la si conosce e qualche volta viene indicata simbolicamente per questo motivo) allora sostituendo si ha z=f(y(t),x) = f(t) .

Quando si scrivono delle equazioni differenziali si scrivono delle equazioni che coinvolgono delle funzioni ( tramite gli operatori di derivazione ), ad esempio si scrive

f'(t)=f(t)   in questo caso si cercano tutte le funzioni f(t) la cui funzione derivata f'(t) ha la stessa forma (grafico, legge) della f(t).

oppure f'(t)= 1/f(t)  in questo caso si cercano tutte le possibili funzioni f(t) la cui funzione derivata ha la forma 1/f(t)  

Metodo di Eulero Risolvere una equazione differenziale con Eulero

Con Eulero possiamo trovare la soluzione a problemi differenziali del tipo y'(t)=f(y(t),t) cioe' ad esempio di problemi del tipo y'(t) = y(t)+ 5-t; stiamo attenti alla lettura della formula scritta y'(t)=f(y(t),t) e' una equazione differenziale perché

• gli elementi che la compongono sono funzioni la prima e' la y'(t) e la seconda e' la f(y(t),t)
• e perché la funzione y'(t) e' calcolata applicando l'operatore di derivazione alla y(t),

della formula y'(t)=y(t) vogliamo trovare tutte le funzioni y(t) le cui funzioni derivate siano uguali alla funzione f , la funzione f dipende dal valore di due variabili , la variabile t e la variabile y=y(t) questa variabile y dipende dal tempo tramite la legge y(t) e non la possiamo scrivere direttamente come espressione del tempo perché la y(t) non e' nota ( e' quello che vogliamo calcolare)

In pratica vogliamo calcolare una funzione y(t) che abbia una particolare derivata y'(t)(funzione pendenza) che dipende secondo la legge f dalla funzione stessa y(t), in questo non c'e' nulla di strano perché la derivata di una funzione dipende sempre dalla forma della funzione stessa .

Ad esempio nel decadimento radioattivo se N(t) esprime il numero di atomi radioattivi del carbonio-14 al variare del tempo quanto veloce e' la trasformazione del carbonio-14 in azoto-14 ? la velocita' cioe' la variazione di N(t) nell'unita' di tempo in un certo istante temporale t0 (cioe' la derivata) e' legata proporzionalmente al numero di atomi radioattivi nell' istante t0 cioe' a N(t0) tramite una costante di proporzionalita' $\lambda$ detta costante di decadimento , con un meno davanti perché N(t) invece di crescere nel tempo diminuisce.

Ora notata questa caratteristica , cioe' il legame fra la grandezza fisica e la sua variazione con Eulero e' possibile anche se non conosciamo la formula analitica di N(t) calcolarne il valore in tanti istanti temporali spaziati fra di loro di una quantita' che decidiamo noi indicata con $\Delta t$ , per risolvere il problema con Eulero bisogna conoscere oltre alla ${\frac {\mathrm {d} }{\mathrm {dt} }}N(t)=-\lambda N(t)$ che e' del tipo y'(t)=f(y(t),t) anche il valore assunto dalla N(t) in un particolare punto t0 cioe' se scegliamo come t0 l'istante iniziale in cui incomincia il decadimento radioattivo bisogna conoscere N(0) che simbolicamente indichiamo con N0

Con Eulero la N(t) viene calcolata negli istanti $t0,t0+\Delta t,t0+2*\Delta t,....,t0+k*\Delta t,t0+(k+1)*\Delta t,...$

che nel caso del decadimento radioattivo del carbonio 14 si semplifica perché t0=0 in $0,\Delta t,2*\Delta t,....,k*\Delta t,(k+1)*\Delta t,...$
Eulero calcola $N(t_{0}+(k+1)*\Delta t)$ in funzione del valore $N(t_{0}+k*\Delta t)$ cioe' quello assunto nel istante temporale precedente tramite la formula

$N(t_{0}+(k+1)*\Delta t)=N(t_{0}+k*\Delta t)+\Delta t*N'(t_{0}+k*\Delta t)$
dove $N'(t_{0}+k*\Delta t)=-\lambda *N(t_{0}+k*\Delta t)$

che possiamo riunire in

$N(t_{0}+(k+1)*\Delta t=N(t_{0}+k*\Delta t)+\Delta t*(-\lambda *N(t_{0}+k*\Delta t))$

ora facciamo 2 calcoli per capire come funziona la tecnica di calcolo che e' iterativa , fissiamo $t_{0}=0$ N0=10000 e $\Delta t=50$ e $\lambda =0.0035$ si parte da N(0)=N0=10000 calcoliamo $N(\Delta t)$

$N(\Delta t)=10000+50*(-0.0035*10000)=8250$

poi da $N(\Delta t)=$ calcoliamo $N(2*\Delta t)=$

$N(2\Delta t)=8250+50*(-0.0035*8250)=6806$

poi da $N(2*\Delta t)=$ calcoliamo $N(3*\Delta t)=$

$N(3\Delta t)=6806+50*(-0.0035*6806)=5614$

Ora il programma in Octave diventa

N0=10000;
delta=50;
lambda=0.0035;
t=0:delta:2000;
y2(1)=N0;
for i=1:40
  y2(i+1)=y2(i)+delta*(-lambda*y2(i));
end
plot(t,y2,"2;eulero;");

La vecchia spiegazione del metodo di eulero era Sia y' = f( y(t),t) e pensiamo di conoscere y(to) vogliamo calcolare y(to+Δt) , y(to+2*Δt) , y(to+3*Δt) ... y(to+n*Δt) Pensiamo che l'asse delle ascisse esprima il tempo ( t al posto di x) Partendo da un generico punto P di una curva f(t) è possibile calcolare (in modo approssimato) un secondo punto y(t+Δt) ( nelle vicinanze di P) sulla curva spostandosi lungo la retta tangente alla curva in P ( la retta che esprime la pendenza per Δt molto piccolo).

Partiamo dalla definizione di derivata di una funzione y'(t) e la approssimiamo :

$y'(t)\approx {\frac {y(t+\Delta t)-y(t)}{\Delta t}}$

che ci permette di calcolare y(t+Δt) come

$y(t+\Delta t)\approx y(t)+\Delta ty'(t)$

visto che y' la possiamo esprimere come funzione del tempo t e della funzione y(t) otteniamo :

$y(t+\Delta t)\approx y(t)+\Delta t*f(t,y(t)).$

Naturalmente calcolato y(to+Δt) e' possibile iterativamente calcolare y(to+2*Δt) y(to+3Δt) e così via (l'errore commesso aumenta a ogni passo rendendo la procedura poco affidabile per predire il valore della y(t) per t distanti da to ) posto:

t₀, t₁ = t₀ + Δt, t₂ = t₀ + 2Δt, …

allora y(t_n). si calcola ricorsivamente con

$y_{n+1}=y_{n}+\Delta t*f(t_{n},y_{n}).$

HEUN

Il metodo Heun e' una versione migliorata del metodo di Eulero è indicato anche come metodo Runge-Kutta di ordine 2 e serve per risolvere delle equazioni differenziali

$y'(t)=f(t,y(t))$ partendo da un punto t0 dove si conosce che $y(t_{0})=y_{0},$

si calcola il punto $y_{i+1}$ in due fasi: prima si calcola una stima di $y_{i+1}$ che indichiamo con il simbolo ${\tilde {y}}_{i+1}$

${\tilde {y}}_{i+1}=y_{i}+\Delta t*f(t_{i},y_{i})$

e poi si calcola $y_{i+1}$ con la

$y_{i+1}=y_{i}+{\frac {\Delta t}{2}}[f(t_{i},y_{i})+f(t_{i+1},{\tilde {y}}_{i+1})],$

dove $\Delta t$ è il passo adottato nella discretizzazione

con $t_{i+1}=t_{i}+\Delta t$.

se usiamo Eulero e la funzione y(t) e' concava attorno a to il calcolo di y(to+Δt) tramite la pendenza della retta tangente è sottostimata e viceversa se e' convessa viene sovrastimata , Heun introduce una compensazione a questo errore.

Runge-Kutta

Il metodo di Runge–Kutta (RK4)sviluppato nel 1900 dai matematici tedeschi C. Runge and M. W. Kutta. migliora ancora la compensazione dell'errore , e risulta il metodo numerico generalmente impiegato . Esegue il calcolo di y(to+Δt) in quattro fasi.Se il problema e' del tipo

$y'(t)=f(t,y(t)),\quad y(t_{0})=y_{0}.$

Dove la y(t) non e' nota se non nel punto to e si conosce y'(t) come funzione di y(t) e t

La soluzione si calcola solo nei punti to , to+Δt, to+2*Δt, etc cioe' con passo Δt. Quando le soluzioni vengono calcolate solo in specifici punti dell'asse temporale si dice che l'asse temporale e' stato discretizzato (con passo Δt)

il valore y(to+(n+1)*Δt) si calcola utilizzando il valore y(to+nΔt) tramite la formula

${\begin{aligned}y_{n+1}&=y_{n}+{\tfrac {\Delta t}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)\\t_{n+1}&=t_{n}+\Delta t\\\end{aligned}}$

dove k1 k2 k3 e k4 sono calcolati come

${\begin{aligned}k_{1}&=f(t_{n},y_{n}),\\k_{2}&=f(t_{n}+{\tfrac {\Delta t}{2}},y_{n}+{\tfrac {\Delta t}{2}}k_{1}),\\k_{3}&=f(t_{n}+{\tfrac {\Delta t}{2}},y_{n}+{\tfrac {\Delta t}{2}}k_{2}),\\k_{4}&=f(t_{n}+\Delta t,y_{n}+\Delta tk_{3}).\end{aligned}}$

Se la f'dipende solo da t allora l'equazione differenziale si riduce al calcolo di un semplice integrale e la RK4 e' equivalente alla regola di Simpson per il calcolo degli integrali definiti.