Riccardo Fazio Dipartimento di Matematica e Informatica Universitá degli Studi di Messina mat521.unime.it/~fazio [email protected] Elementi di Analisi Numerica Elementi di Analisi Numerica c MMXIV Riccardo Fazio Copyright ISBN L’utilizzo di denominazioni generiche, nomi commerciali, marchi registrati, ecc, in quest’opera, anche in assenza di particolare indicazione, non consente di considerare tali denominazioni o marchi liberamente utilizzabili da chiunque ai sensi della legge sul marchio. I programmi di questo libro sono stati aggiunti solo a fini didattici. Nonostante essi siano stati testati, su diverse piattaforme e con diversi linguaggi di programmazione, e si sia cercato in ogni modo di verificare che siano esenti da errori l’autore e l’editore non potranno ritenersi responsabili di eventuali danni derivanti dal loro utilizzo. Riprodotto da copia camera-ready fornita dall’Autore. La figura, in copertina, mostra i polinomi di Taylor per l’approssimazione della funzione sin(x). I diritti di traduzione, di memorizzazione elettronica, di riproduzione e di adattamento anche parziale, con qualsiasi mezzo, sono riservati per tutti i Paesi. I edizione: Settembre 2014 Indice Prefazione x I 1 Elementi di Analisi Numerica 1 Numeri di macchina e errori 1.1 Numeri reali e numeri di macchina . . . . . . . 1.1.1 Cambiamenti di base nella numerazione . 1.2 Numeri di macchina . . . . . . . . . . . . . . . 1.3 Errori . . . . . . . . . . . . . . . . . . . . . . . 1.4 Errori nelle operazioni aritmetiche . . . . . . . . 1.5 Programmi FORTRAN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 4 7 9 12 15 2 Problemi e algoritmi 2.1 Condizionamento di un problema . . 2.2 Proprietà di un algoritmo . . . . . . 2.3 Stabilità di un algoritmo . . . . . . . 2.4 Errore totale . . . . . . . . . . . . . . 2.5 Convergenza e ordine di convergenza 2.6 Programmi FORTRAN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 18 21 24 26 29 30 3 Zeri di funzione 3.1 Il caso lineare . . . . . . . . . 3.2 Il caso non lineare . . . . . . . 3.2.1 Condizionamento . . . 3.2.2 Il metodo di bisezione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 33 33 34 35 vi . . . . . . . . . . . . . . . . vii INDICE 3.2.3 Ordine di convergenza . . . 3.2.4 Iterazioni di punto fisso . . . 3.2.5 Metodo di Newton . . . . . 3.2.6 Analisi di convergenza locale 3.3 Metodi quasi-Newton . . . . . . . . 3.3.1 Zeri multipli . . . . . . . . . 3.4 Programmi FORTRAN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Sistemi lineari: metodi diretti 4.1 Esistenza e unicità . . . . . . . . . . . . . . . . . . 4.2 Matrici diagonali . . . . . . . . . . . . . . . . . . . 4.3 Sistemi triangolari . . . . . . . . . . . . . . . . . . 4.4 Il caso generale . . . . . . . . . . . . . . . . . . . . 4.4.1 Il metodo di Cramer . . . . . . . . . . . . . 4.4.2 Metodo di eliminazione di Gauss . . . . . . 4.5 Esempi di sistema lineari . . . . . . . . . . . . . . . 4.6 Fattorizzazione di Cholesky . . . . . . . . . . . . . 4.6.1 Algoritmo di Cholesky con prodotto interno 4.7 Condizionamento . . . . . . . . . . . . . . . . . . . 4.7.1 Approssimazione nei termini noti . . . . . . 4.7.2 Approssimazione negli elementi di A . . . . 4.7.3 Significato del numero di condizionamento . 4.8 Programmi FORTRAN . . . . . . . . . . . . . . . . 5 Sistemi lineari: metodi iterativi 5.1 Metodi iterativi: generalitá . . . . . . 5.2 Metodo di Jacobi . . . . . . . . . . . 5.3 Metodo di Gauss-Seidel . . . . . . . . 5.4 Metodo di sopra-rilassamento (SOR) 5.5 Programmi FORTRAN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 40 43 51 55 57 60 . . . . . . . . . . . . . . 66 66 68 69 75 75 76 79 85 87 89 90 91 94 95 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 101 103 104 104 105 6 Metodi di interpolazione 6.1 Metodi di interpolazione . . . . . . . . . 6.1.1 Polinomio di Vandermonde . . . . 6.1.2 Polinomi di Lagrange ed esistenza 6.1.3 Polinomi di Newton e unicità . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 112 113 113 114 . . . . . viii INDICE 6.2 6.3 Complessità computazionale . Complessità di Ax = b . . . . 6.3.1 Sistemi triangolari . . 6.3.2 Metodo di eliminazione 6.4 Errore dell’interpolazione . . . 6.5 Interpolazione lineare . . . . . 6.6 Programmi FORTRAN . . . . II . . . . . . . . . . . . . . . . . . di Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Formula di Taylor e applicazioni 120 121 121 123 125 127 132 135 7 Derivazione e integrazione 7.1 Formula di Taylor . . . . . . . . . . . . . . . . . . 7.2 Derivazione numerica . . . . . . . . . . . . . . . . 7.2.1 Ordine di una formula alle differenze finite 7.2.2 Estrapolazione di Richardson . . . . . . . 7.2.3 Differenze finite in precisione finita . . . . 7.3 Integrazione numerica . . . . . . . . . . . . . . . 7.3.1 Metodo di Romberg . . . . . . . . . . . . 7.3.2 Formule di quadratura in precisione finita 7.4 Metodi Monte Carlo . . . . . . . . . . . . . . . . 7.4.1 Calcolo del valore di π . . . . . . . . . . . 7.5 Programmi FORTRAN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 137 138 140 147 150 152 158 160 161 162 163 8 Funzioni di Libreria 8.1 Formula di Taylor . . . . . . 8.2 Errore locale e errore globale 8.3 La funzione esponenziale e le 8.4 Programmi FORTRAN . . . 8.5 Funzioni di libreria e grafica . . . . . . . . . . . . . . . 169 169 171 174 175 177 . . . . 179 179 180 181 184 . . . . . . . . . . . . . . . . . . . . . . formule di Eulero . . . . . . . . . . . . . . . . . . . . . . 9 Sistemi non lineari 9.1 Definiamo il problema . . . . . . . . . 9.2 Il metodo di Newton . . . . . . . . . . 9.2.1 Il metodo di Newton per n = 2 9.3 Un metodo quasi-Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix INDICE 9.3.1 Il problema di Simpson . . . . . . . . . . . . . . 184 9.4 Programmi FORTRAN . . . . . . . . . . . . . . . . . . 186 10 Problemi ai valori iniziali 10.1 Problemi ai valori iniziali . . . . . . 10.2 Metodi del primo e secondo ordine 10.2.1 Risultati numerici . . . . . . 10.3 Teorema di convergenza . . . . . . 10.4 Stabilità assoluta . . . . . . . . . . 10.5 Programmi FORTRAN . . . . . . . A Nozioni preliminari A.1 Il principio di induzione . . . . . A.2 Teoremi per funzioni continue . . A.3 Formula di Taylor . . . . . . . . . A.4 Notazione asintotica dell’errore . A.5 Teoremi fondamentali del calcolo A.6 Matrici, vettori e scalari . . . . . A.7 Norme . . . . . . . . . . . . . . . A.7.1 Norme di vettore . . . . . A.7.2 Norme di matrice . . . . . A.8 Matrici di rango uno . . . . . . . A.9 Partizioni di matrici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 191 193 195 196 199 202 . . . . . . . . . . . 204 205 205 207 209 210 213 217 217 219 220 221 B Introduzione al FORTRAN 77 223 B.1 Convenzioni del FORTRAN . . . . . . . . . . . . . . . 224 B.2 Specifiche per le istruzioni eseguibili . . . . . . . . . . . 226 Bibliografia 234 Elenco delle figure 237 Elenco delle tabelle 238 Indice analitico 241 Capitolo 1 Numeri di macchina e errori The road to wisdom ? Well, it’s plain and simple to express: Err and err and err again but less and less and less. —Piet H. Grooks (1966). Citato in [11, p. 187] 1.1 Numeri reali e numeri di macchina Consideriamo il modo con il quale noi umani rappresentiamo i numeri. Un qualsiasi numero decimale viene rappresentato usando le cifre 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ed una convenzione posizionale per la quale interpretiamo la prima cifra a sinistra del punto decimale come unità, la seconda come decine, la terza come centinaia, ecc.; la prima cifra a destra del punto decimale come decimi, la seconda come centesimi, la terza come millesimi, e così via. La notazione posizionale rappresenta 2 1.1. Numeri reali e numeri di macchina 3 una forma abbreviata di scrittura dei numeri. Ad esempio il numero 123.15 ha il seguente significato: 123.15 = 1 · 102 + 2 · 10 + 3 · 100 + 1 · 10−1 + 5 · 10−2 . Lo stesso numero reale si può rappresentare in infiniti altri modi, basta moltiplicare il numero dato per una potenza intera di dieci e spostare opportunamente il punto decimale. Per 123.15 avremo ad esempio · · · = 1231.5 · 10−1 = 123.15 = 12.315 · 101 = 1.2315 · 102 = · · · . Per avere una rappresentazione univoca possiamo usare la notazione scientifica. Dato un x ∈ IR la sua rappresentazione in notazione scientifica è data da x = ±0.c1 c2 c3 · · · ct ct+1 · · · 10e dove 0 ≤ ci ≤ 9 per i = 1, 2, . . . , con c1 6= 0 ed e ∈ Z (Z insieme degli interi relativi) è un opportuno esponente. La base dieci non è la sola nella quale si possono rappresentare i numeri reali. L’uso della base dieci è dovuta al fatto che noi umani abbiamo dieci dita delle mani e che contare sulle dita è stata una pratica comune a molti popoli. Naturalmente, siamo anche usi a considerare altre basi. Ad esempio per le uova usiamo le dozzine, e questo sembra dovuto al fatto che in alcune comunità agricole le uova venivano contate con la sola mano libera (quella che non prendeva l’uovo) usando il pollice per scorrere le falangi delle altre dita della stessa mano che sono appunto in numero di dodici. Misuriamo poi gli angoli e le ore o i minuti in sessantesimi. La base due è invece la più naturale per scrivere i numeri su un elaboratore elettronico. Il motivo è semplice, con transistori e diodi è immediato realizzare due stati distinti corrispondenti alle due cifre binarie 0, 1. Inoltre, con le due cifre binarie si può costruire un corpo completo, nel quale le operazioni di addizione e moltiplicazioni sono definite dalle tabelle: + 0 1 0 0 1 1 1 10 × 0 1 0 0 0 1 0 1 1.1. Numeri reali e numeri di macchina 4 I numeri binari sono quindi formati con le due cifre binarie e una notazione posizionale che rispetta le regole che indicheremo per i numeri decimali. 1.1.1 Cambiamenti di base nella numerazione Qui ci interessa definire le regole per trasformare un numero dalla base due, o da altre basi di interesse, di norma la base ottale o quella esadecimale, nella base dieci e viceversa. Nella base ottale si usano le cifre 0, 1, 2, 3, 4, 5, 6, 7 e nella base esadecimale i simboli 0 , 1 , 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E, F, dove i simboli da A a F hanno i valori decimali da 10 a 15, rispettivamente. La trasformazione alla base dieci è semplice: basta usare la rappresentazione posizionale nelle diverse basi e moltiplicare e sommare per ottenere il valore in base dieci. Il lettore provi a completare l’esempio seguente. Esempio 1 Desideriamo convertire il numero (101.011)2 dalla base due alla base dieci. Traccia: (101.011)2 = 1 × 22 + 0 × 21 + 1 × 20 + 0 × 2−1 + 1 × 2−2 + 1 × 2−3 . Consideriamo un ulteriore esempio. Esempio 2 Convertiamo il numero esadecimale (D.A2C)16 alla base dieci: (D.A2C)16 = 13 × 160 + 10 × 16−1 + 2 × 16−2 + 12 × 16−3 = 13 + 0.625 + 0.0078125 + 0.0029296875 = (13.6357421875)10 . Trasformare un numero dalla base dieci alle altre basi è più complesso. Iniziamo con il considerare la trasformazione nella base due. Opereremo in due fasi: separando la parte intera dalla parte decimale. Per la parte intera operiamo dividendo ricorsivamente per due e prendendo i resti delle divisioni scrivendoli, dall’ultimo ottenuto al primo, da sinistra a destra. 1.1. Numeri reali e numeri di macchina 5 Esempio 3 Supponiamo di voler trasformare in base due il numero (18)10 . Avremo 18/2 = 9 9/2 = 4 4/2 = 2 2/2 = 1 1/2 = 0 resto resto resto resto resto 0 1 0 0 1 e quindi (18)10 = (10010)2 . Si noti come la prima cifra a sinistra del punto decimale, del numero in base due, sarà sempre zero se il numero in base dieci è pari e uno se dispari. Inoltre, l’ultima cifra del numero in base due è sempre uno. Il programma FORTRAN binary.f, nella sezione 1.5 a pagina 15, implementa l’algoritmo descritto per la conversione di un numero intero in base dieci alla base due. Per la conversione della parte decimale operiamo moltiplicando ricorsivamente per due ed escludendo la parte intera dei risultati che riportiamo da sinistra a destra, dopo il punto decimale, nell’ordine con il quale li otteniamo. Esempio 4 Supponiamo di voler trasformare in base due il numero (0.6)10 . Avremo 0.8 × 2 = 1.6 0.6 × 2 = 1.2 0.2 × 2 = 0.4 0.4 × 2 = 0.8 0.8 × 2 parte parte parte parte intera intera intera intera 1 1 0 0 e a questo punto la sequenza si ripete identica, il che significa che abbiamo trovato un numero in base dieci con un numero finito di cifre decimale che nella base due ha una rappresentazione periodica ma infinita. In definitiva (0.8)10 = (0.1100)2 . 1.1. Numeri reali e numeri di macchina 6 In ultima analisi quello che abbiamo descritto non si può dire essere un algoritmo, a meno che non decidiamo a priori di calcolare un numero finito di cifre significative del corrispondente numero binario. Il programma FORTRAN binary2.f, nella sezione 1.5 a pagina 15, implementa la procedura di conversione della parte decimale di un numero reale dalla base dieci alla base due. A titolo di esempio proviamo a convertire 0.2, otterremo X = 0.200000003 0 RESIDUO = 0.400000006 0 RESIDUO = 0.800000012 1 RESIDUO = 0.600000024 1 RESIDUO = 0.200000048 0 RESIDUO = 0.400000095 0 RESIDUO = 0.800000191 1 RESIDUO = 0.600000381 1 RESIDUO = 0.200000763 0 RESIDUO = 0.400001526 0 RESIDUO = 0.800003052 1 RESIDUO = 0.600006104 1 RESIDUO = 0.200012207 0 RESIDUO = 0.400024414 0 RESIDUO = 0.800048828 1 RESIDUO = 0.600097656 1 RESIDUO = 0.200195312 0 RESIDUO = 0.400390625 0 RESIDUO = 0.80078125 1 RESIDUO = 0.6015625 1 RESIDUO = 0.203125 0 RESIDUO = 0.40625 0 RESIDUO = 0.8125 1 RESIDUO = 0.625 1 RESIDUO = 0.25 0 RESIDUO = 0.5 1 RESIDUO = 0. 0.00110011001100110011001101 La conversione dalla base decimale alla ottale o alla esadecimale segue regole simili a quelle usate nel caso della base due. In partico- 1.2. Numeri di macchina 7 lare, divideremo la parte intera e moltiplicheremo la parte decimale, rispettivamente, per otto o sedici. Notiamo poi che la trasformazione da base due a base otto o sedici è immediata se riflettiamo sul fatti che questi numeri sono potenze di due. Basterà, quindi, raggruppare, a partire dal punto decimale, le cifre binarie a gruppi di tre o di quattro a seconda se vogliamo convertire in base ottale o esadecimale, rispettivamente. Un elaboratore elettronico dispone di una memoria limitata e quindi non tutti i valori reali possono essere rappresentati esattamente. Come abbiamo visto, anche dei numeri che nella base dieci hanno un numero finito di cifre decimali devono essere rappresentati in base due con un numero infinito di valori decimali. Se poi l’elaboratore dispone di una rappresentazione dei numeri in base dieci, tutti i numeri reali √ che hanno infinite cifre decimali, ad esempio 2 e π, non saranno rappresentati esattamente. Si propone quindi l’esigenza di definire come approssimare i numeri reali. 1.2 Numeri di macchina Nel seguito, per semplicità, definiremo troncamento e arrotondamento nel caso di numeri positivi. I primi elaboratori, al fine di risparmiare tempo nella conversione dei dati in numeri di macchina adottavano il troncamento che andiamo a descrivere. Dato un x ∈ IR con x = 0.c1 c2 c3 · · · ct ct+1 · · · B e dove 0 ≤ ci ≤ B − 1 per i = 1, 2, . . . , con c1 6= 0, B ∈ IN − {0, 1}, e è un esponente che verifica le condizioni L ≤ e ≤ U, dove L < 0 e U > 0 sono dei valori fissati, si definisce trt (x) con trt (x) = 0.c1 c2 c3 · · · ct B e , dove la parte decimale è detta mantissa, ed e è denominata caratteristica, del numero di macchina. Nel troncamento si usano quindi solo t cifre significative della mantissa. Nei moderni elaboratori si usa l’arrotondamento, si definisce f lt (x) con f lt (x) = 0.c1 c2 c3 · · · cet B e 1.2. Numeri di macchina 8 dove cet = ( ct ct + 1 se ct+1 < B/2 se ct+1 ≥ B/2 , qui f lt (x) rappresenta l’approssimazione per arrotondamento a t cifre decimali di x. I numeri reali rappresentabili su un elaboratore sono, quindi, in numero finito e appartengono all’insieme dei numeri di macchina che indicheremo con IF. Per la precisione IF risulta univocamente individuato se vengono assegnati i valori di B, t, L, U . Useremo, allora la notazione estesa IF(B, t, L, U ) per indicare particolari insiemi di numeri di macchina. Esercizio 1 Determinare, e rappresentare graficamente (vedi la figura 1.1), i valori che costituiscono l’insieme di numeri di macchina IF(2, 2, −2, 2). Traccia. Per fissare le idee, iniziamo considerando solo i valori positivi. Inoltre fissiamo e = 0. In questo caso, poiché c1 6= 0, i due soli numeri binari che si ottengono sono: 0.10 = 2−1 = 1/2 0.11 = 2−1 + 2−2 = 1/2 + 1/4 = 3/4 . Le distanze fra i diversi numeri di macchina sono visibili nella figura (1.1). 0 −3 −2 0.125 1 1.5 2 3 Figura 1.1: Numeri di macchina di IF(2, 2, −2, 2). Notiamo che il numero zero non appartiene all’insieme IF (infatti c1 6= 0), ma viene aggiunto a priori ai numeri di macchina in ogni compilatore. Inoltre, IF ha un numero finito di elementi. Fissato 1.3. Errori 9 l’esponente c1 può assumere B −1 valori distinti, mentre ogni elemento dopo il primo nella mantissa può assumere B valori distinti e in totale avremo (B − 1) B t−1 numeri positivi per ogni valore dell’esponente. A sua volta l’esponente può assumere U − L + 1 valori distinti. A conti fatti, e tenendo conto anche dei valori negativi, l’insieme dei numeri di macchina IF è formato da 2(U − L + 1)(B − 1) B t−1 elementi distinti. Poiché U e L sono finiti non si possono rappresentare in IF numeri reali in valore assoluto arbitrariamente grandi o piccoli. I valori positivi più piccolo e più grande di IF sono: xmin = B L−1 , e xmax = B U 1 − B −t . L’evenienza di un numero compreso nell’intervallo (−xmin , xmin ) viene segnalata come underflow e trattata dal compilatore in modo speciale o in alcuni casi meno validi sostituito con zero. Un risultato all’esterno dell’intervallo [−xmax , xmax ] viene segnalato come overflow e memorizzato nella variabile Inf o in alcuni casi meno validi il compilatore termina il programma. Come insegna il fallimento del lancio, e l’auto-distruzione, del vettore francese Ariane 5 per la messa in orbita di satelliti nel 1996, l’evenienza di un overflow, nelle applicazioni reali, non deve essere sottovalutata. Si veda nel merito la nota di un Anonimo [2]. Le specifiche dei recenti compilatori utilizzano l’ANSI/IEEE Standard 754 descritto nel documento IEEE [13]. Queste norme prevedono una rappresentazione binaria dei numeri reali, in virgola mobile per L = −U , con t = 23 e U = 7 per la precisione semplice, t = 52 e U = 11 in precisione doppia, e, infine, t = 63 e U = 14 in precisione quadrupla. Nella precisione più estesa i valori positivi rappresentabili appartengono all’intervallo [10−4964 , 104964 ]. 1.3 Errori I seguenti tipi di errori sono presenti nelle applicazioni numeriche: 1. errori sui dati, dovuti, ad esempio, agli errori sistematici compiuti da strumenti di misura; 1.3. Errori 10 2. errori di arrotondamento, causati alla rappresentazione finita dei numeri al calcolatore; 3. errori di discretizzazione (o di troncamento), derivanti dalla sostituzione di un problema discreto ad un problema continuo. Gli errori del primo tipo sono inevitabili e possono essere trattati insieme a quelli del secondo tipo. In questa sezione ci occupiamo di valutare la propagazione degli errori del secondo tipo. Gli errori del terzo tipo rivestono una importanza notevole nell’approssimazione delle funzioni, così come nell’Analisi Numerica dei problemi differenziali. Siamo quindi interessati a dare una valutazione quantitativa dell’errore che commettiamo quando arrotondiamo i valori reali nel modo descritto. A tal fine introduciamo l’errore assoluto e l’errore relativo. Sia x ∈ IR e x∗ una sua approssimazione, chiameremo errore assoluto il valore positivo e dato da e = |x∗ − x| . Se x 6= 0 chiameremo errore relativo il valore positivo ex definito da ex = |x∗ − x| . |x| L’errore relativo è in genere più significativo dell’errore assoluto, in quanto tiene conto anche dell’ordine di grandezza di x. Ad esempio, se consideriamo le approssimazioni x∗ = 2 di x = 1 e x∗ = 1001 di x = 1000, avremo in entrambi i casi un errore assoluto pari a uno, ma l’errore relativo nel primo caso è ancora uno mentre nel secondo è un millesimo. Dal punto di vista dell’approssimazione, è evidente che risulta più accettabile approssimare 1000 con 1001 piuttosto che approssimare 1 con 2. Per quanto concerne l’approssimazione dei numeri reali con i numeri di macchina possiamo sempre usare l’errore relativo in quanto 1.3. Errori 11 nella notazione scientifica avremo x 6= 0. Nel caso del troncamento |0.ct+1 ct+2 ct+3 · · · |B e−t |x − trt (x)| = |x| |0.c1 c2 c3 · · · ct ct+1 · · · |B e B −t ≤ B −1 ≤ B 1−t , dove la prima equazione deriva dal fatto che x e trt (x) hanno le prime t cifre della mantissa uguali, abbiamo quindi semplificato il termine comune B e a numeratore e denominatore, maggiorato il numeratore ponendo ci = B − 1 per i ≥ t + 1 e minorato il denominatore ponendo c1 = 1 e ci = 0 per i ≥ 2. Nel caso dell’arrotondamento l’errore relativo massimo commesso sarà ovviamente un mezzo dell’errore relativo massimo trovato nel caso del troncamento, cioè |x − f lt (x)| B 1−t ≤ . |x| 2 (1.1) Questo errore relativo massimo viene indicato con ǫM e viene denominato precisione di macchina. Poiché ǫM è una maggiorazione dell’errore relativo commesso prescinde dall’ordine di grandezza e quindi dall’esponente. Questo spiega perché la precisione di macchina dipende solo dalla base B e dalla mantissa t, ma non dai limiti per l’esponente L e U . Si può verificare che il valore della precisione di macchina ǫM è la metà della distanza fra uno ed il numero di macchina successivo a uno. In pratica la precisione di macchina si può calcolare su un dato elaboratore che usa un particolare compilatore con un semplice algoritmo, si veda il programma FORTRAN epsm2.f nella sezione 1.5 a pagina 15. Ad esempio su una SUN ultra 5 workstation, che utilizza un compilatore FORTRAN 77, si sono ottenuti i valori di ǫM in precisione semplice (indicata con il simbolo E), doppia (simbolo D) e quadrupla (simbolo Q) riportati nella tabella 1.1. In MATLAB o Octave per ottenere il valore di ǫM basta scrivere sulla linea di comando eps. Nel caso dell’insieme di numeri di 1.4. Errori nelle operazioni aritmetiche 12 Precisione Semplice Doppia Quadrupla REAL REAL*8 REAL*16 ǫM ≤ 6.0E−08 1.2D−16 1.0Q−34 Tabella 1.1: Precisione di macchina in FORTRAN 77. macchina IF(2, 2, −2, 2) la precisione di macchina ǫM ha il valore ǫM = 21−2 1 = , 2 4 proprio metà della distanza fra 3/2 (numero di macchina successivo a uno) e 1. 1.4 Errori nelle operazioni aritmetiche In questa sezione vedremo come gli errori introdotti sui dati si propagano nel calcolo delle operazioni aritmetiche: somma e prodotto. Sottrazione, elevazione a potenza e divisione possono essere trattate insieme alle operazioni di somma e prodotto in quanto opposti, potenze e inversi di elementi di IR sono sempre elementi di IR. Notiamo, in via preliminare, che dato x ∈ IR − {0} e una sua approssimazione con un numero di macchina x∗ , posto ρx = (x−x∗ )/x, si verifica per sostituzione che x∗ si può esprimere come x∗ = x (1 + ρx ) , e inoltre, per le proprietà dei numeri di macchina, |ρx | ≤ ǫM , dove ǫM è la precisione di macchina. Considerando, quindi, le operazioni di somma e prodotto fra due valori approssimati x∗ e y ∗ di IR, avremo x∗ + y ∗ = x (1 + ρx ) + y (1 + ρy ) , x∗ y ∗ = x (1 + ρx ) y (1 + ρy ) , dove, naturalmente, si ha |ρy | ≤ ǫM . A conti fatti vorremmo poter scrivere (x + y) (1 + ρx+y ) e (x y) (1 + ρx y ) con |ρx+y | ≤ ǫM e |ρx y | ≤ 1.4. Errori nelle operazioni aritmetiche 13 ǫM . Con semplici passaggi algebrici, e applicando la diseguaglianza triangolare, otteniamo e y x |ρ | + |ρ | , |ρx+y | ≤ x + y y x + y x |ρx y | ≤ |ρx | + |ρy | + |ρx ρy | . Nel caso del prodotto , poiché |ρx ρy | ≤ ǫM 2 ≪ ǫM otteniamo che |ρx·y | ≤ 2 ǫM . Per la somma , la situazione è più complessa e, dovremo distinguere diversi casi: 1. se x e y sono concordi, cioè hanno lo stesso segno, allora 0 ≤ |x/(x + y)| < 1 e 0 ≤ |y/(x + y)| < 1, e quindi |ρx+y | ≤ 2 ǫM ; 2. se x e y sono discordi e |x| ≪ |y| oppure |y| ≪ |x|, allora i due rapporti |x/(x + y)| e |y/(x + y)| saranno uno circa uno e l’altro prossimo a zero, e quindi |ρx+y | ≤ ǫM ; 3. infine se x e y sono discordi e |x| ≈ |y|, allora i denominatori dei rapporti |x/(x + y)| e |y/(x + y)| sono prossimi a zero, e quindi si potrebbe verificare il caso che |ρx+y | ≫ ǫM . Possiamo affermare che le operazioni aritmetiche propagano l’errore nei dati in modo ragionevole tranne nel caso di somme (sottrazioni) di numeri di segno opposto (dello stesso segno) che siano dello stesso ordine di grandezza. In questo ultimo caso si ha una perdita di accuratezza sui risultati delle elaborazioni e si presenta il rischio di ottenere degli underflow che nei calcoli successivi potrebbero portare a un overflow. Inoltre, alcune delle proprietà verificate dai numeri reali non sono più valide per i numeri di macchina. Ad esempio, i numeri di macchina, a causa della precisione finita nelle operazioni aritmetiche, non verificano la proprietà associativa. Esempio 5 Fissato x = 1 calcolare e stampare i valori ((x+1)−1)/x, x (x + (1 − 1))/x e porre x = 10 fin quando x > 0. 1.4. Errori nelle operazioni aritmetiche 14 x ((x + 1) − 1)/x (x + (1 − 1))/x 1. 0.100000001 0.0100000007 0.00100000005 0.000100000005 1.00000007E − 05 ... 1.00000026E − 11 1.00000032E − 12 1.00000032E − 13 1.00000032E − 14 1.00000032E − 15 1.00000035E − 16 1.0000004E − 17 1.00000046E − 18 1.00000049E − 19 1.00000053E − 20 1.00000057E − 21 ... 9.80908925E − 45 1.40129846E − 45 1. 1. 1. 1. 1. 1. ... 1. 1. 1.0000006 1.00000274 0.99995935 0.999634027 0.997465611 0.9757815 1.08420169 0. 0. ... 0. 0. 1. 1. 1. 1. 1. 1. ... 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. ... 1. 1. Tabella 1.2: Numeri di macchina e proprietà associativa. precisione nelle cifre significative. Perdita di Con un compilatore FORTRAN si ottengono i risultati parzialmente riportati nella tabella 1.2. I risultati ottenuti permettono anche di constatare, nella colonna centrale della tabella 1.2, la perdita di cifre significative quando si effettuano operazioni in virgola mobile. 1.5. Programmi FORTRAN 1.5 15 Programmi FORTRAN Consideriamo, per primo, il programma di conversione da numeri interi in base dieci ai corrispondenti interi in base due: binary.f C C PROGRAM BINARY CALCOLO DEL VALORE BINARIO DI UN INTERO Riccardo Fazio, September 12, 1995 CHARACTER*1 CIFRA CHARACTER*20 BNUMERO BNUMERO = ’ ’ PRINT*, ’DIGITA UN NUMERO INTERO POSITIVO’ READ (5, *) N NN = N DO WHILE (.NOT.(NN.EQ.0)) IF (NN.EQ.2*(NN/2)) THEN CIFRA = ’0’ PRINT*, ’0’ ELSE CIFRA = ’1’ PRINT*, ’1’ ENDIF BNUMERO = CIFRA//BNUMERO NN = NN/2 END DO PRINT*, ’IL BINARIO DI ’, N, ’ E’’ ’, BNUMERO END Per compilare il programma binary.f nell’ambiente cygwin scriveremo sulla linea dei comandi g77 binary.f oppure f77 binary.f. Se la compilazione va a buon fine senza la segnalazione di errori, il compilatore produrrà il file eseguibile a.exe; per l’esecuzione del programma servirà scrivere, sulla linea dei comandi, ./a (senza punto finale). In caso di errori, segnalati dal compilatore, tutti gli errori trovati dovranno necessariamente essere corretti prima di poter eseguire il programma. Riportiamo, quindi, il programma di conversione della parte decimale per numeri reali dalla base dieci alla base due: binary2.f 1.5. Programmi FORTRAN C C 16 PROGRAM BINARY2 CONVERSIONE IN BINARIO DI X CON 0 < X < 1 Riccardo Fazio, March 15, 2011 CHARACTER*1 CIFRA, ZERO, POINT CHARACTER*32 CNUMERO, BNUMERO BNUMERO = ’ ’ CNUMERO = ’ ’ PRINT*, ’DIGITA X CON 0 < X < 1’ READ (5,*) RN PRINT*, RN L = 0 DO WHILE ((RN.GT.0).OR.(L.EQ.30)) RN = 2*RN N = RN IF (N.EQ.1) THEN CIFRA = ’1’ ELSE CIFRA = ’0’ END IF RN = RN-N PRINT*, CIFRA,’ RESIDUO = ’, RN CNUMERO = CIFRA//CNUMERO L = L+1 END DO BNUMERO = CNUMERO(1:1) DO K = 2,L BNUMERO = CNUMERO(K:K)//BNUMERO END DO ZERO = ’0’ POINT = ’.’ BNUMERO = POINT//BNUMERO BNUMERO = ZERO//BNUMERO PRINT*, BNUMERO END Per il calcolo della precisione di macchina ǫM , con un compilatore FORTRAN, possiamo usare il seguente programma epsm2.f: 1.5. Programmi FORTRAN C C C C C PROGRAM EPSM2 CALCOLO DELLA PRECISONE DI MACCHINA Riccardo Fazio, January 14, 1999 revised: April 20, 2000 and March 23, 2007 REAL EPS,TEST REAL*8 EPS,TEST REAL*16 EPS,TEST EPS= 1. TEST = EPS+1. DO WHILE (TEST.GT.1.) EPS = .5*EPS TEST = EPS+1. END DO PRINT*, ’PRECISIONE DI MACCHINA = ’, EPS END Per l’esempio 5 utilizziamo il programma associativa.f: C C C C PROGRAM ASSOCIATIVA RICCARDO FAZIO, 1 MAY 2013 ESEMPIO PER LA PERDITA DI CIFRE SIGNIFICATIVE CALCOLO ((X+1)-1)/X E (X+(1-1))/X REAL X, Y1, Y2 REAL*8 X, Y1, Y2 X = 1. DO WHILE (X.GT.0) Y1 = ((X+1.)-1.)/X Y2 = (X+(1.-1.))/X PRINT*, X, Y1, Y2 X = .1*X END DO END 17 Appendice B Introduzione al FORTRAN 77 Il Fortran (o FORTRAN: il nome in tutte maiuscole è stato sostituito da quello con la sola iniziale maiuscola a partire dal Fortran 90), è uno dei primi linguaggi di programmazione, essendo stato sviluppato, a partire dal 1954 ... —http://it.wikipedia.org/ Per i programmi FORTRAN e gli esempi proposti in questo libro abbiamo usato il compilatore FORTRAN g77, o f77, disponibile gratuitamente dal sito |http://www.cygwin.com|. Il cygwin dispone, alla directory math, anche del software Octave. Come indicato nella prefazione, un’altra distribuzione di Octave è disponibile dal sito web |http://octave.sourceforge.net/|. Il manuale, in inglese, della versione attuale di Octave, la 3.2.3, si trova alla sotto-directory doc/octave/pdf. 223 B.1. Convenzioni del FORTRAN B.1 224 Convenzioni del FORTRAN Un programma FORTRAN deve rispettare alcune regole imposte dal compilatore. La prima regola prevede di specificare al compilatore il nome del file contenente il programma. Ogni codice FORTRAN inizia con la parola riservata, che quindi non può essere usata come nome di variabile, PROGRAM seguita dal nome del programma e si chiude con la parola (riservata) END: PROGRAM <nome> {<istruzioni di specifica>} {<istruzioni eseguibili>} END Ogni linea, di un codice FORTRAN, deve essere scritta rispettando la convenzione sul contenuto delle colonne riassunta nella tabella B.1. Colonna 1 1-5 6 7-72 73-80 Contenuto C or * per una linea di commento Numero della label (intero) & o altro simbolo di continuazione di linea Istruzioni FORTRAN Riservate al compilatore Tabella B.1: Convenzione del Compilatore FORTRAN sul contenuto delle colonne. I sottoprogrammi, SUBROUTINE, del FORTRAN hanno la forma seguente SUBROUTINE <nome sottopro>[<lista variabili locali>] {<istruzioni di specifica>} {<istruzioni eseguibili>} RETURN END B.1. Convenzioni del FORTRAN 225 La parola RETURN, presente almeno una o più volte all’interno della SUBROUTINE, rimanda alla linea successiva alla chiamata del programma principale. Le chiamate dei sottoprogrammi vengono invocate con il comando CALL CALL <nome>[<lista delle variabili>] seguito dal nome del sottoprogramma e da una lista di variabili. Le funzioni FUNCTION [<tipo>] FUNCTION <nome> [(<lista variabili locali>)] {<istruzioni di specifica>} {<istruzioni eseguibili>} RETURN END Variabili. Le variabili del FORTRAN possono essere formate da uno fino a sei caratteri alfanumerici. Il primo carattere deve essere alfabetico. Le variabili che iniziano con I, J, K, L, M, N sono considerate, implicitamente, dal compilatore come intere a meno che non siano dichiarate altrimenti. Implicitamente, il compilatore assume che variabili che iniziano per altre lettere dell’alfabeto siano reali. L’istruzione IMPLICIT NONE, usata dai migliori programmatori e da inserire fra le specifiche, sopprime la definizione implicita delle variabili e costringe il programmatore a dichiarare tutte le variabili usate nel codice. Le istruzioni di specifica, dichiarano il tipo di variabili o parametri di un codice FORTRAN: DOUBLE PRECISION IMPLICIT NONE INTEGER RIGHA,COLONNA REAL AB,BB,MAX,X,Y LOGICAL A,B,C COMPLEX Z,W0 CHARACTER*10 NOME IMPLICIT F,DF PARAMETER (PI=3.1415926535) B.2. Specifiche per le istruzioni eseguibili 226 DATA E/2.178281828/ COMMON AA,BB COMMON /Nome-blocco/ X,Y, PI DIMENSION V(1:50), M(1:10,1:10) REAL V(1:50), M(1:10,1:10) B.2 Specifiche per le istruzioni eseguibili Operazioni aritmetiche: + * / ** Addizione Sottrazione Moltiplicazione Divisione Elevazione a potenza Operatori relazionali: .EQ. .NE. .LT. .GT. .LE. .GE. Uguale a Non uguale a Minore di Maggiore di Minore di o uguale a Maggiore di o uguale a Operatori logici: .NOT. .AND. .OR. Complemento Vero se entrambi gli operandi sono Veri Vero se almeno uno degli operandi è Vero Costanti logiche: .TRUE. .FALSE. Funzioni matematiche: B.2. Specifiche per le istruzioni eseguibili COS(X) SIN(X) TAN(X) EXP(X) ACOS(X) ASIN(X) ATAN(X) ALOG(X) LOG10(X) SQRT(X) ABS(X) INT(X) FLOAT(I) REAL(I) DBLE(X) CMPLX(X) 227 coseno (radianti) seno (radianti) tangente (radianti) esponenziale exp(x) arco-coseno (radianti) arco-seno (radianti) arco-tangente (radianti) logaritmo naturale in base e logaritmo comune in base 10 radice quadrata valore assoluto conversione a intero conversione a reale conversione a reale conversione alla precisione doppia conversione a numero complesso Istruzioni di assegnazione: X = A(1) + 3 Y = A*X**2 + B*X + C Flag = .TRUE. NAME = ’Riccardo’ NB: usare le semplici operazioni elementari con variabili di tipo diverso può portare a delle spiacevoli sorprese, ma può anche essere usato nei nostri algoritmi. Ad esempio, la divisione fra due interi è un intero. Abbiamo usato questa convenzione del FORTRAN nella conversione di un intero decimale nel corrispondente valore binario, infatti il prodotto per due della divisione per due di un intero ci ha permesso di decidere se l’intero è pari o dispari. Al fine di evitare errori si consiglia di usare il punto decimale nell’assegnazione dei valori reali e di usare le funzioni di conversione REAL, INT, DBLE, e CMPLX. Salto incondizionato GOTO 100 Salto condizionato al valore di un intero, ad esempio ISALTO, B.2. Specifiche per le istruzioni eseguibili 228 GOTO(100,200,300,400), ISALTO NB: non usare mai un salto incondizionato GOTO per saltare all’interno del blocco di istruzioni del tipo DO, IF, ELSE, o ELSEIF dall’esterno del blocco. Etichetta (label) per la continuazione o il salto (condizionato o incondizionato) 100 CONTINUE Istruzioni condizionali: • IF (espressione aritmetica) Permette il salto a etichette diverse per valori positivi, negativi o zero dell’(<espressione aritmetica>). Ad esempio: IF (<espressione aritmetica>) 100, 200, 300 • IF (espressione logica) esegue un’istruzione singola se l’(<espressione logica>) è vera. Ad esempio: IF (<espressione-logica>) GOTO 100 IF (<espressione-logica>) WRITE (*,*) ’SI ’ IF (<espressione-logica>) X = A+B • IF (<espressione-logica>) (blocco di istruzioni) esegue il blocco di {<istruzioni eseguibili>} che seguono o trasferisce il controllo alla linea successiva a ELSEIF, ELSE, o ENDIF a seconda del valore dell’(<espressione logica>). Blocco IF, THEN, e ENDIF: IF (<espressione logica>) THEN {<istruzioni eseguibili>} ENDIF Blocco IF, THEN, ELSE, e ENDIF: B.2. Specifiche per le istruzioni eseguibili IF (<espressione logica>) THEN {<istruzioni eseguibili A>} ELSE {<istruzioni eseguibili B>} ENDIF Blocco IF, THEN, ELSEIF, ELSE, e ENDIF: IF (<espressione logica #1>) THEN {<istruzioni eseguibili A>} ELSEIF (<espressione logica #2>) THEN {<istruzioni eseguibili B>} ELSEIF (<espressione logica #3>) THEN {<istruzioni eseguibili C>} ELSE {<istruzioni eseguibili>} ENDIF Esempi. Blocco IF, THEN, e ENDIF: IF (ABS(P3-P1).LT.ABS(P3-P0)) THEN U=P1; P1=P0; P0=U V=Y1; Y1=Y0; Y0=V ENDIF Blocco IF, THEN, ELSE, e ENDIF: IF (Df.EQ.0) THEN Dp=P1-P0 P1=P0 ELSE Dp=Y0/Df P1=P0-Dp ENDIF 229 B.2. Specifiche per le istruzioni eseguibili 230 Blocco IF, THEN, ELSEIF, ELSE, e ENDIF: IF (YC.EQ.0) THEN A=C B=C ELSEIF ((YB*YC).GT.0) THEN B=C YB=YC KR=KR+1 ELSE A=C YA=YC KL=KL+1 ENDIF Istruzioni di iterazione: • DO (Blocco) Per incrementare il contatore di uno, come in K = M1, M1+1, ... , M2, scriveremo DO K = M1, M2 {<istruzioni eseguibili>} END DO NB: M1 e M2 devono essere definiti nel codice che precede l’istruzione DO. Invece, per incrementare il contatore di Mstep, cioè K = M1, M1+Mstep, M1+2 Mstep, ... , M2, porremo DO K = M1, M2, Mstep {<istruzioni eseguibili>} END DO NB: Mstep deve essere definito nel codice che precede l’istruzione DO. Esempi: riportiamo qui di seguito dei semplici codici di esempio. B.2. Specifiche per le istruzioni eseguibili SUM = 0 DO K=0,100,2 SUM = SUM + K END DO SUM = 0 DO K=100,1,-1 SUM = SUM + 1.0/REAL(K) END DO DO J=1,ADIM DO K=1,ADIM A(J,K) = 1.0/FLOAT(1+J+K) END DO END DO SUM = 0 DO K=1,10000 SUM = SUM + 1.0/REAL(K) IF (SUM.GT.5) RETURN END DO • DO WHILE (Blocco) DO WHILE (<espressione logica>) {<istruzioni eseguibili>} END DO Esempio di ciclo DO WHILE. SUM = 0 DO WHILE (K.LT.10000) SUM = SUM + 1.0/REAL(K) IF (SUM.GT.5) EXIT K = K+1 END DO 231 B.2. Specifiche per le istruzioni eseguibili 232 Istruzioni di ingresso (Input) e uscita (Output). READ (con o senza FORMAT) per la lettura: C N READ *,<lista dei dati> READ <formato dei dati>,<lista dei dati> READ <N>,<lista dei dati> dove <N> è l’etichetta di un FORMAT FORMAT(<formato dei dati>) PRINT e WRITE (con o senza FORMAT) per la scrittura: C N PRINT * PRINT *,<lista dei dati> PRINT <formato dei dati>,<lista dei dati> PRINT <N>,<lista dei dati> dove <N> è l’etichetta di un FORMAT FORMAT(<formato dei dati>) C N WRITE (*,*) WRITE (*,*) <lista dei dati> WRITE (*,N) <lista dei dati> dove <N> è l’etichetta di un FORMAT FORMAT(<formato dei dati>) Quando, nella risoluzione di una classe di problemi, i dati possono essere di notevole quantità è bene prevedere la lettura da, o la scrittura su, file esterni in formato ASCII . Ad esempio, nella soluzione dei sistemi lineari, nella sezione 4.8 a pagina 95, abbiamo previsto la lettura dei dati del sistema dal file esterno dati.txt, con le istruzioni FORTRAN OPEN READ READ READ READ READ (7,file=’dati.txt’,status=’old’) (7,*) (7,*) (7,*) N (7,*) (7,*) ((A(J,I), I=1,N), J=1,N) B.2. Specifiche per le istruzioni eseguibili 233 READ (7,*) READ (7,*) (B(J), J=1,N) CLOSE (7) Lo status del file esterno deve essere specificato con old per un file esistente oppure new per un file di risultati che si vuole creare al momento. Istruzione di pausa PAUSE: PAUSE Istruzione di stop STOP: STOP Bibliografia [1] M. Abramowitz and I. A. Stegun, editors. Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables. National Bureau of Standards, Washington D.C., 1964. [2] Anonimo. Inquiry board traces Ariane 5 failure to overflow error. SIAM News, 29(8):1, 12–13, 1996. [3] E. Barbin, J. Borowczyk, J. Chambert, M. Guillemot, A. MichelPajus, L. C. Bernard, A. Djebbar, and J. Martzloff. A History of Algorithms. From the Pebble to the Microchip. Springer, Berlin, 1999. [4] V. Comincioli. Analisi Numerica: Metodi, Modelli, Applicazioni. McGraw-Hill, Milano, 1995. [5] V. Comincioli. Analisi Numerica: Metodi, Modelli, Applicazioni. Apogeo, Milano, 2005. [6] G. Dalquist. A special stability problem for linear multistep methods. BIT, 3:27–43, 1963. [7] C. H. Edwards and D. E. Penney. Differential Equations and Boundary Value Problems. Prentice Hall, Englewood Cliffs, 1996. [8] L. V. Fausett. Applied Numerical Analysis Using MATLAB. Prentice Hall, Upper Saddle River, 1999. [9] W. Gautschi. Numerical Analysis. An Introduction. Birkhauser, Boston, 1997. 234 BIBLIOGRAFIA 235 [10] H. H. Goldstine. A History of Numerical Analysis from the 16th through the 19th Century. Springer-Verlag, New York, 1977. [11] D. J. Higham and N. J. Higham. MATLAB Guide. Philadelphia, 2000. SIAM, [12] N. J. Higham. Handbook of Writing for the Mathematical Sciences. SIAM, Philadelphia, 2a edizione, 1998. [13] IEEE. IEEE standard for binary floating-point arithmetic, ANSI/IEEE standard 754-1985. Institute of Electrical and Electronics Engineers, New York, 1985. [14] Z. Kowalik and T. S. Murty. Numerical Modeling of Ocean Dynamics. World Scientific, Singapore, 1993. [15] J. D. Lambert. Numerical Methods for Ordinary Differential Systems. Wiley, Chichester, 1991. [16] C. R. MacCluer. Industrial Mathematics. Modeling in Industry, Science, and Government. Prentice Hall, Upper Saddle River, 2000. [17] G. Monegato. Fondamenti di Calcolo Numerico. Levrotto & Bella, Torino, 1990. [18] S. G. Nash. A History of Scientific Computing. ACM Press, New York, 1990. [19] G. Prodi. Analisi Matematica. Boringhieri, Torino, 1970. [20] A. Quarteroni, R. Sacco, and F. Saleri. Matematica Numerica. Springer-Verlag Italia, 3a edizione, Milano, 1998. [21] M. Rosati. Lezioni di Geometria. Cortina, Padova, 3a edizione, 1980. [22] L.F. Shampine. Numerical Solution of Ordinary Differential Equations. Chapman & Hall, New York, 1994. BIBLIOGRAFIA [23] G. W. Stewart. Afternotes on Numerical Analysis. Philadelphia, 1996. 236 SIAM, [24] J. Stewart. Calculus. Books/Cole, Pacific Grove, 4a edizione, 1999. [25] C. F. Van Loan. Introduction to Scientific Computing, a MatrixVector Approach Using MATLAB. Prentice Hall, Upper Saddle River, 1997.
© Copyright 2024 ExpyDoc