Date post: | 01-May-2015 |
Category: |
Documents |
Upload: | violetta-bartoli |
View: | 216 times |
Download: | 1 times |
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