Il sistema ideato da Dekker e Brent è considerato uno dei migliori in circolazione e viene utilizzato anche da alcuni ambienti di calcolo conosciuti nell’ambito matematico come ad esempio MATLAB o NAG.

Pseudo codice dell’algoritmo:

a = Penultima approssimazione calcolata

b = Ultima approssimazione

c = punto più recente tale che x cade nell’intervallo [b,c]

Si parte da due valori iniziali a e b tali che sia vero f(a) * f(b) < 0

c = a

bNew = Approssimazione ottenuta con il metodo delle secanti applicato ad a e b se:

Altrimenti

bNew = Approssimazione ottenuta applicando il metodo di bisezione applicato all’intervallo [a,c]

il nuovo intervallo sarà dunque [b,bNew] oppure [bNew,c] a seconda di dove cade lo zero.

codice MATLAB:

%DBMethod Esegue il metodo di Dekker-Brent per il calcolo della radice

%   di una funzione non lineare

%

%   [i,x,tolf,nu]=DBMethod(a,b,f,tolx,iMax)

%

%   I parametri della funzione sono:

%       f -> funzione di cui valutare uno zero

%       a,b -> estremi dell'intervallo in cui si ricerca lo zero di f;

%              si richiede che f(a)f(b)<0

%       tolx -> tolleranza sulla x

%       iMax -> Numero massimo di iterazioni da eseguire per arrivare alla soluzione

%                    se iMax è zero allora viene stimato in maniera autonoma

%

%   I valori di ritorno sono:

%       x -> la soluzione trovata

%       iNumber -> il numero di iterazioni impiegate per ottenere la soluzione

%       tolf -> la tolleranza sulla funzione

%       iFail -> vale uno se i due punti (Il punto A e il punto B) hanno lo stesso segno

%                  rendendo impossibile calcolare un'eventuale zero all'interno dell'intervallo

%                  Vale 2 se iMax è <= di Zero

%

function [iFail,x,tolf,iNumber]=DBMethod(a,b,f,tolx,iMax)

   fa=feval(f,a);

   fb=feval(f,b);   

   iNumber=0;

   if (fa*(fb/abs(fb))>0)

      iFail=1;

      break

   else

      c=a;

      tolf=tolx;

      while (abs(c-b)>tolx) & (abs(fb)>tolf) & (iNumber<iMax)

         iNumber=iNumber+1;

         if (fa==fb)

            %bisezione

            [iFail,bNew,tolf,iNumberTMP]=bisezione(b,c,f,tolx,1);

         else

            %secanti

            [iFail,bNew,tolf,iNumberTMP]=secanti(a,b,f,tolx,1);

            if(bNew<b | bNew>c)

               %bisezione

               [iFail,bNew,tolf,iNumberTMP]=bisezione(b,c,f,tolx,1);

            end

         end

         a=b;

         fa=fb;

         b=bNew;

         fb=feval(f,b);         

         if (fa*(fb/abs(fb))<0)

            c=a;

         end

      end

   end

   x=c;

end