NUMERICI, CALCOLI

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

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 F₁, F₂, ..., F_n, funzionali lineari in E* (il duale di E). Dati gli scalari y₁, y₂, ..., y_n si vuole trovare un n ∈ E tale da aversi F_i(n) = y_i, i ∈ {1, ..., n}. Si dimostra che il problema possiede una e una sola soluzione se e solo se i funzionali F_i 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) = a₀ + a₁ x + .... + a_n-1 x^n-1. Esso è uno spazio vettoriale a n dimensioni e i polinomi 1, x, ... x^n-1 ne costituiscono una base. Prendendo F_i [p(a, x)] = p(a, x_i) per ogni polinomio di E, si vogliono trovare le a_i in modo che p(a, x_i) = y_i, 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 x_i acquista il valore y_i = f(x_i), i ∈ {1, ..., n}.

Il metodo di Lagrange fa uso dei polinomi

vale zero in tutti gli x_j diversi da x_i e vale 1 per x = x_i). 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, ..., (q_i − 1)-esima nel punto x_i (i ∈ {1, ..., n}, si sostituisce a f(x) quel polinomio p(x) di grado ≤ q₁ + q₂ + ... + q_n − 1 (eventualmente ≡ 0) per il quale si ha p^(j)(x_i) = p^(j)(x_i) (j = 1, 2, ..., n_i − 1; i = 1, 2, ..., n).

Posto p_ii(x) per indicare il polinomio di 1° grado che in x_i vale f (x_i) e ha derivata uguale a f′(x_i), p_iii(x) quello di 2° grado che in x_i vale f (x_i) e ha derivate prima e seconda uguali a f′(x_i) e f″(x_i), ..., si può ancora usare il metodo di Neville per ottenere, per es., p₀₀₁(x) a partire da p₀₀(x) e p₀₁(x); p₀₀₁₁(x), noti p₀₀₁(x) e p₀₁₁(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 x₀, x₁, ..., x_n, con x₀ = a, x_n = b, x_i-1 〈 x_i (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 C² [a, b] (possiede cioè derivate prima e seconda continue in [a, b]), vale f (x_i) in x_i, ed è un polinomio di grado ≤ 3 in ogni [x_i-1, x_i] (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 x_k = 2πn/k (k = 0, 1, ..., n − 1) si cerca il polinomio trigonometrico

β_k exp (ikx) (β_k complesso) tale che p(x_k) = f_k (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 x_i, a due a due distinti, e si vuole determinare la funzione razionale

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

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 f_i = f (x_i), (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 x_n dall'ultima equazione, sostituendola nella penultima e calcolando x_n-1 e così via, se tutti gli elementi u_ii 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 P² = 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 l_kk > 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 x⁰ per cui (Ax − B)*(Ax − B) è minimo.

Si trovano una matrice triangolare superiore U con u_ii ≠ o tale che AU = C con C*C = D con D diagonale. Risulta x⁰ = 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⁽¹⁾, x⁽²⁾, ... ottenuta, a partire da un arbitrario x⁽⁰⁾, 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 : Rⁿ → Rⁿ una funzione vettoriale a n componenti delle n variabili x₁, ..., x_n. Si vuole trovare −x tale che f (−x) = 0, oppure, più per esteso, tale che f_i (ù₁, ..., ù_n) = 0 (i = 1, 2, ..., n).

Supposto che le derivate parziali. ∂f_i/∂x_j (i, j = 1, ..., n) esistano in un certo aperto convesso C₀ ⊆ Rⁿ, indicheremo con f′ (x) la matrice (jacobiana) che ha ∂f_i/∂x_j nel posto (i, j). Introdotta una norma ∥•∥ per i vettori e le matrici in oggetto, si può dimostrare quanto segue. Ipotesi: per un certo x⁰ ∈ C₀ e per tre certe costanti α, β, γ risulti: a) ∥ [f′(x₀)]^-1f(x₀) ∥ ≤ α; b) ∥ [f′(x)]^-1 ∥ ≤ β per ogni x ∈ C_0; c) ∥ f′(x) − f′(y) ∥ ≤ γ ∣ x − y ∣ per ogni x ∈ C₀ e y ∈ C₀; d) h = ½ αβγ 〈 1. Tesi: posto r = α(i − h)^-1 e D_r(x⁰) = {x ∣ ∥ x − x⁰ ∥ 〈 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 D_r(x⁰), e converge a un ben determinato ξ ∈ Ã_r(x⁰) per il quale f (ξ) = 0.

Risulta, per ogni k, ∥ x^(k) − ξ ∥ ≤ αh^2k-1 (1 − h^2k-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⁽⁰⁾, x⁽¹⁾, .... Dalla formula di Taylor si ha infatti f (x) ≅ f (x⁰) + f′ (x⁰) (x − x⁰) (qui x − x⁰ va inteso come vettore-colonna), e annullando il 2° membro si ottiene x⁽¹⁾.

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 U^TAU 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 Q_i, R_i tali che A_i = Q_iR_i, Q_i*Q_i = E (l'asterisco indica trasposta coniugata), R_i triangolare superiore, A_i+1 = R_iQ_i. Se gli autovalori di A sono reali e distinti, la successione A_i 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 R²: dato un insieme aperto limitato D ⊆ R², trovare una funzione y: Ã → R (Ã = D ⋃ ∂D, ∂D frontiera di D) tale che − ∂²y/∂x²₁ − ∂²y/∂x²₂ + 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 Δ − ∂²/∂x²₁ + ∂²/∂x²₂ 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 F_h a un numero finito di dimensioni,

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

Per ogni h, si trova u^h ∈ F_h che rende minimo J(u) in F_h. Si dimostra che la successione {u_u} 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 F_h, i cui elementi siano ϕ₁^h, ϕ₂^h, ..., ϕ_nh_h e il calcolo dei coefficienti α₁, ..., α_nh per i quali

rende minimo J(u) in F_h. Ciò porta a risolvere sistemi di equazioni lineari algebriche con matrici (di Gram) di elementi α^h_ij = (Lϕ_i^h, ϕ_j^h) e termini noti (f, ϕ_i^h), 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 Lu^yh − f sia ortogonale a tutti i vettori di F_h. La matrice del sistema di equazioni lineari algebriche che ne deriva ha per coefficienti a_ij = (Lϕ_j^h, ϕ_i^h).

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.