Soluzione di equazioni non lineari

1. Soluzione di equazioni non lineari

2. Introduzione • Una qualsiasi equazione può sempre scriversi nella forma • 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: • in questo caso non è possibile “intrappolare la soluzione” • occorre partire da un valore di prova e arrivare alla soluzione in maniera iterativa

3. 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à

4. Alcuni esempi y f(x) continua f(x) discontinua e limitata f(a) y f(a) b a x b f(b) y a x f(x) discontinua e non limitata f(b) f(a) b a x f(b)

5. 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 x2per cui f(x1)f(x2)<0

6. 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:

7. 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à “floatingpoint” 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

8. Esempi (1) • Consideriamo l’equazione: • 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

9. 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

10. 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

11. 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:

12. 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 2 f(x2) 3 f(x3) x1=a x4 x5 x x3 x2=b f(x5) f(x4) 5 4 f(x1) 1

13. Metodo della falsa posizione • A differenza del metodo della secante, si richiede che la funzione assuma sempre segni opposti negli estremi dell’intervallo y 3 2 f(x3) f(x2) 4 f(x4) 5 f(x5) x1=a x x5 x4 x3 x2=b f(x1) 1

14. Esempio • Consideriamo l’equazione: • 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)

15. 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

16. 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

17. 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:

18. Interpolazione quadratica (2) • Poniamo: • Si ha: dove abbiamo supposto che x=b sia la stima della soluzione ottenuta nell’iterazione precedente

19. Interpolazione quadratica (3) • Poniamo: • Si ha: • Osserviamo che:

20. 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 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

21. 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

22. Metodo di Neton-Raphson (2) y 1 f(x1) 2 f(x2) 3 f(x3) 4 f(x4) 5 f(x5) x x5 x3 x4 x2 x1

23. 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à:

24. 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

25. Esempio • Consideriamo l’equazione: • 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

26. Casi problematici y y • 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 1 2 1 x x 2

27. 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

28. 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 (“tentativeroots” )del polinomio originale, e cercare di migliorarne la precisione (“polishing”)

29. 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: • 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)

30. Metodo di Muller (2) • Si ha: • Utilizziamo le seguenti notazioni: • Effettuando le sostituzioni si ha: • da cui si ottiene:

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

32. 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:

33. 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

34. Metodo di Laguerre (1) • Un polinomio Pn(x) di grado n può essere fattorizzato nel modo seguente: dove x1, x2, ..., xn sono le radici • Osserviamo che: • Calcoliamo le derivate prima e seconda di ln|Pn(x)|:

35. 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:

36. 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:

37. 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)

38. Sistemi di equazioni non lineari • Consideriamo 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 f(x,y)=0 g(x,y)=0

39. 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:

40. 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+x • La ricerca della soluzione viene fatta partendo da una soluzione di prova e procedendo iterativamente