Formulario di EDO e Analisi Numerica

Questo formulario raccoglie gli strumenti principali di equazioni differenziali ordinarie e analisi numerica per ingegneria. Le due parti sono strettamente collegate: le EDO forniscono modelli di evoluzione, mentre l’analisi numerica fornisce metodi affidabili per approssimare soluzioni quando la forma chiusa non è disponibile o non è conveniente.

Le formule vanno lette con le loro ipotesi. Un metodo numerico non è solo una procedura: ha errore, stabilità, costo computazionale e condizioni di applicabilità. Allo stesso modo, una EDO non è solo un’equazione da integrare: è un modello con dati iniziali, parametri, intervallo di validità e possibili fenomeni qualitativi.

1. Linguaggio delle equazioni differenziali ordinarie

Equazione differenziale ordinaria

$$F(t,y,y',\dots,y^{(m)})=0$$

Una EDO lega una funzione incognita $y(t)$ e le sue derivate rispetto a una sola variabile indipendente, di solito il tempo $t$. L’ordine dell’equazione è l’ordine massimo di derivata che compare. In ingegneria le EDO modellano evoluzione di sistemi meccanici, circuiti, reazioni, serbatoi, popolazioni, controlli e dinamiche termiche concentrate.

Forma normale del primo ordine

$$y'=f(t,y)$$

La forma normale esplicita la derivata come funzione del tempo e dello stato. È la forma più adatta per teoria di esistenza, simulazione numerica e interpretazione geometrica tramite campi di direzioni. Ogni punto $(t,y)$ indica una pendenza locale della curva soluzione.

Problema di Cauchy

$$\begin{cases} y'=f(t,y),\\ y(t_0)=y_0. \end{cases}$$

L’equazione differenziale descrive una famiglia di traiettorie; la condizione iniziale sceglie quella che passa per $(t_0,y_0)$. Nel modello fisico, $y_0$ è lo stato misurato o imposto all’istante iniziale. La domanda matematica è se esista una soluzione, se sia unica e per quanto tempo rimanga definita.

Soluzione esplicita e implicita

$$y=\varphi(t),\qquad G(t,y)=C$$

Una soluzione esplicita dà direttamente $y$ come funzione di $t$. Una soluzione implicita descrive invece la traiettoria come curva di livello di una relazione $G(t,y)=C$. Nei problemi reali la forma implicita può essere più naturale o l’unica ottenibile in modo elementare.

Equazione autonoma

$$y'=f(y)$$

Una EDO autonoma non dipende esplicitamente dal tempo. Il comportamento è determinato dallo stato corrente, non dall’istante assoluto. Questo rende possibile l’analisi qualitativa tramite retta di fase, equilibri e stabilità.

Equilibrio

$$f(y_e)=0$$

Un equilibrio è uno stato in cui la derivata si annulla. Se il sistema parte esattamente da $y_e$, resta lì. In applicazioni gli equilibri rappresentano regimi stazionari, posizioni di riposo, concentrazioni di equilibrio o punti operativi.

Ordine e riduzione a sistema

$$y^{(m)}=f(t,y,y',\dots,y^{(m-1)})$$

Una EDO di ordine $m$ può essere trasformata in un sistema del primo ordine introducendo variabili di stato. Questa riduzione è essenziale sia teoricamente sia numericamente, perché molti metodi standard sono formulati per sistemi del primo ordine.

Riduzione di ordine

$$x_1=y,\qquad x_2=y',\qquad \dots,\qquad x_m=y^{(m-1)}$$

Con queste variabili si ottiene $x_1'=x_2$, $x_2'=x_3$ e così via, fino all’equazione per $x_m'$. In questo modo una singola equazione di ordine alto diventa un sistema vettoriale. Il prezzo è aumentare la dimensione dello stato; il vantaggio è uniformare il trattamento.

Linearità

$$a_m(t)y^{(m)}+\cdots+a_1(t)y'+a_0(t)y=g(t)$$

Una EDO è lineare se l’incognita e le sue derivate compaiono solo alla prima potenza e non sono moltiplicate tra loro. La linearità consente sovrapposizione, uso di basi di soluzioni e metodi algebrici. Se $g(t)=0$, l’equazione è omogenea; se $g(t)\ne0$, è non omogenea.

Principio di sovrapposizione

$$L[y_1]=0,\quad L[y_2]=0 \Longrightarrow L[\alpha y_1+\beta y_2]=0$$

Per equazioni lineari omogenee, combinazioni lineari di soluzioni sono ancora soluzioni. Questo principio non vale per equazioni non lineari. È alla base della costruzione della soluzione generale tramite soluzioni fondamentali indipendenti.

2. EDO del primo ordine risolvibili in forma chiusa

Equazione a variabili separabili

$$y'=g(t)h(y)$$

L’equazione è separabile quando la dipendenza da $t$ e quella da $y$ si possono separare come prodotto. Se $h(y)\ne0$, si porta tutto ciò che contiene $y$ da una parte e tutto ciò che contiene $t$ dall’altra. Bisogna però non perdere le soluzioni costanti date dagli zeri di $h$.

Formula di separazione

$$\frac{\mathrm{d}y}{h(y)}=g(t)\,\mathrm{d}t \qquad\Longrightarrow\qquad \int\frac{1}{h(y)}\,\mathrm{d}y=\int g(t)\,\mathrm{d}t+C$$

Questa scrittura è una notazione compatta per integrare lungo una soluzione. Non autorizza manipolazioni prive di ipotesi: la divisione per $h(y)$ è lecita solo dove $h(y)$ non si annulla. Dopo l’integrazione può rimanere una relazione implicita tra $t$ e $y$.

Crescita esponenziale

$$y'=ky \qquad\Longrightarrow\qquad y(t)=Ce^{kt}$$

La variazione è proporzionale alla quantità presente. Se $k>0$ si ha crescita, se $k<0$ decadimento. La costante $C$ si determina con il dato iniziale. Questo modello compare in radioattività, popolazioni elementari, carica/scarica ideale e sistemi lineari del primo ordine.

Equazione logistica

$$y'=ry\left(1-\frac{y}{K}\right)$$

$r$ è il tasso di crescita iniziale e $K$ è la capacità portante. Per $0<y<K$ la soluzione cresce; per $y>K$ decresce; $y=0$ e $y=K$ sono equilibri. Il termine $1-y/K$ introduce saturazione e rende il modello non lineare.

Soluzione logistica

$$y(t)=\frac{K}{1+A e^{-rt}}$$

La costante $A$ dipende dal dato iniziale. La soluzione tende a $K$ per $t\to+\infty$ se $r>0$ e il dato iniziale è positivo. La crescita inizialmente è quasi esponenziale, poi rallenta quando la soluzione si avvicina alla capacità portante.

Equazione lineare del primo ordine

$$y'+a(t)y=b(t)$$

L’equazione è lineare perché $y$ e $y'$ compaiono linearmente. Il coefficiente $a(t)$ descrive decadimento, crescita o accoppiamento proporzionale allo stato; $b(t)$ è una forzante esterna. È uno dei modelli più frequenti in circuiti, bilanci e sistemi del primo ordine.

Fattore integrante

$$\mu(t)=e^{\int a(t)\,\mathrm{d}t}$$

Il fattore integrante trasforma il lato sinistro in una derivata di prodotto. La costante della primitiva dentro l’esponenziale non è importante, perché moltiplica $\mu$ per una costante non nulla che si semplifica nella soluzione finale.

Soluzione dell’equazione lineare

$$(\mu y)'=\mu b \qquad\Longrightarrow\qquad y(t)=\frac{1}{\mu(t)}\left(\int \mu(t)b(t)\,\mathrm{d}t+C\right)$$

La formula mostra la struttura della soluzione: una parte libera, legata alla costante $C$, e una parte forzata, legata a $b(t)$. Nei sistemi fisici la parte libera è il transitorio, mentre la parte forzata è la risposta all’ingresso.

Equazione esatta

$$M(t,y)\,\mathrm{d}t+N(t,y)\,\mathrm{d}y=0$$

L’equazione è esatta se esiste una funzione potenziale $\Phi(t,y)$ tale che $\Phi_t=M$ e $\Phi_y=N$. In quel caso le soluzioni sono curve di livello $\Phi(t,y)=C$. Il problema diventa trovare il potenziale.

Condizione di esattezza

$$\frac{\partial M}{\partial y}=\frac{\partial N}{\partial t}$$

In un dominio semplicemente connesso e con regolarità sufficiente, questa condizione garantisce l’esattezza. Se non vale, talvolta si può moltiplicare l’equazione per un fattore integrante che la renda esatta.

Equazione di Bernoulli

$$y'+a(t)y=b(t)y^\alpha$$

Per $\alpha=0$ o $\alpha=1$ l’equazione è già lineare. Negli altri casi si usa il cambio $z=y^{1-\alpha}$, assumendo che la potenza abbia senso sull’intervallo considerato. Bernoulli è un esempio tipico di non linearità riconducibile a linearità.

Cambio di Bernoulli

$$z=y^{1-\alpha} \qquad\Longrightarrow\qquad z'+(1-\alpha)a(t)z=(1-\alpha)b(t)$$

Il cambio trasforma l’equazione in una lineare del primo ordine per $z$. Dopo aver risolto per $z$, si torna a $y$. Bisogna controllare eventuali soluzioni perse, specialmente $y=0$ quando compare una divisione o una potenza problematica.

Equazione omogenea del primo ordine

$$y'=F\left(\frac{y}{t}\right)$$

L’equazione è omogenea quando il campo di pendenze dipende solo dal rapporto $y/t$. Il cambio naturale è $v=y/t$, cioè $y=vt$. Questo sfrutta l’invarianza per dilatazione delle rette passanti per l’origine.

Cambio per equazioni omogenee

$$y=vt,\qquad y'=v+t v'$$

Sostituendo si ottiene un’equazione separabile per $v(t)$, almeno nei casi regolari. Il cambio richiede attenzione vicino a $t=0$, dove il rapporto $y/t$ non è definito.

3. Esistenza, unicità e comportamento qualitativo

Teorema di esistenza locale

$$f\ \text{continua vicino a }(t_0,y_0) \Longrightarrow \exists\ \text{almeno una soluzione locale}$$

La continuità del campo garantisce l’esistenza di una curva soluzione almeno per un piccolo intervallo attorno a $t_0$. Non garantisce necessariamente unicità. Questo distingue l’esistenza del modello dalla predicibilità univoca.

Condizione di Lipschitz in $y$

$$|f(t,y_1)-f(t,y_2)|\le L|y_1-y_2|$$

La condizione di Lipschitz controlla quanto rapidamente il campo cambia rispetto allo stato. Se vale localmente, insieme alla continuità, garantisce unicità della soluzione del problema di Cauchy. È più forte della semplice continuità, ma spesso segue dalla continuità di $f_y$.

Teorema di Picard-Lindelöf

$$f\ \text{continua},\qquad f\ \text{localmente Lipschitz in }y \Longrightarrow \text{esistenza e unicità locale}$$

Il teorema assicura che il problema di Cauchy è ben posto localmente. Soluzioni con lo stesso dato iniziale non possono biforcarsi. Questa proprietà è fondamentale nella simulazione: senza unicità, il modello non determina un’evoluzione univoca.

Iterazioni di Picard

$$y_{k+1}(t)=y_0+\int_{t_0}^{t}f(s,y_k(s))\,\mathrm{d}s$$

Le iterazioni di Picard costruiscono la soluzione come punto fisso dell’operatore integrale associato alla EDO. La formula è più teorica che computazionale, ma chiarisce il legame tra EDO, integrali e contrazioni.

Prolungamento della soluzione

$$\text{soluzione massimale definita su }(T_-,T_+)$$

Una soluzione locale può essere prolungata finché resta nel dominio del campo e non diverge. Gli estremi dell’intervallo massimale dipendono sia dalla geometria del dominio sia dalla crescita della soluzione. Una soluzione può esplodere in tempo finito.

Campo di direzioni

$$y'=f(t,y)$$

Il campo di direzioni assegna a ogni punto una pendenza. Disegnare questi segmenti aiuta a intuire le traiettorie senza risolvere l’equazione. È particolarmente utile per equazioni non risolvibili in forma chiusa.

Retta di fase per autonome

$$y'=f(y)$$

Sulla retta di fase si segnano gli zeri di $f$ e il segno di $f$ tra gli zeri. Dove $f>0$ le soluzioni crescono; dove $f<0$ decrescono. Questo permette di stabilire stabilità degli equilibri senza integrare esplicitamente.

Stabilità di un equilibrio scalare

$$f(y_e)=0,\qquad f'(y_e)<0 \Longrightarrow y_e\ \text{asintoticamente stabile}$$

Se la derivata del campo in equilibrio è negativa, piccole perturbazioni vengono riportate verso l’equilibrio. Se $f'(y_e)>0$, l’equilibrio è instabile. Se $f'(y_e)=0$, il criterio lineare è inconcludente e bisogna studiare termini di ordine superiore o il segno di $f$.

Linearizzazione scalare

$$y'=f(y),\qquad \eta=y-y_e,\qquad \eta'\approx f'(y_e)\eta$$

Vicino a un equilibrio, la dinamica è approssimata da un’equazione lineare. Questa approssimazione descrive il comportamento locale quando il termine lineare non si annulla. È una delle idee più importanti nel passaggio da modelli non lineari a modelli analizzabili.

4. EDO lineari di ordine superiore

Equazione lineare omogenea di ordine $m$

$$a_m(t)y^{(m)}+\cdots+a_1(t)y'+a_0(t)y=0$$

La soluzione generale, sotto ipotesi regolari e con $a_m(t)\ne0$, è combinazione lineare di $m$ soluzioni indipendenti. La dimensione dello spazio delle soluzioni coincide con l’ordine dell’equazione.

Wronskiano

$$W(y_1,\dots,y_m)(t)= \det \begin{pmatrix} y_1 & y_2 & \cdots & y_m\\ y_1' & y_2' & \cdots & y_m'\\ \vdots & \vdots & \ddots & \vdots\\ y_1^{(m-1)} & y_2^{(m-1)} & \cdots & y_m^{(m-1)} \end{pmatrix}$$

Il Wronskiano misura l’indipendenza delle soluzioni. Se è non nullo in un punto, le funzioni formano un sistema fondamentale di soluzioni. Per EDO lineari omogenee, il Wronskiano o è sempre nullo o non si annulla mai sull’intervallo di regolarità.

Soluzione generale omogenea

$$y_h(t)=C_1y_1(t)+\cdots+C_m y_m(t)$$

Le costanti $C_i$ si determinano imponendo condizioni iniziali o al contorno. La base di soluzioni $y_1,\dots,y_m$ deve essere indipendente; altrimenti non si copre tutto lo spazio delle soluzioni.

Equazione non omogenea

$$L[y]=g(t)$$

La soluzione generale si ottiene sommando una soluzione particolare $y_p$ a tutte le soluzioni dell’omogenea associata. La parte omogenea rappresenta i modi naturali del sistema; la particolare rappresenta la risposta alla forzante.

Struttura della soluzione non omogenea

$$y(t)=y_h(t)+y_p(t)$$

Questa decomposizione vale solo per problemi lineari. La soluzione particolare non è unica: aggiungendo una soluzione omogenea si ottiene un’altra particolare. La somma finale, con le costanti corrette, soddisfa dati iniziali o condizioni al contorno.

Coefficienti costanti

$$a_m y^{(m)}+\cdots+a_1y'+a_0y=0$$

Per coefficienti costanti si cercano soluzioni esponenziali $y=e^{\lambda t}$. Questo funziona perché la derivata di un esponenziale è ancora lo stesso esponenziale moltiplicato per una costante.

Equazione caratteristica

$$a_m\lambda^m+\cdots+a_1\lambda+a_0=0$$

Le radici dell’equazione caratteristica determinano i modi elementari del sistema. Radici reali danno esponenziali reali; radici complesse coniugate danno oscillazioni eventualmente smorzate; radici multiple introducono fattori polinomiali.

Radice reale semplice

$$\lambda\in\mathbb{R} \qquad\Longrightarrow\qquad y(t)=e^{\lambda t}$$

Una radice reale semplice produce un modo esponenziale. Se $\lambda<0$, il modo decade; se $\lambda>0$, cresce; se $\lambda=0$, è costante. In stabilità lineare, il segno delle parti reali delle radici è decisivo.

Radice reale multipla

$$\lambda\ \text{molteplicità }r \qquad\Longrightarrow\qquad e^{\lambda t},\ t e^{\lambda t},\ \dots,\ t^{r-1}e^{\lambda t}$$

La molteplicità richiede più soluzioni indipendenti associate alla stessa radice. I fattori polinomiali compaiono perché una sola esponenziale non basta a generare la dimensione necessaria dello spazio delle soluzioni.

Radici complesse coniugate

$$\lambda=\alpha\pm i\beta \qquad\Longrightarrow\qquad e^{\alpha t}\cos\beta t,\quad e^{\alpha t}\sin\beta t$$

La parte reale $\alpha$ controlla crescita o decadimento dell’ampiezza; la parte immaginaria $\beta$ controlla l’oscillazione. Questo è il linguaggio naturale di vibrazioni, circuiti RLC e sistemi del secondo ordine.

Oscillatore armonico libero

$$y''+\omega_0^2 y=0$$

La soluzione è una combinazione di seno e coseno con pulsazione naturale $\omega_0$. Il sistema oscilla senza smorzamento e conserva energia. È il modello base per vibrazioni ideali.

Oscillatore smorzato

$$y''+2\zeta\omega_0 y'+\omega_0^2y=0$$

$\zeta$ è il rapporto di smorzamento. Se $0<\zeta<1$ il sistema è sottosmorzato e oscilla con ampiezza decrescente; se $\zeta=1$ è criticamente smorzato; se $\zeta>1$ è sovrasmorzato. Questa classificazione è centrale in meccanica e controlli.

Variazione delle costanti

$$y_p=\sum_{i=1}^{m}u_i(t)y_i(t)$$

Il metodo cerca una particolare non omogenea lasciando variare le costanti della soluzione omogenea. Le funzioni $u_i(t)$ si determinano imponendo condizioni ausiliarie che rendono il sistema risolvibile. È generale, ma può essere laborioso.

Metodo dei coefficienti indeterminati

$$L[y]=g(t)$$

Quando i coefficienti sono costanti e $g(t)$ è combinazione di polinomi, esponenziali, seni o coseni, si prova una particolare della stessa famiglia. Se la prova collide con una soluzione omogenea, si moltiplica per una potenza di $t$ sufficiente a ristabilire indipendenza.

5. Sistemi di EDO

Sistema del primo ordine

$$x'=F(t,x),\qquad x(t)\in\mathbb{R}^n$$

Un sistema descrive l’evoluzione simultanea di più variabili di stato. È il formato naturale per modelli fisici con più gradi di libertà, circuiti, reazioni accoppiate e sistemi di controllo. Le traiettorie vivono nello spazio degli stati.

Sistema lineare

$$x'=A(t)x+b(t)$$

La matrice $A(t)$ descrive l’accoppiamento tra variabili di stato; il vettore $b(t)$ rappresenta ingressi o forzanti. Se $b=0$, il sistema è omogeneo. La teoria lineare sfrutta matrici fondamentali, autovalori e sovrapposizione.

Sistema lineare autonomo

$$x'=Ax$$

Con matrice costante, la soluzione è governata dall’esponenziale di matrice. Gli autovalori di $A$ controllano stabilità e oscillazioni. Gli autovettori indicano direzioni caratteristiche nello spazio degli stati.

Esponenziale di matrice

$$e^{At}=\sum_{k=0}^{+\infty}\frac{A^k t^k}{k!}$$

L’esponenziale di matrice generalizza l’esponenziale scalare. La serie converge per ogni matrice quadrata. Non tutte le regole scalari valgono automaticamente, perché le matrici possono non commutare.

Soluzione del sistema lineare autonomo

$$x(t)=e^{A(t-t_0)}x_0$$

La matrice $e^{A(t-t_0)}$ propaga lo stato iniziale al tempo $t$. È detta matrice di transizione. In simulazione e controllo rappresenta l’evoluzione libera del sistema.

Diagonalizzazione

$$A=PDP^{-1} \qquad\Longrightarrow\qquad e^{At}=P e^{Dt}P^{-1}$$

Se $A$ è diagonalizzabile, calcolare l’esponenziale si riduce a esponenziare gli autovalori sulla diagonale. Questo chiarisce la dinamica modo per modo. Se $A$ non è diagonalizzabile, compaiono blocchi di Jordan e fattori polinomiali.

Stabilità lineare

$$\operatorname{Re}\lambda_i<0\quad\forall i \Longrightarrow x=0\ \text{asintoticamente stabile}$$

Per un sistema lineare autonomo, parti reali negative degli autovalori implicano decadimento dei modi. Se almeno un autovalore ha parte reale positiva, l’origine è instabile. Autovalori sull’asse immaginario richiedono analisi più fine.

Sistema non lineare autonomo

$$x'=F(x)$$

La dinamica dipende solo dallo stato. Gli equilibri sono punti in cui il campo si annulla. Per capire il comportamento locale si linearizza attorno agli equilibri, ma il comportamento globale può essere molto più ricco.

Linearizzazione

$$F(x_e)=0,\qquad \eta=x-x_e,\qquad \eta'=J_F(x_e)\eta+o(\lVert\eta\rVert)$$

La Jacobiana nel punto di equilibrio è il modello lineare locale del sistema. Se gli autovalori hanno parti reali non nulle, spesso la linearizzazione determina la stabilità locale. Se compaiono parti reali nulle, servono strumenti non lineari.

Piano delle fasi

$$x'=F(x,y),\qquad y'=G(x,y)$$

Nel piano delle fasi si studiano traiettorie nello spazio degli stati, non grafici temporali separati. Nulcline, equilibri, direzioni del campo e curve integrali permettono una lettura qualitativa del sistema.

Nulcline

$$F(x,y)=0,\qquad G(x,y)=0$$

Le nulcline sono curve dove una componente della velocità si annulla. Le loro intersezioni sono equilibri. Il segno delle componenti nei vari settori aiuta a disegnare il verso del campo e prevedere le traiettorie.

6. Trasformata di Laplace per EDO

Trasformata di Laplace

$$Y(s)=\mathcal{L}\{y(t)\}=\int_0^{+\infty}y(t)e^{-st}\,\mathrm{d}t$$

La trasformata di Laplace converte problemi causali nel dominio del tempo in problemi algebrici nel dominio $s$. È particolarmente adatta a EDO lineari con coefficienti costanti e condizioni iniziali assegnate in $t=0$.

Derivata prima

$$\mathcal{L}\{y'\}=sY(s)-y(0^+)$$

La derivata diventa moltiplicazione per $s$, ma compare il valore iniziale. Per questo Laplace incorpora naturalmente le condizioni iniziali dentro l’equazione trasformata.

Derivata seconda

$$\mathcal{L}\{y''\}=s^2Y(s)-s y(0^+)-y'(0^+)$$

Ogni derivata successiva introduce nuovi dati iniziali. Per EDO del secondo ordine, posizione e velocità iniziale entrano direttamente nella formula. Questo è essenziale per vibrazioni, circuiti e sistemi meccanici.

Convoluzione causale

$$(f*g)(t)=\int_0^t f(\tau)g(t-\tau)\,\mathrm{d}\tau$$

La convoluzione causale combina la storia passata dell’ingresso con la risposta impulsiva del sistema. Un sistema lineare tempo-invariante ha uscita pari alla convoluzione tra ingresso e risposta all’impulso.

Laplace della convoluzione

$$\mathcal{L}\{f*g\}=F(s)G(s)$$

La convoluzione nel tempo diventa prodotto nel dominio trasformato. Questo trasforma equazioni differenziali e integrali in manipolazioni algebriche. È la base della funzione di trasferimento.

Funzione di trasferimento

$$H(s)=\frac{Y(s)}{U(s)}$$

La funzione di trasferimento è il rapporto tra uscita e ingresso con condizioni iniziali nulle. Poli e zeri determinano stabilità, transitorio, risonanza e risposta in frequenza. In controlli e circuiti è una descrizione compatta del sistema.

7. Errori, aritmetica finita e condizionamento

Errore assoluto

$$e_a=|x-\widetilde{x}|$$

L’errore assoluto misura la distanza tra valore esatto $x$ e approssimazione $\widetilde{x}$. È espresso nelle stesse unità della grandezza. È utile quando la scala fisica del problema è fissata.

Errore relativo

$$e_r=\frac{|x-\widetilde{x}|}{|x|}$$

L’errore relativo misura l’errore rispetto alla grandezza del valore esatto. È adimensionale e permette confronti tra quantità di scala diversa. Non è definito se $x=0$, caso in cui serve un criterio assoluto o una scala di riferimento.

Cifre significative

$$e_r\lesssim 10^{-p}$$

Un errore relativo dell’ordine di $10^{-p}$ indica approssimativamente $p$ cifre decimali corrette. La relazione è indicativa, non una definizione rigida. In calcolo scientifico è più importante controllare errori relativi e residui che contare cifre visibili.

Rappresentazione floating point

$$x=\pm m\,\beta^e$$

Un numero floating point è rappresentato da segno, mantissa, base ed esponente. Poiché la mantissa ha lunghezza finita, solo un sottoinsieme discreto dei reali è rappresentabile. Ogni operazione aritmetica può introdurre arrotondamento.

Unità di arrotondamento

$$u\approx \frac{1}{2}\beta^{1-p}$$

L’unità di arrotondamento misura l’errore relativo massimo introdotto arrotondando un numero rappresentabile con $p$ cifre in base $\beta$. È una grandezza di macchina: dipende dal formato numerico. Nei formati binari IEEE è legata alla precisione della mantissa.

Modello di arrotondamento

$$\operatorname{fl}(x\circ y)=(x\circ y)(1+\delta),\qquad |\delta|\le u$$

Per un’operazione elementare $\circ$, il risultato floating point è il risultato esatto perturbato da un piccolo errore relativo. Il modello è valido sotto ipotesi normali, senza overflow, underflow o casi eccezionali. È la base dell’analisi degli errori.

Cancellazione numerica

$$x\approx y \qquad\Longrightarrow\qquad x-y\ \text{può perdere cifre significative}$$

Sottrarre numeri quasi uguali elimina le cifre comuni e amplifica gli errori relativi. La cancellazione non è sempre evitabile, ma va riconosciuta. Molte formule matematicamente corrette sono numericamente instabili proprio per cancellazione.

Condizionamento di un problema

$$\kappa=\frac{\text{errore relativo sull'output}}{\text{errore relativo sull'input}}$$

Il condizionamento misura la sensibilità del problema ai dati. Un problema mal condizionato amplifica inevitabilmente piccoli errori nei dati, indipendentemente dall’algoritmo. Stabilità dell’algoritmo e condizionamento del problema sono concetti diversi.

Stabilità di un algoritmo

$$\text{risultato calcolato}=\text{soluzione esatta di un problema vicino}$$

Un algoritmo stabile produce il risultato esatto di un problema leggermente perturbato. Se il problema è ben condizionato, questo implica un risultato accurato. Se il problema è mal condizionato, anche un algoritmo stabile può produrre risultati sensibili.

Errore in avanti e all’indietro

$$\text{forward error}=|\widetilde{x}-x|,\qquad \text{backward error}=|\Delta d|$$

L’errore in avanti misura quanto la soluzione calcolata differisce da quella esatta. L’errore all’indietro misura quanto bisogna perturbare i dati perché la soluzione calcolata diventi esatta. L’analisi all’indietro è spesso più robusta per algoritmi numerici.

8. Zeri di funzione e metodi iterativi scalari

Problema dello zero

$$f(x)=0$$

Trovare uno zero significa risolvere un’equazione non lineare. In generale non esiste una formula chiusa, quindi si costruisce una successione di approssimazioni. La scelta del metodo dipende da continuità, derivabilità, disponibilità della derivata e intervalli di separazione.

Metodo di bisezione

$$f(a)f(b)<0$$

Se $f$ è continua e cambia segno su $[a,b]$, esiste almeno uno zero nell’intervallo. La bisezione dimezza ripetutamente l’intervallo scegliendo il sottointervallo dove persiste il cambio di segno. È lenta ma molto robusta.

Iterazione di bisezione

$$c_k=\frac{a_k+b_k}{2}$$

Il punto medio divide l’intervallo corrente in due metà. Si valuta $f(c_k)$ e si conserva la metà che contiene un cambio di segno. Il metodo garantisce riduzione deterministica dell’incertezza.

Errore della bisezione

$$|x^\ast-c_k|\le \frac{b-a}{2^{k+1}}$$

Dopo $k$ iterazioni, lo zero è confinato in un intervallo di lunghezza $(b-a)/2^k$. L’errore sul punto medio è al massimo metà della lunghezza dell’intervallo corrente. La convergenza è lineare con fattore $1/2$.

Iterazione di punto fisso

$$x=g(x) \qquad\Longrightarrow\qquad x_{k+1}=g(x_k)$$

Risolvere $f(x)=0$ può essere trasformato in cercare un punto fisso di $g$. La trasformazione non è unica e condiziona la convergenza. Alcune scelte convergono, altre divergono, pur essendo algebricamente equivalenti.

Criterio di contrazione

$$|g'(x)|\le q<1$$

Se $g$ è una contrazione su un intervallo chiuso che manda l’intervallo in sé, l’iterazione converge all’unico punto fisso. Il fattore $q$ controlla la velocità di convergenza. Se $|g'|>1$ vicino al punto fisso, l’iterazione tende a essere instabile.

Metodo di Newton

$$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}$$

Newton sostituisce la funzione con la tangente nel punto corrente e usa lo zero della tangente come nuova approssimazione. È molto rapido vicino a uno zero semplice, ma può fallire se la derivata è piccola, se il punto iniziale è lontano o se la funzione ha geometria sfavorevole.

Convergenza quadratica di Newton

$$e_{k+1}\approx C e_k^2$$

Vicino a uno zero semplice e con regolarità sufficiente, il numero di cifre corrette cresce molto rapidamente. La convergenza quadratica è il grande vantaggio di Newton. Se lo zero è multiplo, la convergenza degrada tipicamente a lineare.

Newton per radice multipla

$$f(x^\ast)=0,\qquad f'(x^\ast)=0$$

Quando lo zero è multiplo, la tangente non taglia in modo efficace l’asse. Se la molteplicità $m$ è nota, una correzione frequente è moltiplicare il passo di Newton per $m$. In pratica, però, la molteplicità raramente è nota con certezza.

Metodo della secante

$$x_{k+1}=x_k-f(x_k)\frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})}$$

La secante approssima la derivata usando due valori della funzione. Evita il calcolo di $f'$, ma richiede due punti iniziali. La convergenza è superlineare, in genere più lenta di Newton ma più economica quando la derivata è costosa.

Criteri di arresto

$$|x_{k+1}-x_k|<\tau_x,\qquad |f(x_k)|<\tau_f$$

Un buon arresto controlla sia la variazione tra iterati sia il residuo. Un residuo piccolo non sempre implica errore piccolo se il problema è mal condizionato; una variazione piccola può indicare stagnazione. Le tolleranze devono essere coerenti con scala e unità del problema.

9. Sistemi non lineari

Sistema non lineare

$$F(x)=0,\qquad F:\mathbb{R}^n\to\mathbb{R}^n$$

Il problema generalizza lo zero scalare. La soluzione è un vettore e il residuo è un vettore. I metodi numerici richiedono norme per misurare errori e residui, e matrici Jacobiane per linearizzare.

Newton multidimensionale

$$J_F(x_k)s_k=-F(x_k),\qquad x_{k+1}=x_k+s_k$$

In più variabili Newton non divide per la derivata: risolve un sistema lineare con la Jacobiana. Il passo $s_k$ è la correzione che annulla il modello lineare di $F$ vicino a $x_k$. Il costo principale è fattorizzare o risolvere il sistema lineare.

Jacobiana

$$J_F(x)=\left(\frac{\partial F_i}{\partial x_j}(x)\right)$$

La Jacobiana raccoglie le sensibilità delle componenti del residuo rispetto alle variabili. Se è singolare o mal condizionata, il passo di Newton può essere instabile o enorme. Nei problemi grandi spesso si usano approssimazioni o metodi quasi-Newton.

Metodo quasi-Newton

$$B_k s_k=-F(x_k)$$

Invece di usare la Jacobiana esatta, si usa una matrice $B_k$ aggiornata iterativamente. L’obiettivo è ridurre il costo evitando derivate esatte o fattorizzazioni complete a ogni passo. La qualità dell’approssimazione determina robustezza e velocità.

Damping del passo

$$x_{k+1}=x_k+\alpha_k s_k,\qquad 0<\alpha_k\le1$$

Ridurre il passo può migliorare la convergenza globale quando il passo di Newton pieno è troppo aggressivo. Il parametro $\alpha_k$ viene scelto con criteri di diminuzione del residuo o funzioni merito. È una protezione pratica contro divergenze.

10. Sistemi lineari numerici

Sistema lineare

$$Ax=b$$

Risolvere un sistema lineare è una delle operazioni fondamentali del calcolo scientifico. Compare direttamente nei modelli lineari e indirettamente dentro Newton, minimi quadrati, discretizzazioni di EDO/EDP e ottimizzazione.

Residuo

$$r=b-A\widetilde{x}$$

Il residuo misura quanto la soluzione approssimata viola l’equazione. Un residuo piccolo è necessario, ma non sempre sufficiente per garantire errore piccolo: se $A$ è mal condizionata, piccoli residui possono corrispondere a grandi errori sulla soluzione.

Condizionamento di matrice

$$\kappa(A)=\lVert A\rVert\,\lVert A^{-1}\rVert$$

Il numero di condizionamento misura quanto la soluzione di $Ax=b$ è sensibile a perturbazioni nei dati. Se $\kappa(A)$ è grande, il problema amplifica errori di arrotondamento e incertezza sui dati. Se $A$ è singolare, $\kappa(A)$ è infinito.

Stima dell’errore relativo

$$\frac{\lVert x-\widetilde{x}\rVert}{\lVert x\rVert} \lesssim \kappa(A)\frac{\lVert r\rVert}{\lVert b\rVert}$$

Questa stima collega errore, residuo e condizionamento. Mostra che il residuo va interpretato alla luce di $\kappa(A)$. Un sistema ben condizionato consente di fidarsi maggiormente del residuo; uno mal condizionato richiede cautela.

Eliminazione di Gauss

$$Ax=b\quad\longrightarrow\quad Ux=c$$

L’eliminazione trasforma il sistema in uno triangolare superiore tramite operazioni elementari sulle righe. Poi si risolve per sostituzione all’indietro. È un metodo diretto: in aritmetica esatta termina dopo un numero finito di operazioni.

Pivoting parziale

$$\text{scambio di righe per scegliere un pivot grande}$$

Il pivoting riduce il rischio di divisioni per numeri piccoli e migliora la stabilità numerica. Senza pivoting, l’eliminazione di Gauss può essere instabile anche per matrici non singolari. In pratica il pivoting parziale è lo standard generale.

Fattorizzazione LU

$$PA=LU$$

Con pivoting, la matrice permutata $PA$ viene fattorizzata in una triangolare inferiore $L$ e una superiore $U$. Risolvere $Ax=b$ diventa risolvere due sistemi triangolari. La fattorizzazione è utile quando si devono risolvere più sistemi con la stessa matrice e termini noti diversi.

Sistemi triangolari

$$Ly=b,\qquad Ux=y$$

Un sistema triangolare inferiore si risolve per sostituzione in avanti; uno superiore per sostituzione all’indietro. Il costo è dell’ordine di $n^2$, molto inferiore alla fattorizzazione, che costa dell’ordine di $n^3$.

Fattorizzazione di Cholesky

$$A=LL^T$$

Se $A$ è simmetrica definita positiva, Cholesky è più efficiente e stabile della LU generale. Compare in problemi di energia, minimi quadrati, elementi finiti e matrici di covarianza. La positività definita è un’ipotesi essenziale.

Matrici sparse

$$\text{molti elementi di }A\ \text{sono nulli}$$

Nei problemi discretizzati grandi, la matrice è spesso sparsa. Conservare e sfruttare la sparsità riduce memoria e costo. Metodi diretti e iterativi devono essere scelti tenendo conto del riempimento e della struttura della matrice.

11. Metodi iterativi per sistemi lineari

Decomposizione iterativa

$$A=M-N,\qquad x^{(k+1)}=M^{-1}Nx^{(k)}+M^{-1}b$$

I metodi iterativi sostituiscono la soluzione diretta con una successione di approssimazioni. La scelta di $M$ deve rendere facile risolvere sistemi con $M$. La convergenza dipende dalla matrice di iterazione $M^{-1}N$.

Criterio spettrale di convergenza

$$\rho(B)<1$$

Per un’iterazione lineare $x^{(k+1)}=Bx^{(k)}+c$, la convergenza per ogni dato iniziale avviene se il raggio spettrale di $B$ è minore di $1$. Il raggio spettrale è il massimo modulo degli autovalori.

Metodo di Jacobi

$$x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}\right)$$

Jacobi aggiorna ogni componente usando solo i valori dell’iterazione precedente. È semplice e parallelizzabile, ma può convergere lentamente. Richiede che gli elementi diagonali usati come divisori siano non nulli.

Metodo di Gauss-Seidel

$$x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum_{j<i}a_{ij}x_j^{(k+1)}-\sum_{j>i}a_{ij}x_j^{(k)}\right)$$

Gauss-Seidel usa subito i valori appena aggiornati. Spesso converge più rapidamente di Jacobi, ma è meno parallelizzabile. L’ordine delle variabili può influenzare la convergenza.

Metodo SOR

$$x^{(k+1)}=(1-\omega)x^{(k)}+\omega x_{\mathrm{GS}}^{(k+1)}$$

SOR introduce un parametro di rilassamento $\omega$. Per $\omega=1$ coincide con Gauss-Seidel. Valori opportuni di $\omega$ possono accelerare molto la convergenza, ma una scelta sbagliata può peggiorarla o far divergere il metodo.

Gradiente coniugato

$$A=A^T>0$$

Il gradiente coniugato risolve sistemi simmetrici definiti positivi sfruttando direzioni coniugate rispetto ad $A$. In aritmetica esatta converge in al più $n$ passi, ma in pratica viene usato come metodo iterativo con arresto anticipato. È centrale per matrici sparse grandi.

Residuo nel gradiente coniugato

$$r_k=b-Ax_k$$

Il residuo guida la costruzione delle direzioni e fornisce un criterio di arresto. Il metodo minimizza l’errore in norma energetica su sottospazi di Krylov. Il precondizionamento è spesso decisivo per renderlo efficace.

Precondizionamento

$$M^{-1}Ax=M^{-1}b$$

Un precondizionatore $M$ approssima $A$ ma deve essere più facile da invertire. L’obiettivo è ridurre il condizionamento effettivo e accelerare la convergenza. Un buon precondizionatore bilancia costo di applicazione e miglioramento spettrale.

12. Interpolazione polinomiale

Problema di interpolazione

$$p_n(x_i)=y_i,\qquad i=0,\dots,n$$

L’interpolazione costruisce un polinomio di grado al più $n$ che passa per $n+1$ dati. Se i nodi $x_i$ sono distinti, il polinomio interpolante esiste ed è unico. L’interpolazione non è sempre approssimazione ottimale: passare per tutti i punti può amplificare oscillazioni.

Base di Lagrange

$$L_i(x)=\prod_{\substack{j=0\\j\ne i}}^{n}\frac{x-x_j}{x_i-x_j}$$

$L_i$ vale $1$ nel nodo $x_i$ e $0$ in tutti gli altri nodi. Questa proprietà rende immediata la costruzione dell’interpolante. La formula è elegante, ma per calcolo numerico diretto può essere meno efficiente di forme alternative.

Polinomio di Lagrange

$$p_n(x)=\sum_{i=0}^{n}y_i L_i(x)$$

Ogni dato $y_i$ pesa la base che vale uno nel proprio nodo. La somma passa automaticamente per tutti i punti assegnati. Se si aggiunge un nuovo nodo, però, la forma di Lagrange va sostanzialmente ricostruita.

Errore di interpolazione

$$f(x)-p_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^{n}(x-x_i)$$

La formula vale se $f$ è sufficientemente derivabile e $\xi$ è un punto dell’intervallo dipendente da $x$. L’errore dipende sia dalla regolarità di $f$ sia dalla disposizione dei nodi. Il prodotto nodale può diventare grande e causare oscillazioni.

Differenze divise

$$f[x_i,x_{i+1}]=\frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}$$

Le differenze divise generalizzano i coefficienti incrementali. Sono la base della forma di Newton dell’interpolante. Permettono di aggiungere nodi aggiornando i coefficienti senza ricominciare da zero.

Forma di Newton

$$p_n(x)=f[x_0]+f[x_0,x_1](x-x_0)+\cdots+f[x_0,\dots,x_n]\prod_{j=0}^{n-1}(x-x_j)$$

La forma di Newton è gerarchica: ogni nuovo nodo aggiunge un termine. È utile quando i dati arrivano progressivamente. Dal punto di vista numerico, la valutazione si fa con uno schema annidato simile a Horner.

Fenomeno di Runge

$$\text{nodi equispaziati ad alto grado}\quad\Longrightarrow\quad\text{oscillazioni agli estremi}$$

Interpolare con polinomi di grado alto su nodi equispaziati può produrre grandi oscillazioni vicino agli estremi. Il fenomeno non è un difetto del software, ma della scelta del modello interpolante. Nodi di Chebyshev o splines riducono il problema.

Nodi di Chebyshev

$$x_k=\cos\left(\frac{2k+1}{2n+2}\pi\right),\qquad k=0,\dots,n$$

I nodi di Chebyshev si addensano vicino agli estremi dell’intervallo. Questa distribuzione riduce il massimo del prodotto nodale e migliora la stabilità dell’interpolazione globale. Sono una scelta naturale per approssimazione polinomiale su intervalli.

Splines cubiche

$$s_i(x)\ \text{cubo su }[x_i,x_{i+1}]$$

Le splines usano polinomi di basso grado a tratti invece di un unico polinomio globale di grado alto. Le cubiche impongono continuità di funzione, derivata prima e derivata seconda nei nodi interni. Sono robuste e molto usate in interpolazione di dati.

Spline naturale

$$s''(x_0)=0,\qquad s''(x_n)=0$$

La spline naturale impone curvatura nulla agli estremi. È una condizione al bordo, non una verità fisica universale. Se si conoscono pendenze o curvature agli estremi, possono essere più appropriate altre condizioni.

13. Approssimazione e minimi quadrati

Approssimazione in norma

$$\min_{p\in V}\lVert f-p\rVert$$

Approssimare significa scegliere, dentro uno spazio semplice $V$, l’elemento più vicino a $f$ secondo una norma. La scelta della norma cambia il significato di “migliore”. In ingegneria si usa spesso la norma quadratica perché porta a problemi lineari e interpreta energia dell’errore.

Minimi quadrati discreti

$$\min_c \sum_{i=1}^{m}\left(y_i-\sum_{j=1}^{n}c_j\phi_j(x_i)\right)^2$$

Si cercano coefficienti $c_j$ che minimizzino la somma dei quadrati dei residui. A differenza dell’interpolazione, non si pretende di passare per tutti i dati. Questo è adatto a misure rumorose o dati sovradeterminati.

Forma matriciale

$$\min_c \lVert Ac-y\rVert_2^2$$

La matrice $A$ contiene i valori delle funzioni base nei punti dati. Il vettore $c$ contiene i coefficienti incogniti. Il problema cerca la proiezione del vettore dati sullo spazio generato dalle colonne di $A$.

Equazioni normali

$$A^TAc=A^Ty$$

Le equazioni normali caratterizzano il minimo quadratico quando $A$ ha rango pieno. Sono semplici, ma possono peggiorare il condizionamento perché coinvolgono $A^TA$. Per problemi mal condizionati sono preferibili QR o SVD.

Decomposizione QR

$$A=QR,\qquad Q^TQ=I$$

La fattorizzazione QR separa una parte ortogonale $Q$ e una triangolare $R$. Le trasformazioni ortogonali preservano la norma e sono numericamente stabili. Per minimi quadrati, QR evita di formare esplicitamente $A^TA$.

Proiezione ortogonale

$$A c^\ast=P_{\operatorname{Col}(A)}y$$

La soluzione dei minimi quadrati produce il vettore nello spazio colonna di $A$ più vicino ai dati. Il residuo $y-Ac^\ast$ è ortogonale a tutte le colonne di $A$. Questa è la geometria dietro le equazioni normali.

14. Derivazione numerica

Differenza in avanti

$$f'(x)\approx \frac{f(x+h)-f(x)}{h}$$

La formula usa un punto a destra di $x$. L’errore di troncamento è dell’ordine di $h$ se $f$ è regolare. Se $h$ è troppo piccolo, però, l’arrotondamento e la cancellazione tra valori vicini possono dominare.

Differenza all’indietro

$$f'(x)\approx \frac{f(x)-f(x-h)}{h}$$

È l’analogo causale della differenza in avanti: usa il valore corrente e uno passato. Nei problemi evolutivi può essere preferibile perché non richiede dati futuri. Anche questa formula è del primo ordine.

Differenza centrata

$$f'(x)\approx \frac{f(x+h)-f(x-h)}{2h}$$

La differenza centrata usa punti simmetrici e cancella il termine di errore di ordine $h$. Per funzioni lisce ha errore di troncamento dell’ordine di $h^2$. È più accurata delle formule unilaterali, ma richiede valori su entrambi i lati.

Seconda derivata centrata

$$f''(x)\approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}$$

Questa formula misura la curvatura tramite la differenza tra il valore centrale e i vicini. È alla base delle discretizzazioni del Laplaciano in una dimensione. L’errore di troncamento è dell’ordine di $h^2$ per funzioni sufficientemente lisce.

Bilancio tra troncamento e arrotondamento

$$E(h)\approx C_1h^p+\frac{C_2u}{h}$$

Ridurre $h$ diminuisce l’errore di troncamento, ma aumenta l’effetto dell’arrotondamento nelle differenze. Esiste quindi un passo ottimo approssimativo. La derivazione numerica è intrinsecamente delicata perché amplifica rumore e errori sui dati.

15. Integrazione numerica

Formula di quadratura

$$\int_a^b f(x)\,\mathrm{d}x \approx \sum_{i=0}^{n}w_i f(x_i)$$

Una formula di quadratura sostituisce l’integrale con una somma pesata di valori della funzione. I nodi $x_i$ e i pesi $w_i$ determinano accuratezza, stabilità e costo. Le quadrature sono indispensabili quando la primitiva non è nota o quando i dati sono campionati.

Formula del punto medio

$$\int_a^b f(x)\,\mathrm{d}x \approx (b-a)f\left(\frac{a+b}{2}\right)$$

La formula usa il valore al centro dell’intervallo. È esatta per funzioni lineari e ha errore legato alla seconda derivata. Spesso è più accurata del rettangolo sinistro o destro perché sfrutta la simmetria.

Formula dei trapezi semplice

$$\int_a^b f(x)\,\mathrm{d}x \approx \frac{b-a}{2}\bigl(f(a)+f(b)\bigr)$$

La funzione viene approssimata con la retta che unisce gli estremi. L’area sotto questa retta è un trapezio. La formula è esatta per polinomi di grado al più uno.

Errore del trapezio semplice

$$E_T=-\frac{(b-a)^3}{12}f''(\xi)$$

L’errore dipende dalla curvatura della funzione. Se $f$ è lineare, la seconda derivata è nulla e la formula è esatta. Il punto $\xi$ è un punto dell’intervallo, generalmente non noto.

Formula di Simpson semplice

$$\int_a^b f(x)\,\mathrm{d}x \approx \frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]$$

Simpson approssima la funzione con una parabola passante per gli estremi e il punto medio. È esatta per polinomi fino al terzo grado. Il peso $4$ al punto medio riflette il ruolo centrale del valore interno.

Errore di Simpson

$$E_S=-\frac{(b-a)^5}{2880}f^{(4)}(\xi)$$

L’errore dipende dalla quarta derivata. Per funzioni molto regolari Simpson converge rapidamente. Se la funzione non è liscia, l’ordine teorico può degradare.

Trapezi composita

$$T_h=h\left[\frac{f(a)+f(b)}{2}+\sum_{i=1}^{n-1}f(a+ih)\right], \qquad h=\frac{b-a}{n}$$

L’intervallo viene diviso in sottointervalli e si applica il trapezio su ciascuno. La formula è semplice e robusta. Per funzioni periodiche lisce può essere sorprendentemente efficace.

Simpson composita

$$S_h=\frac{h}{3}\left[f(x_0)+f(x_n)+4\sum_{\substack{i=1\\i\ \text{dispari}}}^{n-1}f(x_i)+2\sum_{\substack{i=2\\i\ \text{pari}}}^{n-2}f(x_i)\right]$$

La formula richiede un numero pari di sottointervalli. I pesi alternano $4$ e $2$ sui nodi interni. È una delle quadrature composite più usate per funzioni regolari su intervalli finiti.

Quadratura gaussiana

$$\int_a^b f(x)w(x)\,\mathrm{d}x \approx \sum_{i=1}^{n}\omega_i f(x_i)$$

La quadratura gaussiana sceglie nodi e pesi per massimizzare il grado di esattezza polinomiale. I nodi sono legati a polinomi ortogonali rispetto al peso $w$. A parità di valutazioni può essere molto più accurata di formule a nodi equispaziati.

Grado di esattezza

$$\text{formula esatta per ogni polinomio di grado }\le m$$

Il grado di esattezza misura la capacità della quadratura di integrare polinomi. Non garantisce da solo accuratezza su funzioni non polinomiali, ma è un indicatore importante. Per funzioni lisce, maggiore grado di esattezza spesso significa errore più piccolo.

16. Metodi numerici per EDO: concetti base

Griglia temporale

$$t_n=t_0+nh,\qquad y_n\approx y(t_n)$$

La soluzione continua viene approssimata su una griglia discreta. Il passo $h$ controlla risoluzione, costo e stabilità. Un passo più piccolo non garantisce automaticamente un risultato migliore se l’arrotondamento o la stabilità diventano critici.

Metodo one-step

$$y_{n+1}=y_n+h\Phi(t_n,y_n,h)$$

Un metodo a un passo calcola il nuovo valore usando solo il valore corrente. La funzione $\Phi$ approssima la pendenza media sul passo. Euler e Runge-Kutta sono esempi fondamentali.

Errore locale di troncamento

$$\tau_{n+1}=\frac{y(t_{n+1})-y(t_n)}{h}-\Phi(t_n,y(t_n),h)$$

L’errore locale misura quanto il metodo sbaglia in un singolo passo assumendo esatto il valore di partenza. È una proprietà della formula numerica. Un errore locale piccolo non basta: bisogna anche controllare come gli errori si propagano.

Errore globale

$$e_n=y(t_n)-y_n$$

L’errore globale è la differenza tra soluzione esatta e approssimazione dopo molti passi. Include l’accumulo e l’amplificazione degli errori locali. Stabilità e consistenza determinano la convergenza.

Ordine di un metodo

$$\max_n |e_n|=O(h^p)$$

Un metodo è di ordine $p$ se dimezzare il passo riduce l’errore globale di circa $2^p$, nel regime asintotico. L’ordine descrive il comportamento per $h$ piccolo, ma non sostituisce il controllo di stabilità.

Consistenza

$$\tau_{n+1}\to0\qquad(h\to0)$$

Un metodo è consistente se riproduce l’equazione differenziale nel limite di passo nullo. Senza consistenza, anche passi sempre più piccoli non portano alla soluzione corretta. La consistenza è necessaria per la convergenza.

Stabilità numerica

$$\text{piccole perturbazioni non devono amplificarsi incontrollatamente}$$

La stabilità riguarda la propagazione degli errori. Un metodo può essere localmente accurato ma globalmente inutile se amplifica perturbazioni. Nei problemi rigidi la stabilità può imporre passi molto più piccoli di quelli richiesti dall’accuratezza.

17. Metodi a un passo per EDO

Eulero esplicito

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

Eulero esplicito usa la pendenza all’inizio dell’intervallo. È semplice e poco costoso, ma di ordine uno e con stabilità limitata. È utile come introduzione, come predittore o in problemi molto semplici, ma spesso non basta per calcolo affidabile.

Eulero implicito

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

Il nuovo valore compare anche nel lato destro, quindi a ogni passo bisogna risolvere un’equazione, spesso non lineare. Il costo è maggiore, ma la stabilità è molto migliore. È importante per problemi rigidi e dissipativi.

Eulero migliorato o Heun

$$k_1=f(t_n,y_n),\qquad k_2=f(t_n+h,y_n+h k_1)$$ $$y_{n+1}=y_n+\frac{h}{2}(k_1+k_2)$$

Heun usa una pendenza predetta all’inizio e una alla fine del passo, poi ne prende la media. È un metodo esplicito di ordine due. Migliora Eulero senza arrivare al costo del Runge-Kutta classico.

Punto medio esplicito

$$k_1=f(t_n,y_n),\qquad k_2=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1\right)$$ $$y_{n+1}=y_n+h k_2$$

Il metodo stima la pendenza al centro del passo. È anch’esso di ordine due. La simmetria temporale migliora l’accuratezza rispetto all’uso della sola pendenza iniziale.

Runge-Kutta classico del quarto ordine

$$k_1=f(t_n,y_n)$$ $$k_2=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1\right),\qquad k_3=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2\right)$$ $$k_4=f(t_n+h,y_n+h k_3)$$ $$y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$$

RK4 combina quattro valutazioni della pendenza per ottenere ordine quattro. È molto usato perché offre buon equilibrio tra accuratezza, semplicità e costo. Non è però automaticamente stabile per qualunque passo.

Controllo adattivo del passo

$$\text{errore stimato}\le \text{tolleranza}$$

I metodi adattivi scelgono il passo in base a una stima dell’errore locale. Se l’errore è troppo grande, il passo viene ridotto; se è molto piccolo, può essere aumentato. Questo è essenziale quando la soluzione ha zone lente e zone rapide.

Metodo embedded

$$y_{n+1}^{[p]},\qquad y_{n+1}^{[p+1]}$$

Un metodo embedded produce due approssimazioni di ordine diverso usando valutazioni condivise. La differenza tra le due stima l’errore. Le famiglie Runge-Kutta-Fehlberg e Dormand-Prince seguono questa idea.

18. Metodi multistep, rigidità e stabilità

Metodo multistep

$$y_{n+1}=\sum_{j=0}^{k-1}\alpha_j y_{n-j} +h\sum_{j=0}^{k}\beta_j f(t_{n+1-j},y_{n+1-j})$$

I metodi multistep usano più valori precedenti della soluzione e della derivata. Possono essere efficienti perché riutilizzano informazione già calcolata. Richiedono però valori iniziali prodotti da un altro metodo e una gestione attenta della stabilità.

Adams-Bashforth a due passi

$$y_{n+1}=y_n+\frac{h}{2}\left(3f_n-f_{n-1}\right)$$

È un metodo esplicito che extrapola la pendenza usando due valori precedenti. Ha ordine due. Essendo esplicito, è semplice, ma la regione di stabilità è limitata.

Adams-Moulton a due passi

$$y_{n+1}=y_n+\frac{h}{2}\left(f_{n+1}+f_n\right)$$

Questa formula è implicita perché contiene $f_{n+1}$. Coincide con il trapezio applicato all’equazione differenziale. Ha migliore stabilità rispetto a metodi espliciti dello stesso ordine, ma richiede la soluzione di un’equazione a ogni passo.

Predittore-correttore

$$\text{predizione esplicita}\quad\longrightarrow\quad\text{correzione implicita}$$

Si usa un metodo esplicito per stimare $y_{n+1}$ e poi un metodo implicito per correggerlo. Questa strategia combina costo ridotto e maggiore stabilità. Il numero di correzioni influenza accuratezza e costo.

Equazione test di stabilità

$$y'=\lambda y,\qquad \operatorname{Re}\lambda<0$$

Questa equazione rappresenta un modo esponenziale decadente. Un metodo numerico dovrebbe riprodurre il decadimento senza trasformarlo in crescita artificiale. La risposta del metodo a questa equazione definisce la regione di stabilità assoluta.

Stabilità di Eulero esplicito

$$y_{n+1}=(1+h\lambda)y_n$$

La stabilità richiede $|1+h\lambda|<1$. Se $\lambda$ è reale negativo, serve $0<h<2/|\lambda|$. Per problemi con modi molto rapidi, questa restrizione può rendere Eulero esplicito impraticabile.

Rigidità

$$\text{scale temporali molto diverse}$$

Un problema è rigido quando contiene modi rapidi stabili che impongono passi piccolissimi ai metodi espliciti, anche se la soluzione interessante varia lentamente. La rigidità è un problema di stabilità, non solo di accuratezza. Metodi impliciti sono spesso necessari.

A-stabilità

$$\{z\in\mathbb{C}:\operatorname{Re}z<0\}\subset \text{regione di stabilità}$$

Un metodo A-stabile è stabile per l’intero semipiano sinistro sull’equazione test. Questa proprietà è preziosa per problemi rigidi. Eulero implicito è A-stabile; Eulero esplicito no.

L-stabilità

$$R(z)\to0\qquad(z\to-\infty)$$

Un metodo L-stabile non solo è stabile per modi molto rapidi, ma li smorza fortemente. Questo è desiderabile nei problemi dissipativi, dove i transitori veloci devono scomparire senza lasciare oscillazioni numeriche.

19. Differenze finite per problemi al contorno

Problema al contorno lineare

$$-y''(x)=f(x),\qquad y(a)=\alpha,\quad y(b)=\beta$$

Un problema al contorno impone condizioni in punti diversi dell’intervallo. Non si può avanzare semplicemente nel tempo come in un problema di Cauchy. La discretizzazione produce di solito un sistema lineare.

Griglia uniforme

$$x_i=a+ih,\qquad h=\frac{b-a}{N}$$

La griglia divide l’intervallo in sottointervalli uguali. I valori incogniti sono $y_i\approx y(x_i)$ nei nodi interni. Le condizioni al bordo fissano i valori agli estremi.

Approssimazione della seconda derivata

$$y''(x_i)\approx\frac{y_{i-1}-2y_i+y_{i+1}}{h^2}$$

Questa formula centrata è di ordine due. Applicata a ogni nodo interno trasforma l’equazione differenziale in equazioni algebriche. La matrice risultante è tridiagonale nel caso più semplice.

Sistema tridiagonale

$$-y_{i-1}+2y_i-y_{i+1}=h^2 f(x_i)$$

Ogni equazione coinvolge solo il nodo corrente e i due vicini. La struttura tridiagonale consente soluzioni molto efficienti. Nei problemi di diffusione e Poisson unidimensionali questa struttura è ricorrente.

Metodo di shooting

$$y(a)=\alpha,\qquad y'(a)=s$$

Lo shooting trasforma un problema al contorno in problemi di Cauchy parametrizzati dalla pendenza iniziale $s$. Si sceglie $s$ in modo che la soluzione soddisfi la condizione finale. Il metodo può essere instabile se il problema amplifica molto le variazioni di $s$.

Residuo dello shooting

$$\Phi(s)=y_s(b)-\beta$$

Il valore corretto di $s$ è uno zero di $\Phi$. Si possono usare bisezione, secante o Newton. Il costo di valutare $\Phi$ include la soluzione numerica di una EDO.

20. Schemi operativi di risoluzione

Analisi di una EDO

$$\text{ordine},\quad \text{linearità},\quad \text{autonomia},\quad \text{dati},\quad \text{dominio}$$

Prima di risolvere, classificare l’equazione. Riconoscere ordine e linearità suggerisce i metodi disponibili; autonomia e segno del campo suggeriscono analisi qualitativa; dominio e dati iniziali chiariscono dove la soluzione ha senso.

Scelta del metodo analitico

$$\text{separabile},\quad \text{lineare},\quad \text{esatta},\quad \text{coefficienti costanti}$$

Le forme riconoscibili guidano la soluzione. Se l’equazione è separabile, si integra separando le variabili; se è lineare del primo ordine, si usa il fattore integrante; se è esatta, si cerca un potenziale; se ha coefficienti costanti, si usa l’equazione caratteristica.

Scelta del metodo per zeri

$$\text{robustezza}\Rightarrow\text{bisezione},\qquad \text{rapidità locale}\Rightarrow\text{Newton},\qquad \text{derivata non disponibile}\Rightarrow\text{secante}$$

La bisezione è lenta ma garantita se c’è cambio di segno. Newton è veloce ma sensibile al punto iniziale e alla derivata. La secante riduce il costo evitando la derivata, con convergenza di solito intermedia.

Scelta per sistemi lineari

$$\text{denso piccolo}\Rightarrow\text{LU},\qquad \text{simmetrico definito positivo}\Rightarrow\text{Cholesky},\qquad \text{sparso grande}\Rightarrow\text{iterativi}$$

La struttura della matrice decide il metodo. Usare un metodo generale ignorando simmetria, positività o sparsità spreca costo e può peggiorare stabilità. Il precondizionamento è spesso la differenza tra un metodo iterativo utile e uno inutilizzabile.

Scelta per interpolazione

$$\text{pochi nodi}\Rightarrow\text{polinomio},\qquad \text{molti dati}\Rightarrow\text{spline o minimi quadrati}$$

Interpolare molti punti con un unico polinomio di grado alto può essere instabile. Le splines mantengono basso il grado locale; i minimi quadrati accettano residui e sono adatti a dati rumorosi. La scelta dipende dallo scopo: passare esattamente dai dati o approssimare una tendenza.

Scelta per quadratura

$$\text{dati campionati}\Rightarrow\text{trapezi o Simpson},\qquad \text{funzione valutabile}\Rightarrow\text{Gauss o adattiva}$$

Se i dati sono disponibili solo su una griglia, le formule composite sono naturali. Se la funzione può essere valutata liberamente, quadrature gaussiane o adattive possono ridurre molto il numero di valutazioni. Singolarità e discontinuità richiedono spezzare l’intervallo o cambiare variabile.

Scelta per EDO numeriche

$$\text{non rigida}\Rightarrow\text{Runge-Kutta espliciti},\qquad \text{rigida}\Rightarrow\text{metodi impliciti}$$

Per problemi non rigidi, metodi espliciti adattivi sono spesso efficienti. Per problemi rigidi, la stabilità impone metodi impliciti o specializzati. La tolleranza deve essere scelta in base alla scala delle variabili e alla precisione richiesta dal modello.

Controllo finale

$$\text{residuo},\quad \text{errore stimato},\quad \text{condizionamento},\quad \text{stabilità},\quad \text{scala fisica}$$

Un risultato numerico va sempre controllato. Il residuo dice se l’equazione è soddisfatta; l’errore stimato dice quanto è accurata l’approssimazione; il condizionamento dice quanto fidarsi della sensibilità; la stabilità dice se il metodo sta amplificando errori; la scala fisica dice se il risultato ha senso nel modello.