NUMERICI, CALCOLI

Enciclopedia Italiana - IV Appendice (1979)

NUMERICI, CALCOLI (XXV, p. 29; App. III, 11, p. 286)

Enzo Aparo

Introduzione. - La nozione di c. n. si può introdurre, facendo riferimento al termine latino calculus (piccola pietra, pedina), nel modo seguente. Un insieme finito di oggetti, un insieme finito di celle, capaci ciascuna di contenere un oggetto, e un insieme finito di istruzioni, seguendo le quali si perviene in un numero finito di passi da una configurazione iniziale di oggetti nelle celle a una finale, definiscono un calcolo. Si ha un c. n. digitale quando i vari raggruppamenti di oggetti vengono interpretati come numeri, in un certo sistema di numerazione.

Quanto si è detto fa evidente riferimento al calcolatore meccanico o elettronico digitale (più precisamente, alla nozione di macchina di Turing, v. informatica: Informatica matematica, in questa Appendice). Strumenti quali il regolo calcolatore, i planimetri, i calcolatori di tipo analogico trattano invece i numeri come grandezze fisiche continue, e un operatore che esegue due volte la stessa successione di operazioni sugli stessi dati di partenza trova, di regola, risultati molto vicini fra loro ma diversi. Si ha in tal caso un c. n. analogico (in particolare, un calcolo grafico, se le operazioni svolte consistono per la maggior parte nel tracciamento di grafici e nella lettura delle coordinate di certi punti d'intersezione).

Va menzionato che la denominazione di "c. n. o anche c. n. e grafici" tende a scomparire a favore di "matematica numerica" o meglio di "analisi numerica". Viene considerata oggi sotto questo nome quella parte della matematica che si occupa dell'impiego di c. n. nella risoluzione dì problemi relativi alle scienze matematiche, fisiche e naturali e all'ingegneria.

Metodi numerici che forniscono risposte in termini quantitativi a questioni di economia, sociologia, scienze politiche, sono invece per lo più considerati come pertinenti all'analisi dei sistemi o alla ricerca operativa. Diamo ora qualche notizia di problemi e metodi dell'analisi numerica.

Interpolazione.

1) Interpolazione lineare.- Sia E uno spazio vettoriale di dimensione n, ed F1, F2, ..., Fn, funzionali lineari in E* (il duale di E). Dati gli scalari y1, y2, ..., yn si vuole trovare un n ∈ E tale da aversi Fi(n) = yi, i ∈ {1, ..., n}. Si dimostra che il problema possiede una e una sola soluzione se e solo se i funzionali Fi sono linearmente indipendenti. Ne sono casi particolari:

a) Interpolazione polinomiale semplice. - E sia l'insieme formato dai polinomi di grado ≤ n − 1, e da quello identicamente nullo, p(a, x) = a0 + a1 x + .... + an-1 xn-1. Esso è uno spazio vettoriale a n dimensioni e i polinomi 1, x, ... xn-1 ne costituiscono una base. Prendendo Fi [p(a, x)] = p(a, xi) per ogni polinomio di E, si vogliono trovare le ai in modo che p(a, xi) = yi, i ∈ {1, ..., n}. I metodi numerici più in uso sono quelli di Newton, Lagrange e Neville.

Nel metodo di Newton s'introducono le differenze divise:

rispettivamente d'ordine 2, 3, 4, ... relative alla funzione f. vale l'identità:

Se si sopprime al 2° membro l'ultimo termine, si ha il polinomio di grado n − 1 al più, che in xi acquista il valore yi = f(xi), i ∈ {1, ..., n}.

Il metodo di Lagrange fa uso dei polinomi

vale zero in tutti gli xj diversi da xi e vale 1 per x = xi). Il polinomio

è quello cercato.

Nel metodo di Neville si considerano successivamente i polinomi:

sino ad arrivare al polinomio

b) Interpolazione polinomiale di Hermite. - Supposto che f(x) abbia le derivate prima, seconda, ..., (qi − 1)-esima nel punto xi (i ∈ {1, ..., n}, si sostituisce a f(x) quel polinomio p(x) di grado ≤ q1 + q2 + ... + qn − 1 (eventualmente ≡ 0) per il quale si ha p(j)(xi) = p(j)(xi) (j = 1, 2, ..., ni − 1; i = 1, 2, ..., n).

Posto pii(x) per indicare il polinomio di 1° grado che in xi vale f (xi) e ha derivata uguale a f′(xi), piii(x) quello di 2° grado che in xi vale f (xi) e ha derivate prima e seconda uguali a f′(xi) e f″(xi), ..., si può ancora usare il metodo di Neville per ottenere, per es., p001(x) a partire da p00(x) e p01(x); p0011(x), noti p001(x) e p011(x), e così via. I polinomi con indici tutti uguali si scrivono immediatamente mediante la formvula di Taylor.

c) Interpolazione mediante splines. - E stata sviluppata negli ultimi due decenni e la sua applicazione è sempre più diffusa.

Siano x0, x1, ..., xn, con x0 = a, xn = b, xi-1 xi (i = 1, ..., n), n + 1 punti base per l'interpolazione di f (x) in [a, b]. Una funzione s(x) è una spline cubica relativa a f e al sistema di detti punti se è di classe C2 [a, b] (possiede cioè derivate prima e seconda continue in [a, b]), vale f (xi) in xi, ed è un polinomio di grado ≤ 3 in ogni [xi-1, xi] (i = 1, ..., n).

Se ne trovano i coefficienti risolvendo sistemi di equazioni lineari algebriche aventi matrice tridiagonale (cioè con gli elementi, al di fuori della diagonale principale e delle due parallele subito al di sopra e al di sotto, tutti nulli), la cui risoluzione è particolarmente semplice. Una caratteristica notevole delle splines è la loro convergenza a f allorché la massima delle lunghezze degl'intervalli in cui [a, b] è decomposto tende a zero, ciò che non si verifica nel caso dell'interpolazione con un singolo polinomio. Splines di ordine 23 possono definirsi in maniera analoga.

d) Interpolazione trigonometrica. - Dato l'insieme dei punti base xk = 2πn/k (k = 0, 1, ..., n − 1) si cerca il polinomio trigonometrico

βk exp (ikx) (βk complesso) tale che p(xk) = fk (k = 0, 1, ..., n − 1).

Si dimostra che l'unico polinomio trigonometrico che soddisfa le condizioni richieste è quello con

Per calcolare le βk è particolarmente adatto un procedimento recentemente introdotto (1965) da J. W. Cooley e J. W. Tuckey col nome di fast Fourier transform.

2) Interpolazione non lineare. - Interpolazione con funzioni razionali. - È nota f (x) in μ + ν + 1 punti xi, a due a due distinti, e si vuole determinare la funzione razionale

che nei punti xi(i = 0, 1, ..., μ + ν) acquista i valori f(xi).

A differenza dell'interpolazione con polinomi, possono presentarsi casi degeneri nei quali il problema non ha soluzione. Il metodo delle differenze inverse si presta a fornirne, nei casi non degeneri, sotto forma di frazione continua. Prendiamo μ = ν = n e si ponga fi = f (xi), (i = 0, 1, ..., 2n);

Risulta

dove, per es., l'espressione a b + c d va intesa come a/[b + (c/d)]. Allo stesso modo in cui questa formula è analoga all'identità di Newton menzionata in 1a, è possibile sviluppare un procedimento simile a quello di Neville.

Sistemi di equazioni lineari algebriche.

1) Data una matrice A con m righe ed n colonne a elementi reali o complessi e una matrice B con m righe e 1 colonna (vettore-colonna) si cerca un vettore-colonna x a n componenti tale da aversi Ax = B, cioè

simultaneamente. Si dimostra che il problema può avere infinite soluzioni o anche esserne sprovvisto. Quando ha soluzioni, è sempre possibile ridursi al caso di matrici quadrate e perciò considereremo nel seguito solo questo caso.

2) Metodo di Gauss. - È il più usato fra i metodi diretti. Si passa a un sistema equivalente con matrice triangolare superiore U, mediante le operazioni seguenti:

a) scambio fra loro di due equazioni, b) moltiplicazione dei coefficienti e del termine noto di un'equazione per uno stesso numero, c) sottrazione del primo membro di un'equazione dal primo membro di un'altra equazione e del termine noto dell'una dal termine noto dell'altra.

Nel corso di tali trasformazioni può essere necessario scambiare di posto certe incognite. Il sistema Ux - C, al quale si perviene, possiede una e una sola soluzione, facile a ottenersi ricavando xn dall'ultima equazione, sostituendola nella penultima e calcolando xn-1 e così via, se tutti gli elementi uii sulla diagonale principale di U sono diversi da 0, ciò che certamente avviene se il sistema di partenza ha unica soluzione.

Si dimostra che, se A è una matrice quadrata dotata di inversa, le operazioni descritte consentono di scrivere A come prodotto PLU, con L triangolare inferiore, U triangolare superiore, con diagonale principale a elementi ≠ 0, e P matrice permutazione (una matrice con un solo 1 su ogni riga e colonna e zeri altrove). Il sistema PLUx = B viene risolto considerando prima il sistema Ly = PB e poi il sistema Ux = y (si tenga presente che P2 = E, matrice unità). Due modi di organizzare i calcoli recano i nomi di P. D. Crout e di T. Banachiewicz.

Se A è definita positiva, tale cioè che A* = A (l'asterlsco indica l'operazione di trasposizione preceduta o seguita dal passaggio alla matrice coniugata) e x*Ax > 0, per ogni x ≠ 0, sussiste la decomposizione A = LL*, con L triangolare inferiore ed lkk > 0 (h = 1, 2, ..., n). Il calcolo di L, una riga alla volta, costituisce il metodo di Cholesky.

3) Metodo di ortogonalizzazione. - un metodo diretto, applicabile anche nel caso in cui la matrice A sia m × n con m > n e di rango n, se si tratta come "soluzione" del sistema (in generale, incompatibile) Ax = B quel vettore x0 per cui (Ax − B)*(Ax − B) è minimo.

Si trovano una matrice triangolare superiore U con uii ≠ o tale che AU = C con C*C = D con D diagonale. Risulta x0 = UD-1 C*B. Alla determinazione di U e C provvedono i metodi di E. Schmidt o di A. S. Householder.

4) Metodi di Jacobi e di Gauss-Seidel. - Sono i due metodi iterativi più classici. Si decompone A nella somma L + D + U, con D diagonale, L triangolare inferiore e U triangolare superiore.

Nell'ipotesi che D-1 esista, il metodo di Jacobi considera la successione dei vettori x(1), x(2), ... ottenuta, a partire da un arbitrario x(0), tramite lo schema iterativo x(k+1) = D-1 [B − (L + U)]x(k).

Due altri ben noti procedimenti iterativi si basano sulle formule (L + D)x(k+1) + Ux(k) = B, oppure (D + U)x(k+1) + Lx(k) = B (Gauss-Seidel), nell'ipotesi che (L + D)-1, oppure (D + U)-1 esistano. È facile vedere che questi metodi consistono nello scrivere il sistema Ax = B nella forma Cx + (A Cx) = B con C arbitraria matrice non singolare e nell'introdurre lo schema Cx(k+1) + (A C)x(k) = B. Si dimostra che la successione x(k) converge alla soluzione del sistema dato se e solo se gli autovalori della matrice E C-1A (E unità) sono tutti in modulo 〈 1.

5) Metodi Di rilassamento. - Sono una generalizzazione dei metodi iterativi ora detti. Si decompone D nella somma (i/ω) D + [i − (1/ω)] D, e si prende ω in modo tale che il procedimento basato sull'iterazione [L + (i/ω) D]x(k+1) + {[1 − (1/ω)] D + U}x(k)) = B converga il più rapidamente possibile. Si parla di "sovrarilassamento" o "sottorilassamento", rispettivamente, secondo che ω sia > 1 oppure 〈 1. L'impiego di tali metodi è importante specialmente nella risoluzione di equazioni differenziali a derivate parziali, la cui particolare struttura suggerisce volta per volta quale sia il più conveniente valore di ω.

Equazioni e sistemi di equazioni non-lineari.

1) Metodo di Newton. - Sia f : Rn → Rn una funzione vettoriale a n componenti delle n variabili x1, ..., xn. Si vuole trovare −x tale che f (−x) = 0, oppure, più per esteso, tale che fi (ù1, ..., ùn) = 0 (i = 1, 2, ..., n).

Supposto che le derivate parziali. ∂fi/∂xj (i, j = 1, ..., n) esistano in un certo aperto convesso C0 ⊆ Rn, indicheremo con f′ (x) la matrice (jacobiana) che ha ∂fi/∂xj nel posto (i, j). Introdotta una norma ∥•∥ per i vettori e le matrici in oggetto, si può dimostrare quanto segue. Ipotesi: per un certo x0 C0 e per tre certe costanti α, β, γ risulti: a) ∥ [f′(x0)]-1f(x0) ∥ ≤ α; b) ∥ [f′(x)]-1 ∥ ≤ β per ogni x C0; c) ∥ f′(x) − f′(y) ∥ ≤ γ ∣ x − y ∣ per ogni x C0 e y C0; d) h = ½ αβγ 〈 1. Tesi: posto r = α(i h)-1 e Dr(x0) = {x ∣ ∥ x − x0 ∥ 〈 r} la successione definita da x(k+1) = x(k) − [f′ (x)(k)]-1 f (x(k)) (k = o, 1, ...) ha tutti i suoi elementi in Dr(x0), e converge a un ben determinato ξ ∈ Ãr(x0) per il quale f (ξ) = 0.

Risulta, per ogni k, ∥ x(k) − ξ ∥ ≤ αh2k-1 (1 − h2k-1)-1, da cui segue che la convergenza è almeno quadratica.

Il procedimento descritto non è altro che una linearizzazione della f (x) nelle vicinanze dei punti x(0), x(1), .... Dalla formula di Taylor si ha infatti f (x) ≅ f (x0) + f′ (x0) (x x0) (qui x x0 va inteso come vettore-colonna), e annullando il 2° membro si ottiene x(1).

2) Autovalori di una matrice. - Data una matrice A, n × n, reale o complessa, si definiscono suoi autovalori quei numeri λ reali o complessi tali che p(λ) = det (λE A) = 0. Il polinomio p(λ) si dice "polinomio caratteristico" di A. Il calcolo degli zeri di p(A) può essere effettuato direttamente, senza trovarne i coefficienti, mediante vari metodi, il più antico dei quali è quello di Jacobi, valido per matrici reali simmetriche.

Si indichi con U una matrice n × n ottenuta sostituendo nella matrice unità cos θ, sen θ ai due elementi di posto (r, r), (s, s), e − sen θ, cos θ ai due elementi di posto (r, s), (s, r), con r s. Una scelta opportuna di θ fa sì che nella matrice UTAU risultino nulli i due elementi (r, s) ed (s, r). Mediante una successione di moltiplicazioni a sinistra e a destra per matrici del tipo ora descritto si perviene a una matrice B che ha gli stessi autovalori di A e tutti gli elementi fuori della diagonale principale minori, in modulo, di un prefissato numero ε. Gli elementi diagonali forniscono valori approssimati degli autovalori di A. Il metodo attualmente più raccomandato è dovuto a J. F. Francis, ispirato, a sua volta, al cosiddetto metodo LR di H. Rutishauser. Posto A - AI, si costruiscono, in successione, le matrici Qi, Ri tali che Ai = QiRi, Qi*Qi = E (l'asterisco indica trasposta coniugata), Ri triangolare superiore, Ai+1 = RiQi. Se gli autovalori di A sono reali e distinti, la successione Ai tende a una matrice triangolare superiore avente tali autovalori sulla sua diagonale principale. Nell'altro caso, si ha convergenza verso una matrice che differisce da una triangolare superiore per la presenza di sub-matrici quadrate i cui elementi diagonali sono sulla sua diagonale principale. Gli autovalori di A si ricavano dagli autovalori di tali sub-matrici.

3) Zeri reali o complessi di polinomi. - Un modo per calcolare le radici, reali o complesse, di un polinomio, è quello di prendere una matrice che ha tale polinomio come suo polinomio caratteristico (basta bordare opportunamente la matrice unità) e determinarne gli autovalori. Si tende oggi a raccomandarlo perché poco sensibile a perturbazioni sui coefficienti. Nel caso di radici reali, il metodo più semplice è quello di bisezione, consisteme nel dimezzare via via intervalli che sicuramente contengono una certa radice, sino a ottenerne uno di ampiezza arbitrariamente piccola e il cui centro fornisce in pratica la radice cercata. Il procedimento è stato esteso dallo scrivente al calcolo di radici complesse. Per queste ultime, è molto diffusa una variante del metodo di Newton, nota come metodo di Bairstow, che evita d'introdurre calcoli con numeri complessi. Ancora per radici reali, sono impiegati procedimenti basati sulla determinazione di un polinomio di grado ≤ 2 che interpola il polinomio dato nelle vicinanze di uno zero (metodo delle secanti, metodo di Müller).

4) Equazioni differenziali. - Per la difficoltà di esporre succintamente risultati relativi a una vastissima letteratura, si rimanda alla bibliografia, limitandosi a un breve cenno su metodi numerici particolarmente in uso.

Metodi variazionali. - Si abbia, per fissare le idee nel modo più semplice, il problema di Dirichlet in R2: dato un insieme aperto limitato D ⊆ R2, trovare una funzione y: Ã → R (Ã = D ⋃ ∂D, ∂D frontiera di D) tale che − ∂2y/x²1 − ∂2y/x²2 + c(x)y = f(x) per x D, y(x) = 0 per x ∈ ∂D. Le funzioni c(x) ed f (x) si suppongono continue in Ã, e c(x) ≥ 0. Introducendo l'operatore L(u) = − Δu + cu con Δ − ∂2/x²1 + ∂2/∂x²2 e lo spazio vettoriale F di tutte le funzioni continue in D insieme alle loro derivate parziali prime e seconde, e nulle su ∂D, il problema equivale a trovare una soluzione dell'equazione L(u) = f, u F. Si introduce in F un prodotto scalare, indicato con (u, v) e definito da (u, v) = ∉Du(x)v(x) dx. Se D è tale che per esso valgono i teoremi integrali di Gauss e Green, si dimostra che L risulta autoaggiunto e definito positivo.

Nel metodo di Ritz si sostituisce al problema L(u) = f, con L generico operatore autoaggiunto e definito positivo, u F, F spazio di Hilbert, quello consistente nel minimizzare su F il funzionale J(u)=(Lu, u)−2(u, f).

Si considera una successione di sottospazi Fh a un numero finito di dimensioni,

La successione {Fh} si suppone completa in F, nel senso che per ogni u F risulta

Per ogni h, si trova uh Fh che rende minimo J(u) in Fh. Si dimostra che la successione {uu} converge in F e il suo limite è soluzione del problema L(u) = f. In pratica il metodo richiede l'introduzione, per ciascun h, di una base per Fh, i cui elementi siano ϕ1h, ϕ2h, ..., ϕnhh e il calcolo dei coefficienti α1, ..., αnh per i quali

rende minimo J(u) in Fh. Ciò porta a risolvere sistemi di equazioni lineari algebriche con matrici (di Gram) di elementi αhij = (Lϕih, ϕjh) e termini noti (f, ϕih), simmetriche e definite positive, di dimensioni via via crescenti. L'opportunità di avere matrici quanto più possibile sparse (cioè con un'elevata percentuale di elementi uguali a zero) ha portato allo studio di basi formate da funzioni nulle su certi sottoinsiei dell'insieme ove si considera il problema differenziale di partenza (metodo degli elementi finiti). Tali basi trovano altresì impiego nel metodo di Galerkin (proposto anche da S. Faedo), che vale sotto condizioni meno restrittive di quelle di auto-aggiunzione e positività. In esso, posto come prima

s'impone che Luyhf sia ortogonale a tutti i vettori di Fh. La matrice del sistema di equazioni lineari algebriche che ne deriva ha per coefficienti aij = (Lϕjh, ϕih).

Bibl.: J. H. Wilkinson, Rounding errors in algebraic processes, Londra 1963; E. Isaacson, H. B. Keller, Analysis of numerical methods, New York 1966; L. Collatz, Functional analysis and numerical mathematics, ivi 1969; A. Ghizzetti, A. Ossicini, Quadrature formulae, Stoccarda 1970; J. D. Lambert, Computational methods in ordinary differential equations, New York 1973; G. I. Marčuk, Methods of numerical mathematics, New York e Berlino 1975; J. Stoer, Introduzione all'analisi numerica, vol. i, Bologna 1976; J. Stoer, R. Bulirsch, Introduzione all'analisi numerica, vol. 2, ivi 1976.