La seconda rivoluzione scientifica: matematica e logica. L'analisi numerica
La seconda rivoluzione scientifica: matematica e logica. L'analisi numerica
L'analisi numerica
L'analisi numerica moderna comincia a delinearsi verso la metà del XX sec., con le prime ricerche sulla stabilità dei metodi diretti di risoluzione di sistemi lineari, o di inversione di matrici, dovute principalmente a John von Neumann, Herman H. Goldstine e Alan M. Turing. In queste ricerche si affronta in un modo completamente nuovo lo studio degli algoritmi numerici, grazie soprattutto alla costruzione del primo grande calcolatore della storia, l'ENIAC, presso la University of Pennsylvania a Filadelfia intorno al 1945, e ai primi sviluppi di un calcolo scientifico su grande scala, caratterizzato cioè dall'elevata dimensione dei problemi affrontati. L'effetto del calcolatore è stato duplice: da un lato l'accrescimento del potere di risoluzione numerica dei problemi della matematica applicata; dall'altro l'insorgere, proprio a causa di quel nuovo potere, di una serie di difficoltà e di questioni mai prima sollevate. Al posto dell'aritmetica esatta dei numeri reali il calcolatore usa un'aritmetica approssimata, in cui la maggior parte delle proprietà delle operazioni tra numeri cessa di valere. La macchina opera infatti con un insieme finito di numeri finiti e introduce quindi errori di arrotondamento sui dati e sui risultati delle singole operazioni. Le leggi di propagazione di questi errori dipendono di solito dalla dimensione del problema. Per esempio, la sensibilità della soluzione di un sistema di equazioni lineari agli errori sui dati, misurata dall'indice di condizionamento della matrice dei coefficienti, possono dipendere dalla grandezza di questa matrice. La dimensione dei problemi e la mole di calcoli, che la velocità dei processi automatici consente di portare a limiti prima inimmaginabili, possono quindi provocare fenomeni di instabilità numerica, cioè di accrescimento incontrollabile dell'errore. Una difficoltà consiste anche nel fatto che i risultati intermedi di un calcolo sono normalmente nascosti nella memoria della macchina e sconosciuti al programmatore, cosa che rende inevitabile una delega del pensiero umano al processo automatico. Come ha osservato Germund Dahlquist nel 1985, il 'rumore' causato da errori di arrotondamento, da instabilità o da malcondizionamento, deriva più da questa inaccessibilità delle fasi intermedie del calcolo che dall'aumento della mole di operazioni. È questo l'aspetto più tipico, conseguenza della realizzazione del primo calcolatore della storia, di un principio di delega già prefigurato nelle Règles pour la direction de l'esprit (1628) di René Descartes.
La scienza del calcolo: aritmetizzazione e algoritmi efficienti
Il calcolo automatico poteva anche considerarsi una sorta di 'materializzazione' della logica costruttiva, intuizionista, elaborata nei primi decenni del Novecento (la matematica di Brouwer, notava von Neumann in una lettera a Norbert Wiener del 1946, può tradursi in un appropriato meccanismo). Turing, von Neumann, Warren McCulloch e Walter Pitts dimostrarono che la logica intuizionista può essere tradotta integralmente in termini di macchine automatiche e che le proposizioni logiche sono rappresentabili come reti elettroniche o sistemi nervosi idealizzati. Tra logica e rete c'era tuttavia una differenza essenziale: il passaggio dalla formula matematica al processo computazionale rendeva inevitabile l'intromissione di elementi estranei quali lo spazio fisico, necessario a memorizzare dati e risultati delle operazioni, e il tempo di esecuzione. Questa intromissione era comunque intesa come un motivo di arricchimento anziché di snaturamento della formula. Il tempo che intercorre tra il segnale di ingresso e il segnale di uscita, si notava, può infatti contribuire a evitare l'occorrenza di vari circoli viziosi nascosti dovuti a impredicatività, non costruttività o altro. Il tempo di esecuzione dipende inoltre, essenzialmente, dalla complessità dell'algoritmo. Procedure aritmetiche finite e ricorsive, considerate fino a quel momento il nucleo più certo e sicuro della matematica (l'unico possibile fondamento, per David Hilbert, di una teoria coerente dell'infinito), potevano allora rivelarsi troppo complesse (con un numero di operazioni troppo elevato) oltre che numericamente instabili.
Con lo studio sistematico della stabilità e della complessità cominciò a svilupparsi, dalla metà del Novecento, una scienza del calcolo che doveva tener conto di limitazioni fisiche di spazio e di tempo e dimostrare che anche in un sistema di numeri estremamente più povero e meno strutturato di quello dei numeri reali si potevano risolvere i più difficili problemi della matematica applicata. Il termine 'scienza' denota qui un passaggio critico. Prima degli anni Cinquanta, notavano in particolare George Forsythe e Alexander Ostrowski, si poteva parlare soltanto di un'arte del calcolo, cioè di una generica abilità nell'uso di procedure capaci di fornire soluzioni a svariati problemi numerici, ma di tali procedure si capiva spesso relativamente poco. Questa distanza tra prassi e teoria si accompagnava di solito alla propensione per una ricerca della migliore tattica anziché della migliore strategia di risoluzione. Ostrowski osservava che questa maggiore propensione alla tattica finiva per produrre metodi meno efficienti, come per esempio metodi iterativi più lenti, rispetto a quelli che una scienza più sistematica, basata su teoremi di carattere puramente matematico, avrebbe potuto suggerire. La critica era qui rivolta, soprattutto, ai metodi di rilassamento introdotti da Richard V. Southwell intorno al 1946: metodi di tipo variazionale, aciclici, non automatizzabili, la cui esecuzione, passo per passo, era prevalentemente affidata all'estro e all'arbitrio della mente umana, le cui sottigliezze, si diceva, superano regolarmente quelle della macchina.
Teoremi e risultati teorici elaborati negli ultimi decenni dell'Ottocento e nei primi del Novecento favorirono questo passaggio decisivo dall'arte alla scienza del calcolo e ne costituirono anzi il presupposto indispensabile. è sorprendente anche la sopravvivenza, nel calcolo scientifico, di strategie algoritmiche della matematica antica, di cui si è riscoperta l'efficienza nella risoluzione dei più complessi problemi numerici. Il calcolatore esegue infatti, di norma, calcoli di aritmetica elementare e furono appunto i matematici antichi (greci, indiani, cinesi e babilonesi) a introdurre metodi puramente aritmetici, spesso di carattere iterativo, per risolvere equazioni quadratiche con soluzioni irrazionali. L'importanza del calcolo iterativo, puramente aritmetico, nel calcolo scientifico su grande scala era così sottolineata da Goldstine e von Neumann in On the principles of large-scale computing machines (1946): "I nostri problemi sono dati di solito come problemi analitici in termini di variabili continue, spesso di carattere totalmente o parzialmente implicito. Per le esigenze del calcolo digitale questi devono essere sostituiti, o piuttosto approssimati, da procedure puramente aritmetiche, finitiste ed esplicite (per passi successivi o iterative)" (von Neumann 1961-63, V, p. 5).
In altri termini, il programma di sviluppo della scienza del calcolo comportava una sorta di aritmetizzazione della matematica: una riduzione dei modelli differenziali e integrali, utili alla simulazione di fenomeni di varia natura, a processi aritmetici elementari. Questa riduzione doveva poter avvenire per diversi stadi, ciascuno dei quali produceva qualche errore di approssimazione: (a) la discretizzazione di problemi differenziali e integrali ottenuta sostituendo derivate con rapporti incrementali e le variazioni continue di una variabile con le differenze tra i valori assunti dalla variabile nei punti, sufficientemente vicini, di una griglia discreta; (b) la formulazione di un problema aritmetico al posto del modello iniziale, per esempio un sistema di equazioni lineari al posto di un'equazione alle derivate parziali, oppure una funzione algebrica che approssima una funzione trascendente; (c) la scelta di un modello di calcolo e di un algoritmo in base a speciali criteri di efficienza (stabilità, complessità, condizioni di applicabilità, robustezza); (d) l'esecuzione dei calcoli previsti dall'algoritmo codificati in uno speciale linguaggio di programmazione.
I cardini di questo programma erano dunque i metodi di discretizzazione di problemi differenziali e integrali definiti sul continuo, i metodi di approssimazione di una funzione e l'analisi dell'efficienza degli algoritmi diretti e iterativi deputati a risolvere i problemi algebrici o aritmetici ottenuti dalla discretizzazione.
Il processo di aritmetizzazione implicito nei metodi numerici non era naturalmente estraneo alla questione teorica, posta nel corso dell'Ottocento, della possibile aritmetizzazione dell'analisi matematica. In seguito alle ricerche di Augustin-Louis Cauchy, Georg Friedrich Bernhard Riemann, Karl Theodor Wilhelm Weierstrass, Richard Dedekind, Jules-Henri Poincaré e Georg Cantor, per citare solo alcuni nomi, si era formata la convinzione che tutta l'analisi potesse, almeno in linea di principio, fondarsi sulle operazioni dell'aritmetica dei numeri interi e su operazioni di passaggio al limite. Un ruolo importante, per rafforzare questa convinzione, ebbe la teoria delle equazioni integrali di Erik Ivar Fredholm, che prevede l'introduzione di nuclei trattabili come limiti di matrici. In questa teoria era già evidente un rapporto tra 'discreto' e 'continuo' che poteva preludere all'idea di approssimare un'equazione integrale con un sistema di equazioni algebriche. Nei decenni successivi tale idea si consolidò, fino a suggerire nuove strategie computazionali basate su questa approssimazione e su algoritmi efficienti in cui si realizzava un processo di aritmetizzazione non soltanto teorica ma concretamente legata a uno spazio e a un tempo di calcolo, e alla inevitabile intromissione di errori nelle operazioni di macchina. Per esempio, i primi risultati di Wiener sulle serie temporali (1949), elaborati durante la guerra e indipendentemente dalle ricerche di Andrej Nikolaevič Kolmogorov del 1941, contenevano equazioni integrali approssimabili con sistemi di equazioni algebriche lineari. In questi sistemi, come avrebbe dimostrato Norman Levinson nel 1949, la matrice dei coefficienti è di Toeplitz (perché proviene da un nucleo della forma K(t−s)), cioè con una struttura che consente un processo di aritmetizzazione mediante algoritmi particolarmente efficienti. Quindi non si trattava più, propriamente, del classico progetto di fondare l'analisi sul concetto di numero intero, come intendevano i matematici di fine Ottocento e primo Novecento, ma di un programma di riduzione al numero intero per via di un calcolo approssimato effettivo, con i limiti di spazio e di tempo imposti dalla macchina.
Processi di discretizzazione
Il calcolo scientifico su grande scala cominciò a delinearsi verso la metà del XX sec. ma furono orientamenti e teorie precedenti a fornire le basi per il suo pieno sviluppo. La teoria degli spazi di Hilbert, per esempio, avrebbe fornito il necessario presupposto teorico per i processi di approssimazione di funzioni. Il celebre teorema di Weierstrass del 1885 sull'approssimazione di funzioni continue mediante polinomi sarebbe stato dimostrato da Sergej Natanovič Bernštejn nel 1912, utilizzando una classe di polinomi che si sarebbero rivelati, in seguito, di straordinario interesse applicativo (per es., nell'industria delle automobili degli ultimi decenni). I polinomi ortogonali, di centrale importanza nella teoria dell'interpolazione e dell'integrazione numerica (come pure in diversi settori dell'algebra lineare numerica) furono introdotti ben prima della rivoluzione informatica degli anni Cinquanta, tra la fine del XVIII e l'inizio del XIX secolo. Furono poi diversi matematici, come per esempio Pafnutij L′vovič Čebyšev e Charles Hermite, a farne un uso sistematico nel corso dell'Ottocento. Le tabelle di approssimanti di Henri E. Padé risalgono al 1892, anche se analoghi procedimenti di approssimazione razionale di una funzione furono studiati da Carl Gustav Jacob Jacobi e da Ferdinand Georg Frobenius, rispettivamente nel 1846 e nel 1881.
Le tecniche di discretizzazione e diversi algoritmi per la risoluzione numerica di equazioni differenziali e integrali risalgono alla seconda metà dell'Ottocento e alla prima metà del Novecento. Lo stesso vale per diverse procedure iterative applicate a equazioni algebriche, equazioni integrali e alle derivate parziali. I metodi di risoluzione, con formule di quadratura, di equazioni differenziali ordinarie, che perfezionarono tra la fine del XIX sec. e l'inizio del XX l'idea base del metodo 'poligonale' elaborata da Leonhard Euler negli anni 1768-1770, ebbero anche un ruolo insostituibile nello sviluppo delle tecniche di discretizzazione numerica. La letteratura in questo campo nella prima metà del Novecento è enorme e vede impegnati scienziati di spicco, come Richard von Mises e Lothar Collatz, ai quali si deve, per la varietà dei loro contributi, la moderna concezione di una matematica applicabile a tutti i settori della scienza attraverso modelli e tecniche di risoluzione numerica. Von Mises, in particolare, si occupò di analisi numerica (gli si devono importanti contributi sulla stima dell'errore in problemi ai valori iniziali), di geometria, di teoria della probabilità, di matematica statistica, di idrodinamica e di aerodinamica. Uno degli scienziati più impegnati in questo settore fu anche Carl Runge, professore a Gottinga dal 1904 al 1924, al quale si deve la prima idea, nel 1895, del metodo di integrazione numerica noto come metodo di Runge-Kutta.
Analoghe considerazioni valgono per la risoluzione numerica di equazioni alle derivate parziali. I metodi alle differenze applicati a problemi di tipo ellittico fornirono già nei primi decenni del Novecento un procedimento per tradurre il continuo nel discreto (in termini più suggestivi: l'infinito nel finito), approssimando problemi differenziali con sistemi di equazioni algebriche lineari risolvibili con le operazioni dell'ordinaria aritmetica. A equazioni di tipo ellittico si riconducono problemi di torsione, problemi di moto di fluidi viscosi incomprimibili, di flusso stazionario di calore o di elettricità in conduttori omogenei. Le matrici dei sistemi lineari che provengono dalla discretizzazione di queste equazioni hanno una struttura speciale, che semplifica radicalmente il processo di risoluzione numerica.
Nel 1909 Runge propose di risolvere il problema di Dirichlet ∠2u=0 in un dominio regolare Ω del piano (x,y), con condizioni al contorno, introducendo una griglia Ωh di passo h e sostituendo alle derivate di u quozienti alle differenze. Il fatto che l'errore di discretizzazione risultasse dell'ordine di h2 suggeriva che la soluzione Uij del sistema lineare dovesse differire dal corrispondente valore u(xi,yj) assunto da u nel punto (xi,yj) della griglia Ωh per una quantità proporzionale a h2. Un metodo alle differenze finite fu anche l'espediente introdotto da Lewis F. Richardson nel 1910 per una risoluzione numerica di sistemi di equazioni idrodinamiche. Queste equazioni erano in grado di descrivere, come aveva dimostrato lo scienzato norvegese Vilhelm F.K. Bjerknes (1862-1951), complessi fenomeni meteorologici, sulla base di dati osservati sullo stato dell'atmosfera in singoli istanti di tempo. I risultati di Richardson e la descrizione del suo metodo furono raccolti in una monografia pubblicata nel 1922. Nella prefazione Richardson si augurava profeticamente che i calcoli potessero un giorno procedere più velocemente dei mutamenti atmosferici; tuttavia, per le risorse allora disponibili (valutate in 64.000 'automi' umani provvisti di calcolatrici da tavolo), questa eventualità gli appariva un sogno.
La tecnica delle differenze finite fu applicata da Wiener e da Henry B. Phillips alla risoluzione numerica del problema di Dirichlet in un articolo del 1923. L'interesse di Wiener per questa tecnica è spiegabile anche per un'affinità di progetti dell'analisi numerica e dell'informatica nelle prime fasi del loro sviluppo nel XX sec., dovuta in particolare al comune interesse per il concetto di algoritmo e per l'importanza degli algoritmi in ogni tipo di applicazione della matematica, della fisica e dell'ingegneria. La ricerca che Wiener condusse sul problema di Dirichlet è riconducibile anche allo sviluppo della teoria del potenziale, con Pierre-Simon de Laplace, Siméon-Denis Poisson, Carl Friedrich Gauss e George Green, dopo che Joseph-Louis Lagrange aveva introdotto nella meccanica newtoniana il concetto di energia potenziale, dimostrandone l'importanza per la dinamica e la teoria della gravitazione. Di particolare importanza, in questa teoria, era appunto il problema di Dirichlet, che funziona da modello della distribuzione di temperatura in una camera R, dopo che una preassegnata distribuzione di temperatura sulle pareti circostanti ha prodotto un effetto duraturo e uniforme all'interno di R.
Accanto ai metodi di discretizzazione alle differenze finite si cominciava a introdurne altri di tipo variazionale, su cui si fondano anche i più recenti metodi degli elementi finiti. Questi metodi si basavano sulla determinazione di un funzionale F[w] tale che la soluzione y del problema ellittico rendesse minimo F[w] rispetto a tutte le funzioni w sufficientemente regolari ‒ dette funzioni 'ammissibili' ‒ per cui valessero le stesse condizioni al contorno stabilite per y. Sulla base di un'idea di Walter Ritz (1909), sperimentata con successo anche da John W. Strutt nel 1873, per calcolare un'approssimazione discreta di y, le funzioni ammissibili potevano essere vantaggiosamente definite come funzioni continue lineari a tratti dipendenti da un insieme di N parametri reali arbitrari wi, in particolare w(x)=wj+(wj+1−wj)(x−xj)/h, xj≤x≤xj+1,0≤j≤N, essendo h la lunghezza dell'intervallo compreso tra due punti successivi xj e xj+1. In tal caso il funzionale F[w] risultava una forma quadratica nei parametri wi e poteva essere minimizzato come funzione di questi parametri risolvendo le equazioni F[w]/wi=0 per 0≤i≤N. Queste equazioni si traducevano infine in un sistema lineare definito da una matrice la cui struttura avrebbe consentito l'applicazione sistematica di metodi efficienti di risoluzione numerica. Successivamente sono state usate, invece di funzioni lineari a tratti, funzioni regolari, polinomiali a tratti, in particolare funzioni spline. Sia le funzioni spline sia i metodi di tipo variazionale devono la loro scoperta a diverse necessità applicative.
Le funzioni spline erano usate nell'architettura navale (il termine designa le asticelle di legno da unire tra loro per la costruzione di uno scafo). Il metodo degli elementi finiti, in particolare, deve il suo sviluppo negli anni Sessanta a un gruppo di ingegneri orientati alla ricerca. Nell'edizione del 1978 del libro di Olgierd C. Zienkiewicz, The finite element method in engineering science, è inclusa un'ampia rassegna delle possibili applicazioni del metodo.
Ugualmente fondamentale, soprattutto per lo sviluppo dell'algebra lineare numerica nella seconda metà del Novecento, risultò la teoria matematica dei metodi iterativi per risolvere sistemi di equazioni lineari ‒ senza considerare in dettaglio la propagazione dell'errore ‒ che risale all'ultimo Ottocento e primo Novecento, con gli studi effettuati da Friedrich Wilhelm Karl Ernst Schröder (1870), Philipp Ludwig von Seidel (1874), Richardson (1910), Lamberto Cesari (1937), Helmut Wittmeyer (1936), Hilda Geiringer (1949). Gli aspetti matematici più interessanti e più utili di questa teoria dipendono, di norma, dalla speciale 'struttura' che le matrici dei coefficienti dei sistemi lineari ereditano, per via di processi di discretizzazione, dai modelli differenziali (per es., da diversi tipi di problemi al contorno per equazioni alle derivate parziali). è precisamente questa struttura (per es., la distribuzione degli autovalori o l'appartenenza della matrice a un'algebra) che permette di definire algoritmi stabili e veloci. La convergenza dei metodi iterativi per risolvere un sistema lineare dipende in particolare dal comportamento degli autovalori della matrice di iterazione.
Molto importanti, per questo motivo, risultarono gli studi condotti da S. Gerschgorin nel 1931 sulla collocazione degli autovalori di una matrice sul piano complesso e quelli di P. Stein e R.L. Rosenberg del 1948 sulle proprietà spettrali delle matrici di iterazione dei classici metodi iterativi di Jacobi e Gauss-Seidel. In un fondamentale trattato del 1962, Matrix iterative analysis, Richard Varga osservava che ai sistemi lineari che risultavano dalla discretizzazione di problemi ellittici si potevano applicare metodi iterativi definiti da una matrice di iterazione B non negativa (con bij≥0 per ogni i e j). Si poteva allora applicare la teoria di Perron (1907) e di Frobenius (1912) per tali matrici e per ricavare informazioni sugli autovalori di B e sulla convergenza dei metodi. Anche l'elaborazione di metodi variazionali, come il metodo dei minimi quadrati, precedette la costruzione dell'ENIAC e la rivoluzione delle scienze computazionali del secondo Novecento. Il metodo dei minimi quadrati, nel caso più generale non lineare, chiamato metodo del gradiente (o metodo di steepest descent), era stato già proposto da Cauchy nel 1847 e ripreso da diversi autori nei primi decenni del Novecento.
Tecniche del punto fisso
Tra i principali modelli di riferimento per le più recenti tecniche iterative si possono ancora annoverare alcune procedure note fin dalla matematica antica e medievale, come la regula falsi o la procedura che prese il nome di 'algoritmo di Newton-Raphson'. Questo algoritmo, per risolvere numericamente un'equazione f(x)=0, si riassume nella formula iterativa xi+1=g(xi), dove g(x)=x−f(x)/f′(x),i=0,1,2,…, che genera una successione xi convergente al punto fisso α di g(x) (α=g(α)) per un'approssimazione iniziale x0 sufficientemente vicina ad α. Per f(x)=x2−N, ove N è un intero positivo, si tratta di una procedura nota fin dalla più remota antichità, in Mesopotamia, in Cina e in India, facilmente deducibile dalla quarta proposizione del II Libro degli Elementi di Euclide, in cui si fa uso di una costruzione geometrica basata sulla figura dello gnomone quadrato. Il metodo di Newton, brevemente trattato da James M. Ortega e Werner C. Rheinboldt in Iterative solution of nonlinear equations in several variables (1970), risulta il caso particolare di una iterazione xi+1=Gxi, dove G è un operatore che associa elementi di uno spazio metrico X a elementi dello stesso spazio. Se con ∥x−y∥ si indica la distanza tra due elementi x e y di X e se ∥Gx−Gy∥≤λ∥x−y∥, con λ reale e 0〈λ〈1, G si chiama operatore di contrazione. Si dimostra allora (teorema di contrazione) che se X è completo (ogni successione di Cauchy di elementi di X converge a un elemento di X) esiste uno e un solo punto α tale che α=G(α) e α si chiama punto fisso di G. Inoltre il processo iterativo xi+1=Gxi, ove x0 è un'approssimazione iniziale di α, produce una successione xi che converge ad α e l'errore al passo k è dato da ∥xk−α∥≤≲λk/≲1−λ))∥x0−Gx0∥.
Il teorema di contrazione vale in un qualsiasi spazio di Banach (metrico e completo) e può quindi essere usato in una varietà di casi che comprende, oltre a singole equazioni non lineari, sistemi di equazioni lineari e non lineari, equazioni funzionali, equazioni integrali e problemi con valori al contorno. A Stefan Banach si deve una prima formulazione e dimostrazione del teorema di contrazione, per spazi lineari metrici e completi. Il teorema è fondamentale per tutto lo sviluppo dell'analisi numerica del Novecento in quanto comprende, accanto al risultato di esistenza e unicità della soluzione di svariati problemi, anche l'indicazione della tecnica iterativa per costruirla in modo effettivo. Il principio di contrazione, formulato nel 1922, può essere considerato una rielaborazione di un'idea di Charles-émile Picard, che aveva usato, intorno al 1890, un metodo di approssimazioni successive per risolvere numericamente un'equazione differenziale lineare del secondo ordine, riprendendo a sua volta una procedura di Cauchy. Tuttavia il modello di riferimento per le iterazioni basate sulla tecnica del punto fisso resta pur sempre la versione che Joseph Raphson diede degli algoritmi di Viète e Newton verso la fine del XVII sec., versione che Jean-Baptiste-Joseph Fourier e Cauchy collegarono poi alle formule e ai concetti fondamentali dell'analisi matematica.
Diversi studi sul metodo di Newton-Raphson e sulla sua estensione al caso di n equazioni in n incognite si susseguirono nel corso del Novecento: la convergenza del metodo nel caso n-dimensionale fu dimostrata da Henry B. Fine, nel 1916, senza assumere l'esistenza della soluzione. Un trattato sulla risoluzione numerica di equazioni algebriche, in cui è considerato il caso n-dimensionale, si deve a Runge (1921) al quale si deve pure uno studio sui punti di attrazione per equazioni in ℝn pubblicato nel 1899.
Tra il 1933 e il 1936 diverse tecniche di valutazione dell'errore e diversi teoremi di convergenza estesi anche al campo complesso si devono a Ostrowski, che è anche autore di un importante trattato sulla risoluzione numerica di equazioni, Solution of equations and systems of equations (1960), frutto di una serie di lezioni tenute a Washington nel 1952, nel quale egli definisce in modo rigoroso il concetto di punto di attrazione per iterazioni in più variabili e dimostra un teorema che può considerarsi un'estensione al caso n-dimensionale del teorema di contrazione per singole equazioni e la specializzazione del principio di contrazione di Banach per sistemi di equazioni non lineari (ove la costante λ di contrazione è sostituita dal raggio spettrale di una matrice). Più precisamente sia dato il sistema yi=gi(x1,x2,…,xn)≡gi(x) o, in notazione vettoriale, y=G(x), sia J(x) la matrice jacobiana definita da [J(x)]ij=gi(x)/xj. Sia inoltre α un punto fisso di G, cioè α=G(α). Allora, se le funzioni gi sono abbastanza regolari e se il raggio spettrale di J(α) è minore di 1, α è un punto di attrazione e la successione xi generata dalla formula iterativa xi+1=G(xi) converge ad α per i che tende a infinito. A Ostrowski si devono anche rielaborazioni del metodo di falsa posizione (regula falsi) e del metodo di Graeffe, già formulato in modo esplicito da Germinal-Pierre Dandelin nel 1826 e da Nikolaj Ivanovič Lobačevskij nel 1834 (ma tradotto in una procedura pratica ed efficiente soltanto da Karl Heinrich Graeffe nel 1837).
Stabilità e complessità
La risoluzione numerica con metodi iterativi di sistemi di equazioni lineari rientra per un verso nella teoria generale del punto fisso, ma richiede di solito le tecniche più specifiche dell'algebra lineare numerica e del calcolo matriciale. Queste tecniche hanno avuto un ruolo determinante nello sviluppo della scienza del calcolo nel Novecento, dal momento che quasi tutti i problemi dell'analisi numerica possono ricondursi, alla fine, alla risoluzione numerica di un sistema di equazioni lineari. Accanto ai metodi iterativi occorre comunque considerare, per i sistemi lineari, i metodi diretti, cioè i metodi che calcolerebbero, in aritmetica esatta, la soluzione esatta in un numero finito di passi.
Un episodio emblematico, che spiega a quale grado di sviluppo era giunta l'algebra lineare numerica nei primi decenni del Novecento, riguarda precisamente i metodi diretti ed è riportato da James Wilkinson, uno dei principali esponenti della scienza del calcolo del XX secolo. Wilkinson doveva la sua formazione di analista alla scuola di Cambridge, dominata dalle figure di Godfrey H. Hardy e John E. Littlewood, e trascorse gli anni della guerra, dal 1940, presso l'Armament Research Department, occupandosi di problemi di balistica, termodinamica degli esplosivi e, incidentalmente, di problemi numerici. Nel 1946 raggiunse la divisione matematica del National Physical Laboratory di Londra, dove ebbe modo di collaborare con Turing. Fu in particolare all'Armament Research Department che gli fu assegnato un compito all'apparenza molto facile: risolvere un sistema di 12 equazioni lineari in 12 incognite. Wilkinson riferisce del suo tentativo di risolvere il sistema con il metodo di Cramer, consigliato dai manuali dell'epoca ‒ tra questi doveva esserci un libro autorevole di Edmund T. Whittaker e George Robinson, The calculus of observations, del quale si stamparono quattro edizioni, l'ultima nel 1944, che dedica appena sette pagine alla risoluzione numerica di sistemi lineari e prende in considerazione il solo metodo di Cramer ‒, e delle difficoltà impreviste che si trovò ad affrontare con la matrice dei 144 coefficienti del sistema. Pochi anni dopo, nel 1953, Forsythe, un altro protagonista della nuova scienza del calcolo, avrebbe osservato come il metodo di Cramer, basato sul calcolo ricorsivo di n+1 determinanti (per n equazioni in n incognite) richiede almeno n! moltiplicazioni per ciascun determinante e, quindi, un tempo di esecuzione altissimo. Forsythe valutava infatti che, per n=26, il tempo di risoluzione dello SWAC, un calcolatore che impiegava un secondo per 2600 moltiplicazioni, sarebbe stato di circa 1017 anni. Wilkinson scelse di conseguenza un'altra strategia, il metodo (diretto) di eliminazione di Gauss, che richiede un tempo di esecuzione incomparabilmente minore. L'episodio mette in evidenza una delle prime questioni che riguardano il calcolo su grande scala, cioè la complessità degli algoritmi. La propagazione dell'errore non figurava ancora tra le questioni più urgenti: Wilkinson si limitava a titolo di semplice precauzione a eseguire i calcoli con 10 cifre decimali, ma non si sentiva assillato da alcuna pessimistica previsione di instabilità numerica.
Dopo pochi anni, tuttavia, furono pubblicate le prime ricerche sistematiche sul comportamento dell'errore. Una delle prime (1943), dovuta a Harold Hotelling, rafforzò la convinzione pessimistica che la risoluzione numerica di un sistema lineare con un metodo diretto avrebbe richiesto un aumento straordinario della precisione di macchina. Nel 1947 fu pubblicato sul "Bulletin of the American mathematical society" Numerical inverting of matrices of high order, lo storico articolo di von Neumann e Goldstine sui limiti dell'errore nell'inversione di una matrice. Questo articolo è il primo a trattare, in modo esplicito, la stabilità numerica. La stabilità era introdotta con qualche esitazione e usata ancora come sinonimo di continuità. Soltanto in seguito divenne un concetto indipendente, che riguardava la crescita dell'errore algoritmico in un processo effettivo di calcolo. Sulle ricerche di von Neumann e Goldstine ebbe anche una notevole influenza un celebre articolo del 1928 di Richard Courant, Kurt Otto Friedrichs e Hans Lewy, dal quale si desumeva che la stabilità di un'equazione alle derivate parziali non implica necessariamente la stabilità di uno schema alle differenze, indipendentemente dalla precisione con cui questo la approssimi. L'articolo di von Neumann e Goldstine riguardava però soprattutto la stabilità del calcolo aritmetico, cioè il comportamento dell'errore algoritmico dovuto alla sostituzione di operazioni aritmetiche esatte con operazioni aritmetiche approssimate (operazioni di macchina) nel calcolo dell'inversa di una matrice A n×n. Il risultato principale poteva descriversi in questi termini: se A è una matrice definita positiva ed effettivamente invertibile (con gli autovalori sufficientemente distanti da zero) e se X è la sua inversa calcolata con il metodo di eliminazione basato sulla classica fattorizzazione A=LDLT (ove L è triangolare inferiore con lii=1 e D è diagonale), allora l'errore misurato dalla distanza tra AX e la matrice identica (in norma spettrale) è superiormente limitato da 14.24(λ/μ)n2B−s, ove λ e μ sono gli autovalori massimo e minimo, rispettivamente, di A e la macchina esegue le operazioni aritmetiche in base B con s cifre significative. Non si comprese subito tutto il significato e l'importanza del rapporto λ/μ. La valutazione dell'errore fu intesa sulle prime in senso pessimistico e per alcuni anni si affermò la convinzione che i metodi iterativi fossero preferibili a quelli diretti. Wilkinson ne fece poi un'analisi più precisa ed esauriente dimostrando che, con qualche precauzione, il metodo diretto di eliminazione con pivot poteva considerarsi stabile e quindi pienamente affidabile e, a questo fine, sull'esempio di Wallace Givens, sviluppò la tecnica dell'analisi all'indietro dell'errore (Wilkinson 1971a). I metodi diretti riacquistarono credito anche in seguito alle ricerche di Givens e di Alston Householder sulle possibili fattorizzazioni di A per via di trasformazioni unitarie. Nel 1958 Householder propose l'uso di matrici simmetriche unitarie (dette matrici di Householder) nella riduzione di matrici simmetriche a forma tridiagonale e a forma triangolare. Queste matrici erano state già studiate da Herbert W. Turnbull e Alexander Aitken nel 1930.
Il rapporto λ/μ prese in seguito il nome di indice di condizionamento, perché misura anche la sensibilità del problema (inversione di una matrice) agli errori sui dati (gli elementi della matrice), indipendentemente dall'algoritmo di risoluzione. L'indice di condizionamento μ(A) di una matrice A fu definito in modo esplicito come il prodotto delle norme (di Frobenius) di A e di A−1, cioè μ(A)=∥A∥∙∥A−1∥, in un importante articolo di Turing, Rounding-off errors in matrix processes (1948), poco dopo il contributo di von Neumann e Goldstine. In norma euclidea lo stesso prodotto è uguale, per matrici definite positive, precisamente al rapporto λ/μ, per cui si può affermare che la prima apparizione dell'indice di condizionamento di una matrice risale al 1947.
La dimensione elevata dei problemi impone di solito, oltre a un'analisi dell'errore, anche uno studio della complessità algoritmica, e von Neumann fu certo tra i primi ad avvertire la necessità di ridurre il costo, in termini di numero di operazioni, delle procedure aritmetiche. Per questo era opportuno introdurre modelli di calcolo il più possibile somiglianti alla macchina di Turing, sebbene calibrati sulla natura specifica del problema numerico. Prima di von Neumann fu Scholtz a introdurre, in una breve nota del 1937, il concetto di catena additiva per un intero n, un concetto utile per analizzare la complessità di calcolo della semplice funzione xn (la catena additiva è una successione di interi che corrisponde alle potenze di x sufficienti a calcolare xn, per es., 2, 3, 4, 7, 8, 15, 16, 31 per n=31). Nel 1954 Ostrowski si chiese per primo se il metodo di Horner-Ruffini per calcolare un polinomio avesse complessità minima. Negli anni Cinquanta venne quindi introdotto un modello per il calcolo di un polinomio noto come catena polinomiale, esaminato inizialmente da Theodore Motzkin in un articolo del 1955. Nel 1965 James W. Cooley e John W. Tukey introducevano la FFT (fast Fourier transforms), l'algoritmo veloce che calcola la trasformata discreta di Fourier su n punti in sole O(nlogn) operazioni aritmetiche. Nel 1969 Volker Strassen pubblicò un celebre articolo in cui si dimostra che, per risolvere un sistema lineare, si possono definire algoritmi di complessità asintotica minore della classica eliminazione gaussiana. Dal 1970, a cominciare da un importante contributo di Shmuel Winograd (1970), si svilupparono ricerche sistematiche sui limiti inferiori di complessità di ampie classi di funzioni razionali, con sorprendenti risultati che riguardavano il calcolo polinomiale e matriciale e, in genere, le espressioni algebriche o aritmetiche elementari a cui si devono ridurre, per ottenere una plausibile informazione numerica, i modelli della Natura.