IntroduzioneUna qualsiasi equazione può sempre scriversi nella forma

in generale non è detto che la soluzione sia reale e che sia unicanel 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:

in questo caso non è possibile “intrappolare la soluzione”occorre partire da un valore di prova e arrivare alla soluzione in

maniera iterativa

0)( xf

0,,,

0,,,

0,,,

0)(

21

212

211

NN

N

N

xxxf

xxxf

xxxf

xf

Intrappolamento (“bracketing”)Una radice dell’equazione f(x)=0 è

intrappolata (“bracketed”) nell’intervallo [a,b] se f(a) e f(b) hanno segni oppostise 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

x

y

ab

f(b)

f(a)

x

y

ab

f(b)

f(a)

x

y

a b

f(b)

f(a)

f(x) continua f(x) discontinua e limitata

f(x) discontinua e non limitata

IntrappolamentoIl punto di partenza per la soluzione

dell’equazione f(x)=0 è l’intrappolamentoIl problema è quello di cercare due valori a e b

tali che f(a) e f(b) abbiano segno oppostoIn 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)=0f(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’=bse 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:

nn 2

11

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 :

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

0

20 log

2 n

n

Esempi (1)Consideriamo l’equazione:

dove f(x)=x+logx-2osserviamo 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

02log xx

Esempi (2) Consideriamo l’equazione:

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

02 xe x

Esempi (3) Consideriamo ora la funzione:

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

1

52)(

x

xxf

Metodi della secante e della falsa posizioneIn entrambi i metodi si assume che f(x) sia

approssimativamente lineare nella regione di interesseEntrambi i metodi sono validi per funzioni continueSi 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 xnel metodo della secante, ciascuna iterazione riparte dal

punto determinato nella precedente iterazionesi 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ù rapidamentesi dimostra che vale la relazione:

618.1

1lim kkk

c

Metodo della secante

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)

y

x

1

2

x1=ax2=b

f(x1)

f(x2)

x3

f(x3)3

x4

f(x4) 4

x5

f(x5) 5

Metodo della falsa posizione

A differenza del metodo della secante, si richiede che la funzione assuma sempre segni opposti negli estremi dell’intervallo

x3

f(x3)3

y

x

1

2

x1=ax2=b

f(x1)

f(x2)

x4

f(x4)4

x5

f(x5)5

EsempioConsideriamo l’equazione:

dove f(x)=x+logx-2osserviamo 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)

02log xx

Alcune considerazioniIl 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 convergenzay

x

Metodo di BrentIl metodo di Van Wijngaarden-Dekker-Brent

combina la bisezione con l’interpolazione quadraticanell’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)):

Ponendo y=0 si ha:

)()()()(

)()(

)()()()(

)()(

)()()()(

)()()(

afbfcfbf

bafycfy

cfafbfaf

acfybfy

bfcfafcf

cbfyafyyx

)()()()(

)()(

)()()()(

)()(

)()()()(

)()(

afbfcfbf

bafcf

cfafbfaf

acfbf

bfcfafcf

cbfafx

Interpolazione quadratica (2)Poniamo:

Si ha:

dove abbiamo supposto che x=b sia la stima della soluzione ottenuta nell’iterazione precedente

STRcf

afT

af

bfS

cf

bfR

)(

)(

)(

)(

)(

)(

bSR

b

TS

aS

TR

cRTbx

Sbf

Rbf

bSbf

Rbf

TafSaf

aTaf

aSf

TcfRcf

ccRfcTfx

111111

11)(

11)(

)()(

11)(1)(

)()(

1)(1)(

)()(

Interpolazione quadratica (3)Poniamo:

Si ha:

Osserviamo che:

111 TSRQ

Q

P

Q

SRTbRaSScRTbx

111111

TRcSTScSTScRT 1)1( 2

RbSTRbSTTRTbSTRbST

TRTbSRSSSTTb

RSSRTbSRTb

11

111

11111

Interpolazione quadratica (4)Dalle precedenti relazioni segue che:

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:

In altri termini, la nuova soluzione è ottenuta correggendo la vecchia soluzione b di una quantità P/Q che si suppone piccolase 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

abRbcTRTS

RbSTRbSTRaSTRcSTP

1

11

Q

Pbx

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)

x

y

x1

f(x1)1

x2

f(x2)2

x3

f(x3)3

x4

f(x4)4

x5

f(x5)5

Metodo di Newton-Raphson (3)In un intorno di un punto x si può scrivere:

Imponendo la condizione f(x+)=0 si può ricavare :

Dette xi e xi+1 le soluzioni dell’equazione dopo la i-esima e la (i+1)-esima iterazione, si avrà:

)()()( xfxfxf

)(

)(0)(

xf

xfxf

)(

)(1

i

iiiii xf

xfxxx

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:

Posto xn=x+n e xn+1=x+ n+1 si ha:

A differenza del metodo di bisezione, che converge linearmente, il metodo di Newton-Raphson converge quadraticamente

)()()(

)(2

)()(2

)()()(22

xfxfxf

xfxfxfxfxfxf

)(

)(

2

1

)(

)(2

)(

)(

)(

)(

)(

21

2

1

11

xf

xf

xf

xfxf

xf

xf

xf

xfxx

nn

nn

nn

n

nnn

n

nnn

EsempioConsideriamo l’equazione:

dove f(x)=x+logx-2osserviamo 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

02log xx

Casi problematici

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 è infinitouna iterazione riconduce nel punto della iterazione precedente

in questo caso si entra in un ciclo infinito

x

y12

x

y1

2

Radici dei polinomiUn polinomio di grado n ha n radici

le radici possono essere reali o complessenon necessariamente le radici sono distintese 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 generaliradici 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:

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:

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”)

)()( xQrxxP

222 2 baaxxibaxibax

Metodo di Muller (1)Il polinomio P(x) viene interpolato con una

curva di secondo gradoSiano 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:

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)

cxxbxxaxq ii 2)(

Metodo di Muller (2) Si ha:

Utilizziamo le seguenti notazioni:

Effettuando le sostituzioni si ha:

da cui si ottiene:

)()()(

)()()(

)()()(

222

222

112

111

iiiiiii

iiiiiii

iii

xPcxxbxxaxPxq

xPcxxbxxaxPxq

xPcxPxq

1

11

2

22

1122

)()(

)()(

h

xPxP

h

xPxP

xxhxxh

iiii

iiii

222222222

22

111112

1112

1

)()(

)()(

bahhbhahxPxPbhah

bahhbhahxPxPbhah

ii

ii

11

21

212121

ahb

hhahha

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:

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:

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!)

012

1 cxxbxxa iiii

a

acbbxx ii 2

42

1

Metodo di Muller (4)Conviene quindi procedere nel modo seguente:

si pone y=1/(xi+1-xi) e si ha:

le cui soluzioni sono:

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:

000 221

21 abycyc

y

b

y

acxxbxxa iiii

acbbbb

cxx

c

acbbbby ii

4/

2

2

4/21

2

c

acbb

c

acbby

2

4

2

4 22

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

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

acbbbb

cxx ii

4/

221

Metodo di Laguerre (1)Un polinomio Pn(x) di grado n può essere

fattorizzato nel modo seguente:

dove x1, x2, ..., xn sono le radiciOsserviamo che:

Calcoliamo le derivate prima e seconda di ln|Pn(x)|:

nn xxxxxxkxP 21)(

nn xxxxxxkxP lnlnlnln)(ln 21

222

21

2

2

2

2

21

111

)(

)()()()(ln)(

111

)(

)()(ln)(

n

n

nnn

nn

nn

xxxxxx

xP

xPxPxP

dx

xPdxH

xxxxxxxP

xP

dx

xPdxG

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:

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:

nibxx

axx

i ,,3,2 1

22

11)(

11)(

b

n

axH

b

n

axG

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

0)(1)(1

)(21

0)(1

)(1

1

)(2

1

11

1

1)(

1

11

1)(

1)(

1

11

2

22

22

2

2

2

22

xGxnHnxG

na

n

xHnxGnxGxG

a

xHnxGa

xGan

xHn

xG

an

xG

na

axG

nn

axH

axG

nb

Metodo di Laguerre (4)Implementazione del metodo:

si parte da una stima x della radicesi calcola x1=x-a risolvendo il sistema

precedentela nuova stima di x1 è usata come punto di

partenza della iterazione successivala procedura iterativa viene fermata quando a è

sufficientemente piccoloNell’implementazione del metodo bisogna

gestire numeri complessi (le varie quantità possono anche assumere valori complessi)

Sistemi di equazioni non lineariConsideriamo come esempio un sistema di due equazioni in due

incognite:

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

0),(

0),(

yxg

yxf

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:

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:

Nixxxf Ni ,,2,1 0,,, 21

N

jj

j

iii xOx

x

fxfxxf

1

2)(

Metodo di Newton-Raphson generalizzato (2)Trascurando i termini di ordine superiore in x si

ha:

Posto ij=fi/xj e i=-fi(x), il vettore x può essere trovato risolvendo il sistema lineare:

La nuova stima della soluzione è quindi x’=x+xLa ricerca della soluzione viene fatta partendo da

una soluzione di prova e procedendo iterativamente

N

jj

j

iii x

x

fxfxxf

1

0)(0

N

jijij Nix

1

,,2,1