Corso di matematica propedeutica alla fisica: 14

14. Elementi di calcolo numerico

Obiettivo didattico: introdurre i principali metodi numerici per l’approssimazione, la soluzione di equazioni e sistemi lineari/non lineari, la integrazione numerica di ODE/PDE e mostrare applicazioni fisiche rilevanti (fluidodinamica, termodinamica, dinamica molecolare, calcolo orbitale).

1. Concetti fondamentali: errore, ordine, stabilità

Errore assoluto: se $x$ è il valore esatto e $\tilde x$ la stima, l’errore assoluto è $e_a = |x-\tilde x|$ .
Errore relativo: $e_r = \dfrac{|x-\tilde x|}{|x|}$ (se $x\neq0$ ).
Ordine di un metodo: un metodo ha errore locale proporzionale a $h^{p+1}$ e errore globale proporzionale a $h^p$ ; si dice che è di ordine $p$ .
Stabilità numerica: riguarda la propagazione degli errori (es. per integrazione temporale). Un metodo può essere condizionatamente stabile (richiede un passo $h$ abbastanza piccolo) o incondizionatamente stabile.
Condizionamento del problema: anche con algoritmo perfetto, problemi mal condizionati amplificano piccoli errori nei dati.

2. Metodi di approssimazione e interpolazione

2.1 Interpolazione polinomiale (Lagrange)

Dato un insieme di punti $(x_i,y_i)$ , esiste un polinomio unico di grado $\le n$ che interpola. La forma di Lagrange:

P_n(x) = \sum_{i=0}^n y_i \ell_i(x), \qquad \ell_i(x)=\prod_{\substack{j=0\\ j\ne i}}^n \frac{x-x_j}{x_i-x_j}.

Errore di interpolazione (per $f\in C^{n+1}$ ):

f(x)-P_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n (x-x_i),\quad \xi\in\text{conv}(x_0,\dots,x_n).

2.2 Interpolazione spline

Spline cubiche a tratti forniscono approssimazioni più stabili rispetto a polinomi di alto grado; sono due volte differenziabili e usano condizioni di continuità ai nodi.

3. Risoluzione numerica di equazioni non lineari

Obiettivo: trovare radici di $f(x)=0$ .

3.1 Metodo della bisezione (robusto, ordine 1)

Dati $a<b$ con $f(a)f(b)<0$ , iterare:

c=\frac{a+b}{2},\quad \text{se } f(a)f(c)<0 \text{ allora } b=c \text{ else } a=c.

Dopo $n$ iterazioni l’errore è $|b_n-a_n|/2 \le (b-a)/2^{n+1}$ .

Esempio numerico
Risolvere $f(x)=\cos x - x =0$ su $[0,1]$ :

$f(0)=1$ , $f(1)=\cos1-1\approx -0.4597$ → cambia segno.
Primo passo: $c=0.5$ , $f(0.5)\approx 0.3776>0$ → nuovo intervallo $[0.5,1]$ .
Dopo 10 iterazioni l’errore $< (1-0)/2^{11}\approx 0.000488$ .

3.2 Newton–Raphson (veloce, ordine 2)

Iterazione:

x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.

Convergenza quadratica vicino alla radice se $f'(x^*)\neq0$ .

Esempio numerico
Usiamo $f(x)=\cos x - x$ , $f'(x)=-\sin x -1$ . Scegli $x_0=0.7$ :

$x_1 = 0.7 - (\cos0.7-0.7)/(-\sin0.7-1)\approx 0.7391$ .
$x_2\approx 0.739085$ → converge rapidamente alla $x^*\approx 0.7390851332$ .

Nota: richiede derivata e buona inizializzazione; può divergere se $f'$ è piccolo.

3.3 Metodo delle secanti (senza derivata, ordine $\approx1.618$ )

x_{n+1}=x_n - f(x_n)\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}.

4. Sistemi lineari $A\mathbf{x}=\mathbf{b}$

4.1 Eliminazione di Gauss (diretto)

Riduzione a matrice triangolare e sostituzione all’indietro. Complessità $\mathcal{O}(n^3)$ .

Esempio numerico (3×3):

A=\begin{pmatrix}2&1&-1\\-3&-1&2\\-2&1&2\end{pmatrix},\qquad b=\begin{pmatrix}8\\-11\\-3\end{pmatrix}.

Risultato noto: $x=(2,3,-1)^T$ . (Passaggi: pivoting, sottrazioni riga, back-substitution — omessi per brevità ma facilmente riportabili in esercitazione.)

4.2 Fattorizzazione LU

Si scrive $A=LU$ con $L$ unit lower, $U$ upper; poi si risolvono $L\mathbf{y}=\mathbf{b}$ e $U\mathbf{x}=\mathbf{y}$ .

4.3 Metodi iterativi (Jacobi, Gauss–Seidel, SOR)

Per sistemi sparsi o grandi; condizioni di convergenza legate al raggio spettrale della matrice iterativa.

Jacobi: $\mathbf{x}^{(k+1)} = D^{-1}(\mathbf{b} - (L+U)\mathbf{x}^{(k)})$ .
Gauss–Seidel: usa componenti aggiornate immediatamente, solitamente più rapido.
Convergenza: garantita se $A$ è strettamente diagonalmente dominante o simmetrica definita positiva.

5. Integrazione numerica di ODE ordinarie

Problema iniziale: $\dot y = f(t,y),\quad y(t_0)=y_0$ .

5.1 Metodo di Eulero esplicito (ordine 1)

y_{n+1}=y_n + h f(t_n,y_n).

Semplice ma condizionatamente stabile; passo $h$ deve essere piccolo per accuratezza.

Esempio (oscillatore armonico non smorzato):

\ddot x + \omega^2 x = 0 \quad\Rightarrow\quad \begin{cases} \dot x = v,\\ \dot v = -\omega^2 x. \end{cases}

Con $\omega=1$ , $x(0)=1$ , $v(0)=0$ , scelta $h=0.1$ :

$x_1 = x_0 + h v_0 = 1$ .
$v_1 = v_0 + h(-x_0) = -0.1$ .
Questo Eulero esplicito presenta dissipazione/numerical damping non fisico per $h$ non molto piccolo.

5.2 Metodo di Runge–Kutta quarto ordine (RK4, ordine 4)

Passo di aggiornamento:

\begin{aligned} k_1 &= f(t_n,y_n),\\ k_2 &= f\big(t_n+\tfrac h2,\; y_n+\tfrac h2 k_1\big),\\ k_3 &= f\big(t_n+\tfrac h2,\; y_n+\tfrac h2 k_2\big),\\ k_4 &= f(t_n+h,\; y_n + h k_3\big),\\ y_{n+1} &= y_n + \frac{h}{6}(k_1+2k_2+2k_3+k_4). \end{aligned}

Molto usato per la sua efficienza e accuratezza.

Esempio numerico (lo stesso oscillatore): con $h=0.1$ RK4 fornisce approssimazioni molto vicine alla soluzione esatta $x(t)=\cos t$ , con errore dell’ordine $h^4$ .

5.3 Metodi impliciti (Backward Euler, trapezio implicito)

Metodi impliciti sono più stabili per equazioni rigide (stiff): es. Backward Euler

y_{n+1}=y_n + h f(t_{n+1},y_{n+1}),

richiede risoluzione non lineare ad ogni passo (Newton).

6. Metodi per PDE (introduzione)

Consideriamo la equazione del calore (parabolica) in 1D:

\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2},\qquad u(x,0)=u_0(x).

Discretizziamo spazio $x_j=j\Delta x$ , tempo $t^n=n\Delta t$ , $u_j^n\approx u(x_j,t^n)$ .

6.1 Schema FTCS (Forward in time, Centered in space) — esplicito

u_j^{n+1} = u_j^n + r\big(u_{j+1}^n - 2u_j^n + u_{j-1}^n\big),\qquad r=\frac{\alpha\Delta t}{(\Delta x)^2}.

Condizione di stabilità (CFL): per stabilità si richiede $r \le \tfrac{1}{2}$ (in 1D).

6.2 Schema implicito (Backward Euler) — stabile ma richiede soluzione lineare

u_j^{n+1} - r\big(u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}\big) = u_j^n.

Produrrà un sistema tridiagonale risolvibile efficientemente (Thomas algorithm).

6.3 Crank–Nicolson (ordine 2 in tempo e spazio)

Media tra esplicito e implicito:

u_j^{n+1} = u_j^n + \frac{r}{2}\big( u_{j+1}^n - 2u_j^n + u_{j-1}^n + u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1} \big).

Schema non condizionatamente stabile ma numericamente ottimo (A-stable e di ordine 2).

7. Errori, convergenza e analisi di stabilità

Errore locale di truncamento: misura quanto si sbaglia passando dall’equazione continua allo schema discreto in un singolo passo.
Convergenza: uno schema è convergente se l’errore globale tende a zero quando $h\to0$ (o $\Delta x,\Delta t\to0$ ).
Analisi von Neumann: uso di Fourier modes per stabilità lineare di schemi numerici; si ricerca che il fattore di amplificazione $G(k)$ soddisfi $|G(k)|\le1$ .

8. Metodi speciali usati in applicazioni fisiche

8.1 Dinamica molecolare — integratori (Verlet)

Equazioni di Newton $\ddot{\mathbf r} = \mathbf F(\mathbf r)/m$ . Integratore di Verlet (symplectic) molto usato:

\mathbf r(t+\Delta t)=2\mathbf r(t)-\mathbf r(t-\Delta t) + \Delta t^2 \frac{\mathbf F(t)}{m}.

Conserva meglio l’energia su lungo periodo rispetto a metodi non symplectic (importante in simulazioni di sistemi Hamiltoniani).

8.2 Integrazione di orbite — metodi symplectici

Per problemi gravitazionali a lungo termine (es. simulazione del Sistema Solare), si usano integratori symplectici (Leapfrog, Stormer–Verlet, Yoshida schemes) che preservano la struttura geometrica e aumentano stabilità energetica su tempi lunghi.

8.3 Fluidodinamica numerica (CFD)

Discretizzazione: finite difference, finite volume o finite element.
Equazioni: Navier–Stokes (non lineari), richiedono solver robusti (implicit/explicit), trattamenti di pressione (projection methods, SIMPLE).
CFL condition (per schemi espliciti): $\text{CFL} = \dfrac{u\,\Delta t}{\Delta x} \le C_{\max}$ (dipende dallo schema) — garantisce stabilità rispetto all’advection.

9. Esempi numerici svolti

9.1 Newton per radice (già mostrato): convergenza quadratica, pochi iter.

9.2 Sistema lineare semplice con Gauss

Risolvere

\begin{pmatrix} 3 & 2 \\ 1 & 2 \end{pmatrix}\begin{pmatrix} x_1\\x_2\end{pmatrix}=\begin{pmatrix} 5\\5\end{pmatrix}.

Eliminazione:

$R_2 \leftarrow R_2 - \tfrac{1}{3}R_1$ → nuova matrice, poi back-substitution produce $x_1=1$ , $x_2=2$ .

9.3 Integrazione ODE con RK4 (dettagli)

Problema: $\dot y = -2y,\; y(0)=1$ . Soluzione esatta $y(t)=e^{-2t}$ .
Con $h=0.1$ passo singolo RK4:

$k_1=-2y_0=-2$ .
$k_2=-2(y_0 + 0.05 k_1) = -2(1 - 0.1)=-1.8$ .
$k_3=-2(1 - 0.09)=-1.82$ (valori approssimati),
$k_4=-2(1 + 0.1(-1.82))\approx -1.636$ .
$y_1 = 1 + \tfrac{0.1}{6}( -2 - 3.6 - 3.64 - 1.636)\approx 0.8187$ .
Soluzione esatta a $t=0.1$ : $e^{-0.2}\approx 0.81873$ . RK4 mostra ottima accuratezza.

9.4 Heat equation FTCS (esempio numerico)

Take domain $x\in[0,1]$ , $\alpha=1$ , $\Delta x=0.1$ , quindi $r=\Delta t/(\Delta x)^2$ . Per stabilità scegliere $\Delta t \le 0.005$ se vogliamo $r\le 1/2$ . Calcoli iterativi aggiornano $u_j^n$ con formula FTCS; visualizzare profilo temperatura nel tempo dimostra diffusione e smoothing.

10. Linee guida pratiche e suggerimenti di implementazione

Verifica con problemi noti: confronta risultati con soluzioni analitiche quando disponibili.
Convergenza e mesh refinement: esegui studi di convergenza riducendo $h,\Delta x,\Delta t$ e misurando il comportamento dell’errore.
Scelta del metodo: per problemi stiff usare metodi impliciti; per Hamiltoniani usare metodi symplectici.
Prestazioni: sfrutta strutture sparse, fattorizzazioni e solver iterative per grandi sistemi.
Controllo dell’errore adattativo: integra con step adattativi (es. Runge–Kutta–Fehlberg) per mantenere tolleranza desiderata.
Stabilità numerica: controlla CFL e condizioni di stabilità specifiche per lo schema e il problema.
Documentazione e test: mantieni test automatizzati e validazione fisica dei risultati.

11. Bibliografia e risorse consigliate (italiano e internazionale)

A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, (trad. italiana o testi equivalenti) — testo completo su metodi numerici.
R.L. Burden, J.D. Faires, Numerical Analysis — classico introduttivo.
J.W. Thomas, Numerical Partial Differential Equations — finite difference/volume methods.
E. Hairer, S.P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I & II — teoria e metodi (stiff, symplectic).
D. T. Greenwood, Principles of Dynamics — per integrazione di moti e metodi numerici in meccanica (capitoli su integratori).
Le guide/lecture notes di università (appunti su Runge–Kutta, Crank–Nicolson, Thomas algorithm) per esempi implementativi.