Dodicesima lezione

Analizziamo ora il problema di "risolvere", in campo reale, un'equazione ad una incognita, cioè di trovare il passaggio per lo zero di una funzione qualunque.
Come spesso accade in matematica, affrontare il problema per qualunque funzione può essere complicato. Qualunque metodo scegliamo, infatti, può essere messo in crisi scegliendo opportune funzioni cattive a piacere. Nei casi pratici, possiamo spesso immaginare che la funzione in esame sia perlomeno continua. Anche questa richiesta non ci libera tuttavia dalla necessità di trovare un opportuno intervallo iniziale di ricerca. In generale, quando si tratta di algoritmi numerici per la ricerca degli zeri (o dei massimi e dei minimi), la scelta del punto di partenza può sempre influenzare in modo critico il risultato. Questo è vero nel caso ad una dimensione (che stiamo trattando), ma è ancora più vero nel caso a più dimensioni.
Il teorema del valor medio garantisce che una funzione continua che assume valori di segno opposto agli estremi di un intervallo ammetta almeno uno zero nell'intervallo. Se però disgraziatamente la funzione ha più di uno zero nell'intervallo, il punto di partenza e la scelta dell'intervallo possono ancora una volta variare drasticamente il risultato.

In ogni caso, è necessario avere un qualche criterio per definire l'intervallo di ricerca iniziale. Due sono i criteri utilizzati normalmente:

1. Espansione di un intervallo "piccolo" fino a trovare che la funzione cambia di segno.
2. Partizionamento di un intervallo "grande" in intervalli più piccoli, fino a trovare uno (o più) intervalli agli estremi dei quali la funzione cambia di segno.

Nel successivo esempio di codice sono implementati entrambi questi metodi.

Una volta stabilito che un intervallo contiene uno zero, occorre definire la precisione richiesta nella determinazione dello zero stesso. E' molto difficile stabilire questa tolleranza in modo automatico, per la solita ragione che essa dipende dalla precisione disponibile per la rappresentazione dei numeri nella regione dove si trova lo zero. Se cerchiamo uno zero dalle parti di x = 10²³, non ha molto senso chiedere una tolleranza piccola.
Per questa ragione gli algoritmi di ricerca degli zeri si aspettano di ricevere dall'utente l'errore di misura tollerato.

Il criterio di ricerca degli zeri più semplice consiste nel continuare a dividere in due l'intervallo di ricerca iniziato, e continuare a scegliere l'intervallo agli estremi del quale la funzione cambia segno.
Non è un criterio molto efficiente, ma sicuramente converge a qualcosa, anche se la funzione è discontinua o ha punti di singolarità. Se nell'intervallo ci sono più zeri, convergerà ad uno di questi. Se c'è una discontinuità a gradino, convergerà alla discontinuità, e se c'è un punto singolare convergerà al punto singolare.

Vediamo i due criteri di ricerca dell'intervallo presentati sopra ed il metodo di bisezione in azione in questo .

Esistono vari metodi per accelerare il processo di convergenza nella ricerca degli zeri. Anche in questo caso, prendiamo subito in considerazione il metodo "migliore" (se non altro quello più utilizzato). Si tratta del metodo di Newton, detto anche di Newton-Raphson, e si basa sull'assunzione che la derivata prima della funzione in esame possa essere calcolata. Questa la ricetta:

1. Calcolare la derivata in un punto dell'intervallo di ricerca, il che permette di ottenere la retta tangente alla funzione in quel punto.
2. Calcolare quindi il punto in cui la tangente passa per lo zero: se tale punto esce dall'intervallo di ricerca bisogna gettare la spugna.
3. Calcolare di nuovo la derivata in quel punto, e ripetere dal punto (2) fino a convergenza avvenuta.

Un buon approccio di programmazione ci permette però di sommare la rapidità del metodo di Newton-Raphson con la stabilità del metodo di bisezione: basta ricorrere ad un passo di bisezione quando il metodo di Newton tenderebbe ad uscire dall'intervallo oppure converge troppo lentamente. Ecco un che mette in pratica questo metodo. Sottolineaiamo ancora che per applicare il metodo di Newton-Raphson è necessario conoscere la derivata della funzione.

A questo punto, però, una domanda sorge ovvia: non sarebbe possibile calcolare numericamente la derivata della funzione nel punto desiderato?
La risposta a questa domanda richiede una breve riflessione sui particolari del calcolo numerico delle derivate.

Il calcolo della derivata di una funzione sembrerebbe, alla luce di quanto abbiamo fatto finora, la cosa più semplice del mondo. Basta applicare la definizione:

E utilizzare al posto di un incremento "tendente a zero" un incremento "piccolo". In questo modo l'errore nel calcolo della derivata "dovrebbe" essere dato soltanto dai termini successivi al primo dello sviluppo di Taylor:

Questo ovviamente se il calcolatore avesse, come più volte auspicato, precisione infinita. Le cose non stanno così, ed i primi problemi sorgono già nella scelta di un incremento piccolo.
Sappiamo infatti che un incremento, per poter essere significativamente rappresentato dalla macchina, deve avere una qualche relazione con la scala del numero a cui viene aggiunto. Facciamo una piccola prova in singola precisione:

  float x,h;
  float temp;

  x = 10;
  h = 0.0000005;

  temp = x + h;

  printf("x = %.8f h = %.8f x+h = %.8f x+h-x = %.8f\n",
         x,h,temp,temp-x);

Questo il risultato:

x = 10.00000000 h = 0.00000050 x+h = 10.00000095 x+h-x = 0.00000095

L'incremento che abbiamo effettivamente aggiunto è quasi il doppio di quello previsto; di conseguenza dobbiamo in primo luogo accertarci che l'incremento piccolo, possa essere esattamente rappresentato dalla macchina. Si può usare proprio un trucco di questo tipo:

  float x,h;
  volatile float temp;

  x=...
  h=...

  temp = x+h;
  h = x-temp;

dove la parola chiave volatile istruisce il compilatore a non pensare di avere completo controllo sulla variabile (un conto così potrebbe infatti sembrare stupido anche al compilatore stesso, che potrebbe tentare di eliminarlo in un passo di ottimizzazione).

Dopo aver determinato esattamente il valore dell'incremento, restano due contributi all'indeterminazione nel calcolo della derivata:

1. L'errore di arrotondamento nel calcolo del rapporto incrementale, dovuto alla precisione finita della macchina. Tolta l'indeterminazione sull'incremento, esso vale all'incirca kf(x)/h, dove k è l'errore relativo di arrotondamento per il valore di di f(x).
2. L'approssimazione nel calcolo della derivata che, come abbiamo visto sopra, è dell'ordine di hf"(x).

Se ora scegliamo il valore di h che minimizza queste indeterminazioni (lasciamo perdere i conti), scopriamo che in questo, che è il migliore dei casi, l'errore relativo che grava sul calcolo della derivata è pari alla radice quadrata dell'indeterminazione relativa intrinseca k, dove k è un numero piccolo... Questa non è certo una buona notizia!

Possiamo allora ricorrere ad uno sporco trucco:

In questo modo l'errore nel calcolo della derivata utilizzando una differenza finita dipende dalla derivata terza e da h². Il miglioramento ottenuto, tuttavia (anche qui lasciamo perdere i conti), è solo tale da portare l'errore relativo sulla derivata a k^2/3.

Se ci abbandonassero su un'isola deserta, forse quest'ultimo sistema potrebbe anche andare bene, ma siccome abitiamo una ricca comunità umana, ringraziamo C.J.F. Ridders, che ci ha offerto (1982) un algoritmo che applica un metodo di estrapolazione polinomiale (di cui non ci occupiamo) ad una successione di valori, che sono i rapporti incrementali finiti, dove l'incremento viene di volta in volta diminuito di un valore fissato.

Lasciamo a questo punto la parola a Numerical recipes in C, che ci offre questa implementazione.

Esercizio: provare a disegnare qualche derivata di funzione, calcolata con l'algoritmo di Ridders.