Calcolo numerico/Metodo di Brent

Struttura del metodo

È una combinazione del metodo di bisezione e delle secanti, in cui si sfruttano la sicurezza di convergenza del metodo di bisezione e la velocità di convergenza di quello delle secanti.

Il metodo richiede come input una bracket $(b,c)$ con $f(b)\ast f(c)<0$. Restituisce come output lo zero della funzione nella variabile chiamata $b$. Nell'algoritmo, $b$ rappresenta l'approssimazione della radice cercata. Il metodo si basa sull'idea che se

$|f(b)|<|f(c)|$

allora $b$ è più vicino alla radice di quanto lo sia $c$. Questo in realtà non è sempre vero, infatti se consideriamo una funzione che interseca l'asse $x$ in $(0,\alpha )$ e poi decresce lentamente dopo tale punto, si ha che, pur scegliendo $b$ lontano dalla radice, segue che $|f(b)|<|f(c)|$ ma $c$ è più vicina alla radice di quanto lo sia $b$.

Nell'algoritmo, $a,b$ rappresentano le due iterate con cui il metodo delle secanti calcola l'iterata successiva, cioè $b=x_k$, $a=x_{k-1}$, mentre $b$ e $c$ sono gli estremi della bracket.

Algoritmo

1. Definizione di $a,b,c$: in base all'idea precedente, se $|f(b)|>|f(c)|$ allora scambio fra loro $b$ e $c$. Poniamo poi $a=c$ per far partire il metodo delle secanti.

2. Test d'arresto sull'ampiezza della bracket: finché $|b-c|>\epsilon$ le istruzioni vengono eseguite.

3. Nuovo passo: ci si avvicina alla radice attraverso l'assegnamento

$b=b+d_d$

dove $d_d$ può essere posto uguale al passo del metodo delle secanti oppure uguale al passo del metodo di bisezione, e questa scelta viene valutata a ogni iterazione.

4. Valutazione dei due passi possibili: per scegliere quali dei due metodi usare, si valuta di quanto ci si sposta da $b$ in entrambi i metodi:

$d_m=(c-b)/2\text{passo di bisezione}$

Prima di calcolare $d_s$ (passo del metodo delle secanti), si valuta $d_f=f(b)-f(a)$. Se $d_f=0$, si pone $d_d=d_m$, perché $d_s$ non è applicabile. Altrimenti

$d_s=a-\frac{b-a}{f(b)-f(a)}\ast f(b)$

5. Scelta del passo: se entrambi i passi sono applicabili, si eseguono le seguenti istruzioni:

$\mathrm{i}\mathrm{f}\mathrm{(}\mathrm{s}\mathrm{g}\mathrm{n}{\mathrm{d}}_{\mathrm{m}}\mathrm{=}\mathrm{=}\mathrm{s}\mathrm{g}\mathrm{n}{\mathrm{d}}_{\mathrm{s}}\mathrm{\wedge }\mathrm{|}{\mathrm{d}}_{\mathrm{s}}\mathrm{|}<\mathrm{=}\mathrm{|}{\mathrm{d}}_{\mathrm{m}}\mathrm{|}\mathrm{)}$

$d_d=d_s$

$\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}$

$d_d=d_m$

$\mathrm{e}\mathrm{n}\mathrm{d}$

Nota: queste istruzioni si giustificano nel seguente modo. Il passo di bisezione rimane sempre nella bracket per definizione del metodo, invece, se il passo di secanti ha direzione opposta, si esce dalla bracket. Una volta assicurato che entrambi i metodi stanno dentro la bracket, si sceglie il passo di ampiezza più piccola, che fa allontanare di meno dall'approssimazione corrente della radice.

6. Controllo

7. Applicazione del passo scelto: $d=b+d_d$

8. Preparazione della terna di punti per un eventuale passo successivo:

$a_1=a$

$b_1=d$

$\mathrm{i}\mathrm{f}\mathrm{(}\mathrm{s}\mathrm{g}\mathrm{n}\mathrm{f}\mathrm{(}\mathrm{d}\mathrm{)}\mathrm{\ne }\mathrm{s}\mathrm{g}\mathrm{n}\mathrm{f}\mathrm{(}\mathrm{c}\mathrm{)}\mathrm{)}$

$c_1=c$

$\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}$

$c_1=b$

$\mathrm{e}\mathrm{n}\mathrm{d}$

Nota: se $\mathrm{s}\mathrm{g}\mathrm{n}f(d)\ne \mathrm{s}\mathrm{g}\mathrm{n}f(c)$, allora segue necessariamente che $\mathrm{s}\mathrm{g}\mathrm{n}f(d)\ne \mathrm{s}\mathrm{g}\mathrm{n}f(b)$, perché $b,c$ erano una bracket)

9. Assegnamento:

$b=b_1$

$c=c_1$

$\mathrm{i}\mathrm{f}\mathrm{(}\mathrm{|}\mathrm{f}\mathrm{(}\mathrm{b}\mathrm{)}\mathrm{|}>\mathrm{|}\mathrm{f}\mathrm{(}\mathrm{c}\mathrm{)}\mathrm{|}\mathrm{)}$

$temp=b;$

$b=c;$

$c=temp;$

$\mathrm{e}\mathrm{n}\mathrm{d}$

$a_1=c$

Metodo implementato:

function [b,it] = Mybrent(f,a,b,c)
if(abs(feval(f,b))>=abs(feval(f,c)))
temp = b;
b = c;
c = temp;
end
while(abs(b-c)>eps)
dm = (c-b)/2;
q = feval(f,a)-feval(f,b);
if(q==0)
ds = dm;
else
ds = a-(b-a)/(feval(f,b)-feval(f,a))*feval(f,b);
end
if(sign(ds)==sign(dm)&&abs(ds)<=abs(dm))
dd = ds;
else
dd = dm;
end
if(dd<eps)
d = eps/2*sign(dd);
end
d = b+dd;
b1 = d;
if(sign(b1)==sign(b))
c1 = c;
else
c1 = b;
end
b = b1;
c = c1;
if(abs(feval(f,c))<=abs(feval(f,b)))
temp = b;
b = c;
c = temp;
end
end