SOLUZIONE NUMERICA DELLE EQUAZIONI DIFFERENZIALI

Le equazioni differenziali sono equazioni che legano una funzione e le sue derivate. Soluzione delle equazioni differenziali è una funzione detta funzione integrale dell'equazione differenziale data. La forma generale delle equazioni differenziali è:

ϕ(x,y,y,y,.....,yn)=0ϕ\left(x,y,y',y'',.....,y^n\right)=0
e l'ordine maggiore della derivata presente dà l'ordine alla equazione differenziale.
Le equazioni differenziali ordinarie (ODE) sono quelle in cui appare una sola variabile indipendente e quindi solo derivate totali. Se le variabili indipendenti coinvolte sono più di una, e quindi appaiono anche derivate parziali, si hanno le equazioni differenziali parziali (PDE).
I metodi analitici per la soluzione delle ODE permettono di ricavare la funzione y= f(x) se sono note n costanti di integrazione, ovvero i valori iniziali della funzione incognita y0 = f(x0) e delle sue derivate y'0 = f'(x0), y''0 = f''(x0), ....., yn-10 = fn-1(x0) (i metodi analitici sono detti  "problema di Cauchy").
I metodi numerici forniscono invece solo la conoscenza di un insieme finito di yi in corrispondenza ad un insieme finito di xi .  Quindi la soluzione numerica sarà tanto più vicina alla soluzione analitica quanto più ampio sarà l'insieme delle xi . L’esigenza di avere dei metodi numerici per la risoluzione delle ODE e PDE è  motivata dal fatto che solo una piccola parte delle ODE e PDE  possono essere risolte attraverso l’analisi matematica classica, anche se in realtà non è l’unica ragione infatti è importante sottolineare che spesso anche la valutazione di una soluzione analitica richiede una approssimazione numerica e quindi affetta da errore, e anche ottenere questa approssimazione può avere un elevato costo computazionale.

Una prima importante considerazione è che "tutte le ODE di ordine n si possono ridurre a un sistema di n ODE di primo ordine".
Per esempio per una ODE del 2° ordine, introducendo la variabile intermedia y' si ottiene:

d2ydx2=f(x,y){dydx=y(x,y)dydx=f(x,y)\frac{d^2y}{dx^2}= f(x,y)⇒\begin{cases} \frac{dy}{dx}=y'\left(x,y\right) \\ \frac{dy'}{dx} = f\left(x,y\right)\end{cases}
Con ODE di ordine superiore possono essere introdotte nuove variabili per le derivate seconde, terze ecc.
Conseguentemente  i metodi che funzionano per le ODE di primo ordine possono essere applicati per risolvere ODE di qualsiasi ordine.
Per cui d'ora in poi si parlerà di ODE del primo ordine della forma:
y(x,y)=dydx=f(x,y)y'(x,y)=\frac{dy}{dx}=f(x,y)
I metodi numerici per la soluzione delle ODE si possono dividere in due diversi tipi: metodi one-step e metodi multi-step.
Nei metodi one-step il valore yi viene calcolato utilizzando unicamente l’approssimazione precedente yi-1; nei metodi multi-step invece il valore yi viene calcolato utilizzando più approssimazioni precedenti.
Inoltre i metodi possono essere di carattere esplicito o implicito. Si dice che il metodo è esplicito se il calcolo di yi non è funzione di yi-1, altrimenti si dice implicito.

METODO DI EULERO
Il metodo di Eulero è il più semplice metodo one-step esplicito per risolvere una ODE. Questo metodo assume che, entro un piccolo incremento della variabile indipendente dal punto iniziale x0 , la derivata y' si mantenga costante al valore f(x0,y0) :
f(x0,y0)=[dydx]x=x0f(x_0,y_0)= \left[\frac{dy}{dx}\right]_{x=x_0}e,nota la y', :

ΔyΔx=f(x0,y0)Δy=f(x0,y0)Δxy1-y0=f(x0,y0)Δx\frac{Δy}{Δx}=f(x_0,y_0)⇒Δy=f(x_0,y_0)\cdot\Delta x⇒ y_1-y_0=f(x_0,y_0) \cdot \Delta x⇒y1=y0+f(y0,x0)Δx⇒y_1=y_0+f(y_0,x_0) \cdot \Delta

Il processo può essere ripetuto usando (x1, y1) come punto di partenza in modo da calcolare la y'(y1, x1) e quindi una nuova y2 (è uso scrivere h al posto di Δx):

y2=y1+hf(x1,y1)y_2= y_1 + h \cdot f(x_1,y_1)
E così via. La formula generale del metodo di Eulero è:
yi+1=yi+hf(xi,yi)y_{i+1}=y_i+ h \cdot f(x_i,y_i)
Si può osservare come il metodo di Eulero in fondo corrisponde all'espansione di Taylor troncata al primo ordine :
yi+1=yi+hyi+h22!y(ξ)y_{i+1}= y_i + h \cdot y'_i+\frac{h^2}{2!} y''(ξ)
l'errore di troncamento è quindi dell'ordine di O(h²).
Per migliorare l'approssimazione si può diminuire il valore di h e in questo caso diminuisce anche l'errore di troncamento in modo proporzionale ad h².  Tuttavia, poiché, per il dato insieme di x entro il quale l'integrazione è eseguita, il numero di passi è proporzionale a 1/h (n=xf-xihn=\frac{x_f-x_i}{h}) allora l'errore totale diventa del tipo O(h) ovvero proporzionale ad h.
Si verifica spesso che il metodo di Eulero non fornisce un’approssimazione accurata della soluzione e l’errore cresce allontanandosi dal punto x1 in cui è assegnata la condizione iniziale. Bisogna scegliere un passo h molto piccolo per avere un errore globale accettabile.

METODO DI HEUN
Il difetto principale del metodo di Eulero è che la pendenza all'inizio dell'intervallo f(xi,yi) non  è rappresentativa del suo valore sull'intero intervallo. Una soluzione è  sostituirla con la media aritmetica tra il valore all'inizio e quello alla fine :
yi+1=yi+h2[f(xi,yi)+f(xi+1,yi+1)]y_{i+1}=y_i+\frac{h}{2}\left[f\left(x_i,y_i\right)+f(x_{i+1},y_{i+1})\right]
Tuttavia non conosciamo la yi+1 ( infatti le derivate f(x,y) in generale sono una funzione delle x e delle y).
Una prima stima di yi+1 può essere fornita dal metodo di Eulero:
yi+10=yi+hf(xi,yi)y^0_{i+1}=y_i + h \cdot f(x_i, y_

e con questa yi+10y^0_{i+1} (predictor) calcoliamo la pendenza f(xi+1,yi+10) f(x_{i+1},y^0_{i+1}) (corrector). Adesso abbiamo le due pendenze iniziale e finale con costruire la pendenza media. La formula finale corretta del  metodo di Heun è:yi+11=yi+hf(xi,yi)+f(xx+1,yi+10)2y_{i+1}^1=y_i+h \cdot \frac{f(x_i,y_i)+f(x_{x+1},y^0_{i+1})}{2}
In questo caso l'errore locale è di tipo O(h³) per cui l'errore totale è del tipo O(h²). È anche possibile riutilizzare in modo iterativo il valore ottenuto per una seconda stima della pendenza f(xi+1,yi+11) f(x_{i+1},y^1_{i+1}) e ricalcolare la pendenza media per ottenere una nuova stima di yi+1 e così via.

METODO DEI POLIGONI O DEL PUNTO MEDIO
Nel metodo dei poligoni o del punto medio  si migliora la tecnica per calcolare il  predictor  .
Infatti si utilizza il valore della pendenza a metà intervallo utilizzando il metodo di Eulero:

yi+1/20=yi+h2f(xi,yi)y^0_{i+1/2}=y_i+\frac{h}{2} \cdot f(x_i,y_i)
e con questo predictor si calcola la pendenza f(xi+1/2,yi+1/20)f(x_{i+1/2},y^0_{i+1/2}) corretta.

La formula finale del metodo dei poligoni è:
yi+1=yi+hf(xi+1/2,yi+1/20)=yi+hf(xi+h2,yi+h2f(xi,yi))y_{i+1}=y_i+h \cdot f(x_{i+1/2},y^0_{i+1/2})= y_i + h \cdot f(x_i + \frac{h}{2}, y_i + \frac{h}{2} \cdot f(x_i,y_i))
Anche in questo caso l'errore locale è di tipo O(h³) per cui l'errore totale dovuto all'applicazione su più intervallini è O(h²).
METODI DI RUNGE-KUTTA
I metodi Runge-Kutta (in breve, RK) sono ancora metodi ad un passo come i precedenti, ma che tuttavia comportano diverse valutazioni della funzione f(x,y) in ciascun intervallo [xi,xi+1]. Possiamo ipotizzare un metodo di Eulero generalizzato della forma:

yi+1=yi+hϕ(xi,yi,h))y_{i+1}=y_i+h \cdot ϕ\left(x_i,y_i,h)\right)
dove ϕ(xi,yi,h)ϕ(x_i,y_i,h)   è una funzione  incremento che rappresenta una pendenza media pesata entro l'intervallo ed ha la forma generale:

ϕ=a1k1+a2k2+......+anknϕ= a_1k_1 + a_2k_2 + ...... + a_nk_n
le ai sono costanti e le ki sono pendenze valutate in corrispondenza di punti scelti.
La forma delle pendenze ki è:
{k1=f(xi,yi)k2=f(xi+hp1,yi+hq11k1)k3=f(xi+hp2,yi+h(q21k1+q22k2))....kn=f(xi+hpn-1,yi+h(qn-1,1k1+qn-1,2k2+...+qn-1,nkn-1))\begin{cases} k_1=f(x_i,y_i) \\ k_2 = f(x_i+h p_1, y_i + h q_{11}k_1) \\ k_3 = f(x_i+h p_2, y_i + h (q_{21}k_1+q_{22}k_2))\\ ....\\ k_n = f(x_i+h p_{n-1}, y_i + h (q_{n-1,1}k_1+q_{n-1,2}k_2+ ... +q_{n-1,n}k_{n-1} ) )\end{cases}
La particolare scelta dei coefficienti ai e delle pendenze ki fornisce un particolare metodi di Runge-Kutta (RK)
RK al primo ordine. Se a1=1a_1=1   e a2=a3=.....=an=0a_2=a_3= .....=a_n=0 si ha :  yi+1=yi+hf(xi,yi)y_{i+1}=y_i + h \cdot f(x_i,y_i) e si ottiene il metodo di Eulero ad un passo (ovvero l'espansione di Taylor al primo ordine) . 
RK al secondo ordine. Cerchiamo le costanti ai e le pendenze ki in modo da riprodurre anche l'espansione di Taylor al secondo ordine.
L'espansione di Taylor attorno il punto xi al secondo ordine è:yi+1=yi+hf(xi,yi)+h22f(xi,yi)=yi+hf(xi,yi)+h22df(xi,yi)dxi=yi+hf(xi,yi)+h22(f(xi,yi)xi+f(xi,yi)yidyidxi)y_{i+1}=y_i + h \cdot f(x_i,y_i) + \frac{h^2}{2} \cdot f'(x_i,y_i)=y_i + h \cdot f(x_i,y_i) + \frac{h^2}{2}\cdot \frac{df(x_i,y_i)}{dx_i}=y_i + h \cdot f(x_i,y_i) + \frac{h^2}{2}\left(\frac{∂f(x_i,y_i)}{∂x_i}+\frac{∂f(x_i,y_i)}{∂y_i}\cdot \frac{dy_i}{dx_i}\right)
 la formula di Runge-Kutta (RK) al secondo ordine invece è: yi+1=yi+h(a1k1+a2k2)y_{i+1}=y_i + h \cdot (a_1k_1 + a_2k_2)
con k1=f(xi,yi)k_1=f(x_i,y_i)k2=f(xi+p1h,yi+q11k1h)=f(xi,yi)+p1hf(xi,yi)xi+q11k1hf(xi,yi)yi+O(h2)k_2=f(x_i+p_1h,y_i+q_{11}k_1h)=f(x_i,y_i)+ p_1 h \frac{∂f(x_i,y_i)}{∂x_i}+q_{11}k_1h\frac{\partialf(x_i,y_i)}{∂y_i}+O(h^2) 
(ricordare che g(x+Δx,y+Δy)=g(x,y)+g(x,y)xΔx+g(x,y)yΔyg(x+Δx,y+Δy)=g(x,y)+\frac{∂g(x,y)}{∂x} Δx+\frac{∂g(x,y)}{∂y} Δy)
Quindi, sostituendo k1 e k2 in Runge-Kutta, :yi+1=yi+h(a1k1+a2k2)=yi+ha1f(xi,yi)+ha2f(xi,yi)+ha2p1hf(xi,yi)xi+ha2q11k1hf(xi,yi)yi+O(h3)=y_{i+1}=y_i + h \cdot (a_1k_1 + a_2k_2)=y_i+ha_1f(x_i,y_i)+ha_2f(x_i,y_i)+ha_2 \cdot p_1 h \frac{∂f(x_i,y_i)}{∂x_i}+ha_2\cdot q_{11}k_1h\frac{f(x_i,y_i)}{∂y_i}+O(h^3)==yi+h[a1+a2]f(xi,yi)+h2[a2p1f(xi,yi)xi+a2q11f(xi,yi)f(xi,yi)yi]+O(h3)=y_i+h\left[a_1+a_2\right]f(x_i,y_i)+h^2 \cdot\left[a_2 p_1 \frac{∂f(x_i,y_i)}{∂x_i}+ a_2q_{11}f(x_i,y_i)\frac{f(x_i,y_i)}{∂y_i}\right]+O(h^3)Adesso confrontiamo l'espansione ottenuta con l'espansione di Taylor al secondo ordine:

yi+1=yi+hf(xi,yi)+h22(f(xi,yi)xi+f(xi,yi)yidyidxi)+O(h3)y_{i+1}=y_i + h \cdot f(x_i,y_i) + \frac{h^2}{2}\left(\frac{∂f(x_i,y_i)}{∂x_i}+\frac{∂f(x_i,y_i)}{∂y_i}\cdot \frac{dy_i}{dx_i}\right)+O(h^3)
Dal confronto  si ricava che deve essere:

{a1+a2=1a2p1=12a2q11=12\begin{cases} a_1+a_2=1 \\ a_2p_1=\frac{1}{2}\\ a_2q_11=\frac{1}{2}\end{cases}
Si ottiene comunque un sistema di 3 equazioni a 4 incognite che non è determinato. Si può scegliere  il coefficiente a2, che dipende dal punto dove viene calcolata la pendenza k2.
Alcune possibili scelte:
  • Se a2=1/2 allora p1= 1 e si ottiene il metodo di Heun.

  • Se a2=1 allora p1=1/2 e si ottiene il metodo del punto medio.
  • Se a2= 2/3 allora p1= 3/4 e si ottiene il metodo detto metodo di Ralston

Runge-Kutta genera il metodo del punto medio o dei poligoni

Runge-Kutta genera il metodo di Ralston
Questi comunque sono tutti sviluppi di Taylor fino al secondo ordine.

RK al terzo ordine
. Un metodo di RK al terzo ordine genera un errore  locale è di tipo O(h4) e un  l'errore totale dovuto all'applicazione su più intervalli di tipo  O(h³).  Nel RK al terzo ordine la pendenza media pesata è:

Φ(xi,yi,h)=a1k1+a2k2+a3k3Φ(x_i,y_i,h)= a_1 k_1 + a_2 k_2 + a_3 k_3dove : {k1=f(xi,yi)k2=f(xi+p1h,yi+q11k1h)k3=f(xi+p2h,yi+q21k1h+q22k2h)\begin{cases} k_1= f(x_i,y_i) \\ k_2=f(x_i+p_1h, y_i+q_{11}k_1h)\\k_3=f(x_i+p_2h,y_i+q_{21}k_1h+q_{22}k_2h) \end{case

Adesso ci sono due pendenze da scegliere (a parte la prima che è nota) per risolvere il sistema che viene generato dal confronto con lo sviluppo di Taylor fino al terzo ordine.
La versione più utilizzata è: yi+1=yi+h[16(k1+4k2+k3)]y_{i+1}=y_i+h \cdot \left[\frac{1}{6} \left(k_1+4k_2+k_3\right)\right]dove : {k1=f(xi,yi)k2=f(xi+h2,yi+h2k1)k3=f(xi+h,yi-k1h+2k2h)\begin{cases} k_1= f(x_i,y_i) \\ k_2=f(x_i+\frac{h}{2},y_i+\frac{h}{2}k_1) \\ k_3= f(x_i+h,y_i-k_1h+2k_2h) \end{cases}
Questa versione si riduce alla regola di Simpson 1/3 se la derivata è funzione solo di x.

RK al quarto  ordine. La versione del metodo di RK del quarto ordine è data da:

yi+1=yi+h[16(k1+2k2+2k3+k4)]y_{i+1}= y_i + h \cdot \left[\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\right]
Le quattro pendenze sono così determinate: {k1=f(xi,yi)k2=f(xi+h2,yi+h2k1)k3=f(xi+h2,yi+h2k2)k4=f(xi+h,yi+k3h)\begin{cases} k_1= f(x_i,y_i) \\ k_2=f(x_i+\frac{h}{2},y_i+\frac{h}{2}k_1) \\ k_3=f(x_i+\frac{h}{2},y_i+\frac{h}{2}k_2) \\ k_4= f(x_i+h,y_i+k_3h) \end{cases}In questo caso l'errore locale è del tipo O(h5), per cui l'errore totale dovuto all'applicazione su più intervalli è O(h4)


EQUAZIONI DIFFERENZIALI ALLE DERIVATE PARZIALI

Una equazione differenziale alle derivate parziali (PDE) presenta derivate di funzioni rispetto a 2 o più variabili indipendenti.
In particolare le PDE lineari del secondo ordine di due variabili indipendenti hanno la forma generale:

A2ux2+B2uxy+C2uy2+Dux+Euy+Fu+G=0A\frac{∂^2u}{∂x^2}+B\frac{∂^2u}{∂x∂y}+C\frac{∂^2u}{∂y^2}+D\frac{∂u}{∂x}+E\frac{∂u}{∂y}+F\,u+G=0
dove i coefficienti A, B, ..., G possono essere costanti o funzioni di x e y. In effetti il comportamento delle soluzioni di questa PDE dipende dai coefficienti e in particolare dal segno del discriminante Δ=B2-4ACΔ=B^2-4\,A\,C:
PDE ELLITTICHE. Consideriamo l'equazione di Laplace nel caso della distribuzione di temperatura in una lastra conduttrice dovuta alla conduzione del calore:

2Tx2+2T2y=0\frac{∂^2T}{∂x^2}+\frac{∂^2T}{∂^2y}=0
Si vuole conoscere la temperatura in una serie di punti (nodi) equispaziati  che formano un reticolo (la variabile indipendente in più genera un piano a cui punti è associato un valore di temperatura. Una prima approssimazione è rendere discreto l'insieme di punti). In questo schema la coordinata x è la stessa per tutti i nodi che hanno lo stesso indice i, anche se hanno diverso indice j.  Allo stesso modo, tutti i nodi che hanno lo stesso indice j hanno la stessa coordinata y.  Questo rende estremamente semplice la costruzione di schemi alle differenze finite per il calcolo PDE rispetto a x o y.
Ci sono N=(xN-x1)/ΔxN=(x_N-x_1)/Δx  intervalli nelle ascisse positive e M=(yM-y1)/ΔyM=(y_M-y_1)/Δy  intervalli ngeNΔx=Δy=hnelle ordinate positive.

Una first central approximation di 2T/x2∂^2T/∂x^2 nel punto (i,j) è:

[2Tx2]i,j=Ti+1,j-2Ti,j+Ti-1,jΔx2+O(Δx2)\left[\frac{∂^2T}{∂x^2}\right] _{i,j}=\frac{T_{i+1,j}-2T_{i,j}+T_{i-1,j}}{Δx^2}+O(Δx^2)
in modo simile si ottiene la 2T/y2∂^2T/∂y^2:

[2Ty2]i,j=Ti,j+1-2Ti,j+Ti,j-1Δy2+O(Δy2)\left[\frac{∂^2T}{∂y^2}\right] _{i,j}=\frac{T_{i,j+1}-2T_{i,j}+T_{i,j-1}}{Δy^2}+O(Δy^2)
Per semplificare ancora assumiamo che Δx=Δy=hΔx=Δy=h, sostituendo le precedenti differenze finite nell'equazione di Laplace otteniamo l'equazione numerica:

Ti+1,j+Ti-1,j+Ti,j+1+Ti,j-1-4Ti,j=0T_{i+1,j}+T_{i-1,j}+T_{i,j+1}+T_{i,j-1}-4T_{i,j}=0
Con questa tecnica possiamo scrivere un'equazione numerica per ogni nodo in modo da ottenere un sistema lineare la cui soluzione fornisce una soluzione approssimata della PDE (anche perché gli errori di troncamento O(h) sono stati trascurati). Se h è sufficientemente piccolo la soluzione approssimata sarà sempre più accurata (la limitazione è data dal tempo macchina disponibile e dalla memoria del computer).
Tuttavia il sistema ancora non è determinato se non vengono definite le costanti di integrazioni ovvero le condizioni al contorno. Ci possono essere due possibili scelte:
  • Con le condizioni al contorno di Dirichlet sono note le temperature nei punti sul bordo della lastra ovvero:

    T1,1,T1,2,,T1,N,T2.1,T3,1,,TM,1,TM,2,,TM,N,T2,N,T3,N,,TM-1,NT_{1,1},T_{1,2},⋯,T_{1,N},\quadT_{2.1},T_{3,1},⋯,T_{M,1},\quadT_{M,2},⋯,T_{M,N},\quadT_{2,N},T_{3,N},⋯,T_{M-1,N}
    e occorre solo da calcolare le temperature nei N·M nodi interni (risolvendo il sistema lineare).
  • Con le condizioni al contorno di Neumann è noto il calore (flussi temici) nel bordo della lastra ovvero le derivate parziali :

T1.1y,T1,2y,,T1,My,T1.1x,T2,1x,,TN,1x,TN,1y,TN,2y,TN,My,T1,Mx,T2,Mx,,TN,Mx\frac{∂T_{1.1}}{∂y},\frac{∂T_{1,2}}{∂y},⋯,\frac{∂T_{1,M}}{∂y},\quad\frac{∂T_{1.1}}{∂x},\frac{∂T_{2,1}}{∂x},⋯,\frac{∂T_{N,1}}{∂x},\quad\frac{∂T_{N,1}}{∂y},\frac{∂T_{N,2}}{∂y}⋯\frac{,∂T_{N,M}}{∂y},\quad\frac{∂T_{1,M}}{∂x},\frac{∂T_{2,M}}{∂x},⋯,\frac{∂T_{N,M}}{∂x}




Avendo le derivate nei bordi abbiamo anche delle relazioni fra i valori di temperatura in nodi contigui ai bordi. Per esempio per la prima colonna individuata da j=1 possiamo scrivere:

F1,1=T1,1x=T2,1-T0,12ΔxF_{1,1}=\frac{∂T_{1,1}}{∂x}=\frac{T_{2,1}-T_{0,1}}{2Δx}
ma dobbiamo introdurre un bordo immaginario aggiuntivo. Le temperature nei nodi aggiuntivi si possono però ricavare noti i flussi:T0,1=T2,1-2F1,1ΔxT_{0,1}=T_{2,1}-2\,F_{1,1}Δx
Così si ritorna ad avere un sistema di (N+1)·(M+1) equazioni con altrettante incognite.

PDE PARABOLICHE. Consideriamo ora il flusso di calore non stazionario e unidimensionale attraverso una barra di un materiale con una conducibilità termica k costante. La PDE che descrive la distribuzione di temperatura T(x,t) nello spazio e nel tempo è: Tt=k2Tx2\frac{∂T}{∂t}=k\;\frac{∂^2T}{∂x^2}Qui occorre fissare le condizioni al contorno : {T(0,t)=T1T(L,t)=T2\begin{cases} T(0,t)=T_1 \\ T(L,t)=T_2\end{cases} (valori di temperatura agli estremi della barra dove sono posti i termostati) che non dipendono dal tempo, e le condizioni iniziali T(x,0)=Ti(x)T(x,0) = T_i(x) all'istante t=0 (distribuzione della temperatura nella barra all'istante iniziale)
Dividiamo la barra in N= L/Δx nodi e consideriamo un'evoluzione temporale discreta, ad intervalli Δt.
Ad un certo tempo tlt_l  la derivata seconda  spaziale può essere calcolata applicando il metodo esplicito (riprendendo la first central approximation delle derivate seconde qui chiamata CFDA, Central Finite Difference Approximation ) :

2T(xi,tl)x2=Ti+1l-2Til+Ti-1lΔx2\frac{∂^2T(x_i,t_l)}{∂x^2}=\frac{T^l_{i+1}-2\;T_i^l+T^l_{i-1}}{Δx^2}
mentre la derivata prima temporale tra l'istante tlt_l e tl+1t_{l+1} nel nodo xix_i è : Tt=Til+1-TilΔt\frac{∂T}{∂t}=\frac{T_i^{l+1}-T_i^l}{Δt}Sostituendo nella PDE parabolica:kΔx2(Ti+1l-2Til+Ti-1l)=Til+1-TilΔt\frac{k}{Δx^2}\cdot\left(T^l_{i+1}-2\;T^l_i+T^l_{i-1}\right)=\frac{T^{l+1}_i-T^l_i}{Δt}
Da cui esplicitando la temperatura nel nodo xix_i  all'istante di tempo tl+1t_{l+1} si trova la temperatura in ogni nodo nell'istante di tempo successivo:

Til+1=Til+kΔtΔx2(Ti+1l-2Til+Ti-1l)T_i^{l+1}=T_i^l+\frac{kΔt}{Δx^2}\left(T_{i+1}^l-2\,T_i^l+T_{i-1}^l\right)
Si può dimostrare che il metodo esplicito fornisce soluzioni stabili solo se kΔtΔx212\frac{kΔt}{Δx^2}\le \frac{1}{2}.
Ora  consideriamo la PDE parabolica con la derivata seconda spaziale al tempo tl+1t_{l+1}:kΔx2(Ti+1l+1-2Til+1+Ti-1l+1)=Til+1-TilΔt\frac{k}{Δx^2}\cdot\left(T^{l+1}_{i+1}-2\;T^{l+1}_i+T^{l+1}_{i-1}\right)=\frac{T^{l+1}_i-T^l_i}{Δt}
ponendo λ=kΔtΔx2\lambda=\frac{kΔt}{Δx^2}, si può scrivere:-λTi+1l+1+(1+2λ)Til+1-λTi-1l+1=Til-\lambda\;T^{l+1}_{i+1}+\left(1+2\;\lambda\right)T^{l+1}_i-\lambda\;T^{l+1}_{i-1}=T^l_i
E si possono ottenere, in modo implicito,  le temperature nei nodi al tempo  tl+1t_{l+1} note quelle al tempo tlt_l risolvendo un sistema lineare algebrico triangolare  senza che ci siano problemi di stabilità.
Un miglioramento ulteriore al metodo implicito è quello di Crank-Nicolson che utilizza come derivata spaziale la media delle derivate seconde spaziali al tempo tlt_l e tl+1t_{l+1}: 2Tx2=12Δx2[Ti+1l+1-2Til+1+Ti-1l+1+Ti+1l-2Til+Ti-1l]\frac{∂^2T}{∂x^2}=\frac{1}{2Δx^2}\left[T^{l+1}_{i+1}-2\;T^{l+1}_i+T^{l+1}_{i-1}+T^{l}_{i+1}-2\;T^{l}_i+T^{l}_{i-1}\right]
Da cui la PDE :k2Δx2[Ti+1l+1-2Til+1+Ti-1l+1+Ti+1l-2Til+Ti-1l]=Til+1-TilΔt\frac{k}{2Δx^2}\left[T^{l+1}_{i+1}-2\;T^{l+1}_i+T^{l+1}_{i-1}+T^{l}_{i+1}-2\;T^{l}_i+T^{l}_{i-1}\right]=\frac{T^{l+1}_i-T^l_i}{Δt}
e  l'equazione generatrice del sistema triangolare algebrico:
 -λTi+1l+1+2(1+λ)Til+1-λTi-1l+1=λTi+1l+2(1-λ)Til+λTi-1l-\lambda\;T^{l+1}_{i+1}+2\;\left(1+\;\lambda\right)T^{l+1}_i-\lambda\;T^{l+1}_{i-1}=\lambda\;T^l_{i+1}+2\left(1-\lambda\right)T^l_i+\lambda\;T^l_{i-1}
PDE paraboliche bidimensionali. Ora se si ha  un flusso di calore non stazionario bidimensionale attraverso una lastra  di un materiale con una conducibilità termica k costante, la PDE è :
Tt=k(2Tx2+2Ty2)\frac{∂T}{∂t}=k\left(\frac{∂^2T}{∂x^2}+\frac{∂^2T}{∂y^2}\right)
Anche in questo caso è possibile applicare la CFDA nel nodo xi,yix_i, y_i  :

k[Ti+1,jl-2Ti,jl+Ti-1,jlΔx2+Ti+1,jl-2Ti,jl+Ti-1,jlΔy2]=Ti,jl+1-Ti,jlΔtk\cdot\left[\frac{T^{l}_{i+1,j}-2\;T^{l}_{i,j}+T^{l}_{i-1,j}}{Δx^2}+\frac{T^{l}_{i+1,j}-2\;T^{l}_{i,j}+T^{l}_{i-1,j}}{Δy^2}\right]=\frac{T^{l+1}_{i,j}-T^l_{i,j}}{Δt}
e risolvere l'equazione trovando in modo esplicito la temperatura Ti,jl+1T^{l+1}_{i,j} nel nodo xi,yix_i, y_i  all'istante di tempo tl+1t_{l+1}.   Ma  anche qui il metodo esplicito fornisce soluzioni stabili solo se kΔtΔx2+kΔtΔy212\frac{kΔt}{Δx^2}+\frac{kΔt}{Δy^2}\le \frac{1}{2}.
La formulazione implicita del metodo è:

k[δx2(Tl+1)Δx2+δy2(Tl+1)Δy2]=Til+1-TilΔtk\cdot\left[\frac{\delta^2_x\left(T^{l+1}\right)}{Δx^2}+\frac{\delta^2_y\left(T^{l+1}\right)}{Δy^2}\right]=\frac{T^{l+1}_i-T^l_i}{Δt}

con : {δx2(Tl+1)=Ti-1,jl+1-2Ti,jl+1+Ti+1,jl+1δy2(Tl+1)=Ti,j-1l+1-2Ti,jl+1+Ti,j+1l+1\begin{cases} \delta^2_x\left(T^{l+1}\right)=T^{l+1}_{i-1,j}-2\;T^{l+1}_{i,j}+T^{l+1}_{i+1,j} \\ \delta^2_y\left(T^{l+1}\right)=T^{l+1}_{i,j-1}-2\;T^{l+1}_{i,j}+T^{l+1}_{i,j+1} \end{cases}
E quella di Cranck-Nicolson è:

k[δx2(Tl+1)+δx2(Tl)2Δx2+δy2(Tl+1)+δy2(Tl)2Δy2]=Til+1-TilΔtk\cdot\left[\frac{\delta^2_x\left(T^{l+1}\right)+\delta^2_x\left(T^{l}\right)}{2Δx^2}+\frac{\delta^2_y\left(T^{l+1}\right)+\delta^2_y\left(T^{l}\right)}{2Δy^2}\right]=\frac{T^{l+1}_i-T^l_i}{Δt}con : {δx2(Tl)=Ti-1,jl-2Ti,jl+Ti+1,jlδy2(Tl)=Ti,j-1l-2Ti,jl+Ti,j+1l{δx2(Tl+1)=Ti-1,jl+1-2Ti,jl+1+Ti+1,jl+1δy2(Tl+1)=Ti,j-1l+1-2Ti,jl+1+Ti,j+1l+1\begin{cases} \delta^2_x\left(T^l\right)=T^l_{i-1,j}-2\;T^l_{i,j}+T^l_{i+1,j} \\ \delta^2_y\left(T^l\right)=T^l_{i,j-1}-2\;T^l_{i,j}+T^l_{i,j+1} \end{cases}\qquad\qquad\begin{cases} \delta^2_x\left(T^{l+1}\right)=T^{l+1}_{i-1,j}-2\;T^{l+1}_{i,j}+T^{l+1}_{i+1,j} \\ \delta^2_y\left(T^{l+1}\right)=T^{l+1}_{i,j-1}-2\;T^{l+1}_{i,j}+T^{l+1}_{i,j+1} \end{case
Sia il metodo implicito che quello di Cranck-Nicolson sono stabili ma possono portare alla risoluzione di sistemi algebrici con un numero di equazioni molto elevato. In figura è mostrata la matrice di coefficienti per l'applicazione del metodo di Crank-Nicolson ad una griglia di 4x4 nodi.

Col metodo ADI (Alternating Direct Implicit) si effettuano due mezzi passi temporali uno in ogni direzione. Uno in  direzione y:

k[Ti+1,jl-2Ti,jl+Ti-1,jlΔx2+Ti+1,jl+1/2-2Ti,jl+1/2+Ti-1,jl+1/2Δy2]=Ti,jl+1/2-Ti,jlΔt/2k\cdot\left[\frac{T^{l}_{i+1,j}-2\;T^{l}_{i,j}+T^{l}_{i-1,j}}{Δx^2}+\frac{T^{l+1/2}_{i+1,j}-2\;T^{l+1/2}_{i,j}+T^{l+1/2}_{i-1,j}}{Δy^2}\right]=\frac{T^{l+1/2}_{i,j}-T^l_{i,j}}{Δt/2}e, da quell'istante di tempo tl+1/2,,,,t_{l+1/2} , uno in direzione x :

k[Ti+1,jl+1-2Ti,jl+1+Ti-1,jl+1Δx2+Ti+1,jl+1/2-2Ti,jl+1/2+Ti-1,jl+1/2Δy2]=Ti,jl+1-Ti,jl+1/2Δt/2k\cdot\left[\frac{T^{l+1}_{i+1,j}-2\;T^{l+1}_{i,j}+T^{l+1}_{i-1,j}}{Δx^2}+\frac{T^{l+1/2}_{i+1,j}-2\;T^{l+1/2}_{i,j}+T^{l+1/2}_{i-1,j}}{Δy^2}\right]=\frac{T^{l+1}_{i,j}-T^{l+1/2}_{i,j}}{Δt/2}Posto  λx=Δt2Δx2\lambda_x=\frac{Δt}{2Δx^2 }  e λy=Δt2Δy2\lambda_y=\frac{Δt}{2Δy^2} si ottiene la relazione relativa al primo mezzo passo temporale:

-λxTi+1,jl+1/2+(1+2λx)Ti,jl+1/2-λxTi-1,jl+1/2=-λyTi+1,jl+(1-2λy)Ti,jl+λyTi-1,jl-\lambda_x\;T^{l+1/2}_{i+1,j}+\left(1+2\;\lambda_x\right)T^{l+1/2}_{i,j}-\lambda_x\;T^{l+1/2}_{i-1,j}=-\lambda_y\;T^{l}_{i+1,j}+\left(1-2\;\lambda_y\right)T^{l}_{i,j}+\lambda_y\;T^{l}_{i-1,j}
In modo simile si ricava l’espressione relativa al secondo mezzo passo temporale:

-λyT1,j+1l+1+(1+2λy)Ti,jl+1-λyTi,j-1l+1=λxTi+1,jl+1/2+(1-2λx)Ti,jl+1/2+λxTi-1,jl+1/2-\lambda_y\;T^{l+1}_{1,j+1}+\left(1+2\;\lambda_y\right)T^{l+1}_{i,j}-\lambda_y\;T^{l+1}_{i,j-1}=\lambda_x\;T^{l+1/2}_{i+1,j}+\left(1-2\;\lambda_x\right)T^{l+1/2}_{i,j}+\lambda_x\;T^{l+1/2}_{i-1,j}
Queste espressioni generano un sistema lineare più semplice perché puramente tridiagonale che può essere risolto con l'algoritmo di Thomas.