EANweb - Università degli Studi di Messina

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.