Computazionali, metodi

Enciclopedia della Scienza e della Tecnica (2007)

Computazionali, metodi

Alfio Quarteroni

I metodi computazionali permettono di risolvere con i computer, nell'ambito delle scienze applicate, problemi complessi formulabili tramite il linguaggio della matematica e raramente risolubili per via analitica. In generale, la loro soluzione non ammette una rappresentazione in forma esplicita: basti pensare a equazioni o sistemi non lineari (algebrici, differenziali o integrali), di cui non siano note formule risolutive. In altri casi la forma esplicita, pur nota, non è praticamente utilizzabile per determinare valori quantitativi della soluzione stessa, o perché essa richiederebbe una quantità proibitiva di operazioni, o perché espressa a sua volta attraverso altri enti matematici (quali integrali o sviluppi in serie), di cui sarebbe difficile dare una valutazione quantitativa. La risoluzione di problemi di interesse applicativo comporta diverse fasi: (a) la comprensione qualitativa e quantitativa del problema; (b) la sua modellazione attraverso equazioni matematiche che possono essere algebriche, funzionali, differenziali o integrali (modello matematico); (c) l'individuazione di metodi della matematica numerica che siano idonei ad approssimare tale modello matematico (discretizzazione); (d) l'implementazione di tali metodi numerici di approssimazione su calcolatori.

I metodi computazionali comprendono ed esauriscono le fasi (c) e (d) del processo. La scelta di tali metodi non può prescindere da una conoscenza adeguata delle proprietà qualitative della soluzione del modello matematico e da una buona percezione della fenomenologia del problema originario. È pertanto giustificata la terminologia di modelli computazionali che viene usata in alternativa a quella di metodi computazionali. La modellistica computazionale è un settore di ricerca interdisciplinare, che si trova alla confluenza della matematica, dell'informatica e delle scienze applicate. La sua evoluzione è stata molto rapida negli anni recenti, grazie al vigoroso sviluppo dei calcolatori e degli algoritmi, nonché a una maggior consapevolezza del ruolo che il calcolo scientifico può rivestire nella simulazione di complessi problemi di interesse scientifico/teorico, industriale, ambientale e sociale.

fig. 1

Si può riassumere in forma schematica il processo (a)-(d) con l'ausilio di un diagramma (fig. 1). Indichiamo con uf la soluzione del problema originario, qui detto per brevità problema fisico. Il modello matematico che descrive il problema fisico è posto nella forma F(u,d)=0, dove si suppone che d descriva l'insieme dei dati, u la soluzione, e F la relazione funzionale che lega fra loro dati e soluzione. Vari e diversificati sono i problemi matematici che, in astratto, possono essere formulati in questo modo. Si pensi all'integrazione definita di una funzione, alla ricerca delle radici di equazioni algebriche, alla soluzione di problemi differenziali (ordinari o alle derivate parziali).

Svariati altri modelli matematici sono oggetto d'interesse per la modellistica computazionale. Si pensi, per esempio, al problema della rappresentazione di curve e superfici, al trattamento numerico di immagini e segnali, alla minimizzazione di funzionali liberi o vincolati, alla regolarizzazione numerica e risoluzione di problemi inversi, ai problemi integro-differenziali, nonché al controllo di sistemi differenziali, integrali o di natura stocastica. La rilevanza di tali modelli nell'analisi e simulazione di processi meccanici, fisici, economici, biologici ecc. è notevole e destinata a crescere nel tempo.

Convergenza e stabilità dei modelli numerici

La prima fase della modellistica computazionale consiste nella sostituzione del modello matematico F(u,d)=0 con uno numerico Fn(un,dn)=0, dipendente da un parametro n il cui significato varia a seconda della famiglia di problemi considerata, ma che possiamo sempre pensare riferito alla dimensione del modello numerico medesimo. In tale contesto, dn è un'approssimazione di d, Fn una discretizzazione della relazione funzionale F e un la soluzione numerica. La proprietà fondamentale che si richiede a un modello numerico è la convergenza, ovvero che (in un'opportuna topologia) limn→∞un=u. Condizione necessaria perché ciò avvenga è che Fn traduca correttamente la legge funzionale F, ovvero che Fn(u,d)→0 per n→∞, essendo d un dato ammissibile per il modello matematico e u la corrispondente soluzione. Tale proprietà, detta consistenza, non è peraltro sufficiente a garantire la convergenza della soluzione numerica a quella esatta: a tale scopo serve che il modello numerico sia stabile. Al fine di definire cosa si intende con stabilità, ricordiamo che un modello matematico si dice ben posto se ammette una sola soluzione e se, a piccole perturbazioni dei dati, corrispondono variazioni controllabili della soluzione. In altri termini, se supponiamo che d e d* rappresentino due insiemi di dati, δd=d*−d la loro differenza, u e u* le soluzioni corrispondenti e δu=u*−u la loro differenza, allora il problema sarà ben posto in corrispondenza del dato d se per ogni ε>0 esiste K(ε,d) tale che se ∥δd〈ε allora ∥δu∥≤K,d)∥δd∥. Abbiamo indicato con lo stesso simbolo ∥∙∥ delle norme opportune (non necessariamente le stesse) per l'insieme dei dati e per quello delle soluzioni. Se d≠0, il numero di condizionamento K(d) del modello è definito come:

[1] formula

formula

dove l'estremo superiore viene preso sull'insieme D di tutte le perturbazioni (non nulle) ammissibili dei dati.

Il numero K esprime dunque il massimo valore possibile del rapporto fra l'errore relativo sulla soluzione e l'errore relativo sui dati. I problemi per i quali K è dell'ordine dell'unità si dicono ben condizionati, se invece K è molto grande si diranno mal condizionati. Supponiamo che la soluzione u del problema modello si possa rappresentare come u=f(d) e analogamente la soluzione del modello numerico si possa ottenere come un=fn(dn), essendo f e fn due opportune leggi di corrispondenza. Allora, se f è differenziabile, l'uso del teorema del valor medio consente di ottenere per K(d) la seguente stima: K(d)=∥f′'(d)∥∙∥d∥/∥f(d)∥ (avendo trascurato infinitesimi di ordine superiore a ∥δd∥). Con ragionamenti del tutto analoghi possiamo associare al modello numerico un numero di condizionamento Kn(dn)=∥f′n(dn)∥∙∥dn∥/∥fn(dn)∥. Se l'insieme dei dati {dn} coincide con {d}, possiamo definire K* (d)=lim→∞supn≥Kn(d). Il numero K* viene detto numero di condizionamento asintotico del modello numerico.

Siamo ora in grado di precisare il concetto di stabilità numerica. Dato un modello matematico ben condizionato con numero di condizionamento K, diremo che la famiglia (dipendente dal parametro n) di modelli numerici associati è stabile se K* è dello stesso ordine di grandezza di K. In caso contrario, si parlerà di modello numerico instabile. Si noti che, in particolare, la stabilità assicura la limitatezza uniforme (rispetto a n) delle soluzioni numeriche. A titolo d'esempio, consideriamo il seguente modello matematico elementare: trovare x tale che x22px+1=0, dove p è il dato e le soluzioni sono le due radici x±=F(p)=p±√(p2−1). Per p>1, il numero di condizionamento corrispondente è K=∣p∣/(p2−1). Esso degenera per p→1 in quanto le due radici sono coincidenti per p=1; tuttavia, anche per p=1 il problema della ricerca delle radici è ben posto, come si vede trasformandolo nel problema equivalente F(x,t)=x2−[(1+t2)/t]x+1=0, con t=p+(p2−1), le cui radici x=t e x+=1/t risultano coincidenti per t=1.

Il cambio di parametro rimuove la singolarità presente nella precedente rappresentazione delle radici in funzione di p e il calcolo del condizionamento restituisce K≃1 per ogni valore di t. Supponendo p≫1 in modo che le radici siano ben separate, il problema è ben condizionato, ma il calcolo di x attraverso la formula esplicita è instabile per via della cancellazione di cifre significative, tipica dell'aritmetica finita del calcolatore. Se si usa invece il processo iterativo di Newton (adottato da quasi tutte le calcolatrici tascabili per il calcolo delle radici quadrate), ovvero xn=xn−(x2−1n−2pxn+1)/(2xn−2p)≡≡fn(p), n≥1, si ottiene Kn=∣p∣/∣xn−p∣. Se x0 è scelto opportunamente allora xn, per n→∞, converge a una delle radici e si vede che K*=∣p∣/∣x±p∣=K. Si tratta dunque di un modello numerico stabile. Questo elementare esempio evidenzia come problemi ben condizionati possano essere affrontati con metodi numerici stabili o instabili.

Errore computazionale e modelli adattivi

Il modello matematico mira a legare fra loro quantità di interesse fisico attraverso relazioni matematiche, spesso semplificate rispetto alla complessità del problema originario (per le nostre notazioni faremo costante riferimento al diagramma riportato nella fig. 1). La soluzione u del modello matematico sarà di conseguenza diversa da quella uf del problema fisico in esame, e tale differenza em=uf−u renderà conto sia della significatività (e ragionevolezza) delle ipotesi semplificative adottate, sia dell'accuratezza con la quale le quantità fisiche reali sono state rappresentate dall'insieme dei dati del modello. L'errore em è pertanto intrinseco al processo. Il modello computazionale creerà un altro errore, indicato con ec, che esprime la distanza fra la soluzione effettivamente calcolata ūn e u. L'errore globale e, differenza fra la reale soluzione fisica uf e quella calcolata ūn, risulterà dalla sovrapposizione dei due errori em e ec.

L'obiettivo della modellistica computazionale è di garantire che l'errore ec sia piccolo e controllabile (affidabilità) con il minor costo computazionale possibile (efficienza). Per costo computazionale si intende la quantità di risorse (tempo di calcolo, occupazione di memoria) necessaria per calcolare la soluzione ūn. L'affidabilità è un requisito cruciale di un modello computazionale: l'analisi mira a trovare stime dell'errore ec (in funzione dei dati del problema e dei parametri della discretizzazione) che consentano di garantire che esso stia al di sotto di una certa soglia di precisione fissata a priori. A questo scopo, è utile ricorrere ad algoritmi adattivi che adottano una procedura di retroazione (o feedback) a partire dai risultati già ottenuti per modificare i parametri della discretizzazione e migliorare la qualità della soluzione. I metodi numerici che implementano una strategia adattiva per il controllo dell'errore sono detti modelli computazionali adattivi.

Vediamo come viene generato l'errore computazionale. Il modello numerico deve essere implementato su calcolatore con l'ausilio di un algoritmo opportuno. Si generano di conseguenza nuovi errori, alcuni intrinseci all'algoritmo, altri attribuibili all'aritmetica finita di cui fa uso ogni calcolatore per rappresentare i numeri reali ed eseguire operazioni algebriche elementari. Per esempio, se il modello matematico è rappresentato dal problema differenziale u′=G(u) e per la discretizzazione della derivata si usano rapporti incrementali all'indietro, si ottiene come modello numerico un sistema di equazioni algebriche (lineari o non, a seconda della linearità o meno di G rispetto all'incognita u). Nel caso lineare possiamo rappresentare il sistema nella forma generale Ax=b, con A una matrice n×n, b un termine noto, x un vettore incognito (avendo indicato con n l'insieme dei nodi in cui si cerca la soluzione numerica). Anche nel caso di problemi differenziali non lineari, si perviene a un sistema lineare (previa opportuna linearizzazione) utilizzando tecniche à la Newton; si ottengono sistemi lineari in innumerevoli altre circostanze, in particolare nell'approssimazione numerica di problemi differenziali alle derivate parziali.

Nell'ipotesi in cui A sia non singolare, la soluzione del sistema lineare x=A−1 è nota in forma esplicita attraverso la formula di Cramer, che peraltro richiede un numero di operazioni elementari superiore a n!. Tale numero è proibitivo anche se n non è troppo grande (nelle applicazioni ai problemi differenziali, n tipicamente è dell'ordine di 10p per 2≤p≤6), anche per calcolatori, come gli attuali, in grado di effettuare diverse centinaia di miliardi di operazioni al secondo. Algoritmi numerici efficienti, quali l'algoritmo di fattorizzazione di Gauss, o algoritmi iterativi in sottospazi di Krylov, multigriglia o multilivello, consentono di ridurre il costo sino alla soglia ottimale Nn, introducendo tuttavia un ulteriore errore. La combinazione di tale errore con gli errori di arrotondamento dell'aritmetica finita (quest'ultimo decresce esponenzialmente con −t, dove t è il numero di cifre significative della macchina) fornisce l'errore algoritmico ea=unūn. L'errore del modello computazionale (ec=uūn) è la somma dell'errore numerico (en=uun) e di quello algoritmico.

L'obiettivo ultimo dell'analisi dell'errore ec del modello computazionale è quello di dimostrare che esso tende a zero per n→∞. Più precisamente, si cercano stime del tipo

[2] formula

formula

per un'opportuna norma ∣∣∙∣∣ che dipenderà dal problema in esame. L'esponente p (ordine di infinitesimo dell'errore rispetto a 1/n) è l'ordine del modello numerico. Nell'ipotesi ottimale in cui l'errore algoritmico ea sia controllabile con la tolleranza desiderata, p è lo stesso ordine di infinitesimo con cui en tende a zero rispetto a 1/n. La tecnica dimostrativa che consente di ottenere stime dell'errore del tipo indicato è basata sull'analisi di 'buona posizione' del modello matematico e su quella di stabilità del modello numerico. Essa è detta analisi a priori, in quanto può operarsi prima di implementare effettivamente il modello numerico e dipende in maniera implicita dal numero di condizionamento matematico K e da quello numerico K*. La costante C(u,d) dipende dai dati del problema, nonché da u e dalle sue derivate (supposto che esistano) fino a un ordine opportuno, legato al valore di p. La sua quantificazione numerica è generalmente assai difficile. Di conseguenza, la stima sopra riportata è incompleta: consente di prevedere la tendenza di decadimento dell'errore per n crescenti, ma non di fornirne una quantificazione numerica accurata per n fissato.

È compito dell'analisi a posteriori completare questa lacuna. Essa, sulla base della soluzione ūn effettivamente calcolata in corrispondenza di un certo valore del parametro di discretizzazione n e del valore del residuo rn=F(ūn,d) (una misura indiretta dello scostamento di ūn da u), consente di stabilire se l'errore sia o meno al di sotto di una prestabilita tolleranza. Se ciò non è vero, un algoritmo adattivo indica come migliorare la discretizzazione (fornendo un criterio per la modifica di ūn) al fine di garantire una riduzione dell'errore. Nel caso elementare di un sistema lineare, se x−indica una soluzione calcolata (con un qualunque algoritmo numerico) e r=bAx− è il residuo a essa associato, si può facilmente verificare che ∥x−x−∥/∥x∥≤K(A)∥r∥/∥b∥, dove ∥∙∥ indica la norma euclidea e K(A) è il numero di condizionamento della matrice (se A è una matrice simmetrica con autovalori reali positivi K(A)=λmaxmin, dove λmax e λmin sono rispettivamente il massimo e il minimo autovalore di A); dunque l'errore computazionale (relativo, o percentuale) è maggiorabile con il residuo (relativo) a meno del fattore K(A) che, in questo contesto, viene detto fattore di stabilità. Per problemi più complessi lo spirito è il medesimo, pur essendo in generale più difficile determinare i fattori di stabilità. I modelli computazionali adattivi combinano pertanto l'analisi a priori con quella a posteriori al fine di ottenere stime quantitativamente accurate dell'errore computazionale.

I metodi numerici

Intrinseco al concetto di modello computazionale è quello di metodo numerico per la realizzazione della fase (c) del processo che stiamo esaminando (fig. 1). Le componenti essenziali di un metodo numerico si possono classificare in rapporto al loro scopo, che può essere quello di: (1) rappresentare in modo approssimato una funzione (o una successione finita di dati sperimentali); (2) approssimare il valore numerico di un integrale; (3) approssimare (localmente o globalmente) derivate di ordine arbitrario; (4) approssimare soluzioni di relazioni funzionali (lineari e non) attraverso successioni il cui termine generale sia facile da ottenere e la cui velocità di convergenza sia la più elevata possibile; (5) decomporre matrici di struttura generale nel prodotto di matrici più semplici, per esempio triangolari, allo scopo di facilitare la risoluzione di sistemi lineari di grandi dimensioni.

Approssimazione di funzioni

Se f:[a,b]→ℝ è una funzione reale nota, definita sull'intervallo [a,b] della retta reale, ci si pone spesso il problema di darne una rappresentazione approssimata, per esempio tramite un polinomio o un polinomio a tratti, ovvero una funzione continua che sia un polinomio su ogni sottointervallo [ak,bk] di una certa partizione di [a,b]. Le ragioni per le quali si cercano approssimanti polinomiali sono molteplici. Innanzitutto, i polinomi sono semplici da trattare: è immediato derivarli e integrarli, ed è facile, grazie all'algoritmo di Horner, valutarli in punti prestabiliti. È inoltre noto dal teorema di Weierstrass che ogni funzione continua può approssimarsi con ordine arbitrario di accuratezza con un polinomio di grado opportuno. Precisamente, per ogni fC([a,b]) e ε>0 esistono n=n(ε) e pnn tale che ∥f−pnε, ove qui e nel seguito si indica con ∥g∥ il massimo valore di ∣g(x)∣ nell'intervallo [a,b], g essendo una generica funzione continua. Abbiamo denotato con ℙn l'insieme dei polinomi algebrici di grado inferiore o uguale a n.

Il teorema di Weierstrass, tuttavia, non è costruttivo: esso non permette di conoscere n(ε), né di costruire praticamente pn. Esempi di risultati costruttivi sono invece forniti sia dal teorema di Taylor sia dalla teoria dei polinomi ortogonali. Il primo assicura che se fCn([a,b]), il polinomio

[3] formula

formula

di grado n è un'eccellente approssimazione di f e delle sue derivate in un intorno del punto x0∈(a,b), ma sfortunatamente l'accuratezza diminuisce se x si allontana da x0. Si parla in tal caso di approssimazione locale.

Esempi di approssimazione globale sono forniti dalle serie troncate di sviluppi di f rispetto a basi di polinomi ortogonali. La famiglia {φk, k=0,1,…} di polinomi ortogonali in [a,b] rispetto alla funzione peso w(x)>0 è definita come segue:

formula

se

formula

Se a=−1 e b=1, per w(x)=1 si ottengono i polinomi di Legendre e per w(x)=1/√(1−x2) si ottengono quelli di Cebycev. Su un intervallo arbitrario tali definizioni si estendono in modo ovvio. Ora,

formula

è il coefficiente k-esimo di Fourier (generalizzato) di f; se la serie

formula

converge, allora si ha

formula

Indicando dunque con ∥f∥ = =(f,f)1/2 la cosiddetta norma L2 di f in (a,b) e, per ogni n, con

formula

la serie troncata di f (si noti che fn∈ℙn), abbiamo

[4] formula

formula

o, equivalentemente,

[5] formula.

formula

La convergenza è tanto più rapida quanto più f è regolare, per n→∞ si ha che ∣fn∣ tende a zero come (1/n)p se fCp([a,b]) per un certo p≥1. Tuttavia, la costruzione dei coefficienti {fn, k=0,…,n} è assai problematica e deve essere fatta in modo approssimato ricorrendo a formule di integrazione numerica gaussiane, in modo da preservare l'ordine di convergenza della serie.

Un approccio più semplice consiste nel sostituire f con un polinomio ∏nf di grado n che interpoli f in n+1 nodi distinti x0,…,xn di [a,b], ovvero ∏nf(xj)=f(xj) dove j=0,…,n. Tale polinomio si può rappresentare come

formula

essendo l'unico polinomio di grado n tale che Lj(xi)=δij. L'errore che si genera è

[6] formula

formula

essendo ξ un punto opportuno di [a,b] e

[7] formula

formula

il cosiddetto polinomio nodale. Se i nodi sono equispaziati, si possono fornire esempi di funzioni per le quali la proprietà di convergenza uniforme

formula

non è verificata. A tal fine è possibile ricorrere a nodi opportunamente distribuiti, detti di tipo gaussiano: si considerano gli zeri del polinomio ortogonale φn+1 oppure le radici del polinomio (x−a)(x−b)φ′n(x). Con tali scelte si ottiene, oltre alla convergenza uniforme, una stima dell'errore del tipo ∥En≤Cn−p, nell'ipotesi che f∈Cp([a,b]) (qui e nel seguito, C denota una costante positiva indipendente da n). Un'altra strategia conduce alla cosiddetta interpolazione composita: si suddivide dapprima [a,b] in sottointervalli disgiunti, di ampiezza h, e in ognuno si approssima f con una funzione continua hf,r la cui restrizione su ogni intervallo è un polinomio interpolatore di grado basso, per esempio r=1,2. In tal caso l'errore Ehr generato verifica una disuguaglianza del tipo ∥Ehr (k)∥≤Chr+1−k con k=0,1, avendo supposto f sufficientemente regolare, per esempio f∈Cr+1([a,b]). L'idea dell'interpolazione composita si estende senza difficoltà al caso multidimensionale ed è alla base del processo di discretizzazione degli elementi finiti, nonché dell'approssimazione con funzioni spline, nel qual caso si cercano raccordi anche per la derivata prima agli estremi dell'intervallo di integrazione.

Approssimazione di integrali

Per l'approssimazione del valore

formula

un'idea naturale è quella di calcolare I(fn), dove fn è un'approssimazione polinomiale di f. Prendendo per esempio fn=nf, si otterrà

formula

(dove formula),

formula

ovvero una formula di integrazione numerica in cui i punti {xj} sono ancora detti nodi, mentre i coefficienti {αj} sono detti pesi. Nel caso di interpolazione con nodi equispaziati, si ottengono le cosiddette formule di Newton-Cotes. Otterremo invece formule di integrazione gaussiana (assai più precise nell'integrazione di funzioni polinomiali) nel caso in cui si usino nodi di interpolazione gaussiana. La valutazione dell'errore di integrazione I(f)−I(fn) viene ereditata dalla stima del corrispondente errore En. Una naturale evoluzione in questa direzione si basa sull'uso di polinomi di interpolazione composita ∏hf,r, nel qual caso l'integrale su [a,b] verrà approssimato da una somma di integrali approssimati su ogni sottointervallo.

Un approccio diverso dai precedenti è alla base dei metodi di tipo Monte Carlo, dove i nodi sono scelti in modo statistico in funzione dei valori assunti da variabili casuali aventi una distribuzione di probabilità nota. Nella pratica computazionale i nodi sono ottenuti tramite un generatore di numeri casuali, e le formule che ne derivano si dicono di integrazione numerica pseudo-casuali. L'errore di integrazione su regioni dello spazio ℝn decresce come 1/√N (N essendo il numero di nodi), indipendentemente dalla dimensione n.

Approssimazione delle derivate di una funzione

Per approssimare in opportuni nodi i valori della derivata di una funzione f, una via naturale è ricorrere alla definizione stessa di derivata e sostituire f′(xk) con opportuni rapporti incrementali nei nodi adiacenti. Nel caso di nodi equispaziati, indicando con h la distanza fra due nodi consecutivi e con fk il valore f(xk), si può approssimare f′(xk) con (fk+1−fk−1)/h (ma anche con (fk−fk−1)/h) a meno di infinitesimi di ordine superiore a h. Un'accuratezza di ordine superiore si ottiene con il rapporto incrementale centrato (fk+1−fk−1)/2h, come si vede applicando la formula di Taylor. Naturalmente, iterando i precedenti rapporti incrementali si ottengono formule di approssimazione per derivate di ordine superiore. Questo processo, di tipo locale, è alla base dell'approssimazione alle differenze finite di problemi di Cauchy. Un approccio globale si ottiene invece approssimando f′(x) su un intero intervallo (a,b) con la derivata esatta di un polinomio approssimante f in (a,b), ottenuto tramite uno dei metodi precedentemente illustrati. Utilizzando per esempio il polinomio ∏nf interpolante f in n+1 nodi di tipo gaussiano, si ottiene dall'analisi precedente che ∥f′−(∏nf)′∥Cn−(p−1), avendo ancora assunto fCp([a,b]). Questo processo è alla base dell'approssimazione di problemi ai limiti con metodi spettrali.

Iterazioni funzionali

Capita assai di frequente di dover trovare per via approssimata la soluzione di equazioni funzionali della forma F(x)=0, dove F è un sistema di equazioni, lineari o non. Una metodologia generale assai naturale consiste nel realizzare un processo iterativo della seguente forma: assegnato x0, per ogni k≥0 si ottiene xk +1

da xk risolvendo il problema Gk(xk+1−xk)=−F(xk), dove Gk è un'approssimazione opportuna della matrice jacobiana F′ ( xk) di F nel punto xk. Chiameremo la matrice Gk operatore di precondizionamento e la supporremo invertibile con inversa Hk. Nel caso in cui F sia una trasformazione affine F(x)=Ax−b, la matrice Gk costituirà un'opportuna approssimazione della matrice A, e il processo iterativo ricomprenderà i classici metodi iterativi dell'algebra lineare numerica (di Jacobi, di Gauss-Seidel, di Richardson, del gradiente, e così via). Indicato con ek=x−xk l'errore all'iterazione k-esima e con I la matrice identità, si ha ek +1=(I−HkA)ek e per la convergenza del processo iterativo si richiederà che ϱk=ϱ(I−HkA) sia minore di 1 (ϱk è il raggio spettrale della matrice I−HkA). Nel caso in cui F(x)=0 sia un'equazione (o un sistema di equazioni) non lineare, il processo descritto è il paradigma di numerosi metodi classici, quali il metodo delle corde, di Newton e del gradiente.

L'analisi dell'errore si può fare analogamente al caso lineare. Supponendo infatti F differenziabile in un intorno della soluzione, si avrà ek+1=(I−HkF′ (yk))ek, dove yk è compreso fra x e xk. Il processo iterativo precedente comporta a ogni iterazione la risoluzione di un sistema lineare nella matrice Gk. Per la sua efficienza si richiede un compromesso fra la semplicità computazionale e la velocità di convergenza: infatti, Gk dovrebbe essere al tempo stesso di forma semplice (diagonale, triangolare, a banda, ecc.) e tale che lo spettro di HkF′(yk) sia il più compatto possibile, in modo da garantire che ϱk (o un'opportuna norma di I−HkF′(yk) nel caso non lineare), sia più piccolo di 1 e più vicino possibile a 0.

Decomposizione di matrici

Per una matrice quadrata A non singolare di n righe si possono in generale trovare due fattori: una matrice L triangolare inferiore (ovvero con elementi nulli al di sopra della diagonale principale) e una U triangolare superiore (con elementi nulli al di sotto della diagonale principale), tali che A=LU. Questa decomposizione (che esiste a meno di possibili permutazioni di righe o colonne di A) è unica se, per esempio, si fissano uguali a 1 gli elementi diagonali di L (ottenendo in tal modo la cosiddetta decomposizione di Gauss) o quelli di U (decomposizione di Crout). Note le matrici L e U, la risoluzione di un sistema Ax=b può ricondursi a quella di due sistemi triangolari Ly=b e Ux=y (in sequenza); ognuno di essi può essere risolto con il metodo delle sostituzioni in avanti per L, all'indietro per U, con un costo computazionale di n2 operazioni. Nel caso in cui A sia simmetrica con autovalori positivi, essa ammette la decomposizione (detta di Cholesky) A=HHt, dove H è triangolare inferiore e Ht ne è la trasposta, che naturalmente risulta triangolare superiore.

Conviene inoltre citare la cosiddetta decomposizione QR per la risoluzione di sistemi lineari sovradeterminati Ax=b, dove b è un vettore di m componenti, x un vettore di n componenti e A una matrice di m righe ed n colonne, con m>n. Non esistendo una soluzione in senso classico, se ne può cercare una x*, detta soluzione nel senso dei minimi quadrati, risolvendo il sistema n×n: At(Ax*−b)=0. Il classico problema della determinazione della retta di regressione lineare, che minimizza la somma dei quadrati delle distanze da una distribuzione di osservate {bi,i=1,…,m}, si formula esattamente in questo modo (in questo caso n=2 è il numero di coefficienti in grado di individuare la retta di regressione). Siano Q,R due matrici rispettivamente m×n ed n×n, l'una ortogonale e l'altra triangolare superiore, per le quali si abbia A=QR. Indicata con D la matrice diagonale QtQ, il sistema trasformato si può riscrivere nel modo seguente: RtQt(QRx*−b)=0, ovvero Rt(DRx*−Qtb)=0. Si è dunque condotti alla risoluzione in sequenza di due sistemi di struttura elementare: uno diagonale, Dy=Qtb, l'altro triangolare, Rx*=y.

Modellistica computazionale per problemi differenziali

La complessità dei problemi differenziali, ordinari e alle derivate parziali, e la loro rilevanza nella modellistica fisico-matematica hanno stimolato negli ultimi decenni un enorme sviluppo di metodi numerici e computazionali adatti alla loro risoluzione.

Equazioni differenziali ordinarie

Consideriamo l'equazione differenziale ordinaria al prim'ordine y′(x)=f (x,y (x)), per x≥x0, con un dato iniziale y(x0)=y0. Supporremo f lipschitziana rispetto all'argomento y in un intorno di (x0,y0), ovvero che esista una costante L>0 tale che ∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣ per ogni y1,y2. Ciò assicura che esiste una soluzione del problema, quantomeno in un intorno I=[x0,b] del punto iniziale x0. Consideriamone un'approssimazione basata sul metodo delle differenze finite. Passo preliminare è la suddivisione di I in n sottointervalli di ampiezza h, i cui estremi xj=x0+jh, con j=0,…,n, sono detti nodi. Il modello numerico genererà dei valori {uj} approssimanti la soluzione esatta y nei nodi {xj}. Una famiglia di metodi numerici è basata sull'idea di soddisfare l'equazione differenziale solo nei nodi, dopo aver sostituito la derivata prima con una sua conveniente approssimazione che faccia uso dei valori nodali. Come approssimante si può scegliere sia un rapporto incrementale sia la derivata prima del polinomio interpolatore (in quest'ultimo caso i nodi dovranno essere non più equispaziati, ma di tipo gaussiano e h denoterà la massima distanza fra due nodi consecutivi). Un'altra strategia si basa sull'uso della formula fondamentale del calcolo integrale

formula

,

dove l'integrale stesso è sostituito con un'opportuna formula interpolatoria in alcuni dei nodi {x0,x1,…,xj +1}.

Astrattamente, in entrambi i casi il problema numerico assumerà la forma seguente:

formula

,

dove i coefficienti {αk} e {βk} (non tutti necessariamente diversi da zero) caratterizzano lo schema. Possiamo ora esaminare il significato dei concetti generali di convergenza, consistenza e stabilità, precedentemente enunciati, in questo contesto. Ponendo yj=y(xj) per ogni j=0,…,n, lo schema è detto consistente se, posto

formato

,

si ha che τ(h)=maxj∣τj (h )∣→0 per h→0. Esempi significativi sono forniti dal metodo di Euler in avanti (EA), uj +1−uj=hf (xj, uj), e dal metodo dei trapezi o di Crank-Nicolson (CN), uj+1−uj=(h/2)[f (xj,uj)f (xj+1,uj+1)]. Per l'analisi dell'errore di discretizzazione locale τj (h) si ricorre al teorema di Taylor: poiché f(xj,yj)=y′ (xj), si ottiene che τ(h)=O(h) per EA, mentre τ(h)=O(h2) per CN. Diremo che il problema numerico è convergente, se la soluzione numerica {uj} converge a quella esatta nei nodi, ovvero se ∣y(xj)−uj∣≤C(h) per ogni j, C(h) essendo un infinitesimo quando h→0. Nel caso in cui un metodo converga, generalmente si ha che C(h)=Chp, dove p è l'ordine di infinitesimo di τ(h).

Pertanto, l'analisi dell'ordine di accuratezza di un metodo numerico si può ricondurre all'analisi dell'errore di troncamento locale. Sappiamo che la consistenza non è sufficiente a garantire la convergenza, richiedendosi anche la stabilità. Quest'ultima (detta anche zero-stabilità) si traduce nella richiesta seguente: per ogni b>x0, h>0, sia N(h) il numero di nodi in [x0,b] e cioè xj=x0+jh, con xN≤b. Se {uj, j=1,…,N(h)} e {vj, j=1,…,N(h)} sono le successioni generate dal metodo numerico a partire da due diversi valori iniziali differenti, u0 e v0, allora esiste una costante K, dipendente eventualmente da b, taleche per ogni h sufficientemente piccolo si abbia ∣uj−vj∣≤K∣u0−v0∣ per ogni j=1,…,N(h). A titolo esemplificativo, consideriamo il metodo EA. Vale l'identità yj−uj=(yj−y*j)+(y*j−uj), dove yj* è definito dalla relazione yj*−yj=hf(xj,−1, yj−1) ovvero è la soluzione di EA supponendo di conoscere in xj−1 la soluzione esatta. Dalla definizione di τj(h) discende che yj*−yj−1=hτj(h), grandezza che esprime l'errore introdotto al passo j-esimo. Evidentemente la quantità globale di errore introdot-ta in [x0,b] non potrà superare N(h)[hτ(h)], ovvero (b−x0)τ(h). Peraltro, yj*−uj=h[f(xj−1,yj−1)−f(xj,−1,uj−1)] e, essendo f lipschitziana con costante L rispetto al secondo argomento, si avrà ∣yj*−uj∣≤hL∣yj−1−uj−1∣ (una relazione che può considerarsi inerente alla stabilità del metodo, in quanto esprime la variazione della soluzione numerica in rapporto a quella sui dati al generico passo j-esimo). Posto allora ej=∣yj−uj∣, si ha la relazione ricorsiva ej〈hτ(h)+hLej−1 per ogni j, la quale consente di concludere che ej≤Cτ(h) per ogni j=1,…,N(h) con C=[exp(b)−1]/L, evidenziando come anticipato il contributo sinergico di stabilità e consistenza di un metodo numerico.

Il limite intrinseco del precedente concetto di zero-stabilità sta nella possibile dipendenza esponenziale della costante K da b; questa dipendenza rende tale proprietà inadeguata a governare comportamenti su intervalli di integrazione grandi o addirittura illimitati. In tal caso sarà opportuno richiedere che il metodo sia assolutamente stabile, ovvero che la soluzione {uj} del problema modello y′(x)=λy(x) per x>0, con y(0)=1 e Reλ〈0, tenda a zero per xj→∞. Nel caso di EA si ha uj+1=(1+hλ)uj e pertanto uj=(1+hλ)j; la condizione di stabilità assoluta diviene ∣1+hλ∣〈1, ovvero hλ deve appartenere al cerchio di raggio unitario e centro (−1,0) nel piano complesso. Ciò induce una limitazione su h del tipo h〈2∣Reλ∣/∣λ∣2 e il metodo EA verrà detto condizionatamente assolutamente stabile. Al contrario, il metodo CN soddisfa la relazione uj=[(1+hλ)/(1−hλ)]j per ogni j, e pertanto uj→0 per xj→+∞ incondizionatamente rispetto ad h. I metodi che godono di questa proprietà vengono detti A-stabili e sono appropriati per integrare equazioni differenziali su intervalli illimitati e anche sistemi differenziali di tipo stiff, ovvero sistemi la cui soluzione ha componenti che decadono nel tempo con dinamiche molto diverse.

Equazioni alle derivate parziali

La modellistica numerica e computazionale è di interesse primario in questo settore, la cui vastità non consente una sintesi soddisfacente in poche righe. Faremo un cenno alle più comuni metodologie di approssimazione. Se Ω⊂ℝd è una regione limitata del piano (d=2) o dello spazio (d=3), molti problemi fisici si possono modellare con il seguente problema matematico: trovare una funzione u (scalare o vettoriale) dipendente dal tempo e dallo spazio tale che ogni x=(x1,..., xd)∈Ω e t ›0 valga

[8] formula

formula

con u=u0 assegnata al tempo iniziale t=0. Il simbolo ∂/∂t esprime la derivata parziale rispetto alla variabile t (con x fissato), mentre per ogni funzione vettoriale w=w(x,t) la divergenza è definita da

formato

(dove ∂/∂xi è la derivata parziale lungo la direzione xi). L'identità precedente esprime una legge di conservazione (o un sistema di leggi di conservazione nel caso vettoriale). In effetti, integrando la [8] su ogni sottoregione T di Ω, grazie al teorema della divergenza di Gauss si ottiene:

[9] formula

formula

avendo indicato con Γ il bordo di T e, per ogni x∈Γ, con n=n(x) il versore di componenti ni esterno a Γ (ovvero il vettore di modulo unitario perpendicolare al piano tangente a Γ nel punto x). Infine, f(u) è una funzione vettoriale, a d componenti f1(u),…, fd(u), detta flusso di u, e

formula

La relazione precedente esprime la conservazione di u nella sottoregione T; questo significa che la variazione temporale del suo integrale ∫Tudx è compensata dal bilancio di flusso attraverso la frontiera.

Sia {T} una partizione di Ω di elementi poligonali (triangoli o tetraedri) non sovrapponentisi. Più precisamente, due elementi distinti di {T} possono avere in comune solo vertici, lati e facce, e la loro unione ricopre interamente Ω. Il metodo ai volumi finiti è basato sulla proprietà di conservazione, imposta su ogni sottoregione T (detta volume finito). Per ogni coppia di volumi finiti adiacenti Ti e Tj, detta Γij l'interfaccia fra Ti e Tj, i flussi normali su Γij sono deputati a stabilire le relazioni di interdipendenza fra i valori delle incognite ui=u∣Ti e uj=u∣Tj (generalmente supponendo ui costante su Ti, per ogni i). Tipicamente avremo

[10] formula

formula

dove la somma è estesa a tutti gli indici j per i quali Tj è adiacente a Ti, mentre fap(ui,uj) è una conveniente approssimazione del flusso normale f(u)∙n su Γij.

Il metodo delle differenze finite si basa invece sulla sostituzione delle derivate spaziali ∂/∂xi in [8] con opportuni rapporti incrementali in nodi allineati lungo la i-esima direzione coordinata (analogamente a quanto fatto precedentemente per approssimare altre derivate). Infine, il metodo degli elementi finiti viene derivato dalla forma [8], supponendo che u sia approssimata con una funzione uh la cui restrizione a ogni T (stavolta detto elemento finito) sia un polinomio di grado 1 (o più elevato), continuo sulle interfacce e soluzione del problema

[11] formula

formula

dove vh è una funzione (della sola x) polinomiale su ogni T dello stesso grado di uh, fh(uh) un'opportuna approssimazione di f(uh) e g esprime l'eventuale contributo del flusso normale (qualora venga assegnato) sulla frontiera fisica ∂Ω della regione Ω. Per flussi fh(uh) di forma complessa, gli integrali richiederebbero di essere approssimati numericamente. Nelle notazioni precedenti, h indica la massima lunghezza dei lati dei triangoli T, così che uh→u per h→0. L'approssimazione agli elementi finiti si ispira alla cosiddetta formulazione debole del problema [8], formulazione che (nel caso di problemi stazionari e simmetrici) esprime a livello discreto il principio dei lavori virtuali: le vh sono quindi le approssimanti dei cosiddetti spostamenti virtuali. Risulta evidente che, per l'effettiva risoluzione dei problemi sopra riportati, converrà discretizzare anche la derivata temporale, utilizzando per esempio schemi alle differenze finite (esplicite o implicite) come visto sopra nel trattamento di equazioni differenziali ordinarie del primo ordine, oppure lo stesso metodo agli elementi finiti. Nel caso implicito, l'equazione [11] genererà per ogni livello temporale un sistema di equazioni per la cui soluzione si possono usare le tecniche numeriche illustrate prima.

Strategie multi-livello

Come si è visto in molte situazioni, per l'approssimazione di problemi differenziali (siano essi ordinari o alle derivate parziali) si è ricondotti alla soluzione di sistemi lineari della forma Ahxh=fh, dove xh indica il vettore delle incognite nodali di una griglia di calcolo la cui ampiezza caratteristica è h, mentre fh e Ah indicano rispettivamente il termine noto e la matrice associati al sistema. Se si usa un metodo iterativo (quale il metodo di Jacobi o quello di Gauss-Seidel), si osserva che poche iterazioni sono sufficienti per abbattere le alte frequenze dell'errore ek=xkh−xh, mentre per la riduzione delle basse frequenze servirebbero moltissime iterazioni (xkh indica l'approssimazione di xh dopo k  iterazioni). Allo scopo di accelerare la convergenza si può ricorrere a un problema ausiliario AHxH=fH, che possiamo immaginare associato alla discretizzazione dello stesso problema, ma su una griglia più rada di ampiezza caratteristica H>h.

A ogni vettore vH, composto da valori nodali sulla griglia rada, possiamo associare un vettore vh=IhvH composto da valori nodali sulla griglia fine originaria. Tipicamente, Ih è un processo di interpolazione rappresentabile con una matrice rettangolare. Indichiamo con Ith la matrice trasposta di Ih, che rappresenta dunque il processo di restrizione dalla griglia fine alla griglia rada. Indichiamo inoltre con QH la matrice la cui inversa è

formula

e con Qh la matrice Ah o un suo conveniente precondizionatore. Per la risoluzione del sistema originario si può usare il metodo iterativo che illustriamo di seguito. Se ukh è nota (per k≥0), si definisca uhk tramite la soluzione del sistema lineare (sulla griglia rada)

[12] formula

formula

e quindi si definisca uhk come soluzione del sistema (sulla griglia fine)

[13] formula.

formula

Combinando i due processi si trova che

[14] formula,

formula

ovvero uhk è il risultato di una singola iterazione di Richardson sul problema originario, dove si introduce il precondizionatore P tale che P−1=QH−1+Qh−1−Qh−1AhQH−1. Si ritrova in questo modo il metodo multigriglia nella sua versione più semplice (a due livelli). I metodi a più livelli sono ottimali dal punto di vista della complessità computazionale, in quanto per la risoluzione del sistema originario richiedono un numero di operazioni proporzionale al numero di incognite (ovvero al numero di righe della matrice Ah). Una versione ridotta del metodo multilivello più adatta a calcolatori di tipo multi-processore si ricava sostituendo P con Ppar, dove Ppar−1=QH−1+Qh−1 si ottiene trascurando l'interazione non lineare fra QH−1 e Qh−1. I cosiddetti metodi di decomposizione per sottodomini (o sottoregioni), introdotti a metà degli anni Ottanta del Novecento per la risoluzione parallela di problemi differenziali, possono essere formulati come esempi notevoli di tali procedure.

Bibliografia

Atkinson 1989: Atkinson, Kendall E., An introduction to numerical analysis, 2. ed., New York, Wiley, 1989.

Eriksson 1996: Eriksson, Kenneth e altri, Computational differential equations, Cambridge-New York, Cambridge University Press, 1996.

Golub, Van Loan 1989: Golub, Gene H. - van Loan, Charles F., Matrix-computations, 2. ed., Baltimore (Md.), Johns Hopkins University Press, 1989.

Lambert 1991: Lambert, John D., Numerical methods for ordinary differential systems. Initial value problem, Chichester, Wiley, 1991.

Quarteroni 2002: Quarteroni, Alfio - Sacco, Riccardo - Saleri, Fausto, Matematica numerica, 2. ed., Milano, Springer, 2002.

Quarteroni 2006: Quarteroni, Alfio, Modellistica numerica per problemi differenziali, 3. ed., Milano, Springer, 2006.

Quarteroni, Saleri 2006: Quarteroni, Alfio - Saleri, Fausto, Introduzione al calcolo scientifico, 3. ed., Milano, Springer, 2006.

Quarteroni, Valli 1994: Quarteroni, Alfio - Valli, Alberto, Numerical approximation of partial differential equations, Berlin-New York, Springer, 1994.

Quarteroni, Valli 1999: Quarteroni, Alfio - Valli, Alberto, Domain decomposition methods for partial differential equations, Oxford-New York, Clarendon, 1999.

Saad 1996: Saad, Yousef, Iterative methods for sparse linear systems, Boston (Mass.)-London, PWS, 1996.

Smith 1996: Smith, Barry - Bjørstad, Petter - Gropp, William, Domain decomposition, Cambridge-New York, Cambridge University Press, 1996.

Van der Vorst 2003: van der Vorst, Henk A., Iterative Krylov methods for large linear systems, Cambridge, Cambridge University Press, 2003.

CATEGORIE
TAG

Formula fondamentale del calcolo integrale

Equazione differenziale ordinaria

Equazioni alle derivate parziali

Sistema di equazioni, lineari

Metodo degli elementi finiti