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

[vpb]

(Не совсем понятно, куда написать, так что напишу сюда.)

Уважаемые коллеги, все мы знаем метод половинного деления.

Предположим, что $a<b$, $f(a)<0$, $f(b)>0$. Дальше на Сях напишу:

x1=a, x2=b;
while (x2-x1>=delta) 
  { 
     x_mid=(x1+x2)/2;
     if (f(x_mid)<0) 
        x1=x_mid;
     else 
        x2=x_mid;
   }

Вопрос такой: а приходилось ли лично вам его использовать, и из каких соображений выбирается delta ?

[Someone]

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

из каких соображений выбирается delta ?

А с какой погрешностью Вам нужно получить результат? Вот её и возьмите в качестве $\delta$.

[Dmitriy40]

Приходилось.
Из двух соображений, когда одно из, когда оба вместе: 1) из физических соображений (например для моделирования падения в обычных условиях микросекундного шага достаточно с большим запасом), 2) когда при уменьшении дельты скажем вдвое перестаёт сильно меняться результат (т.е. расти точность). Иногда добавляю условие на время счёта, чтобы не ждать сутками промежуточных/отладочных результатов, а прогнать и отладить всё за минуты/часы и лишь потом запускать с мелким шагом (малой дельта) на просчёт финального результата.
Универсального математического метода наверняка нет, погрешности всегда же разные допустимы.

[Утундрий]

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

[vpb]

Спасибо за ответы.
Dmitriy40
То есть Вы назначали $\delta$ исходя из нужд конкретной прикладной задачи, по обстоятельствам, что называется ad hoc. Примерно то же, что и коллега Someone советует.

Утундрий в сообщении #1466436 писал(а):

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

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

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

[Dmitriy40]

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

Dmitriy40
То есть Вы назначали $\delta$ исходя из нужд конкретной прикладной задачи, по обстоятельствам, что называется ad hoc. Примерно то же, что и коллега Someone советует.

Да. Только не всегда оценивал заранее, часто уже по итогам прогона, вижу что долго считает — увеличу дельту (уменьшу точность), вижу что быстро — увеличу точность, вижу что при уменьшении дельты точность итогового результата перестала расти — дальше нет смысла уменьшать дельту. Как-то так, эмпирически. Далеко не всегда можно (и стоит этим заниматься) оценить влияние всех дельт в программе на итоговый результат.
А иногда наоборот, из физических соображений очевидно какую дельту можно выбрать "чтобы с гарантией" (например координаты автомобиля знать с субмикронной точностью — излишество).

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

Так как из школьной формулы могут возникнуть ошибки.

При математически корректных вычислениях — не могут. А если вычисления не абсолютно точные, то это уже совсем другой вопрос, о точности компьютерных вычислений и там о-о-очень много других подводных камней кроме выбора дельты. И далеко не всегда метод половинного деления будет лучше других.

[vpb]

Dmitriy40 в сообщении #1466450 писал(а):

При математически корректных вычислениях — не могут. А если вычисления не абсолютно точные, то это уже совсем другой вопрос, о точности компьютерных вычислений и там о-о-очень много других подводных камней кроме выбора дельты.

Несомненно, вычислительная математика --- достаточно сложная наука. Я немного ее изучал (не считая того, что в университете было, но то уже в основном забылось). Вот, кстати, про квадратное уравнение такой известный математик Дж. Форсайт еще в 1966 г. написал статью с интригующим названием "How do you solve a quadratic equation ?". 19 стр., не хухры-мухры.

[ozheredov]

Утундрий в сообщении #1466436 писал(а):

о самой функции ничего толком неизвестно

Она как минимум должна быть монотонна, не? )

[Dan B-Yallay]

ozheredov в сообщении #1466454 писал(а):

Она как минимум должна быть монотонна, не? )

Для половинного деления - необязательно. Если корней несколько - значит найдётся какой-то один из них.

[slavav]

Про выбор дельты есть забавное соображение: каждая итерация цикла добавляет ещё один бит к точности результата. А таких битов в большинстве случаев не больше 54 (стандартный double). После того как все они будут найдены переменные x1 и x2 сравняются. Следовательно можно написать заголовок цикла как

while (x2 > x1)

. С точки зрения математики - некорректно, с точки зрения программирования - всё работает и дельту определять не нужно.

[Sender]

slavav в сообщении #1466520 писал(а):

Про выбор дельты есть забавное соображение: каждая итерация цикла добавляет ещё один бит к точности результата. А таких битов в большинстве случаев не больше 54 (стандартный double).

Нет ли тут опасности, что при неравных x1 и x2 отличие их полусуммы от одного из них уже выйдет за пределы машинной точности, приведя к зацикливанию?

-- Вт июн 02, 2020 12:28:18 --

Конечно, именно в данном случае double в интеловской архитектуре - это ещё не машинная точность, но...

[vpb]

slavav в сообщении #1466520 писал(а):

всё работает и дельту определять не нужно.

И Вы сами так делали ?

[Dmitriy40]

Я только что проверил, да, на 53..59 итерации (в зависимости от исходных значений) double переменные становятся равными. Но есть и исключение: если корень близко к нулю (если не ошибаюсь, примерно $|x|<10^{-308}$), то итераций может быть до 1080. Да, денормализованные числа однако. Они кстати кажется ещё и страшно медленно обрабатываются.

[vpb]

Dmitriy40 в сообщении #1466557 писал(а):

Я только что проверил, да, на 53..59 итерации (в зависимости от исходных значений) double переменные становятся равными.

Гм. Так быть не должно. На всякий случай поставил простейший эксперимент.

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>

double f(double);


double f(double x)
  { return x*x-2; }


main()
 {
    int count, maxcount; 
    double a,b,x1,x2,x3;


    a=1.1, b=1.8;
    x1=a, x2=b, count=0, maxcount=5000;
    while(x1<x2 && count < maxcount) {
      x3=(x1+x2)/2;
      if(f(x3)<0) x1=x3;
      else x2=x3;
      count++;
    } 

    printf("count=%d\n", count);
    printf("left=%f,  right=%f, difference=%e\n", x1,x2, x2-x1);

    return;
 }

Выдача: count=5000, left=1.14214, right=1.14214, difference=2.220446e-016.
То есть концы интервала устанавливаются в два соседних числа с плавающей точкой, и дальше зацикливание. Т.е. как и должно быть.

А можно на Ваш код взглянуть ?

[Dmitriy40]

vpb
Я попеременно писал новое значение то в левый, то в правый предел, без всякой функции вообще, видимо поэтому и не зациклилось. Если сделать нормально, с функцией, то да, зацикливается на соседних представимых числах (из-за округления полусуммы всегда к одному из них). Код банальный:

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

uses    SysUtils;
var     i: integer;
        a, b, m: double;
begin
        i := 0; a := StrToFloat(ParamStr(1)); b := StrToFloat(ParamStr(2));
        while b > a do begin
                Inc(i);
                m := (a + b) / 2;
//              if f(m) < 0 then a := m else b := m;//Правильно
                if (i and 1) > 0 then a := m else b := m;//Ошибочно
                Writeln(i, ':', a:20:20, '..', b:20:20, '=', b-a);
        end;
end.

С другой стороны, можно выбирать дельту как $10^{-15}$ от величины чисел (скажем середины интервала) и точность всегда будет близка к максимально возможной.

А ещё лучше проверить что новая середина стала равна любому из пределов и тогда дальше не считать, это исключает зацикливание в любом случае, оставляя точность вычислений максимально возможной для чисел в любом формате single/double (с той же возможностью уйти в денормализованные числа около нуля):

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

for (;;) {
        m=(a+b)/2;
        if ((a==m) || (m==b)) break;
        if (f(m)<0) a=m; else b=m;
}

Если 53-54 (или 64-65 если вдруг оптимизатор всё оставит в регистрах FPU или до 1075 около нуля) итерации допустимы всегда, то это видимо лучшее решение.