1. Soluzione numerica di equazioni non-lineari




Data una funzione y=f(x), in generale nonlineare nella variabile x, trovare i valori di x per cui f(x)=0 significa trovare gli zeri di f(x). Questa operazione, assai diffusa nel calcolo numerico, richiede in generale l’utilizzo di un metodo di approssimazione della soluzione di tipo iterativo. Infatti, anche nel caso in cui la funzione sia un polinomio algebrico, è ben noto che una formula risolutiva che richieda un numero finito di operazioni, come ad esempio la formula delle radici di un polinomio del secondo grado, non esiste già a partire dal grado 5.

Trattandosi di un metodo iterativo, ci sono due aspetti specifici di cui è opportuno tenere conto nell’analisi e nell’applicazione del metodo: la garanzia di convergenza alla soluzione e la velocità con cui questo accade. Questa velocità è tipicamente rappresentata dall'ordine di convergenza ed è definita nel modo seguente:   si dice che la successione {x(i)} converge ad alfa con ordine p, se esiste C>0 e finito tale che:

(   abs( x(i+1) - alfa )  )   /   (  abs( x(i) - alfa )^p  )  ≤  C

per ogni  i>K  , dove "K"  è un intero opportuno.    Per valutare sperimentalmente l’ordine di convergenza di un metodo iterativo è sufficiente, una volta finita l’esecuzione del metodo, calcolare il rapporto definito qui sopra, ponendo p=1, lungo la successione {x(i)} calcolata: i metodi a convergenza lineare mantengono valori costanti, mentre se i valori decrescono significa che la convergenza è più che lineare.

Vediamo dunque di costruire e di sperimentare qualche metodo iterativo per il calcolo degli zeri di f(x):




1.    Metodo di Bisezione :

Se la funzione f(x) è continua in [a,b] e se f(a)*f(b)<0 , f(x) contiene sicuramente almeno uno zero in [a,b]. Per semplicità supponiamo che ne contenga uno solo.
Il metodo di bisezione suddivide ad ogni passo l’intervallo, scegliendo poi di proseguire nella metà in cui i valori di f(x) agli estremi “nuovo_a” e “nuovo_b”  verifichino:  f(nuova_a)*f(nuova_b)<0  .
Si ottiene un metodo di convergenza certa anche se lenta (il metodo ha ordine di convergenza pari a 1).

Il comportamento del metodo è visualizzato dal seguente m-file :     presenta_bisezione.m

Notare che il numero di iterazioni può essere impostato a-priori, in base alla tolleranza d’errore ammessa , toll , dall’ampiezza dell’intervallo iniziale,  b-a  . Ciò è possibile perchè:
| I(k) | = (1/2)k | b-a |

| e(k) |  <  ½ | I(k) |

nitermin > log2( (b-a)/toll ) -1

Una possibile implementazione del metodo è contenuta nel seguente m-file :     metodo_di_bisezione.m

Esercizio:  

      provare il metodo nei due casi seguenti :
  1. f(x) = 5-3*x                   
  2. f(x) = log(x)-0.001        
     avendo prima determinato (ad es. per via grafica) un intervallo che contenga sicuramente una e una sola soluzione.
      NOTA:  è tipico dei problemi non-lineari il richiedere una certa conoscenza a-priori sulla localizzazione della soluzione; un approccio puramente numerico a                         partire da una aprossimazione iniziale qualunque incorre generalmente in problemi di minimi locali oppure di convergenza verso una soluzione diversa                     da quella desiderata.




2.    Metodo di Newton :

Il metodo di Newton sfrutta più informazioni su f(x) di quanto faccia il metodo di bisezione, ed infatti la sua velocità di convergenza è maggiore (nel caso di una radice semplice esso ha ordine pari a 2); però la convergenza in generale non è garantita per ogni scelta di x(0) : si dice che il metodo non è globalmente convergente .
Supponiamo che f(x) sia derivabile, e scriviamo l’equazione della retta tangente nel generico punto x(k):

y(x) = f(x(k)) + f ’(x(k))*(x-x(k))

Se scegliamo x(k+1) come quel valore di x tale per cui y(x(k+1))=0, che equivale ad approssimare f(x) localmente con la retta tangente in x(k) ed a trovare lo zero di questa , otteniamo la formula seguente per x(k+1) :

x(k+1) = x(k)  –  f(x(k))  / f ’ (x(k))

Ovviamente deve essere f ‘ (x(k))  ≠  0 .
E’ necessario sapere quando arrestare le iterazioni; un metodo che dà buoni risultati in questo caso è il controllo dell’incremento : le iterazioni si arrestano quando è verificata la

| x(k) – x(k-1) | < toll

Notare che anche per il metodo di Newton, nel caso in cui ci sia convergenza, il numero di iterazioni dipende da toll , ma anche dall’andamento di f(x), a differenza di quanto accadeva per il metodo di bisezione.

Il funzionamento del metodo è esposto nel seguente m-file :    presenta_Newton.m

Una possibile implementazione del metodo è la seguente:        metodo_di_Newton.m

Esercizio:  individuare per quale classe di funzioni f(x) il metodo di Newton converge in una sola iterazione.





Esercizio:  combinazione dei due metodi

provare a combinare i due metodi, bisezione e Newton, in modo da ottenere un metodo avente convergenza certa ed una velocità superiore alla semplice bisezione.
Suggerimento:    partendo da un intervallo [a,b] ed un valore iniziale, es x(0)=a, ripetere iterativamente la seguente procedura:
Confrontare su un paio di esempi gli algoritmi di bisezione, di Newton, quello combinato e la routine "fsolve" di Matlab/Octave (che implementa l'algoritmo di Brent, attualmente l'algoritmo più accreditato per tovare lo zero di una funzione non-lineare in un intervallo).