Calcolo Numerico e programmazione - Appunti

Funzionamento del calcolatore

L'elemento alla base del funzionamento di un calcolatore è la CPU – Central Processing Unit – che lavora secondo un ciclo Fetch – Decode – Execute:

• Fetch: la CPU preleva da una memoria i dati da elaborare
• Decode: la CPU decodifica il dato
• Execute: la CPU esegue, se possibile, l'operazione richiesta

Esempio: Si vuole eseguire il calcolo “2+3=5”.

Esempio di ciclo CPU

1° Ciclo:

• Fetch: la CPU preleva il dato “2”
• Decode: la CPU traduce il dato nel numero “2”
• Execute: non si hanno tutti i dati disponibili, quindi non avviene la fase Execute

2° Ciclo:

• Fetch: la CPU preleva il dato “3”
• Decode: la CPU traduce il dato nel numero “3”
• Execute: non si hanno tutti i dati disponibili, quindi non avviene la fase Execute

3° Ciclo:

• Fetch: la CPU preleva l'istruzione “+”
• Decode: la CPU traduce il dato nel termine “+”
• Execute: la CPU esegue l'istruzione e fornisce il risultato “5”

Il ciclo FDE è denominato “ciclo di clock”. Il ciclo di clock è un'operazione continua, perciò, se non avviene una delle fasi del ciclo, questo ricomincia.

La velocità del ciclo di clock è la velocità del processore (es. 1 GHz = 10 cicli/secondo). Per consentire il salvataggio dei risultati delle operazioni nella memoria, la fase di Fetch è bidirezionale. Per memorizzare i dati, la CPU ricorre a sequenze di numeri in linguaggio binario: una cifra binaria è detta bit.

Calcolo della memoria

• 8 bit = 1 byte (B)
• 1024 B = 1 KB
• 1024 KB = 1 MB
• 1024 MB = 1 GB
• 1024 GB = 1 TB

I calcolatori lavorano su “pacchetti di bit” di dimensioni variabili – es. i calcolatori più vecchi lavorano a pacchetti da 32 bit, i moderni calcolatori a pacchetti da 64 bit, i telefoni cellulari lavorano a pacchetti tra i 12 e i 24 bit.

Memoria

La prima memoria è detta RAM – Random Access Memory – e la CPU interagisce con essa tramite i cosiddetti “bus”.

Problemi della RAM

Problema 1: i bus lavorano a una frequenza compresa tra i 200 e 400 MHz, perciò la CPU non può lavorare alla sua velocità effettiva. Inoltre, produrre bus capaci di lavorare alla velocità della CPU richiede costi di produzione non indifferenti. Per risolvere il problema in maniera più economica, alla memoria RAM si affianca un'ulteriore memoria, di tipo RAM, più piccola ma più veloce, detta “memoria cache”, dotata di algoritmi di previsione delle operazioni: attraverso i buffer della memoria cache, la CPU può lavorare ad una velocità più prossima alla sua velocità effettiva.

Problema 2: la memoria RAM è una componente elettronica, quindi non è capace di mantenere i dati senza una continua alimentazione. Per motivi pratici, quindi, alla memoria di tipo RAM si affianca una seconda memoria, più lenta ma indipendente dall'elettricità, poiché magnetica: questa memoria è detta Hard Disk.

L'architettura del calcolatore che si delinea è detta Architettura di Van Neumann.

Memorizzazione dei numeri

A 32 bit

Ipotesi 1: il pacchetto da 32 bit viene diviso in due blocchi, ossia un blocco da 1 bit per il segno e un blocco da 31 bit per il numero. In questo modo, il numero più alto descrivibile è: N_max = 2³¹ - 1 = 2,147,483,647.

A 64 bit

Ipotesi 1: il pacchetto da 64 bit viene diviso in due blocchi, ossia un blocco da 1 bit per il segno e un blocco da 63 bit per il numero. In questo modo, il numero più alto descrivibile è: N_max = 2⁶³ - 1 = 9,223,372,036,854,775,807.

Problemi di memorizzazione

Problema 1: se si lavora con numeri prossimi a N_max, il calcolatore è soggetto al cosiddetto “Integer Overflow”.

Esempio: a 32 bit, sommiamo N_max e N = 1: il risultato vero sarebbe N_max + 1, ma per il calcolatore il risultato è 0.

Problema 2: memorizzando ogni cifra, compresi gli zeri prima della prima cifra non nulla, si memorizzano informazioni inutili, con conseguente spreco di memoria.

Esempio: (-1)₃₂ = 10000000000000000000000000000001

Si utilizza, quindi, la Notazione in Virgola Mobile Normalizzata o “floating point”: in base 2, f(b) = ±1.m·2^f, dove m è detto “mantissa” e 1 è il cosiddetto “hidden bit”.

A 32 bit, i dati sono disposti in tre blocchi, di cui uno da 1 bit per il segno, uno da 8 bit per l'esponente f, e uno da 23 bit per la descrizione del numero 1.m. A 64 bit, i dati sono disposti in tre blocchi, di cui uno da 1 bit per il segno, uno da 11 bit per l'esponente f, e uno da 52 bit per la descrizione del numero 1.m.

Problema: f può essere negativo, ma tra i bit dedicati all'esponente non vi è il bit del segno. Si utilizza quindi il “biased exponent”: poiché a 32 bit f è compreso tra 0 e 255, e a 64 bit f è compreso tra 0 e 2047, si pone f = 127 trasposto, f_{32 reale} = 0, f_{64 reale} = 1023.

A 32 bit, il numero massimo e il numero minimo descrivibili diventano quindi: N_max ≈ 1.7·10³⁸, N_min ≈ 10^-38.

A 64 bit, il numero massimo e il numero minimo descrivibili diventano: N_max ≈ 3·10³⁰⁸, N_min ≈ 10^-308.

Struttura dell'elaboratore

Cos'è un file?

Un file è un insieme di indirizzi di celle di memoria, che si presenta nella forma: nome.estensione dove “.estensione” fa riferimento a un preciso programma di lettura.

Esempio:

• nome.doc → Word
• nome.pdf → Adobe Reader

Cos'è una directory?

Una directory è una collezione di file, ossia un file che contiene una lista di altri file.

Struttura delle directory – File System

La struttura delle directory è una struttura ramificata detta File System, leggibile tramite appositi programmi di lettura.

Esempio:

• Microsoft Windows → Gestione Risorse
• Fedora → Comandi al Terminale

(In Fedora) i simboli “.” e “..” sono interpretati dal calcolatore come “directory corrente” e “directory precedente”.

Compilazione in Fedora

Un file del tipo: nome.f è un codice scritto in linguaggio FORTRAN tramite un apposito editor di testo (Gedit).

Esempio di file: program Salutami
write(*,*) 'Buongiorno!'
stop
end

Un file nome.f non è eseguibile: per questo, deve essere convertito in linguaggio binario da un compilatore, ossia un programma che analizza il file nome.f, lo converte e lo trascrive in un nuovo file eseguibile del tipo: nome.out.

In Fedora, la compilazione viene eseguita tramite il comando:

gfortran file.f -o file.out

Successivamente, si potrà eseguire il file compilato tramite il comando:

./file.out

Attraverso il quale la CPU preleverà file.out (Fetch), lo tradurrà in formato ASCII (Decode) e mostrerà i risultati dell'esecuzione sullo schermo (Execute).

Esempio:

~]$ ./file.out
Buongiorno!

Per compilare file.f, il compilatore lavora in due fasi:

• Compilazione: viene generato il file file.o, denominato “object file”
• Link: collega file.o e file.f, prelevando le istruzioni necessarie dal file libfortran.a per produrre l'eseguibile file.out

Precisione di macchina

Nel momento in cui si descrive un numero nella forma esponenziale f(b) = ±1.m·2^f, si commette un errore in m, poiché il calcolatore lavora a precisione finita: il modulo dell'errore è:

ε = 2-(s+1), s ∈ ℕ

Dove s è il numero di cifre significative.

A 32 bit si ha ε₃₂ = 10^-8, mentre a 64 bit si ha ε₆₄ = 10^-16.

Risoluzione di equazioni non lineari – Introduzione

Quando si risolve un problema, poiché ε ≠ 0, è necessario valutare se l'errore è tale da distruggere l'accuratezza del conto, o se invece è possibile considerare i risultati esatti entro il limite dell'errore.

Stabilità e condizionamento

Definizione: uno schema/algoritmo si definisce stabile se la propagazione dell'errore non è distruttiva dell'accuratezza del calcolo.

Definizione: uno schema/algoritmo si definisce ben condizionato se piccole variazioni dei dati non determinano grandi variazioni dei risultati.

Esempio 1: Intersezione tra due rette r e s

In questo caso la variazione del risultato è comparabile con la variazione del dato, quindi per piccole variazioni del dato si hanno piccole variazioni dei risultati: il problema è ben condizionato.

Esempio 2: Intersezione tra due rette r e t

In questo caso la variazione del risultato è molto maggiore della variazione del dato, quindi per piccole variazioni del dato NON si hanno piccole variazioni dei risultati: il problema è MAL condizionato.

Esempio 3: Analizzare l'integrale

I(x) = ∫01 enx x dx = [ex - n·e]01 - n·In-1 = 1 - n·In-1

Tramite l'integrazione per parti si ha la formula ricorsiva. Introduco uno schema: una volta calcolato I_n=0, ottengo I_n+1 in maniera ricorsiva.

Parto quindi dal dato iniziale: I₀ = 1 - e^-1.

Problema: I₀ ha infinite cifre decimali a causa del termine e^-1. Approssimo quindi il termine e^-1, e pongo:

I0 ≈ I0~ = 1 - e-1

Poiché I₀ ≠ I_0~, vi è un errore ε₀, da cui:

I = I0~ + ε0

Si ha quindi: ε_n = (-n)·ε_n-1 = ... = (-1)ⁿ·n!·ε₀

Problema: L'errore cresce in modo esponenziale, quindi, per n sufficientemente grande, l'accuratezza del calcolo è nulla, quindi lo schema è instabile!

Adotto, quindi, un altro metodo: poiché lim_n→+∞ I_n = 0, scelgo un termine sufficientemente grande e pongo I_N ≈ 0. Ottengo quindi: se vale la relazione |I_n - I_N| < τ, l'ipotesi è buona; in caso contrario, sarà necessario trovare un termine N maggiore sufficientemente grande e ripetere lo schema.

Per N sufficientemente grande, varrà: ε_N = ε₀·N!·(-1)^N

Si ha, quindi, un errore che si riduce progressivamente: lo schema è stabile.

Problema del punto fisso

La risoluzione di equazioni non lineari è correlata al classico problema: “Data la funzione f, trovare ψ ∈ ℝ | f(ψ) = 0”.

Esempio: f(x) = x² + 3x - 4

Tale problema può dare una o più soluzioni in numero finito (es. 2⁴ - x), infinite soluzioni (es. sen(x)) o nessuna soluzione (es. x³ + 3x).

Per valutare se la mia soluzione tende al valore vero in modo sufficientemente accurato, valuto f(x): se ottengo |f(x)| < τ, allora si ha ψ ≈ x.

ATTENZIONE: Il metodo di verifica non è sempre valido!

Esempio: f(x) = x² - τ·10^-8

Ponendo x = 0, si ha f(0) = τ·10^-8 ≪ τ, ma il problema originario non ha soluzioni!

Esercizio: Trovare ψ ∈ ℝ dove f(ψ) = 0, f(x) = x² - 3x + 2

La risoluzione analitica del problema fornisce i seguenti risultati:

Δ = 12 - 2 ⇒ x = 3 ± √Δ / 2

Problema: Un elaboratore non lavora analiticamente, ma solo numericamente! Costruisco, quindi, una successione tale che {x_k} → ψ con k → +∞. La successione è data dall'equazione:

xk+1 = g(xk)

Tale equazione è un problema di punto fisso, e rappresenta l'intersezione tra la curva y = g(x) e la retta y = x: la funzione g si ottiene esplicitando la variabile x nell'equazione f(x) = 0.

Esempio: x + log(x) = 0 → x = -log(x) o x = e^-x

Definizione: uno schema si definisce fortemente consistente se, sostituendo il valore vero della soluzione, si ottiene un'identità. Lo schema del punto fisso, o di Picard, è fortemente consistente per definizione: (ψ = g(ψ)) → (x → ψ).

Nel nostro caso, f(x) = x² - 3x + 2, potrò quindi esplicitare x in due modi diversi, ottenendo dunque due funzioni g:

g1(x) = 3 - √(x2 + 2)
g2(x) = 3 - 2/x

Definiamo, inoltre, una terza funzione, di cui si discuterà successivamente:

g3(x) = (x2 - 2)/3

Compongo, dunque, lo schema per la prima funzione:

x0 = 3, x1 = g1(x0) ≈ 2.645
x2 = g1(x1) ≈ 2.436
...

Lo schema mi dà, effettivamente, una delle radici di f(x). Componendo lo schema con la seconda funzione di punto fisso, si ha x_k → +∞: lo schema è dunque inefficace!

Componendo lo schema con la terza funzione, si ha:

x0 = 3, x1 = g3(x0) ≈ 2.333
x2 = g3(x1) ≈ 2.067
x3 = g3(x2) ≈ 2.004
...

Già alla terza iterazione, lo schema fornisce una buona approssimazione della soluzione x = 2, migliore del valore fornito alla terza iterazione dalla prima funzione: la terza funzione risulta, quindi, MOLTO più veloce della prima funzione!

Definendo l'errore ε_k e studiando il rapporto |ε_k+1/ε_k|, si osserva che esso, per g₁(x), assume un valore M pressoché costante approssimabile a 0.73.

Si ha quindi: ε_k+1 ≈ 0.73·ε_k → 0 con k → +∞.

Studiando il rapporto ε_k+1/ε_k per g₂(x), si osserva che esso assume un valore M costante e maggiore di 1: per questo motivo, poiché vale ancora la definizione ε_k+1 ≈ M·ε_k, si deduce che lo schema è instabile, poiché ε_k+1 → +∞ con k → +∞.

Per g₃(x), si osserva che il rapporto ε_k+1/ε_k assume un valore prossimo allo 0.

Si può dimostrare che i valori di ε_k corrispondono ai valori delle derivate prime delle funzioni g₁, g₂, e g₃ per x = ψ.

Si osserva, quindi, che lo schema più stabile e consistente, g₃(x), è quello per cui si ha il valore minore della derivata prima per x = ψ.

Definizione: convergenza dello schema

Cosa significa che uno schema funziona “fino a convergenza”? A che valore dell'iterazione è possibile fermarsi?

Si ha convergenza quando: lim_k→+∞ |x_k - ψ| = 0 = g(x_k) - x_k.

Osservo: Se si ha una variazione, o scarto, minore della tolleranza, calcoli ulteriori non sono utili, poiché i risultati non compaiono tra le cifre significative di x_k!

NB: Sarà necessario verificare SE e QUANDO vale la relazione: |ε_k+1| ≤ |x_k+1 - ψ|·C ≤ τ.

Implementazione dello schema di Picard all'elaboratore

Per inserire questo algoritmo in un codice FORTRAN, sarà sufficiente inserire le seguenti istruzioni:

program Picard
implicit none
integer iter, itmax
real*8 x0, xk, xkp1, scarto, toll

Dove il comando implicit none indica che ogni variabile DEVE essere dichiarata: in caso non si imponga l'implicit none, il formato della variabile è legato alla iniziale del suo nome (es. i = intero).

I comandi integer e real*8 indicano che le variabili elencate dopo il rispettivo comando sono rispettivamente intere e reali a 16 cifre decimali – a “doppia precisione”. Esiste anche il comando real*4, che indica le variabili intere a 8 cifre decimali – a “precisione singola”.

ATTENZIONE: Imponendo l'implicit none è OBBLIGATORIO descrivere il formato di OGNI variabile utilizzata!

Si impone il valore iniziale di x, il numero limite di iterazioni, la tolleranza: è possibile imporre dei valori nel codice o eseguendo una lettura da terminale o una lettura da file.

Esempio: nel codice

...
x0 = 1.0d0
itmax = 100
toll = 1e-10
...

da terminale

write(*,*) 'inserire il valore di x0, itmax, toll'
read(5,*) x0, itmax, toll

da file

...