In questo articolo andremo a parlare di differenze finite. Questo sarà un articolo introduttivo all’argomento.
Oltre alla descrizione del metodo vedremo un paio di esempi molto semplici scritti con Matlab, dove andremo a risolvere l’equazione di Poisson su un intervallo $I\subset\mathbb{R}$ e una sua variante.
Se vuoi vedere anche un video che ho fatto su questo argomento lo trovi sul canale Youtube 😉
Di sicuro ti è stato detto o comunque hai studiato e letto da qualche parte che è davvero piccolo l’insieme di equazioni differenziali risolvibili in maniera analitica ed esatta. Molto poche ammettono una soluzione esprimibile tramite una funzione che ha una sua espressione precisa. Descrivibile in forma chiusa.
Per questo motivo è necessario trovare un’alternativa alla procedura analitica. La procedura esatta che ci permette di arrivare ad una soluzione delle equazioni è infatti spesso limitante.
Ok, è importante saper risolvere gli esercizi in cui viene chiesto di trovare un integrale generale di un’equazione differenziale, ma questi sono appunto esercizi. Spesso le equazioni che definiscono un modello matematico, una volta che si riesce a mostrare che una soluzione esista, sono trattati in termini numerici.
Infatti tutto quello che è presente nel mondo, nella realtà, è descritto da una variazione di certe quantità, di certe proprietà mentre il tempo scorre liberamente.
Quindi come possiamo descrivere tutti questi fenomeni? Beh, intanto dobbiamo necessariamente coinvolgere delle equazioni differenziali. Quindi chiaramente non possiamo fermarci davanti al fatto che non sia possibile risolvere un’equazione di questo tipo in maniera esatta, in forma chiusa.
Se a noi interessa fare previsioni su qualche modello, su qualche fenomeno, dobbiamo trovare comunque un modo per ottenere informazioni sulla soluzione. A questo punto si aprono due strade molto interessanti:
L’analisi qualitativa (di cui magari ci occuperemo in altri articoli e puoi trovare già un esempio in questo articolo https://www.mathone.it/pendolo-semplice/) ma puoi già trovare un mio video sull’argomento qui di seguito:
e l’approssimazione numerica della soluzione, argomento di cui inizieremo ad occuparci proprio in questo articolo.
Per questa prima introduzione parleremo di equazioni differenziali ordinarie, quindi del caso in cui compaiono solo derivate ordinarie e c’è una sola variabile. Non andiamo quindi a coinvolgere le equazioni alle derivate parziali anche perché in quel caso il metodo alle differenze finite è abbastanza limitante perché non è comodo per lavorare con domini di dimensioni di forma particolari perché è necessario avere delle griglie fatte in un certo modo (spesso in quel caso si usa il metodo degli elementi finiti).
Comunque di sicuro porterò qualche esempio riguardo al metodo applicato al caso delle derivate parziali perchè permette di analizzare, senza andare troppo nel complesso, sistemi che evolvono in spazio e tempo, ampliando così di molto la classe dei modelli che potremo analizzare.
Ma torniamo alle equazioni differenziali ordinarie. Questa tipologia di equazioni solitamente ci interessa risolverle in un certo dominio. Per poterle risolvere numericamente dobbiamo imporre un’importante condizione su questo dominio: deve essere limitato.
Numericamente infatti non possiamo direttamente risolvere un’equazione differenziale su tutta la retta reale, ma dobbiamo considerarne un sottodominio compatto della forma $[0,L]$ con $L<+\infty$.
Le differenze finite si prestano molto bene a risolvere equazioni differenziali ordinarie che descrivono fenomeni stazionari, ovvero nel caso le quantità non varino nel tempo ma nello spazio. Spesso ci si riferisce ad essi come problemi al bordo (boundary value problems o BVP). Per esempio parliamo dell’equazione di Poisson $$-\frac{d^2u(x)}{dx^2}= f(x)$$ con $x\in[0,L]$ ed $f:\mathbb{R}\rightarrow \mathbb{R}$ una funzione qualsiasi, con delle opportune condizioni al bordo $u(0)=a$ e $u(L)=b$.
Chiaramente vediamo subito che serve un minimo di regolarità per la funzione $f$ per poter dire di avere una soluzione classica in questo esempio, ovvero siccome dobbiamo calcolare la derivata seconda di $u$ si può richiedere di avere $u\in\mathcal{C}^2(\mathbb{R},\mathbb{R})$ e quindi segue naturalmente $f\in\mathcal{C}^0(\mathbb{R},\mathbb{R})$. Se hai già sentito parlare di soluzioni deboli sai che in realtà si può chiedere meno regolarità per $f$ in generale, ma non preoccupiamocene per questo articolo.
Andiamo quindi ad introdurci alle tecniche approsimative che ci porteranno a definire uno schema alle differenze finite per risolvere l’equazione differenziale precedente, che possiamo per esempio complicare anche passando a questa dove comprare anche la derivata prima :$$\frac{d^2u(x)}{dx^2}+\frac{du(x)}{dx}=f(x)$$ per ogni $x\in[0,L]$ e ancora delle buone condizioni al bordo.
La prima idea che dobbiamo avere per approcciare l’approssimazione di una derivata, perché è questo che vogliamo fare, con delle strategie alternative è quello di definire una discretizzazione del nostro dominio.
Cos’è una discretizzazione? Semplicemente vogliamo dividere il nostro intervallo $[0,L]$ in intervallini sufficientemente piccoli, la cui unione restituisce esattamente l’intero dominio:
Definiamo la discretizzazione $$\tau = \{x_1=0<x_2<…<x_{N-1}<x_N=L\}$$ in modo dale che $$\cup_{i=1}^{N-1} [x_i,x_{i+1}] = [0,L].$$
Questa discretizzazione può essere fatta in maniera uniforme o non uniforme nel senso che le distanze $$\Delta x_k = x_{k+1}-x_k$$ possono essere rispettivamente tutte uguali o diverse.
Per semplicità d’ora in poi nella trattazione e anche nel codice supporremo tale discretizzazione uniforme e chiamiamo quindi $\Delta x = x_{n+1}-x_n$.
Benissimo ora siamo pronti a fare lo step fondamentale dietro l’idea delle differenze finite.
Se ti ricordi un po’ come ti è stato introdotto il concetto di derivata, ti ricorderai senz’altro che è coinvolto il limite del rapporto incrementale
Infatti la derivata di una funzione $u:\mathbb{R}\rightarrow\mathbb{R}$ è definita come il limite del rapporto incrementale, qualora esso sia finito:
$$ u'(x) = \lim_{\Delta x\to 0}\frac{u(x+\Delta x)-u(x)}{\Delta x}. $$
Allo stesso modo si può definire anche la derivata seconda:
$$ u”(x) = (u'(x))’ = \lim_{\Delta x\to 0}\frac{u'(x+\Delta x)-u'(x)}{\Delta x}.$$
Introduciamo quindi l’ultima notazione: $u(x_k) \approx u_k$, ovvero noi quello che andremo a calcolare sarà $u_k$ che è l’approssimazione numerica della soluzione esatta in $x_k$.
Ottimo, ora possiamo finalmente fornire un’approssimazione alle differenze finite di queste due derivate. Per farlo basta la semplice idea: invece di passare al limite su $\Delta x$, definiamo una discretizzazione sufficientemente raffinata del dominio $[0,L]$, ovvero tale che gli elementi $x_k$ e $x_{k+1}$ distano sufficientemente poco.
Ecco quindi una prima approssimazione della derivata
$$ u'(x_k)\approx \frac{u_{k+1}-u_k}{\Delta x}.$$
Questa è però una stima un po’ rozza infatti si può mostrare, espandendo con i polinomi di Taylor, che $$|u'(x_k)-\text{questa approssimazione}|$$ va a zero con la stessa velocità con cui ci va $\Delta x$, quindi è un’approssimazione di ordine 1:
$$ u(x_{k+1}) = u(x_k) + u'(x_k)(x_{k+1}-x_k) + o(\Delta x^2) $$
$$\frac{1}{\Delta x}(u(x_{k+1})-u(x_k)) = \frac{1}{\Delta x} (u(x_k)+u'(x_k)\Delta x + o(\Delta x^2)-u(x_k))$$
$$ = u'(x_k) + o(\Delta x).$$
Un modo per ottenere un’approssimazione più precisa, del secondo ordine, è quello di procedere con una differenza finita centrata invece che in avanti come abbiamo visto prima. Ti chiedo di provare a verificare da solo che la prossima approzione è precisa al secondo ordine 😉
$$ u'(x_k) \approx \frac{u_{k+1}-u_{k-1}}{2\Delta x}. $$
Similmente, partendo dalla differenza finita in avanti vista prima, si può ottenere un’approssimazione accurata al secondo ordine della derivata seconda come segue:
$$ u”(x_k) \approx \frac{\frac{u_{k+1}-u_k}{\Delta x} – \frac{u_k-u_{k-1}}{\Delta x}}{\Delta x} $$
$$ = \frac{u_{k+1}-2u_k+u_{k-1}}{\Delta x^2}.$$
Ottimo, direi che con la “teoria” siamo a posto. Vediamo di applicare questi risultati alle due equazioni prima introdotte. Prima però è importante rimarcare il fatto che il risultato del metodo delle differenze finite sarà un vettore che corrisponde alle approssimazioni della soluzione dell’equazione analizzata nei nodi della discretizzazione. Otterremo quindi un vettore $\vec{u}\in\mathbb{R}^N$ definito come segue:
$$ \vec{u} \approx \begin{bmatrix} u(x_1) \\ u(x_2) \\ . \\ .\\ . \\ u(x_N) \end{bmatrix} $$
e solitamente quando rappresenteremo graficamente la soluzione ottenuta si costruisce un’interpolazione lineare di tali valori, ovvero negli intervalli $(x_k,x_{k+1})$ si congiungono i punti $(x_k,u_k)$ e $(x_{k+1},u_{k+1})$ con un segmento come puoi vedere qui sotto:
Ottimo direi che la parte introduttiva può dirsi chiusa, qui di seguito oltre ai codici che ho scritto per Matlab e che puoi scaricare cliccando sul link di GitHub, ti riporto l’idea in breve dietro l’implementazione. La cosa interessante da precisare per il codice è che l’ho scritto in forma matriciale. Ho definito quindi due matrici $D1$ e $D2$ in modo tale che la loro azione sul vettore $u$ permetta di ottenere rispettivamente l’approssimazione della derivata prima e della derivata seconda. Ecco qui cosa intendo:
$$D1\, u = \begin{bmatrix}
0& 0&0& \dots &\dots \\
-1/(2\Delta x)& 0& 1/(2\Delta x)&0&\dots\\
0& -1/(2\Delta x)& 0& 1/(2\Delta x)&0\\ 0&\ddots&\ddots&\ddots&\dots\\
\dots& \dots& \dots& 0& 0
\end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ . \\ .\\ . \\ u_N \end{bmatrix} = \begin{bmatrix} 0 \\ \frac{u_3-u_1}{2\Delta x} \\ . \\ .\\ \frac{u_N-u_{N-2}}{2\Delta x}\\ 0 \end{bmatrix} $$
$$D2\, u = \frac{1}{\Delta x^2}\begin{bmatrix}
1& 0&0& \dots &\dots&\dots&\dots \\
1& -2& 1&0&0&0&\dots\\
0& 1& -2&1&0&\dots&\dots\\ 0&\ddots&\ddots&\ddots&\dots&\dots&\dots\\
\dots& \dots&\dots&\dots& \dots& 0& 1
\end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ . \\ .\\ . \\ u_N \end{bmatrix} = \begin{bmatrix} \frac{u_1}{\Delta x^2} \\ \frac{u_3-2u_2+u_1}{\Delta x^2} \\ . \\ .\\ \frac{u_N-2u_{N-1}+u_{N-2}}{\Delta x^2}\\ \frac{u_N}{\Delta x^2} \end{bmatrix} $$
E quindi i due problemi si ridurranno semplicemente a risolvere un sistema lineare. Il problema
$$\begin{cases}-u”(x) = 1,\quad x\in(0,1) \\ u(0)=u(1)=0\end{cases} $$
diventa quindi
$$ -D2\, u = \begin{bmatrix} 0\\ 1 \\ . \\ .\\ 1\\ 0 \end{bmatrix} .$$
Ecco quindi cosa otteniamo, dove la soluzione analitica con cui ho comparato quella numerica è la seguente parabola $$u_{\text{esatta}} (x)= -\frac{1}{2}x(1-x).$$
Invece il secondo problema
$$\begin{cases}u”(x)+u'(x) = 0,\quad x\in(0,1) \\ u(0)=0,u(1)=1\end{cases} $$
diventa quindi
$$ (D2+D1)\, u = \begin{bmatrix} 0\\ 0 \\ . \\ .\\ 0\\ \frac{1}{\Delta x^2} \end{bmatrix} .$$
Ecco quindi cosa otteniamo, dove la soluzione analitica con cui ho comparato quella numerica è la seguente $$u_{\text{esatta}} (x)= \frac{e-e^{1-x}}{e-1}.$$
Lascia un commento