Решение уравнения F(x) = 0 методом Брента

Метод Брента - метод поиска корней, окруженных на интервале [a, b] (это значит, что F(a)F(b) < 0). Он сочетает в себе высокую скорость сходимости метода Риддлера с надежностью дихотомии. В то время, как метод Риддлера быстр, но может "споткнуться" и очень сильно замедлиться на какой-нибудь сложной функции, скорость работы дихотомии не столь высока, но не зависит от обрабатываемой функции. Метод Брента сочетает в себе одновременно два подхода - на "хороших" функциях используется трехточечная интерполяционная схема, дающая квадратичную скорость сходимости, а в случае слишком низкой скорости работы "продвинутого" подхода (или расходимости) используется деление отрезка на две части.

Если нашли ошибку в алгоритме - сообщите!

Реализации алгоритма на различных языках:

Реализация алгоритма на AlgoPascal:

unit BrentUnit;

declaration
    function F(X:Real):Real;

interface SearchRootBrent;

implementation

const
    //максимальное число итераций
    ItMax : Integer = 100;

(*************************************************************************
Поиск корня методом Брента

Если на интервале [x1,x2] непрерывная функция F(x)  удовлетворяет  условию
F(x1)F(x2)<0, то корень уравнения F(x)=0 может быть найден методом Брента.

Процесс    продолжается,  пока  корень функции не будет найден с точностью
epsilon или число итераций не станет слишком велико.

HasRoot содержит   True, если останов  произошел  из-за  нахождения корня,
(в переменной x содержится корень) и False,  если число   итераций слишком
велико или функция не принимает противоположные знаки на концах.
*************************************************************************)
procedure SearchRootBrent(
    const   x1      :   Real;
    const   x2      :   Real;
    const   Epsilon :   Real;
    out     HasRoot :   Boolean;
    out     x       :   Real);
var
   a    :   Real;
   b    :   Real;
   c    :   Real;
   d    :   Real;
   e    :   Real;
   min1 :   Real;
   min2 :   Real;
   min0 :   Real;
   fa   :   Real;
   fb   :   Real;
   fc   :   Real;
   p    :   Real;
   q    :   Real;
   r    :   Real;
   s    :   Real;
   tol1 :   Real;
   xm   :   Real;
   iter :   Integer;
begin
    a:=x1;
    b:=x2;
    fa:=F(a);
    fb:=F(b);
    fc:=fb;
    for Iter:=1 to ITMax do
    begin
        if fb*fc>0 then
        begin
            c := a;
            fc := fa;
            d := b-a;
            e := d
        end;
        if AbsReal(fc)<AbsReal(fb) then
        begin
            a:=b;
            b:=c;
            c:=a;
            fa:=fb;
            fb:=fc;
            fc:=fa
        end;
        tol1:=2*MachineEpsilon*AbsReal(b)+0.5*Epsilon;
        xm:=0.5*(c-b);
        if (AbsReal(xm)<=tol1) or (fb=0) then
        begin
            x:=b;
            HasRoot:=True;
            Exit;
        end;
        if (AbsReal(e)>=tol1) and (AbsReal(fa)>AbsReal(fb)) then
        begin
            s:=fb/fa;
            if a=c then
            begin
                p:=2.0*xm*s;
                q:=1.0-s
            end
            else
            begin
                q:=fa/fc;
                r:=fb/fc;
                p:=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
                q:=(q-1.0)*(r-1.0)*(s-1.0)
            end;
            if p>0 then
                q:=-q;
            p:=AbsReal(p);
            min1:=3*xm*q-AbsReal(tol1*q);
            min2:=AbsReal(e*q);
            if min1<min2 then
                min0:=min1
            else
                min0:=min2;
            if 2*p<min0 then
            begin
                e:=d;
                d:=p/q
            end
            else
            begin
                d:=xm;
                e:=d;
            end;
        end
        else
        begin
            d:=xm;
            e:=d;
        end;
        a:=b;
        fa:=fb;
        if AbsReal(d)>tol1 then
        begin
            b:=b+d;
        end
        else
        begin
            if xm>0 then
            begin
                b:=b+AbsReal(tol1);
            end
            else
            begin
                b:=b-AbsReal(tol1);
            end;
        end;
        fb:=F(b);
    end;
    HasRoot:=False;
    x:=b;
end;

end.

Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU