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 :
- f(x) = 5-3*x
- 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:
- calcolare il punto successivo x(i+1) con il metodo di Newton;
- se tale punto appartiene ad [a,b],
- prendere tale punto come punto di suddivisione di [a,b] e
proseguire con l’iterata successiva di Newton, aggiornando però
[a,b] con la regola della bisezione, ovvero diventa il sottointervallo
scelto in base
alla regola f(a)*f(b)<0;
- altrimenti,
- non considerare tale punto e procedere con un normale passo di
bisezione;
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).