Решение уравнения F(x) = 0 "безопасным" методом Ньютона (гарантированная сходимость)

Этот метод является аналогом метода Брента для случая, когда известна производная функции F. Сочетая в себе гарантированную сходимость дихотомии с высокой скоростью сходимости обычного метода Ньютона, он требует не только знания производной функции, но и "окружения" корня, т.е. знания отрезка, на концах которого функция имеет разные знаки.

Сначала метод делает шаг по правилу Ньютона. Если сделанный шаг не укладывается в указанный отрезок или же сходимость слишком низкая, то делается один шаг метода дихотомии. Затем процесс повторяется.

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

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

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

unit SafeNewtonUnit;

declaration
    procedure F(X:Real; out F: Real; out dF: Real);
    
interface
    SearchRootSafeNewton;

implementation

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

(*************************************************************************
Поиск корня комбинированным методом Ньютона (подобен   методу  Брента,  но
использует информацию о производной функции).

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

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

HasRoot содержит   True, если останов  произошел  из-за  нахождения корня,
(в      переменной Root содержится корень) и False,  если число   итераций
слишком велико или функция не принимает противоположные знаки на концах.
*************************************************************************)
procedure SearchRootSafeNewton(
        x1      :   Real;
        x2      :   Real;
        epsilon :   Real;
    out HasRoot :   Boolean;
    out Root    :   Real);
var
   df   :   Real;
   dx   :   Real;
   dxold:   Real;
   f0   :   Real;
   fh   :   Real;
   fl   :   Real;
   swap :   Real;
   temp :   Real;
   xh   :   Real;
   xl   :   Real;
   rts  :   Real;
   j    :   Integer;
begin
    HasRoot:=True;
    F(x1,fl,df);
    F(x2,fh,df);
    if fl*fh>=0 then
    begin
        HasRoot:=False;
        Exit;
    end;
    if fl<0 then
    begin
        xl:=x1;
        xh:=x2
    end
    else
    begin
        xh:=x1;
        xl:=x2;
        swap:=fl;
        fl:=fh;
        fh:=swap
    end;
    rts:=0.5*(x1+x2);
    dxold:=AbsReal(x2-x1);
    dx:=dxold;
    F(rts,f0,df);
    for j := 1 to ItMax do
    begin
        if (((rts-xh)*df-f0)*((rts-xl)*df-f0)>=0)
            or (AbsReal(2*f0)>AbsReal(dxold*df)) then
        begin
            dxold:=dx;
            dx:=0.5*(xh-xl);
            rts:=xl+dx;
            if xl=rts then
            begin
                Root:=rts;
                HasRoot:=True;
                Exit;
            end;
        end
        else
        begin
            dxold := dx;
            dx := f0/df;
            temp := rts;
            rts := rts-dx;
            if temp=rts then
            begin
                Root:=rts;
                HasRoot:=True;
                Exit;
            end;
        end;
        if AbsReal(dx)<Epsilon then
        begin
            Root:=rts;
            HasRoot:=True;
            Exit;
        end;
        F(rts,f0,df);
        if f0<0 then
        begin
            xl:=rts;
            fl:=f0;
        end
        else
        begin
            xh:=rts;
            fh:=f0;
        end;
    end;
    Root:=rts;
    HasRoot:=False;
end;

end.

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