Libero

Adavescienzeapproxcalc_apx.html

Segnala abuso

Metodi numerici per il calcolo approssimato

Il calcolo approssimato è un approccio utilissimo a certi problemi di difficile soluzione, come il calcolo delle aree descritte da una curva complicata oppure trovare lo zero di una equazione di grado n. In questi casi una risoluzione algebrica può essere molto difficile, ma se il risultato che si cerca non è teorico (pi^2/3) ma pratico (2.145...) ci si può affidare al calcolatore. Il valore che si ottiene è approssimato con precisione a piacere (compatibilmente con i limiti del calcolatore utilizzato).
In questa pagina si esaminano vari metodi di approssimazione numerica al calcolatore con relativa traduzione in due linguaggi di programmazione: il linguaggio "universale e sempreverde" C, e il linguaggio Pascal, ancora diffuso nella scuola media superiore.

Sommario

Le approssimazioni al calcolatore
Zeri di una funzione
Integrali
1. Metodo dei trapezi (o di Bézout)
2. Metodo delle parabole (o di Cavalieri-Simpson)
Derivate
1. Metodo dei tre punti

Fonti

Appunti universitari e di liceo scientifico
http://www.liceofoscarini.it/didattic/numerici/index.html
Liceo Classico Statale Foscarini (VE) - Materiale didattico, Metodi numerici
http://www.fisica.uniud.it/%7Eercolessi/MQ/mq/mq.html
Università di Udine - Appunti per il corso di Meccanica quantistica

Strumenti

Grafici tracciati con Mathplot per Linux; equazioni scritte con MS Equation Editor (parte di Office per Windows); diagrammi di flusso realizzati con Kivio (parte di Koffice per Linux).

Le approssimazioni al calcolatore

Realizzare delle approssimazioni al calcolatore significa prima di tutto riuscire a scomporre il problema in parti più semplici, ed identificare un algoritmo (una "strategia", o, meglio, una sequenza di operazioni) che lo risolva; poi si tratta di tradurre l'algoritmo in un programma che faccia i calcoli. Esistono alcuni algoritmi standard, e sono quelli presentati di seguito; altri possono essere inventati con un po' di lavoro teorico.

Per sapere se un algoritmo è "buono" oppure no ci sono alcuni parametri: velocità, leggerezza, flessibilità. La velocità è un indice del numero di passaggi necessari per ottenere un certo grado di approssimazione; la leggerezza è legata al numero di operazioni che devono essere eseguite ad ogni passaggio; la flessibilità coinvolge le limitazioni, cioè le eventuali condizioni che devono verificarsi perchè l'algoritmo possa essere applicato. Un algoritmo che impiega troppo tempo, che genera troppo carico sul processore o può essere applicato solamente in un numero ristretto di casi è da scartare.

Una volta individuato l'algoritmo, è necessario implementarlo, cioè realizzare un programma che lo esegua. Gli algoritmi qui presentati sono relativamente semplici, e coinvolgono conoscenze di base sulla programmazione e sulle funzioni. Non esiste un linguaggio particolarmente adatto alla loro implementazione; ne riporto la trascrizione in linguaggio C, il linguaggio standard del mondo Unix / Linux, e in Turbo Pascal, utilizzato nei corsi sperimentali dei Licei Scientifici. Le descrizioni e i diagrammi di flusso permettono comunque una facile trasposizione in altri linguaggi.

programma in C

programma in Pascal

Nota tecnica: i programmi presentati aderiscono strettamente all'algoritmo di risoluzione; ad esempio a volte si utilizzano due diverse variabili quando sarebbe indifferente per il calcolatore utilizzarne una sola. In questo senso si è privilegiata la chiarezza sintattica dei listati piuttosto che l'ottimizzazione del programma risultante.

Tutti gli algoritmi per il calcolo approssimato, proprio perchè si tratta di approssimazioni successive, consistono quasi unicamente di un ciclo, caratterizzato da delle condizioni di uscita, e spesso da una coppia di valori che contengono il calcolo precedente e quello attuale. I parametri iniziali dei programmi, l'uso delle variabili e le condizioni di uscita sono descritti volta per volta.

Problemi risolvibili con metodi numerici:

Calcolo degli zeri di una funzione

Calcolo delle aree (Integrali)

Zeri di una funzione

Il primo problema che affrontiamo è quello di calcolare in modo approssimato gli zeri di una funzione, cioè i valori di x per cui si annulla. Abbiamo bisogno di alcune ipotesi importanti, cioè che la funzione sia continua e ammetta uno zero nell'intervallo [a,b]. Questo perchè gli algoritmi qui proposti si basano sul progressivo restringimento dell'intervallo che contiene lo zero; ecco quindi che dobbiamo già avere un'idea approssimativa di dove la funzione si annulla. Se la funzione non è continua nell'intervallo i metodi proposti cessano di avere fondamento teorico, anche se è possibile in alcuni casi ottenere ugualmente qualche risultato.
Inoltre, condizione indispensabile per poter applicare questi metodi è che il segno della f(x) a destra dello zero sia opposto al segno della f(x) a sinistra dello zero, cioè la funzione deve essere secante all'asse delle x e non tangente. La necessità di questa ipotesi è palese considerando gli esempi grafici in questa pagina.

Per ottenere l'intervallo di partenza [a,b] è sufficiente fare qualche prova, cercando due valori di x per cui le due y=f(x) abbiano segno opposto, oppure si può tracciare un grafico sommario a mano o con qualche programma apposito (come l'ottimo Mathplot, per Linux) e prendere due valori di x vicini a x₀.

Metodi per calcolare gli zeri di una funzione:

Metodo di bisezione

Metodo delle tangenti

Metodo delle secanti

Un confronto tra i metodi

Il metodo di bisezione

Il metodo di bisezione è il più semplice: si considera solamente l'intervallo [a,c] con c punto medio fra a e b, e si calcola il valore della funzione nei due estremi. Se i due valori hanno segno opposto (il prodotto è negativo) allora la funzione si annulla in [a,c]; si pone b = c, e si ripete la procedura. Se invece il prodotto è positivo, la funzione si annulla in [c,b]; si pone a = c e si ripete il ciclo. Il programma si ferma quando l'intervallo (e quindi l'errore su x₀) è diventato più piccolo di un "epsilon", valore a scelta; si può scegliere |a-b|>e (errore assoluto), oppure |a-b|>e·a (errore relativo).

Nei seguenti programmi a è l'estremo sinistro dell'intervallo, b è l'estremo destro, c è il punto medio degli intervalli successivi, e è l'incertezza desiderata; alla fine del ciclo si può prendere indifferentemente a o b come approssimazione dello zero.
La funzione presa in considerazione è f(x)=1/x-ln(x) nell'intervallo [1,2].

#include <stdio.h>
#include <math.h>

float f(float x) {
float y;
y = 1/x - log(x);
return y;
}
main (void) {
float a, b, c, e;
a = 1;
b = 2;
e = 0.0001;
while (b-a > e) {
   c = (b+a)/2;
   if ( f(a)*f(c) < 0) { b = c }
   else { a = c }
   }
printf("Lo zero si ha in x = %f", a);
}

program bisezione;

var a, b, c, e :real;

function f(x :real) :real;
begin
f := 1/x - ln(x);
end;

begin
a := 1;
b := 2;
e := 0.0001;
while (b-a > e) do begin
   c := (b+a)/2;
   if ( f(a)*f(c) < 0) then b := c
   else a := c;
   end;
write('Lo zero si ha in x = ', a);
end.

NOTA: in alcuni casi sono possibili alcuni problemi con la condizione di arresto:

considerando |a-b|>e (errore assoluto), con un x₀ molto grande gli errori di arrotondamento del calcolatore potrebbero essere maggiori di e e l'algoritmo si fermerebbe troppo presto;
considerando |a-b|>e·a (errore relativo) ci possono essere problemi per x prossimo allo 0;
se la pendenza di f(x) in prossimità dello zero è molto piccola, ci può essere un intero intervallo in cui i valori di f(x) si equivalgono nell'approssimazione del calcolatore.

Il metodo delle tangenti

Il metodo delle tangenti è chiaramente illustrato dalla figura a lato: partendo da un punto prossimo allo zero, l'intersezione dell'asse delle x con la tangente alla funzione in quel punto determina la successiva approssimazione dello zero cercato. In formule, questo si traduce con

dove il primo valore di x_n è uno qualsiasi dei due estremi dell'intervallo che contiene lo zero. La formula converge verso lo zero, a patto che la funzione sia monotona fra l'estremo scelto e lo zero: la tangente non deve mai cambiare di segno durante il processo di avvicinamento, altrimenti l'algoritmo "torna indietro". Il difetto è facilmente verificabile considerando la funzione x²-8x+16+tan(x) e un intervallo di ampiezza 4.

Dimostrazione

La funzione presa in esame negli esempi è 1/x-ln(x); a e b sono gli estremi dell'intervallo considerato, e l'algoritmo comincia con x₀ = a. L'esecuzione termina quando due approssimazioni consecutive differiscono di meno di un epsilon (e) dato all'inizio del programma (si utilizza il valore assoluto per comprendere anche il caso in cui le approssimazioni siano sia a destra e sinistra dello zero; può capitare con particolari funzioni senza che la validità dell'algoritmo sia compromessa).
x2 e x1 rappresentano l'approssimazione precedente e quella corrente. Inizialmente a queste variabili vengono dati i valori degli estremi dell'intervallo, in modo da poter effettuare da subito (e con esito positivo) il confronto che rappresenta la condizione di continuazione del ciclo. Al primo passaggio del ciclo l'estremo sinistro viene assegnato ad x2, mentre x1 viene calcolato, e così via.

#include <stdio.h>
#include <math.h>

float f(float x) {
float y;
y = 1/x - log(x);
return y;
}
float f1(float x) {
float y;
y = -1/(x*x) - 1/x;
return y;
}
main (void) {
float a, b, e, x1, x2;
a = 1;
b = 2;
e = 0.0001;
x1 = b;
x2 = a;
while (fabs(x1-x2)>e) {
   x2 = x1;
   x1 = x2 - f(x2)/f1(x2);
   }
printf("Lo zero si ha in x = %f", x1);
}

program tangenti;

var a, b, e, x1, x2 :real;

function f(x :real) :real;
begin
f := 1/x - ln(x);
end;

function f1(x :real) :real;
begin
f1 := -1/(x*x) - 1/x;
end;

begin
a := 1;
b := 2;
e := 0.0001;
x1 := b;
x2 := a;
while (abs(x1-x2)>e) do begin
   x2 := x1;
   x1 := x2 - f(x2)/f1(x2);
   end;
write('Lo zero si ha in x = ', x1);
end.

Il metodo delle secanti

Questo metodo lavora "dall'altra parte" rispetto al metodo della tangente: infatti utilizza l'intersezione di una secante con l'asse delle x per determinare l'approssimazione successiva. Rispetto al metodo della tangente non richiede la disponibilità della derivata, ma necessita di un punto fisso c da utilizzare come secondo punto per definire le secanti. La formula che ne risulta è

la cui dimostrazione è riassunta a lato.
L'algoritmo converge nella stragrande maggioranza dei casi, anche se è possibile che in qualche "caso patologico" possa "perdere" lo zero; ad ogni modo la sua efficacia è meno legata alle particolari caratteristiche della funzione rispetto al metodo della tangente, anche se è un po' più lento come numero di iterazioni necessarie.
Per il resto, la struttura ciclica e le condizioni di uscita dei due metodi sono identiche.

Dimostrazione

Il programma è costruito esattamente come quello del metodo delle tangenti: stessa funzione di esempio, stesse variabili, stesso ciclo, con l'unica variazione della diversa formula utilizzata per il calcolo del termine successivo.

#include <stdio.h>
#include <math.h>

float f(float x) {
float y;
y = 1/x - log(x);
return y;
}
main (void) {
float a, b, e, x1, x2;
a = 1;
b = 2;
e = 0.0001;
x1 = b;
x2 = a + e/2;   /* ! vedi nota ! */
while (fabs(x1-x2)>e) {
   x2 = x1;
   x1 = x2 - (f(x2)*(x2-a))/(f(x2)-f(a));
   }
printf("Lo zero si ha in x = %f", x1);
}

program secanti;

var a, b, e, x1, x2 :real;

function f(x :real) :real;
begin
f := 1/x - ln(x);
end;

begin
a := 1;
b := 2;
e := 0.0001;
x1 := b;
x2 := a + e/2; (* ! vedi nota ! *)
while (abs(x1-x2)>e) do begin
   x2 := x1;
   x1 := x2 - (f(x2)*(x2-a))/(f(x2)-f(a));
   end;
write('Lo zero si ha in x = ', x1);
end.

NOTA: la prima assegnazione di x2 non è a, bensì a+e/2. Questo perchè la formula utilizza al numeratore (x-c), con c estremo sinistro, nel nostro caso ancora a. Sostituendo, si avrebbe (a-a), che dà zero e fa terminare precocemente l'algoritmo. Quindi il primo valore di x2 non può essere a, ed a+e/2 è un buon valore perchè è comunque molto vicino ad a, e vi si discosta di un valore minore della precisione desiderata per il risultato. In questo modo non dovrebbero crearsi problemi neppure in caso di intervalli molto stretti.

Un confronto tra i metodi

Come si vede dalla tabella sottostante, il metodo migliore sembra quello della secante; ciò non toglie che con particolari macchine o particolari funzioni convenga utilizzare gli altri due metodi. A volte per risolvere con buona precisione un determinato problema è meglio addirittura combinare due metodi, applicandone uno e passando ad un altro da un certo momento in poi (il primo dei due sarà quasi certamente l'algoritmo di bisezione).

	Bisezione	Tangenti	Secanti
Velocità	lento	molto veloce	veloce
Leggerezza	molto buona	media (deve calcolare anche la derivata prima)	buona
Flessibilità	si ferma troppo presto per funzioni con pendenza quasi nulla allo zero	medio-bassa: può uscire dall'intervallo dato, ma pone alcune condizioni sulla derivata	buona: può uscire dall'intervallo, se un estremo è sbagliato; non richiede la derivata

Integrali

Il calcolo degli integrali (l'area compresa fra la curva di una funzione e l'asse delle x) è un altro dei problemi classici che non sempre trova facile soluzione per via algebrica.

Entrambi i metodi considerati dividono l'intervallo di integrazione in intervalli più piccoli, e all'interno degli intervalli la funzione data viene approssimata in qualcosa di geometricamente più semplice; i diversi modi di approssimare la funzione caratterizzano i due metodi: il metodo dei trapezi utilizza delle rette, mentre il metodo delle parabole utilizza, appunto, delle parabole. Viene quindi calcolata la somma delle aree di tutti i sottintervalli.

Va da sè che quanti più sono gli intervalli considerati, tanto minore sarà la loro ampiezza e tanto più l'approssimazione si avvicinerà alla funzione data, ottenendo un risultato più preciso.

Nota: mentre i metodi per il calcolo degli zeri di una funzione sono necessariamente iterativi, cioè ogni passaggio si basa sul precedente e ne migliora la precisione, per approssimare un integrale ci si può accontentare un calcolo solo con un determinato numero di intervalli, scelto alto a piacere. Un programma di ricerca dell'area potrebbe prevedere un ciclo che rimpicciolisca la dimensione degli intervalli ad ogni passaggio e si fermi quando due calcoli successivi dell'area differiscono di meno di una costante data (questa sarà quindi l'incertezza sull'area); lo lasciamo come esercizio.

Metodi per calcolare gli integrali:

Metodo dei trapezi

Metodo delle parabole

Il metodo dei trapezi

L'integrale di una funzione secondo Riemann è definito come il limite delle due aree delimitate da due funzioni a scalino, una che approssima la funzione data per eccesso e una che la approssima per difetto, per n numero di intervalli tendente ad infinito (una funzione si dice a scalino se il suo dominio può essere diviso in intervalli nei quali la funzione è costante; il grafico risultante è proprio "a scalini").
Le due funzioni sono rappresentate in rosso e in blu nella figura in alto; ogni intervallo ha uguale ampiezza h = (b-a/n).

con h ampiezza degli n intervalli calcolata come (b-a)/n.

La formula che ne risulta sommando le aree di tutti gli intervalli e raccogliendo è

come si può verificare con pochi passaggi partendo dalla formula dell'area di singolo intervallo.

Il programma seguente calcola l'integrale di e^x fra 0 e 1, che noi sappiamo essere uguale a e-1 = 1.7... dividendo l'intervallo di integrazione in intervalli più piccoli, e applicando la semplificazione dei trapezi in ognuno di questi. Nei listati, a e b sono gli estremi di integrazione, n il numero di intervalli (più grande è n, più piccolo è h, più preciso è il valore calcolato), h l'ampiezza di un intervallo, x1 il valore della x di volta in volta considerato e area la variabile contenente la somma delle aree fino ad ora calcolate.

#include <stdio.h>
#include <math.h>

float f(float x) {
float y;
y = exp(x);
return y;
}
main (void) {
int n;
float a, b, h, x1, area;
a = 0;
b = 1;
n = 100;
h = (b-a)/n;
area = 0;
x1 = a;
while (x1<b) {
   area += f(x1) + f(x1+h);
   x1 += h;
   }
area += f(b);
area = area*h/2;
printf("L'integrale vale %f", area);
}

program trapezi;

var a, b, h, x1, area :real;
     n :integer;

function f(x :real) :real;
begin
f := exp(x);
end;

begin
a := 0;
b := 1;
n := 100;
h := (b-a)/n;
area := 0;
x1 := a;
while (x1<b) do begin
   area := area + f(x1) + f(x1+h);
   x1 := x1 + h;
   end;
area := area + f(b);
area := area*h/2;
write('L''integrale vale ', area);
end.

Il metodo delle parabole

Il metodo delle parabole è simile al metodo dei trapezi ma utilizza una serie di parabole come approssimazione della funzione data; ogni parabola interpola tre valori di f(x), così da occupare due intervalli; è quindi necessario che il numero di intervalli sia pari.
La formula finale, la cui dimostrazione è relativamente lunga, è:

Il programma seguente calcola l'integrale di e^x fra 0 e 1, che noi sappiamo essere uguale a e-1 = 1.7... dividendo l'intervallo di integrazione in intervalli più piccoli, e applicando la semplificazione delle parabole in ognuno di questi. Nei listati, a e b sono gli estremi di integrazione, n il numero di intervalli (più grande è n, più piccolo è h, più preciso è il valore calcolato - n deve essere pari!), h l'ampiezza di un intervallo, x1 il valore della x di volta in volta considerato e area la variabile contenente la somma delle aree fino ad ora calcolate.

#include <stdio.h>
#include <math.h>

float f(float x) {
float y;
y = exp(x);
return y;
}
main (void) {
int n;
float a, b, dx, x1, area;
a = 0;
b = 1;
n = 100; /* n dev'essere pari */
dx = (b-a)/n;
area = 0;
x1 = a;
while (x1+h<b) {
   area += f(x1) + 4*f(x1+dx) + f(x1+2*dx);
   x1 += 2*dx;
   }
area = area*dx/3;
printf("L'integrale vale %f", area);
}

program parabole;

var a, b, h, x1, area :real;
     n :integer;

function f(x :real) :real;
begin
f := exp(x);
end;

begin
a := 0;
b := 1;
n := 100; (* n dev'essere pari *)
h := (b-a)/n;
area := 0;
x1 := a;
while (x1+h<b) do begin
   area := area + f(x1) + 4*f(x1+h) + f(x1+2*h);
   x1 := x1 + 2*h;
   end;
area := area*h/3;
write('L''integrale vale ', area);
end.

Derivate

Per approssimare il valore della derivata di una funzione in un punto si parte dalla definizione di derivata, ossia come limite del rapporto incrementale. La formula utilizzata è proprio il rapporto incrementale con dx piccolo a piacere: più è piccolo, più preciso sarà il valore trovato.

Ci sono vari modi di utilizzare la formula (rapporto incrementale destro/sinistro, misto...), e se ne possono utilizzare anche due o tre per lo stesso problema e poi fare una media.

Nota: mentre i metodi per il calcolo degli zeri di una funzione sono necessariamente iterativi, cioè ogni passaggio si basa sul precedente e ne migliora la precisione, per approssimare la derivata ci si può accontentare un calcolo solo che contenga subito in dx l'approssimazione voluta. Ovviamente in questo modo si determina l'incertezza sulla x e non sulla derivata, ma con valori di dx molto piccoli la differenza non è rilevante. Un programma con tutti i crismi dovrebbe invece prevedere un ciclo che rimpicciolisca dx ad ogni passaggio e si fermi quando due calcoli successivi della derivata differiscono di meno di una costante data (questa sarà quindi la vera incertezza della derivata); lo lasciamo come esercizio.

Metodi per calcolare le derivate:

Metodo dei tre punti

Il metodo dei tre punti

Il metodo dei tre punti è il meno banale degli altri ma solamente perchè li comprende un po' tutti facendo poi una media (der) dei vari risultati. Un risultato (der0) è il valore della pendenza del segmento che unisce i punti (x-dx,f(x-dx)) e (x+dx,f(x+dx)), mentre gli altri due risultati sono il rapporto incrementale destro e sinistro (der1 e der2).

#include <stdio.h>
#include <math.h>

float f(float x) {
float y;
y = sin(x);
return y;
}
int main(void) {
float x, dx, der, der0, der1, der2;
x = 3.14;
dx = 0.001;
der0 = (f(x+dx)-f(x-dx))/(2*dx);
der1 = (f(x+dx)-f(x))/dx;
der2 = (f(x)-f(x-dx))/dx;
der = (der0+der1+der2)/3;
printf("La derivata in %f vale %f",x,der);
}

program derivata;

var x, dx, der0, der, der1, der2 :real;

function f(x :real) :real;
begin
f := sin(x);
end;

begin
x := 3.14;
dx := 0.001;
der0 := (f(x+dx)-f(x-dx))/(2*dx);
der1 := (f(x+dx)-f(x))/dx;
der2 := (f(x)-f(x-dx))/dx;
der := (der0+der1+der2)/3;
write('La derivata in ',x,' vale ', der);
end.

2003/2004, Cester Davide