Поиск приближения нуля функции на отрезке модифицированным методом Ньютона

newton_search_root.pas

program NewtonSearchRoot;

type
    FunctionType = function(x: double): double;

const
    maxSteps = 100;         { Максимально допустимое количество шагов }
    eps = 0.000000000001;   { Требуемая точность вычислений }
    printInterRes = false;  { Печатать ли все промежуточные приближения } 

{ ********* Секция пользовательских изменений ********* }

{ Исследуемая функция }
function f(x: double): double;
begin
    f := sqr(x) * x - 2 * sqr(x) - x + 1;
end;

{ Первая ее производная }
function df(x: double): double;
begin
    df := 3 * sqr(x) - 4 * x - 1;
end;

{ Вторая производная }
function d2f(x: double): double;
begin
    d2f := 6 * x - 4;
end;

{ ----------------------------------------------------- }

{ Исследуемая функция }
function f2(x: double): double;
begin
    f2 := x * exp(-x) + 2;
end;

{ Первая ее производная }
function df2(x: double): double;
begin
    df2 := exp(-x) * (1 - x);
end;

{ Вторая производная }
function d2f2(x: double): double;
begin
    d2f2 := exp(-x) * (x - 2);
end;

{ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ }


{ f - исследуемая функция,
  a, b, c - коэффициенты уравнения параболической касательной в точке prev,
  left, right - границы отрезка,
  n, m - счетчики числа шагов параболического и линейного типа соответственно }
function getNextApproximation(f: FunctionType; a, b, c, prev, left, right: double; var n, m: integer): double;
var
    d, x1, x2: double;
begin
    if a = 0 then begin                 { При нулевой второй производной исходной функции параболическая касательная вырождается }
        getNextApproximation := -c / b; { в линейную - возьмем в качестве следующего приближения ноль этой прямой                }
        inc(m);
        exit;
    end;

    d := sqr(b) - 4*a*c;
    if d < 0 then begin                 { Если параболическая касательная не пересекает ось Оx, возьмем в качестве приближения   }
        getNextApproximation := prev - f(prev) / df(prev);      { ноль линейной касательной в этой же точке }
        inc(m);
        exit;
    end;

    x1 := 0.5 * (-b - sqrt(d)) / a;
    x2 := 0.5 * (-b + sqrt(d)) / a;
                                        { Из двух нулей параболической касательной выберем ближайший к истинному нулю            }
    if (x1 >= left) and (x1 <= right) then begin            { и принадлежащий отрезку поиска            }
        if (x2 >= left) and (x2 <= right) then begin
            if abs(f(x1)) < abs(f(x2)) then
                getNextApproximation := x1
            else
                getNextApproximation := x2;
        end
        else
            getNextApproximation := x1;
    end
    else
        getNextApproximation := x2;
    inc(n);
end;


procedure newtonSearchRoot(f, df, d2f: FunctionType; left, right, eps: double);
var
    prev, next, a, b, c: double;
    n, m: integer;
begin
    if f(left) * d2f(left) >= 0 then  { При ином выборе начально приближения можно вылететь за отрезок }
        prev := left
    else
        prev := right;
    
    n := 0;
    m := 0;
    while true do begin
        if printInterRes then
            write(prev:16:16, ' ');

        a := 0.5 * d2f(prev);         { Расчет коэффициентов параболической касательной }
        b := df(prev) - 2 * a * prev;
        c := f(prev) - a * sqr(prev) - b * prev;

        next := getNextApproximation(f, a, b, c, prev, left, right, n, m);

        if abs(next - prev) < eps then begin                        { Достигнута необходимая точность }
            writeln;
            writeln('Segment [', left:5:2, ', ', right:5:2, ' ]');
            writeln('    has root: ', next:16:16);
            writeln('(Calculated in ', n, ' parabolic and ', m, ' linear steps with esp = ', eps:16:16, ')');
            exit;
        end;

        if (n + m) >= maxSteps then begin                           { Превышен лимит вмремени работы }
            writeln('The required accuracy is not reached in ', maxSteps, ' steps');
            exit;
        end;

        prev := next;
    end;
end;


begin
    writeln('[ Function f  ]');
    NewtonSearchRoot(@f, @df, @d2f, 1.6, 2.6, eps);
    writeln;

    writeln('[ Function f2 ]');
    NewtonSearchRoot(@f2, @df2, @d2f2, -2.0, 0.0, eps);
    readln;
end.

Успешно компилируется Free Pascal Compiler version 2.4.4