Introduzione
 Una qualsiasi equazione può sempre scriversi nella forma
f ( x)  0
 in generale non è detto che la soluzione sia reale e che sia unica
 nel caso unidimensionale è possibile effettuare la ricerca della soluzione
“intrappolandola” in un intervallo e procedendo per approssimazioni
successive
 Se il problema è multidimensionale, avremo un sistema di equazioni:
 f1 x1 , x2 , , xN   0

 f 2 x1 , x2 ,, xN   0
f ( x)  0  



 f N x1 , x2 ,, xN   0
 in questo caso non è possibile “intrappolare la soluzione”
 occorre partire da un valore di prova e arrivare alla soluzione in maniera
iterativa
Intrappolamento (“bracketing”)
 Una radice dell’equazione f(x)=0 è intrappolata
(“bracketed”) nell’intervallo [a,b] se f(a) e f(b) hanno
segni opposti
 se f(x) è continua, allora nell’intervallo (a,b) c’è almeno
una radice dell’equazione (teorema del valor medio)
 se f(x) è discontinua ma limitata, invece di una radice
potrebbe esserci una discontinuità a gradino attraverso
lo zero
 se f(x) è discontinua e non limitata, potrebbe esserci una
singolarità

in questo caso, un algoritmo di ricerca delle radici potrebbe
convergere proprio alla singolarità
Alcuni esempi
y
f(a)
y
f(a)
f(x) continua
f(x) discontinua
e limitata
b
b
a
f(b)
x
a
x
f(b)
y
f(x) discontinua
e non limitata
f(a)
a
f(b)
b
x
Intrappolamento
 Il punto di partenza per la soluzione dell’equazione
f(x)=0 è l’intrappolamento
 Il problema è quello di cercare due valori a e b tali che
f(a) e f(b) abbiano segno opposto
 In alcuni casi tali valori di a e b non sono noti, ma
vanno ricercati
 si può procedere per estrapolazione
 dati due valori di partenza x1 e x2, l’intervallo [x1,x2] viene
esteso in maniera geometrica finché non si arriva a ottenere
un intervallo per cui f(x1)f(x2)<0
 alternativamente, si può procedere per interpolazione
 dati due valori di partenza x1 e x2, l’intervallo [x1,x2] viene
suddiviso in N parti e si cerca una coppia di valori x1 e x2 per
cui f(x1)f(x2)<0
Metodo della bisezione (1)
 Si parte da un intervallo [a,b] che intrappoli la radice
dell’equazione f(x)=0
 f(a)f(b)<0
 Si calcola il valore della funzione nel punto medio xm
dell’intervallo [a,b]
 se f(xm) ha lo stesso segno di f(a) si pone a’=xm e b’=b
 se f(xm) ha lo stesso segno di f(b) si pone a’=a e b’=xm
 Si ripete la procedura, dimezzando l’intervallo [a’,b’]
 Dopo ogni iterazione l’ampiezza dell’intervallo
contenente la radice viene dimezzata:
1
 n 1   n
2
Metodo della bisezione (2)
 Detta 0=b-a l’ampiezza dell’intervallo di partenza, e detta 
la precisione richiesta per la soluzione, si può calcolare il
numero n di iterazioni necessarie ad ottenere la soluzione
con la precisione :

0
2n
 n  log 2
0

 Con quale precisione è possibile ottenere la soluzione?
 Nella modalità “floating point” il computer usa sempre un
numero fisso di cifre per rappresentare i numeri


per esempio, richiedere =10-6 può aver senso se la radice dell’equazione
è un numero dell’ordine di 1, ma non ha senso se è un numero
dell’ordine di 1026!
si potrebbe usare un criterio basato sul rapporto (xm+1-xm)/xm, ma tale
criterio non avrebbe senso per soluzioni vicine allo zero
 solitamente si sceglie una precisione pari al prodotto della
precisione del computer per la semiampiezza dell’intervallo di
partenza
Esempi (1)
 Consideriamo l’equazione:
x  log x  2  0
 dove f(x)=x+logx-2
 osserviamo che f(1)=-1<0 e f(2)=log2>0
 Cerchiamo la soluzione con una precisione =10-5
 la procedura di bisezione converge dopo 18 iterazioni e il risultato è x=1.55715
Esempi (2)
 Consideriamo l’equazione:
2  e x  x  0
 dove f(x)=2-e-x-x1/2
 osserviamo che f(0)=+1>0 e f(4)=-e-4<0
 Cerchiamo la soluzione con una precisione =10-5
 la procedura di bisezione converge dopo 20 iterazioni e il risultato è x=3.92112
Esempi (3)
 Consideriamo ora la funzione:
f ( x) 
2x  5
x 1
e applichiamo il metodo di bisezione all’intervallo [0,2]
 f(0)=5>0; f(2)=-1<0
 f(x) ha una singolarità in x=1
 Il metodo di bisezione in questo caso restituisce la singolarità
 se si cerca la soluzione dell’equazione f(x)=0 con una precisione =10-5, il metodo
di bisezione restituisce il valore x=0.99998 dopo 19 iterazioni
Metodi della secante e della falsa posizione
 In entrambi i metodi si assume che f(x) sia
approssimativamente lineare nella regione di interesse
 Entrambi i metodi sono validi per funzioni continue
 Si parte da due punti a e b con f(a) e f(b) di segni opposti
 si traccia la secante che passa per il punto [a,f(a)] e per il punto
[b,f(b)]
 si considera il punto in cui la secante interseca l’asse delle x
 nel metodo della secante, ciascuna iterazione riparte dal punto
determinato nella precedente iterazione
si considera l’intervallo [xn-1,xn] senza curarsi dei segni di f(xn-1) e f(xn)
 nel metodo della falsa posizione, si richiede che la radice
dell’equazione sia sempre “intrappolata” nell’intervallo in
esame
 il metodo della secante converge più rapidamente
 si dimostra che vale la relazione:
lim  k 1  c  k
k 
1.618
Metodo della secante
y
2
f(x2)
3
f(x3)
x1=a
f(x5) f(x )
4
f(x1)
x4
4
x5
5
x3
x2=b
x
1
 Si procede iterativamente considerando le secanti tra il
punto n-esimo e quello (n+1)-esimo
 possono esserci problemi se f(xn+1)=f(xn)
Metodo della falsa posizione
y
3
f(x3)
f(x2)
f(x4)
f(x5)
2
4
5
x1=a
x5
f(x1)
x4
x3
x2=b
x
1
 A differenza del metodo della secante, si richiede che la
funzione assuma sempre segni opposti negli estremi
dell’intervallo
Esempio
 Consideriamo l’equazione:
x  log x  2  0
 dove f(x)=x+logx-2
 osserviamo che f(1)=-1<0 e f(2)=log2>0
 Cerchiamo la soluzione con una precisione =10-5
 la procedura della falsa posizione converge dopo 6 iterazioni e il risultato è
x=1.55715 (contro le 18 iterazioni richieste dal metodo di bisezione)
Alcune considerazioni
 Il metodo della secante e quello della falsa posizione in genere
convergono più rapidamente di quello di bisezione
 Ci sono però particolari tipi di funzione per cui entrambi i
metodi possono richiedere molte iterazioni per raggiungere la
convergenza
y
x
Metodo di Brent
 Il metodo di Van Wijngaarden-Dekker-Brent combina
la bisezione con l’interpolazione quadratica
 nell’interpolazione quadratica si usano tre punti per
approssimare la funzione f(x) con una funzione
quadratica inversa (si esprime x come funzione
quadratica di y)
 il valore di x calcolato per y=0 viene preso come stima
della radice dell’equazione f(x)=0
 nel metodo di Brent si verifica se la radice così calcolata
cade all’interno o all’esterno dell’intervallo di partenza

se la radice calcolata con l’interpolazione quadratica cade al di
fuori dell’intervallo di partenza, si procede per bisezione
Interpolazione quadratica (1)
 Applichiamo la formula di Lagrange e calcoliamo il
polinomio di secondo grado x(y) che passa per i tre
punti (a,f(a)), (b,f(b)), (c,f(c)):
y  f (a)y  f (b)c  y  f (b)y  f (c)a
 f (c)  f (a) f (c)  f (b)  f (a)  f (b) f (a)  f (c)
y  f (c)y  f (a)b

 f (b)  f (c) f (b)  f (a)
x( y ) 
 Ponendo y=0 si ha:
f (a) f (b)c
f (b) f (c)a

 f (c)  f (a) f (c)  f (b)  f (a)  f (b) f (a)  f (c)
f (c) f (a)b

 f (b)  f (c) f (b)  f (a)
x
Interpolazione quadratica (2)
 Poniamo:
f (b)
R
f (c )
f (b)
S
f (a)
f (a)
T
f (c )
R  ST
 Si ha:
f (b) f (b)
f (a)
b
a
Tf (c) Rf (c)c
R
S
T
x


f (c)1  R  f (c)1  T 
 1
 1
 1
f (a)1  S  f (a) 1   f (b) 1   f (b) 1  
 T
 R
 S
cRT
aS
b
 x b 


b
R  1T  1 S  1T  1 R  1S  1
Sf (a)
dove abbiamo supposto che x=b sia la stima della
soluzione ottenuta nell’iterazione precedente
Interpolazione quadratica (3)
 Poniamo:
Q  R 1S 1T 1
 Si ha:
x b 
cRT S  1  aS R  1  bT  11  R  1S  1 P

Q
Q
 Osserviamo che:
cRT (S  1)  cST 2 S  1  cST R  T 
bT  11  R  1S  1  bT  1R  S  RS  
bT  1ST  S  RS   bS T  11  R  T  
 bST R  T   bS T  1  R  T   bST R  T   bS 1  R 
Interpolazione quadratica (4)
 Dalle precedenti relazioni segue che:
P  cST R  T   aS 1  R   bST R  T   bS 1  R 
 S T R  T c  b   1  R b  a 
 Le operazioni algebriche su P servono ad esprimere tale quantità in
termini delle differenze (c-b) e (b-a)
 Questa espressione di P risulta utile quando deve essere
implementata in una routine
 In conclusione si ha:
x b
P
Q
 In altri termini, la nuova soluzione è ottenuta correggendo la
vecchia soluzione b di una quantità P/Q che si suppone piccola
 se la correzione P/Q non è piccola, o se l’intervallo ottenuto con la
nuova soluzione non si restringe abbastanza rapidamente,
l’algoritmo di Brent effettua una bisezione
Metodo di Newton-Raphson (1)
 A differenza dei metodi visti in precedenza, il metodo
di Newton richiede il calcolo sia della funzione f(x) che
della sua derivata f’(x)
 Dato un punto P(x,f(x)), si considera la tangente al
grafico della funzione nel punto P e la si prolunga fino
ad incontrare l’asse delle ascisse
 Il punto di intersezione della tangente con l’asse x
rappresenta la nuova stima della soluzione
Metodo di Neton-Raphson (2)
y
1
f(x1)
2
f(x2)
3
f(x3)
f(x4)
f(x5)
4
5
x5 x4 x3 x2 x1
x
Metodo di Newton-Raphson (3)
 In un intorno di un punto x si può scrivere:
f ( x   )  f ( x)  f ( x)
 Imponendo la condizione f(x+)=0 si può ricavare :
f ( x)
f (x   )  0    
f ( x)
 Dette xi e xi+1 le soluzioni dell’equazione dopo la i-esima
e la (i+1)-esima iterazione, si avrà:
f ( xi )
xi 1  xi   i  xi 
f ( xi )
Metodo di Newton-Raphson (4)
 Sviluppando la funzione f(x) in serie di Taylor in un intorno
del punto x in cui f(x)=0 si ha:
f ( x   )  f ( x)  f ( x) 
2
f ( x)    f ( x) 
2
f ( x   )  f ( x)  f ( x)  
2
2
f ( x)  
 Posto xn=x+n e xn+1=x+ n+1 si ha:
xn 1  xn 
 n 1   n 
f ( xn )
f ( xn )
  n 1   n 

f ( xn )
f ( xn )
 n f ( x) 
 n2
2
f ( x)
f ( x)
1
f ( x)
  n 1    n2
2
f ( x)
 A differenza del metodo di bisezione, che converge
linearmente, il metodo di Newton-Raphson converge
quadraticamente
Esempio
 Consideriamo l’equazione:
x  log x  2  0
 dove f(x)=x+logx-2
 osserviamo che f(1)=-1<0 e f(2)=log2>0
 Cerchiamo la soluzione con una precisione =10-5
 la procedura della falsa posizione converge dopo 6 iterazioni e
il risultato è x=1.55715 (contro le 18 iterazioni richieste dal
metodo di bisezione)
 partendo da x=2, la procedura di Newton-Raphson converge
dopo 4 iterazioni e il risultato è x=1.55715
 la rapidità della convergenza dipende dal punto di partenza:



se si parte da x=1.5 sono sufficienti 3 iterazioni
se si parte da x=10 occorrono 6 iterazioni
se si parte da x=20 occorrono 8 iterazioni
Casi problematici
y
y
2
1
1
x
x
2
 La procedura di Newton-Raphson può fallire quando:
 in una iterazione si incontra un massimo (o un minimo) locale di f(x)

in questo caso la tangente alla curva è orizzontale ed il valore di x determinato
nell’iterazione successiva è infinito
 una iterazione riconduce nel punto della iterazione precedente

in questo caso si entra in un ciclo infinito
Radici dei polinomi
 Un polinomio di grado n ha n radici
 le radici possono essere reali o complesse
 non necessariamente le radici sono distinte
 se i coefficienti del polinomio sono reali, le radici complesse
si presentano in coppie coniugate:

se a+ib è una radice del polinomio, lo è anche a-ib
 se i coefficienti del polinomio sono complessi, non è detto che
vi siano relazioni tra le radici
 Per cercare le radici di un polinomio si possono usare gli
algoritmi generali
 radici multiple (o radici fra loro molto vicine) introducono
delle difficoltà

per esempio P(x)=(x-a)2 ha una radice doppia in x=a e risulta P(x)>0
ovunque eccetto che in x=a
 in questo caso non si può intrappolare la radice
 neanche il metodo di Newton-Raphson funziona bene, perché in
un intorno di x=a sia la funzione P(x) che la sua derivata tendono
ad annullarsi
Fattorizzazione (“deflation”)
 Una volta trovata una radice r del polinomio, questo può
essere fattorizzato nella forma:
P( x)  x  r Q( x)
 Le radici di Q(x) sono le radici restanti di P(x)
 In caso di radici complesse del tipo a+ib e a-ib, il polinomio
P(x) va diviso per il fattore quadratico:
x  a  ib x  a  ib   x 2  2ax  a 2  b2 
 La fattorizzazione va usata con attenzione:
 gli errori di calcolo su una radice possono dar luogo ad errori
nella determinazione dei coefficienti del polinomio Q(x), e
quindi nella determinazione delle altre radici
 conviene trattare le radici trovate di volta in volta come radici
approssimate (“tentative roots” )del polinomio originale, e
cercare di migliorarne la precisione (“polishing”)
Metodo di Muller (1)
 Il polinomio P(x) viene interpolato con una curva di
secondo grado
 Siano xi-2, xi-1 e xi le ultime tre stime di una radice di
P(x)
 La curva quadratica che passa per i punti (xi-2, P(xi-2)),
(xi-1, P(xi-1)) e (xi, P(xi)) avrà equazione:
q( x)  ax  xi   bx  xi   c
2
 I coefficienti a,b e c si trovano imponendo che sia
q(xi)=P(xi), q(xi-1)=P(xi-1), q(xi-2)=P(xi-2)
Metodo di Muller (2)
 Si ha:
q ( xi )  P( xi )  c  P( xi )
q ( xi 1 )  P( xi 1 )  a xi 1  xi   b xi 1  xi   c  P( xi 1 )
2
q ( xi  2 )  P( xi  2 )  axi  2  xi   bxi  2  xi   c  P( xi  2 )
2
 Utilizziamo le seguenti notazioni:
h2  xi  2  xi
2 
h1  xi 1  xi
P ( xi  2 )  P( xi )
h2
1 
P( xi 1 )  P ( xi )
h1
 Effettuando le sostituzioni si ha:
ah12  bh1  P( xi )  P( xi 1 )  ah12  bh1  1h1  ah1  b  1
ah22  bh2  P( xi )  P( xi 2 )  ah22  bh2   2 h2  ah2  b   2
 da cui si ottiene:
a h1  h2   1   2  a 
b  1  ah1
1   2
h1  h2
Metodo di Muller (3)
 Una volta calcolati a,b e c si può trovare la stima successiva della
radice, xi+1, imponendo che q(xi+1)=0:
axi 1  xi   bxi 1  xi   c  0
2
 La precedente è un’equazione di secondo grado nell’incognita
(xi+1-xi)
 Occorre scegliere la soluzione più piccola in modulo, poiché la
stima (i+1)-esima della radice deve essere vicina alla stima
precedente
 l’equazione potrebbe essere risolta in modo tradizionale ricavando:
 b  b 2  4ac
xi 1  xi 
2a
 la soluzione più piccola in modulo è quella col segno + se b>0, ed è
quella col segno – se b<0
 in ogni caso, si deve fare una sottrazione che potrebbe portare ad un
risultato nullo (per via della precisione della macchina!)
Metodo di Muller (4)
 Conviene quindi procedere nel modo seguente:
 si pone y=1/(xi+1-xi) e si ha:
axi 1  xi   bxi 1  xi   c  0 
2
a b
2


c

0

cy
 by  a  0
2
y
y
 le cui soluzioni sono:
 b  b 2  4ac
b  b 2  4ac
y

2c
2c
 poiché (xi+1-xi) deve essere un numero di modulo piccolo, y
deve essere un numero di modulo grande
 la soluzione di modulo più grande si ottiene prendendo il
segno + se b>0, il segno – se b<0 nell’ultima formula

in ogni caso, si deve fare una somma di due quantità con lo stesso
segno, evitando problemi di cancellazione:
y
b  b / b  b 2  4ac
2c
 xi 1  xi 
2c
b  b / b  b 2  4ac
Metodo di Muller (5)
 Si parte da tre valori arbitrari x1, x2 e x3
 Si determina il valore successivo di x sfruttando la relazione
di ricorrenza
xi 1  xi 
2c
b  b / b  b 2  4ac
 Si elimina il valore di x più “vecchio” e si procede con una
nuova iterazione
 La procedura iterativa viene terminata quando le differenze
h1 e h2 sono sufficientemente piccole
 Nell’implementazione del metodo bisogna ricordare che c’è
la possibilità che le varie quantità assumano valori
complessi
 Il metodo di Muller può essere usato anche per funzioni
che non siano polinomi
Metodo di Laguerre (1)
 Un polinomio Pn(x) di grado n può essere fattorizzato
nel modo seguente:
Pn ( x)  k x  x1 x  x2 x  xn 
dove x1, x2, ..., xn sono le radici
 Osserviamo che:
ln Pn ( x)  ln k  ln x  x1  ln x  x2    ln x  xn
 Calcoliamo le derivate prima e seconda di ln|Pn(x)|:
G ( x) 
d ln Pn ( x)
H ( x)  
dx
Pn( x)
1
1
1




Pn ( x) x  x1 x  x2
x  xn
d 2 ln Pn ( x)
dx 2
P( x) P( x)  Pn( x)
 n

2
Pn ( x)
1
1
1




x  x1 2 x  x2 2
x  xn 2
2
Metodo di Laguerre (2)
 Sia x una stima approssimata della radice x1
 Si assume che la radice x1 sia ad una distanza a da x, e che le
altre radici x2,x3,...,xn siano tutte equidistanti da x, alla stessa
distanza b:
x  x1  a
x  xi  b
i  2,3,, n
 questa assunzione ha senso per valori di x effettivamente
prossimi a x1:

se |x-xi|>>|x-x1|, le distanze |x-xi| sono in pratica infinite e possono
essere considerate tutte uguali a un numero molto grande b
 Con le assunzioni di Laguerre si ha:
1 n 1
G ( x)  
a
b
1 n 1
H ( x)  2  2
a
b
Metodo di Laguerre (3)
 Per ricavare la stima di x1 occorre calcolare il valore di a
dal precedente sistema e ricavare x1=x-a. Si ha:
1
1 
1

G
(
x
)

b n  1 
a 
1
1
H ( x)  2  n  1
a
n  12
2
1

G ( x)  a  
1  2G ( x) 1 G 2 ( x)

1  n  1  n  1 a  n  1  H ( x)  0 


1
1
n 2  2G ( x)  G 2 ( x)  n  1H ( x)  0 
a
a
1
a2


1 G ( x)  G 2 ( x)  n G 2 ( x)  n  1H ( x)


a
n
n
a
G ( x)  n  1 nH ( x)  G 2 ( x)


Metodo di Laguerre (4)
 Implementazione del metodo:
 si parte da una stima x della radice
 si calcola x1=x-a risolvendo il sistema precedente
 la nuova stima di x1 è usata come punto di partenza della
iterazione successiva
 la procedura iterativa viene fermata quando a è
sufficientemente piccolo
 Nell’implementazione del metodo bisogna gestire
numeri complessi (le varie quantità possono anche
assumere valori complessi)
Sistemi di equazioni non lineari
 Consideriamo come esempio un sistema di due equazioni in due incognite:
 f ( x, y )  0

 g ( x, y )  0
 Risolvere il sistema significa trovare le intersezioni tra le linee di contorno
f(x,y)=0 e g(x,y)=0
 in generale non è possibile dire nulla sull’esistenza di soluzioni del sistema
e sulla loro eventuale molteplicità
 per un sistema di N equazioni la soluzione è l’intersezione di N iperpiani di
contorno, ciascuno di dimensione N-1
f(x,y)=0
g(x,y)=0
Metodo di Newton-Raphson
generalizzato (1)
 Sia dato un sistema di N equazioni non lineari nelle N
variabili x1,x2,...,xN:
f i x1 , x2 ,, xN   0
i  1,2,, N
 Indicato con x=(x1,x2,...,xN) il vettore delle variabili che
rappresenta una soluzione approssimata del sistema, e
con x=(x1,x2,...,xN) una piccola perturbazione di x,
in un intorno di x si ha:
f i
x j  O x 2
j 1 x j
N
f i x  x   f i ( x)  
 
Metodo di Newton-Raphson
generalizzato (2)
 Trascurando i termini di ordine superiore in x si ha:
f i
x j  0
j 1 x j
N
f i x  x   0  f i ( x)  
 Posto ij=fi/xj e i=-fi(x), il vettore x può essere
trovato risolvendo il sistema lineare:
N
 x
j 1
ij
j
 i
i  1,2,, N
 La nuova stima della soluzione è quindi x’=x+x
 La ricerca della soluzione viene fatta partendo da una
soluzione di prova e procedendo iterativamente