NUMERICI CALCOLI

NUMERICI CALCOLI (XXV, p. 29)

Generalità. - Il concetto di calcolo numerico si può introdurre da un punto di vista generale, come segue. 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 disposizione iniziale degli oggetti nelle celle ad una finale, definiscono un calcolo. La successione delle istruzioni da seguire si dirà un programma. Si ha un calcolo numerico allorché gli oggetti sui quali si opera sono interpretati come cifre, relative ad una certa base di numerazione. L'analisi numerica può dirsi quella parte della matematica che si occupa dell'impiego dei calcoli numerici nella risoluzione di problemi (intendendo come problema la ricerca, in un insieme S, di un sottoinsieme di elementi che godono di una certa proprietà, p). Non è usualmente compito dell'analisi numerica appurare se un tale insieme è vuoto, o se è costituito da un solo elemento, di pervenire, cioè, alla dimostrazione di teoremi di esistenza e di unicità.

Discretizzazione. - Consiste nel rappresentare l'insieme di partenza, S, su un altro, S′, i cui elementi sono raggruppamenti finiti di cifre, e nel sostituire la proprietà p richiesta per le soluzioni in S, con un'altra proprietà, p′, il cui verificarsi o meno è riconoscibile, in S′, mediante un programma. Detto {x′} l'insieme di tutti gli x′ soluzione del problema (S′,p′), si prenderà in considerazione, se {x′} non è vuoto, l'insieme P {x′} di tutti quegli x di S, ai quali corrispondono tali x′. È di grande importanza trovare una maggiorazione del diametro di Γ {x′} (relativo a una certa metrica o pseudo-metrica da introdursi in S): Di questa ricerca fa parte lo studio degli errori di troncamento (sostituzione di algoritmi finiti ad algoritmi infiniti) e di arrotondamento (sostituzione di numeri con un numero finito di cifre ai numeri reali).

Approssimazione. - Un passo intermedio nella discretizzazione è spesso costituito dalla rappresentazione di S su un opportuno sottoinsieme di S stesso. Ad esempio, quella dello spazio delle funzioni reali, continue in un intervallo (a, b) dell'asse reale, sul sottoinsieme di tale spazio formato dalle funzioni razionali intere in una variabile reale (polinomî).

Prenderemo ora in esame alcuni tipici problemi di calcolo numerico, indicando, per ciascuno di essi, uno o più metodi di risoluzione.

Principali problemi di calcolo numerico.

1) Interpolazione secondo Newton e secondo Hermite. - Il problema è la determinazione di quel polinomio g(x), di grado ≤ n, che in n + 1 punti x₀, x₁ ,....., x_n a due a due distinti dell'intervallo (a, b) dell'asse reale acquista certi n + 1 valori assegnati y₀, y₁---, y_n.

Indicando con f_i0, i₁...., i_r, quel polinomio di grado ≤ r che in x_ik vale y_ik, si ha:

in cui va inteso f₀ = y₀, f₁ = y₁, ... Troviamo dapprima

poi

e così via (procedimento di neville). Il polinomio g(x) = f_{0, 1, ..., n}, è quello cercato. La formula [1] si presta bene anche alla ricerca di quel polinomio f(x), di grado ≤ n, che soddisfa le condizioni

in certi s + 1 punti a due a due distinti dell'asse reale. Sarà indicato con

quel polinomio che ha grado ≤ ν₀ + ν₁ + ... + ν_r − 1 e soddisfa le:

Si applicherà la l1] se il polinomio a primo membro non ha tutti gli indici uguali. Nell'altro caso, tale polinomio si ottiene immediatamente mediante la formula di Taylor.

2) Sistemi di equazioni lineari algebriche. - Indichiamo con A una matrice m × n (cioè ad m righe ed n colonne) ad elementi reali; con b una matrice m × 1 (vettore-colonna ad m componenti). L'equazione Ax = b, nel vettore-colonna incognito di componenti x_I, ..., x_n, costituisce un sistema di equazioni lineari algebriche, che può anche scriversi nella forma

Poiché, com'è noto, la risoluzione di un sistema compatibile a matrice m × n può essere sempre ricondotta a quella di un sistema a matrice quadrata, supponiamo senz'altro m = n. Si fa distinzione fra due tipi di metodi di risoluzione: a) metodi diretti, b) metodi iterativi. Fra i metodi del tipo a) citiamo i seguenti:

Metodo di eliminazione (o di Gauss). - Sia a₁₁ ≠ o. Sottraendo alla h^a equazione la 1^a moltiplicata per a_h1/a₁₁ (h = 2,..., n) si passa ad un sistema equivalente al dato (cioè avente la stessa soluzione), con matrice della forma

in cui A₂₂ è una matrice (n − 1) × (n − 1) e O indica una colonna di n − 1 zeri. Siano b′₁, b′₂, ..., b′_n i nuovi termini noti. Si opera nello stesso modo sul sistema A₂₂ x⁽¹⁾ = b, in cui x⁽¹⁾ indica il vettore di componenti x₂, ..., x_n, e b⁽¹⁾ quello di componenti b′₂, ..., b′_n. Ripetendo tale procedimento, si perviene ad un sistema formato da una singola equazione in una incognita. Ottenuta questa ultima, si sostituisce in una delle due equazioni del sistema immediatamente precedente, e così via.

Metodo dei gradienti coniugati. - Si parte da un arbitrario vettore x⁽⁰⁾, e si calcola in successione, r⁽⁰⁾ = b − Ax⁽⁰⁾, p⁽⁰⁾ = A^Tr⁽⁰⁾,

in cui A^T denota la matrice trasposta di A. Il vettore x(ⁿ) fornisce la soluzione cercata.

Metodo di ortogonalizzazione. - Nell'ipotesi che A sia una matrice m × n, con m ≥ n ed a colonne linearmente indipendenti, questo metodo fornisce il vettore x per cui è minimo (Ax − b)^T (Ax − b), e quindi, in particolare, la soluzione di Ax = b se m = n. si trovano le due matrici C e P con P triangolare superiore (in cui cioè sono tutti nulli gli elementi al disotto della diagonale principale), avente tutti 1 sulla diagonale principale, tali che AP = C, C^T C = D, con D diagonale. Risulta allora x =PD^-1 C^Tb, in cui D^-1 denota la matrice diagonale avente per elementi i reciproci di quelli di D.

Fra i metodi del tipo b) citiamo i seguenti:

Metodo diagonale. - Si decomponga A nella somma L + D + U, in cui L è triangolare inferiore con tutti 0 sulla diagonale principale; D è diagonale (la diagonale principale di A) e U è triangolare superiore. Partendo da x⁽⁰⁾, vettore arbitrario, si calcolano in successione i vettori x⁽ⁱ⁾ (in base alla formula: Dx^(i+I) = b − (L + U)x⁽ⁱ⁾. Risulta

soluzione di Ax = b, se e solo se [D^-1 (L + U)]ⁱ → O per i → + ∞.

Metodo di Seidel. - Usando i simboli del metodo precedente, la successione x⁽ⁱ⁾ (i = o, 1, ...) è definita da: (L + D)x^{(i + 1)} = b − Ux⁽ⁱ⁾. La condizione [(L + D)^-1 U]ⁱ →- O per i→ + ∞ è necessaria e sufficiente per la convergenza di x⁽ⁱ⁾ ad x. Vanno anche tenute presenti, come sufficienti, l'una o l'altra delle due condizioni: 1) A matrice simmetrica, definita positiva; 2)

Metodo dei gradienti. - Si assuma a piacere una matrice B simmetrica, definita positiva. Partendo da x⁽⁰⁾ arbitrario si calcolano in successione i vettori x(ⁱ) secondo lo schema:

3) Polinomio caratteristico e radici caratteristiche di una matrice. - Sia : A una matrice n × n, E la matrice unitaria dello stesso ordine. Il determinante di λE-A, sviluppato secondo le potenze di λ, è il polinomio caratteristico di A e le sue radici si dicono radici caratteristiche di A. Ricordiamo ancora che il simbolo trA (leggi "traccia" di A) denota la somma degli elementi sulla diagonale principale di A.

Tra i metodi per calcolare il polinomio caratteristico, o, direttamente, le radici caratteristiche di A, ricordiamo:

Metodo di Soitriaa-Faddeiev. - Posto A₁ = A, si calcola p₁ = tr A₁; B₁ = A₁ − p₁E; A₂ = AB;

B₂ = A₂ − p₂E, e così via sino ad A_n = AB_n-1,

B_n = A_n − p_nE. Il polinomio

è il polinomio caratteristico di A. Si dimostra che B_n = O, ciò che costituisce un efficace controllo dei calcoli.

Metodo di Danilewski. - Se A ha la forma

con E₂₁ matrice unitaria, d'ordine n − 1, ed O matrice (n − 1) × 1 con tutti gli elementi nulli, il suo polinomio caratteristico si scrive subito, avendo come coefficienti p₁, ..., p_n, ordinatamente, gli elementi situati sulla prima riga di A. Se A ha la forma

con O matrice q × (n − q) ad elementi tutti nulli, il suo polinomio caratteristico si spezza nel prodotto dei due polinomî caratteristici di A₁₁ e di A₂₂. È noto che, se B è una matrice n × n a determinante ??? 0 e B^-1 è la sua inversa, la matrice B^-1AB ha come polinomio caratteristico quello stesso di A. Il metodo di Danilewski consiste nel sostituire successivamente ad una matrice A assegnata, delle matrici A₁ = B¹₁ AB₁, A₂ = B^-¹₂ A₁ B₂, ecc., così da pervenire, in un numero finito di passi, alla forma [1] o [2].

Metodo di Givens. - Indichiamo con U_ij(θ) la matrice n × n che si ottiene dalla matrice unità ponendo cos θ nel posto (i, i), sin θ nel posto (i, i), sin θ nel posto (i, j) cos θ in quello (j, j). Data una matrice A n × n simmetrica, la matrice A⁽¹⁾ = UT_ij(θ) A U_ij(θ), indicando con UT_ij(θ) la trasposta di Uij(θ), ha, per un opportuno θ, nulli i due elementi nei posti (i,j) e (j, i). Mediante successive trasformazioni del tipo ora indicato, in numero finito, si riduce la matrice A ad una T di forma tridiagonale, cioè con elementi tutti nulli nei posti (i, j) per cui ∣ i − j ∣ > 1. Il polinomio caratteristico di T si trova facilmente, osservando che, detta T_i la matrice formata con le prime i righe e colonne di T, ed f_i(λ) il suo polinomio caratteristico, risulta f₁(λ) = λ − t₁₁, f₂(λ) = (λ − t₂₂) f₁(λ) − t₁₂t₂₁,

Il polinomio f_n(λ) è il polinomio caratteristico di A.

Metodo di Jacobi. - Nell'ipotesi che A sia simmetrica, si passa da A ad A⁽¹⁾ come nel metodo precedente, con l'avvertenza di scegliere per i e j gli indici corrispondenti a un elemento avente massimo modulo. Operando nello stesso modo a partire da A^(I), e così via, si costruisce una successione di matrici, convergente ad una matrice diagonale, la quale ha come elementi le radici caratteristiche di A.

4) Equazioni non Lineari. - Per il metodo di Newton e il metodo di Gräffe nella risoluzione delle equazioni non lineari rinviamo alla voce Approssimazione nel vol. III, rispettivamente alle pp. 770, 771.

Metodo di Lin-Bairstow. - Sia

un polinomio a coefficienti reali. Presi due numeri reali r⁽⁰⁾, s⁽⁰⁾ a piacere, si costruiscono le due successioni definite da

e si risolve il sistema nelle due incognite Δr⁽⁰⁾, Δs⁽⁰⁾

Si itera il procedimento a partire da r⁽¹⁾ = r⁽⁰⁹ + Δr⁽⁰⁾, s⁽¹⁾ = s⁽⁰⁾ + Δs⁽⁰⁾ in luogo di r⁽⁰⁾, s⁽⁰⁾ e così via. Le r⁽ⁱ⁾, s⁽ⁱ⁾ convergono a due numeri r, s, tali che z² + rz + s risulta un divisore di f(z).

Metodo di Bernoulli. - Sia

un polinomio a coefficienti reali, e sia a^o = 1. Distinguiamo due casi: nel caso 1), dette α_i (i = 1,..., r) le radici di f(z) = o, a due a due distinte, con ∣ α₁ ∣ ≥ ∣ α₂ ∣ ... ≥ ∣ α_r ∣ risulti ∣ α₁ ∣ > ∣ α_i ∣ (i = 2, ..., r); nel caso 2) invece risulti ∣ α₁ ∣ ≥ ∣ α₂ ∣ > ∣ α_i ∣ (i = 3, ..., r).

Si ponga

Nel caso 1) è b_h1 definitivamente ??? 0, e ciascuna successione

converge, per h → + ∞. Detto b_k il suo limite, il polinomio

è un divisore di f(z). Nel caso 2) si costruiscano come prima le b_hk e si ponga c_hk = b_h1 b_h+1,k − b_hk b_h+1,1, (k n= 2, ..., n). È c_h2 definitivamente ??? 0, e ciascuna successione c_hk/c_h2 (k = 2,..., n) converge, per h → + ∞. Detto c_k il suo limite, il polinomio

è un divisore di f(z).

Metodo dicotomico nel campo complesso. - Sia

a coefficienti reali, e a_n =1; si ponga 2μ_h = n − h 0 n − h − i a seconda che n − h sia pari o no (h = 0,1, .., n), e si definiscano i numeri a_i0, a_i1 (i = 0,1 ..., μ₀ e 0,1, ...μ₁ rispettivamente) come segue:

Successivamente, si calcolino i numeri α_ij definiti da α_ji = 0 per i > μ_j, α_ij = − α_o,j-2 α_{i+1, j-1} − α_{0, j-1} α_{i+1, j-2} per i = 9,1, ..., μ_j, j = 0,1, ..., n. La condizione α_ij > 0 (in = 0, ..., μ_j, j = 0, ..., n) è necessaria e sufficiente perché f(z) sia un polinomio di Hurwitz (abbia cioè, tutie le sue radici a sinistra dell'asse immaginario). Si supponga che, dall'esame delle α_ij, si sia riconosciuto che f(z) è di Hurwitz. Si cerca allora un h 〈 0, tale che, trovati i coefficienti di f (z + h) e costruite le corrispondenti α_ij, almeno una di esse sia ≤ 0. Nella striscia compresa fra le rette x = − h ed x = 0 cadono le radici di f(z) = 0 aventi parte reale massima. Operando nello stesso modo su f(z + h/2) si riduce a metà la larghezza della striscia, e così via, sino a trovarne una, sulla cui mediana, x = c si riterrà che cadano tutte le radici di f(z) = 0 aventi parte reale massima. Nel caso che f(z) non sia di Hurtwitz, si prenderà in esame f(z + h) con h > 0 e si procederà come prima. Supposto f(c) ??? 0, i coefficienti dell'immaginario delle radici complesse di f(z) situate sulla retta x = c si trovano esaminando le varie colonne di numeri a_ij (j = 0,1, ...) costruiti per f(z + c). Si constata che per un certo p ≥ 2 è α_oo > 0, α₀₁ > 0, ..., α_{0, n-p} > 0, mentre α_{0, n-p+1} = 0ao. Posto ν = n − p, l'equazione in y:

ha come sue radici, tutte reali, i detti coefficienti.

5) Sistemi di equazioni non lineari. - Metodo dei gradienti. - Sia A un dominio rettangolare dello spazio reale euclideo Rⁿ; f_i (x₁, ..., x_n) (i = 1, .., n) n funzioni reali continue in A ed F(x₁, ..., x_n) una funzione reale, definita in A, due volte parzialmente derivabile rispetto alle x_i in A con derivate ivi continue. Inoltre il sistema f_i (x₁, ......, x_n) = 0 (i = 1, ..., n) abbia una soluzione, x1*, ..., xn*; in tal punto F(x₁, ..., x_n) abbia un minimo proprio. Se il punto x0₁, ..., x0_n è interno a un dominio contenente P* racchiuso da una superficie di livello F(x₁, ..., x_n) = c, e nel quale l'unico a F punto ove le derivate parziali

sono simultaneamente nulle è P*, il problema di Cauchy:

ha come soluzione un vettore x₁(t), ..., x_n(t), le cui componenti convergono a x1*, ..., xn* per t→ + ∞.

Metodo delle funzioni implicite. - Si sostituisce al sistema 1) f_i(x₁, ..., x_n) = o (i = 1, ..., n) l'altro: 2) g_i(x₁, ..., x_n, λ) = 0 (i = 1, ..., n), con le g_i continue e parzialmente derivabili in A × I (I intervallo (0,1)), scelte in modo da ridursi alle f_i per λ =1, e tali che il sistema

sappia risolversi senza difficoltà. Sia (x0₁, ..., x0_n, 0) soluzione di quest'ultimo. Se il sistema 2) e le condizioni g_i(x0₁, ..., x0_n, 0) = 0 (i = 1, ..., n) definiscono implicitamente n funzioni x₁(λ), ..., x_n(λ) in tutto (0,1), la n − upla x_i(1) (i = 1, ..., n) fornisce una soluzione di 1).

6) Risoluzione numerica del problema di Cauchy per sistemi di equazioni differenziali ordinarie, del 1° ordine. - Sia A un intervallo aperto dell'asse reale, B un insieme aperto dello spazio reale euclideo R^m, f(x, y) una funzione vettoriale ad m componenti (reali) continua in A × B, e (x₀, y₀) un punto fissato di A × B. (A × B denota l'insieme delle coppie ordinate x, y con x in A ed y in B). Il problema è quello di trovare un intervallo aperto I ⊂ A, e contenente x₀, ed una funzione y(x) ad m componenti (reali) verificante in I il sistema di m equazioni differenziali ordinarie in m incognite

e la condizione y(x₀) = y₀. Supposti esistenti I ed y(x), e y(x) unica, siano x₀, x₁, ..., x_N punti di I, con x_n 〈 x_n+1 (n = 0,1..., N − 1). Si vuole costruire una tabella di vettori y₀, y₁, .., y_N ad m componenti approssimanti i vettori y(x₀), y(x₀), ..., y(x_N), in cui y(x) è la soluzione "esatta". Di solito x_n+1 − x_n (n = 0,1, ..., N − 1) è una costante, h.

Metodi a passo singolo. - Si calcola y_n+1, una volta noto y_n, mediante una formula del tipo y_n+1 = y_n + hϕ(x_n, y_n, h), con ϕ(x, y, h) opportuna funzione vettoriale verificante la condizione

in A × B. Si dirà che il metodo è d'ordine p, se, posto

risulta r(x, h) = O (h^p+1) per h → 0. Nelle tabelle alla pag. precedente sono riportati alcuni casi particolari più in uso.; per ciascuno di essi è indicato il modo di calcolare la funzione ϕ(x, y, h) a partire dalla funzione nota f (x, y) e l'ordine del metodo.

Metodi a passi multipli. - Intervengono nel calcolo di y_n+1, i k precedenti valori y_n, y_n-1, ..., y_n-k+1(k > 1). Posto f_n = f(x_n, y_n) si ía uso di una formula del tipo:

La [2] fornisce esplicitamente y_n+1 se b₀ = 0, implicitamente se b₀ ??? 0. Spesso una formula con b_0??? 0 è più accurata di una con b₀ = 0; si fa allora uso della seconda per "predire" y_n+1, calcolare f_n+1 e della prima per "correggere" y_n+1, sostituendovi a secondo membro il valore f_n+1 precedentemente trovato. Esempio (metodo di Milne)

Bibl.: a) Opere a carattere generale: D. R. Hartree, Numerical analysis, Oxford 1952; F. B. Hildebrand, Introduction to numerical analysis, New York 1956; A. S. Householder, Principles of numerical analysis, ivi 1953; L. V. Kantorovic e V. I. Krylov, Näherungsmethode der höheren Analysis, Berlino 1957; E. Grabbe, S. Ramo, D. Wooldridge, Handbook of automation, computation and control, vol. I e II, 1959; A. Korganoff, Méthodes de calcul numérique, vol. I, Parigi 1961.

b) Memorie (Il numero progressivo si riferisce al paragrafo in questo articolo). 1. E. H. Neville, Iterative interpolation, in Journ. Indiana Math. Soc., XX (1933), pp. 87-120; 2. M. R. Hestenes e E. Stiefel, Methods of conjugate gradients for solving linear systems, in Journ Res. Nat. Bur. Standards, IL (1952), pp. 409-436; 3. A. M. Danilewski, Sulla risoluzione numerica delle equazioni secolari, in Mat. Sbornik 2, XLIV (1937), pp. 169-172.