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

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

Пусть в точках a, b функция имеет противоположные знаки. Пусть x₁  = a, x₂  = b, x₃  = (a+b)/2, f₁  = f(x₁ ), f₂  = f(x₂ ), f₃  = f(x₃ ).

Тогда x₄  = x₃ +(x₃ -x₁ )*sign(f₁ -f₂ )*f₃ /sqrt(f₃ *f₃ -f₁ *f₂ ), f₄  = f(x₄ ). Перейдем к следующему шагу метода. Если f₃ *f₄  < 0, то полагаем a = x₃ , b = x₄ . В противном случае за [a, b] берем тот из отрезков [x₁ , x₄ ], [x₄ , x₂ ], на котором функция меняет знак.

Алгоритм останавливает работу, если достигнута заданная точность или выполнено указанное число итераций.

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

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

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

unit SearchRootRiddlerUnit;

declaration
    function F(X:Real):Real;

interface
    SearchRootRiddler;
implementation

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

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

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

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

HasRoot содержит   True, если останов  произошел  из-за  нахождения корня,
(в переменной x содержится корень) и False,  если число   итераций слишком
велико или функция не принимает противоположные знаки на концах.
*************************************************************************)
procedure SearchRootRiddler(
            x1      :   Real;
            x2      :   Real;
    const   Epsilon :   Real;
    out     HasRoot :   Boolean;
    out     x       :   Real);
var
    I   :   Integer;
    x3  :   Real;
    x4  :   Real;
    F1  :   Real;
    F2  :   Real;
    F3  :   Real;
    F4  :   Real;
begin
    HasRoot := False;
    F1:=F(x1);
    F2:=F(x2);
    if F1*F2>=0 then
    begin
        HasRoot:=False;
        Exit;
    end;
    
    for I:=1 to ItMax do
    begin
        //проверка значений функции
        if AbsReal(x2-x1)<Epsilon then
        begin
            if AbsReal(F1)>AbsReal(F2) then
                x := x1
            else
                x := x2;
            HasRoot:=True;
            Exit;
        end;

        //шаг метода
        x3:=0.5*(x1+x2);
        F3:=F(x3);
        x4:=x3 + (x3-x1)*Sign(F1-F2)*F3/Sqrt(F3*F3-F1*F2);
        F4:=F(x4);
        
        //выбор новой пары точек
        if F3*F4<=0 then
        begin
            X1:=X3;
            F1:=F3;
            X2:=X4;
            F2:=F4;
        end
        else
        begin
            if F1*F4<=0 then
            begin
                X2:=X4;
                F2:=F4;
            end
            else
            begin
                X1:=X4;
                F1:=F4;
            end;
        end;
    end;
end;


end.

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