Алгоритм поиска интервала нахождения корня нелинейной функции

Выяснение интервала, на котором корень содержится является важной проблемой поиска корня нелинейной функции действительной переменной. Здесь приведен алгоритм поиска такого интервала и ограничения на его применение. Примем, что корень функции f(x) окружен на интервале [a,b], если f(a) и f(b) имеют противоположные знаки. Чтобы окруженный согласно этому определению корень действительно существовал на этом интервале, достаточно непрерывности f(x), а для его единственности - еще и монотонности. При невыполнении этих свойств возможно отсутствие корня на [a,b] или неопределенность его позиции.

При использовании компьютера, всегда имеем дело с дискретным набором возможных представлений чисел (хотя и достаточно плотным). Кроме того, монотонность вычисленной функции может быть слегка нарушена в пределах точности ее вычисления. Это в ряде случаев усложняет вычисление окруженных корней функции, если к их точности предъявляются завышенные требования.

Окружение корня функции при гарантии ее определения на неограниченном интервале, производится по следующему итерационному алгоритму.
1. Для начального приближения x₀, найти f₀=f(x₀), задать начальный интервал поиска D и его инкремент d>1.
2. Вычислить a=x₀-D, b=x₀+D; f_a=f(a), f_b=f(b).
3. Увеличить интервал поиска: D=D*d. Если интервал превысил некоторый заданный предел, то выйти с индикацией ошибки.
4a. Если знаки f_a и f₀ отличаются, то считать корень окруженным на [a,x₀] -> выход.
4b. Если знаки f_b и f₀ отличаются, то считать корень окруженным на [x₀,b] -> выход.
4c. Если f₀>0 (случай меньше нуля делается аналогично) алгоритм продолжается:
5. Проверяется, какое из f_a или f_b наименьшее. Если оба одинаковы, то переходим к 6a (двусторонний поиск), если f_b -- производим поиск вправо 6b, иначе -- поиск влево 6c.
6a. Находим a=a-D, b=b+D, f_a=f(a), f_b=f(b), идем к пункту 3.
6b. Присваиваем последовательно a=x₀, x₀=b, f_a=f₀, f₀=f_b; находим b=b+D, f_b=f(b), идем к пункту 3.
6c. Аналогично 6b, только направление поиска -- влево.

Так как интервал поиска постоянно расширяется, то в конце концов используя указанный алгоритм корень будет окружен. Возможны модификации алгоритма в двух направлениях:
1) Увеличивать интервал не в геометрической прогрессии, а в арифметической либо по заданному сценарию;
2) Если область определения функции заведомо ограничена, то расширение интервала поиска также следует ограничивать имеющимися пределами, либо доопределять функцию там, где ее оригинал не определен.

 

Программа окружения корня нелинейной функции

Ниже расположена программа окружения корня нелинейной функции, которая реализует данный алгоритм.

/* Окружение корня функции. 
Предполагается, что функция существует на бесконечном и непрерывном промежутке.

   int BracketRoot(double x0,double *a,double *b,double d0,
                    double di, double dmax,double (*fun)(double));
   Параметры:
    x0 - начальное приближение;
    a  - левая граница;
    b  - правая граница;
    d0 - начальный интервал поиска;
    di - интервал расширения (в геометрической прогрессии);
    dmax - максимальный интервал;
    fun - указатель на функцию.
   Возвращает: 
    1 - если корень окружен;
    0 - сбой
*/

int BracketRoot(double x0,double *a,double *b,double d0,
                 double di, double dmax,double (*fun)(double)) {
  double fa,fb,f0;
  f0=(*fun)(x0); *a=x0-d0; *b=x0+d0; fa=(*fun)(*a); fb=(*fun)(*b);
  while((d0*=di)<dmax) {
    if(f0>=0.) {
      if(fa<0.) {*b=x0;return(1);}
      if(fb<0.) {*a=x0;return(1);}
      if(fa>fb) {*a=x0; x0=(*b); *b+=d0; fa=f0; f0=fb; fb=(*fun)(*b);}
      else if(fa<fb) {*b=x0; x0=(*a); *a-=d0; fb=f0; f0=fa; fa=(*fun)(*a);}
      else {*a-=d0; *b+=d0; fa=(*fun(*a);fb=(*fun)(*b);}
    }
    else if(f0<0.) {
      if(fa>=0.) {*b=x0;return(1);}
      else if(fb>=0.) {*a=x0;return(1);}
      if(fa<fb) {*a=x0; x0=(*b); *b+=d0; fa=f0; f0=fb; fb=(*fun)(*b);}
      else if(fa>fb) {*b=x0; x0=(*a); *a-=d0; fb=f0; f0=fa; fa=(*fun)(*a);}
      else {*a-=d0; *b+=d0; fa=(*fun(*a);fb=(*fun)(*b);}
    }
  }
  return(0);
}