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^[1] e costruito su un precedente algoritmo di Theodorus Dekker.^[2] Di conseguenza, l'algoritmo è conosciuto anche come metodo di Brent–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 ${\displaystyle f(x)=0}$. Come per il metodo della bisezione, serve inizializzare il metodo di Dekker con due punti, ${\displaystyle a_{0}}$ e ${\displaystyle b_{0}}$, tali che ${\displaystyle f(a_{0})}$ e ${\displaystyle f(b_{0})}$ abbiano segno opposto. Se ${\displaystyle f(x)}$ è continua in ${\displaystyle [a_{0},b_{0}]}$, il teorema di Bolzano garantisce l'esistenza di una radice fra ${\displaystyle a_{0}}$ e ${\displaystyle b_{0}}$. In ogni iterazione sono coinvolti tre punti:

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

${\displaystyle 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

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

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

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

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

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

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

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 ${\displaystyle \delta }$, se il passo precedente ha utilizzato il metodo della bisezione, deve valere entrambe la disuguaglianza

${\displaystyle |\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

${\displaystyle |\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 ${\displaystyle k}$-esima iterazione, il metodo della bisezione sarà utilizzato in al massimo ${\displaystyle 2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )}$ successive iterazioni, poiché le precedenti condizioni obbligano le ${\displaystyle |b_{k}-b_{k-1}|}$ a dimezzarsi ogni due iterazioni. Così dopo al massimo ${\displaystyle 2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )}$ iterazioni, ${\displaystyle |b_{k}-b_{k-1}|}$ sarà più piccolo di ${\displaystyle \delta }$ e quindi viene eseguita una bisezione. Brent dimostrò che questo metodo richiede al massimo ${\displaystyle N^{2}}$ iterazioni, dove ${\displaystyle N}$ indica il numero di volte in cui metodo di bisezione è stato utilizzato. Se la funzione ${\displaystyle 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 ${\displaystyle f(b_{k})}$, ${\displaystyle f(a_{k})}$ e ${\displaystyle f(b_{k-1})}$ sono distinti, l'efficienza aumenta leggermente. Di conseguenza, la condizione per accettare ${\displaystyle s}$ (il valore proposto sia dall'interpolazione lineare o quadratica) deve essere cambiata: ${\displaystyle s}$ deve trovarsi tra ${\displaystyle (3a_{k}+b_{k})/4}$ e ${\displaystyle b_{k}}$.

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

    ${\displaystyle 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

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

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

    ${\displaystyle 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)

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

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

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

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

• 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

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

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

Portale Matematica: accedi alle voci di Wikipedia che trattano di matematica