Метод половинного деления

[GAA]

На мой взгляд, если пытаться выжать всё возможное из соответствующего типа чисел с плавающей точкой, то проще выйти на зацикливание, а затем уже выполнять дополнительные проверки и возвращать соответствующий результат:

если среднее не совпадает с концами отрезка, то возвращать среднее.
если совпадает, то возвращать конец отрезка с меньшим абсолютным значением функции.

Пусть функция $f(x) = x^{33}$.

Используется синтаксис Delphi

program Dichotomy;

type Float = Double;

function f(x: Float): Float;
begin
 f:= sqr(sqr(sqr(sqr(sqr(x)))))*x;
end;

var
 a, b, x_c, f_c, f_prev : Float;
 count, maxcount: integer;

begin
 a:= -0.5; b := 1; count:= 0;
 f_c := f(a);
 repeat
      f_prev:= f_c;
      count:= count+1;
      x_c:= (a+b)/2;
      f_c:= f(x_c);
      if f_c < 0
       then a := x_c
       else if f_c > 0
             then b:= x_c;
 until f_prev = f_c;
 if (x_c <> a) and (x_c <> b)
  then x_c := (a+b)/2
  else if abs(f(a)) < abs(f(b))
        then x_c := a
        else x_c := b;
 writeln('left= ', a,  ' right= ', b, ' x_c =', x_c);
 writeln('difference= ', b-a, ' count=', count);
 writeln('f(a)= ', f(a), ' f(b)= ', f(b), ' f(x_c)= ', f(x_c));
end.

Выдача

left= -4.65661287307739E-0010 right=  2.32830643653870E-0010 x_c =-1.16415321826935E-0010
difference=  6.98491930961609E-0010 count=33
f(a)= -1.11253692925360E-0308 f(b)=  1.29516344663408E-0318 f(x_c)= -0.00000000000000E+0000 

Не очень текст обдумывал. Просто для примера. Но вообще, когда в давние времена смотрел чужие программы, то там проверки на значения функции присутствовали.
И конечно, в реальности погрешности вычисления функции могут быть устроены очень сложно. Скорее всего, нужно будет допиливать под конкретную функцию.

Upd else break в тексте потерял, но сути это не меняет.

[vpb]

Dmitriy40, GAA
Спасибо за ответы. И вообще спасибо всем принявшим участие в теме. (Я сам не специалист по вычислительным вещам, в некоторый момент понял, что не вполне понимаю про этот метод, решил узнать, что люди, больше меня занимающиеся вычислениями, думают. В общем, представление себе составил.)

[B@R5uk]

vpb в сообщении #1466443 писал(а):

Мой весьма скромный опыт говорит, что иногда и квадратное уравнение лучше этим методом решать. Так как из школьной формулы могут возникнуть ошибки.

Оооо! Это вы зря, очень зря! Численное нахождение корня (или корней) вблизи экстремума функции (минимум, максимум) или даже вблизи просто критической точки — это настоящий ад! Если есть честная формула для обратной функции, то всегда лучше воспользоваться ей. Разумеется, если в формуле есть разности, то её наверняка придётся доработать напильником. Но в любом случае, формула предоставляет возможность доработки, при численном счёте, как правило, такой роскоши нет.

Проблема с численным обращением функции вблизи нуля производной заключается в погрешностях округления чисел с плавающей запятой. Вот как, например, выглядит функция $f(x)=x^3-5.88x+5.488$ вблизи "корня" $1.4$:



Пример утрированный, но даёт понять проблему.

[Утундрий]

Чисметы для классического математика вообще довольно мутная тема, уже в силу того, что кажущиеся ему эквивалентными преобразования внезапно перестают быть таковыми.

[EzikBro]

Использовать такую программу, конечно, можно, но гораздо лучше использовать проверенный временем вещественный бинарный поиск с условием "найти минимальный x такой, что $f(x) \geqslant 0$".
Тогда программа будет выглядеть примерно так:

Используется синтаксис C++

double xl = a, xr = b;
for (int i = 0; i < 100; i++)
{
double m = (xl + xr) / 2;
if (f(m) >= 0)
xr = m;
else
xl = m;
}
 

При этом количество итераций выбирается обычно наобум, но считается, что 100 хватает в любом случае.

[venco]

Ста итераций может сильно не хватить, если начальные границы сильно отличаются по порядку от корня. Например, для универсального поиска одного корня начиная с $a=-DBL\_MAX$, $b=DBL\_MAX$, и если корень окажется чуть больше нуля, то после ста итераций диапазон уменьшится всего лишь до $a=0$, $b=2.8\cdot10^{278}$.

[Dmitriy40]

О, а об этом я не подумал, проверял лишь на близких начальных значениях.
Много итераций понадобится и в случае решения $|x|<10^{-16}$ из-за возможности работы с денормализованными числами (с потерей относительной точности).

GAA:

Обсуждение IEEE 754 выделено в самостоятельную ветку «IEEE 754»