Metodo di Brent

Nell'analisi numerica, il metodo di Brent è un algoritmo che permette il calcolo numerico della radice di una funzione combinando il metodo della bisezione, il metodo delle secanti e l'interpolazione quadratica inversa. L'algoritmo ha la robustezza della bisezione ma può essere veloce come altri metodi meno sicuri. Il metodo prova ad usare il potenzialmente veloce metodo delle secanti o dell'interpolazione quadratica inversa se possibile, ma può ripiegare sulla più robusta bisezione se necessario. Il metodo di Brent, come suggerisce il nome, è dovuto a Richard Brent e costruito su un precedente algoritmo di Theodorus Dekker. Di conseguenza, l'algoritmo è conosciuto anche come metodo di Brent–Dekker.

Metodo di Dekker

L'idea di combinare il metodo della bisezione con quello delle secanti risale a Dekker (1969).

Si suppoga che si voglia risolvere l'equazione $f(x)=0$. Come per il metodo della bisezione, serve inizializzare il metodo di Dekker con due punti, $a_{0}$ e $b_{0}$, tali che $f(a_{0})$ e $f(b_{0})$ abbiano segno opposto. Se $f(x)$ è continua in $[a_{0},b_{0}]$, il teorema di Bolzano garantisce l'esistenza di una radice fra $a_{0}$ e $b_{0}$. In ogni iterazione sono coinvolti tre punti:

• $b_{k}$ è la iterata corrente, cioè l'approssimazione temporanea dello zero di $f$.
• $a_{k}$ è il "contro-punto", cioè un punto tale che $f(a_{k})$ e $f(b_{k})$ abbiano segno opposto, cosicché l'intervallo$[a_{0},b_{0}]$ contenga la soluzione. Oltretutto, $|f(b_{k})|$ dovrebbe essere minore o uguale di $|f(a_{k})|$, in modo che $b_{k}$ sia una migliore approssimazione della soluzione rispetto a $a_{k}$.
• $b_{k-1}$ è la precedente iterata (per la prima iterazione, si pone $b_{k-1}=a_{0}$).

Ad ogni iterazione vengono calcolati due valori provvisori. Il primo è dato dall'interpolazione lineare, anche conosciuto come il metodo delle secanti:

$s={\begin{cases}b_{k}-{\frac {b_{k}-b_{k-1}}{f(b_{k})-f(b_{k-1})}}f(b_{k}),&{\mbox{se }}f(b_{k})\neq f(b_{k-1})\\m&{\mbox{altrimenti }}\end{cases}}$

e il secondo è dato dal metodo della bisezione

$m={\frac {a_{k}+b_{k}}{2}}.$

Se il risultato delle secanti, $s$, giace strettamente tra $b_{k}$ e $m$, allora diventa la prossima iterata ($b_{k+1}=s$), altrimenti viene preso il punto medio ($b_{k+1}=m$).

Dopo, si sceglie il valore del nuovo contro-punto tale che $f(a_{k+1})$ e $f(b_{k+1})$ abbiano segno diverso. Se quest'ultimi hanno segni opposti, allora il contro-punto rimane lo stesso: $a_{k+1}=a_{k}$. Altrimenti, $f(b_{k+1})$ e $f(b_{k})$ hanno segno diverso e perciò il nuovo contro-punto diventa $a_{k+1}=b_{k}$.

Infine, se $|f(a_{k+1})|<|f(b_{k+1})|$, allora $a_{k+1}$ è probabilmente una migliore approssimazione della soluzione rispetto a $b_{k+1}$, e quindi vengono scambiati i valori di $a_{k+1}$ e $b_{k+1}$.

Questo conclude la descrizione di una singola iterazione del metodo di Dekker.

Il metodo di Dekker funziona bene se $f$ si "comporta" ragionevolmente bene. Tuttavia, ci sono delle circostanze in cui ogni iterazione utilizza il metodo delle secanti, ma $b_{k}$ converge molto lentamente (in particolare, $|b_{k}-b_{k-1}|$ può essere arbitrariamente piccolo). In questo caso, il metodo di Dekker richiede molte più iterazione dell'algoritmo di bisezione.

Metodo di Brent

Brent (1973) propose una piccola modifica per evitare questo problema, infatti inserì un test aggiuntivo che deve essere soddisfatto prima che il risultato del metodo delle secanti sia accettato per l'iterazione successiva. Devono essere soddisfatte contemporaneamente due disuguaglianze: Data una specifica tolleranza numerica $\delta$, se il passo precedente ha utilizzato il metodo della bisezione, deve valere entrambe la disuguaglianza

$|\delta |<|b_{k}-b_{k-1}|,\quad |\delta |<|s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k}-b_{k-1}|$

per eseguire l'interpolazione, altrimenti viene usata nuovamente la bisezione.

Se il passo precedente ha usato l'interpolazione, invece si usano

$|\delta |<|b_{k-1}-b_{k-2}|\quad |s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k-1}-b_{k-2}|$

per decidere se eseguire l'interpolazione (quando le disuguaglianze sono entrambe soddisfatte) o la bisezione (in caso contrario).

Questa modifica assicura che alla $k$-esima iterazione, il metodo della bisezione sarà utilizzato in al massimo $2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )$ successive iterazioni, poiché le precedenti condizioni obbligano le $|b_{k}-b_{k-1}|$ a dimezzarsi ogni due iterazioni. Così dopo al massimo $2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )$ iterazioni, $|b_{k}-b_{k-1}|$ sarà più piccolo di $\delta$ e quindi viene eseguita una bisezione. Brent dimostrò che questo metodo richiede al massimo $N^{2}$ iterazioni, dove $N$ indica il numero di volte in cui metodo di bisezione è stato utilizzato. Se la funzione $f$ non è , allora il metodo di Brent viene usato con sia la interpolazione lineare che quella quadratica inversa, e in tal caso l'algoritmo avrà una velocità di convergenza superlineare.

Oltretutto, il metodo di Brent usa l'interpolazione quadratica inversa invece di quella lineare (come usato nel metodo delle secanti). Se $f(b_{k})$, $f(a_{k})$ e $f(b_{k-1})$ sono distinti, l'efficienza aumenta leggermente. Di conseguenza, la condizione per accettare $s$ (il valore proposto sia dall'interpolazione lineare o quadratica) deve essere cambiata: $s$ deve trovarsi tra $(3a_{k}+b_{k})/4$ e $b_{k}$.

Algoritmo

Di seguito si fornisce lo pseudocodice di questo metodo:

input a, b, e (un puntatore a) una funzione f
calcola f(a)
calcola f(b)
if f(a)f(b) ≥ 0 then esci dalla funzione perché la radice potrebbe non essere di (a,b).
if |f(a)| < |f(b)| then scambia (a,b) end if
c := a
set mflag
repeat until f(b o s) = 0 or |b − a| è abbastanza piccolo (convergenza)
  if f(a) ≠ f(c) and f(b) ≠ f(c) then

    $s:={\frac {af(b)f(c)}{(f(a)-f(b))(f(a)-f(c))}}+{\frac {bf(a)f(c)}{(f(b)-f(a))(f(b)-f(c))}}+{\frac {cf(a)f(b)}{(f(c)-f(a))(f(c)-f(b))}}$  (interpolazione quadratica inversa)

  else

    $s:=b-f(b){\frac {b-a}{f(b)-f(a)}}$ (metodo delle secanti)

  end if
  if (condizione 1) s non è tra ${\frac {3a+b}{4}}$ e b or
     (condizione 2) (mflag è impostata and $|s-b|\geq |b-c|/2$) or
     (condizione 3) (mflag è libera and $|s-b|\geq |c-d|/2$) or
     (condizione 4) (mflag è impostata and $|b-c|<|\delta |$) or
     (condizione 5) (mflag è libera and $|c-d|<|\delta |$)
  then

    $s:={\frac {a+b}{2}}$ (metodo della bisezione)

    set mflag
  else
    clear mflag
  end if
  calcola f(s)
  d := c  (d qui è stata assegnata per la prima volta; non sarà usata sopra nella prima iterazione perché mflag è impostata)
  c := b
  if f(a)f(s) < 0 then b := s else a := s end if
  if |f(a)| < |f(b)| then scambia (a,b) end if
end repeat
output b o s (restituisce la radice)

Esempio

Si supponga che si stia cercando uno zero della funzione $f(x)=(x+3)(x-1)^{2}$.

Si prende $[a_{0},b_{0}]=[-4,4/3]$ come intervallo iniziale.

Si ha $f(a_{0})=-25$ e $f(b_{0})=0.48148$ (tutti i numeri in questa sezione sono arrotondati), perciò le condizioni $f(a_{0})f(b_{0})<0$ e $|f(b_{0})|\leq |f(a_{0})|$ sono soddisfatte.

Grafico di $f(x)=(x+3)(x-1)^{2}$

1. Nella prima iterazione, si utilizza l'interpolazione lineare tra $(b_{-1},f(b_{-1}))=(a_{0},f(a_{0}))=(-4,-25)$ e $(b_{0},f(b_{0}))=(1.33333,0.48148)$, che produce $s=1.23256$. Poiché giace tra $(3a_{0}+b_{0})/4$ e $b_{0}$, il valore è accettabile. Oltretutto, $f(1.23256)=0.22891$, così si imposta $a_{1}=a_{0}$ e $b_{1}=s=1.23256$.
2. Nella seconda iterazione, si utilizza l'interpolazione quadratica inversa tra $(a_{1},f(a_{1}))=(-4,-25)$, $(b_{0},f(b_{0}))=(1.33333,0.48148)$ e $(b_{1},f(b_{1}))=(1.23256,0.22891)$. Questo porta a 1.14205, che è compreso tra $(3a_{1}+b_{1})/4$ e $b_{1}$. Inoltre, la disuguaglianza $|1.14205-b_{1}|\leq |b_{0}-b_{1}|/2$ è soddisfatta, e quindi questo valore è accettabile. Oltretutto, $f(1.14205)=0.083582$, perciò si imposta $a_{2}=a_{1}$ e $b_{2}=1.14205$.
3. Al terzo passo, si utilizza l'interpolazione quadratica inversa tra $(a_{2},f(a_{2}))=(-4,-25)$, $(b_{1},f(b_{1}))=(1.23256,0.22891)$ e $(b_{2},f(b_{2}))=(1.14205,0.083582)$. Si ottiene così 1.09032, che giace tra $(3a_{2}+b_{2})/4$ e $b_{2}$. Ma qui si deve tenere conto della condizione aggiuntiva di Brent: la disuguaglianza $|1.09032-b_{2}|\leq |b_{1}-b_{0}|/2$ non è soddisfatta, perciò il valore viene scartato. Invece si considera il punto medio $m=-1.42897$ dell'intervallo $[a_{2},b_{2}]$, e si ha $f(m)=9.26891$, dunque si mette $a_{3}=a_{2}$ e $b_{3}=-1.42897$.
4. Nella quarta iterazione, si utilizza l'interpolazione quadratica inversa tra $(a_{3},f(a_{3}))=(-4,-25)$, $(b_{2},f(b_{2}))=(1.14205,0.083582)$ e $(b_{3},f(b_{3}))=(-1.42897,9.26891)$, restituendo $1.15448$. Quest'ultimo non giace tra $(3a_{3}+b_{3})/4$ e $b_{3}$, quindi è sostituito dal punto medio $m=-2.71449$. Di conseguenza si ha $f(m)=3.93934$, quindi si imposta $a_{4}=a_{3}$ e $b_{4}=-2.71449$.
5. Nella quinta interazione, usando l'interpolazione quadratica inversa, si ottiene $-3.45500$, che giace nell'intervallo richiesto. Tuttavia, la precedente iterazione ha usato la bisezione, perciò deve essere soddisfatta la disuguaglianza $|3.45500-b_{4}|\leq |b_{4}-b_{3}|/2$. Quest'ultima è falsa, quindi si usa il punto medio $m=-3.35724$. Si ricava $f(m)=-6.78239$, allora $m$ diventa il nuovo contro-punto ($a_{5}=-3.35724$) e l'iterata rimane la stessa ($b_{5}=b_{4}$).
6. Nella sesta iterazione, non si può usare l'interpolazione quadratica perché $b_{5}=b_{4}$. Quindi, si usa l'interpolazione lineare tra $(a_{5},f(a_{5}))=(-3.35724,-6.78239)$ e $(b_{5},f(b_{5}))=(-2.71449,3.93934)$. Il risultato è $s=-2.95064$, che soddisfa tutte le condizioni. Ma poiché l'iterata non è cambiata nel passo precedente, si rifiuta il risultato e si ritorna alla bisezione. Si aggiorna così $s=-3.03587$, e si ha $f(s)=-0.58418$.
7. Nella settima iterazione, si può nuovamente usare l'interpolazione quadratica inversa. Il risultato è $s=-3.00219$, che soddisfa tutte le condizioni. Ora, $f(s)=-0.03515$, quindi si imposta $a_{7}=b_{6}$ e $b_{7}=-3.00219$ ($a_{7}$ e $b_{6}$ sono scambiati in modo che la condizione $|f(b_{7})|\leq |f(a_{7})|$ sia soddisfatta).
8. Nell'ottavo passo del metodo, no si può usare l'interpolazione quadratica perché $a_{7}=b_{6}$. Dall'interpolazione lineare si ricava $s=-2.99994$, che è accettabile.
9. Nelle successive iterazioni, ci si avvina molto rapidamente alla radice $x=-3$: $b_{9}=-3+6\cdot 10^{-8}$ e $b_{1}0=-3-3\cdot 10^{-15}$ (E si ha alla 9ª iterazione $f(s)=-1.4\cdot 10^{-7}$ e $f(s)=6.96\cdot 10^{-12}$ in quella successiva).

Implementazioni

• Brent (1973) pubblicò una implementazione in Algol 60.
• Netlib contiene una traduzione in Fortran di questa implementazione con leggere modifiche.
• Il metodo solve di PARI/GP implementa il metodo di Brent.
• Altre implementazione dell'algoritmo (in C++, C, e Fortran) si trovano nella serie di libri Numerical Recipes.
• La libreria matematica di Apache Commons implementa l'algoritmo in Java.
• Il modulo di ottimizzazione SciPy implementa l'algoritmo in Python.
• La libreria standard di Modelica implementa il metodo in Modelica.
• La funzione optimize implementa l'algoritmo in R (software).
• Le librerie boost implementano l'algoritmo in C++ nel toolkit di matematica ("Locating function minima").
• Il pacchetto Optim.jl implementa l'algoritmo in Julia

Bibliografia

• Atkinson, Kendall E. (1989). "Sezione 2.8.". An Introduction to Numerical Analysis (2ª ed.). John Wiley and Sons. ISBN 0-471-50023-2.
• Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. (2007). "Section 9.3. Van Wijngaarden–Dekker–Brent Method" Archiviato l'11 agosto 2011 in Internet Archive.. Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8.

Collegamenti esterni

• zeroin.f in Netlib.
• modulo Brent in C++ (anche C, Fortran, Matlab) Archiviato il 5 aprile 2018 in Internet Archive. di John Burkardt
• GSL implementazione.
• Boost C++ implementazione.
• Metodo di brent Archiviato il 21 aprile 2018 in Internet Archive.