метод Брента

Информацию об алгоритме обнаружения цикла Брента см. в разделе « Обнаружение цикла § Алгоритм Брента» .

В численном анализе метод Брента представляет собой гибридный алгоритм поиска корня, сочетающий в себе метод деления пополам , метод секущего и обратную квадратичную интерполяцию . Он надежен, как разделение пополам, но может быть таким же быстрым, как и некоторые менее надежные методы. Алгоритм пытается использовать потенциально быстро сходящийся метод секущего или обратную квадратичную интерполяцию, если это возможно, но при необходимости он возвращается к более надежному методу деления пополам. Метод Брента принадлежит Ричарду Бренту. и основан на более раннем алгоритме Теодора Деккера . Следовательно, этот метод также известен как метод Брента – Деккера .

Современные усовершенствования метода Брента включают метод Чандрупатлы, который проще и быстрее для функций, плоских относительно своих корней; Метод Риддерса , который выполняет экспоненциальную интерполяцию вместо квадратичной, обеспечивая более простую замкнутую формулу для итераций; и метод ITP , который представляет собой гибрид правила-ложного и деления пополам, который обеспечивает оптимальные гарантии наихудшего случая и асимптотические гарантии.

метод Деккера

Идея объединить метод деления пополам с методом секущих восходит к Деккеру (1969) .

Предположим, что мы хотим решить уравнение f ( x ) = 0. Как и в случае с методом деления пополам, нам нужно инициализировать метод Деккера двумя точками, скажем a ₀ и b ₀ , такими, что f ( a ₀ ) и f ( b ₀ ) имеют противоположные знаки. Если f непрерывно на [ a ₀ , b ₀ ], теорема о промежуточном значении гарантирует существование решения между a ₀ и b ₀ .

В каждой итерации участвуют три точки:

• b _k — текущая итерация, т. е. текущее предположение для корня f .
• a _k — это «контраточка», то есть точка такая, что ( ak ) _и f ( b k ₎ имеют противоположные знаки, поэтому интервал [ ak _f , b _k ] содержит решение. Кроме того, | ж ( б _k )| должно быть меньше или равно | f ( a _k )|, так что b _k является лучшим предположением для неизвестного решения, чем a _k .
• b _{k −1} — предыдущая итерация (для первой итерации мы устанавливаем b _{k −1} = a ₀ ).

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

$s={\begin{cases}b_{k}-{\frac {b_{k}-b_{k-1}}{f(b_{k})-f(b_{k-1})}}f(b_{k}),&{\mbox{if }}f(b_{k})\neq f(b_{k-1})\\m&{\mbox{otherwise }}\end{cases}}$

а второй дается методом деления пополам

$m={\frac {a_{k}+b_{k}}{2}}.$

Если результат метода секущего s лежит строго между b _k и m , то он становится следующей итерацией ( b _{k +1} = s ), в противном случае используется средняя точка ( b _{k +1} = m ).

Затем значение новой контраточки выбирается такое, чтобы f ( a _{k +1} ) и f ( b _{k +1} ) имели противоположные знаки. Если f ( a _k ) и f ( b _{k +1} ) имеют противоположные знаки, то контрапункт остается прежним: a _{k +1} = a _k . В противном случае f ( b _{k +1} ) и f ( b _k ) имеют противоположные знаки, поэтому новая контрапункт становится a _{k +1} = b _k .

Наконец, если | ж ( а _{k +1} )| < | f ( b _{k +1} )|, то a _{k +1,} вероятно, является лучшим предположением для решения, чем b _{k +1} , и, следовательно, значения a _{k +1} и b _{k +1} меняются местами.

На этом заканчивается описание одной итерации метода Деккера.

Метод Деккера работает хорошо, если функция f ведет себя достаточно хорошо. Однако бывают случаи, когда на каждой итерации используется метод секущего, но итерации b _k сходятся очень медленно (в частности, | b _k − b _{k −1} | может быть сколь угодно малым). В этом случае метод Деккера требует гораздо больше итераций, чем метод деления пополам.

метод Брента

Брент (1973) предложил небольшую модификацию, чтобы избежать проблем с методом Деккера. Он вставляет дополнительный тест, который должен быть выполнен, прежде чем результат метода секущего будет принят в качестве следующей итерации. Одновременно должны выполняться два неравенства:

Учитывая конкретный числовой допуск $\delta$, если на предыдущем шаге использовался метод деления пополам, неравенство ${\textstyle |\delta |<|b_{k}-b_{k-1}|}$ должен удерживаться для выполнения интерполяции, в противном случае выполняется метод деления пополам, и его результат используется для следующей итерации.

Если на предыдущем шаге выполнялась интерполяция, то неравенство ${\textstyle |\delta |<|b_{k-1}-b_{k-2}|}$ вместо этого используется для выполнения следующего действия (выбора) интерполяции (когда неравенство истинно) или метода деления пополам (когда неравенство неверно).

Также, если на предыдущем шаге использовался метод деления пополам, неравенство ${\textstyle |s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k}-b_{k-1}|}$ должно выполняться, в противном случае выполняется метод деления пополам, и его результат используется для следующей итерации. Если на предыдущем шаге выполнялась интерполяция, то неравенство ${\textstyle |s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k-1}-b_{k-2}|}$ вместо этого используется.

Эта модификация гарантирует, что на k-й итерации шаг деления пополам будет выполнен не более чем за $2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )$ дополнительных итераций, поскольку вышеуказанные условия вынуждают размеры последовательных шагов интерполяции уменьшаться вдвое каждые две итерации, и после не более $2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )$ итераций, размер шага будет меньше, чем $\delta$, который вызывает шаг деления пополам. Брент доказал, что его метод требует не более N ² итераций, где N обозначает количество итераций метода деления пополам. Если функция f ведет себя хорошо, то метод Брента обычно будет использовать либо обратную квадратичную, либо линейную интерполяцию, и в этом случае он будет сходиться сверхлинейно .

Кроме того, метод Брента использует обратную квадратичную интерполяцию вместо линейной интерполяции (как используется в методе секущего). Если f ( b _k ), f ( a _k ) и f ( b _{k −1} ) различны, это немного увеличивает эффективность. Как следствие, условие принятия s (значение, предложенное либо линейной интерполяцией, либо обратной квадратичной интерполяцией) должно быть изменено: s должно лежать между (3 a _k + b _k ) / 4 и b _k .

Алгоритм

input a, b, and (a pointer to) a function for f
calculate f(a)
calculate f(b)
if f(a)f(b) ≥ 0 then 
    exit function because the root is not bracketed.
end if
if |f(a)| < |f(b)| then
    swap (a,b)
end if
c := a
set mflag
repeat until f(b or s) = 0 or |b − a| is small enough (convergence)
    if f(a) ≠ f(c) and f(b) ≠ f(c) then
        ${\textstyle s:={\frac {af(b)f(c)}{(f(a)-f(b))(f(a)-f(c))}}+{\frac {bf(a)f(c)}{(f(b)-f(a))(f(b)-f(c))}}+{\frac {cf(a)f(b)}{(f(c)-f(a))(f(c)-f(b))}}}$ (inverse quadratic interpolation)
    else
        ${\textstyle s:=b-f(b){\frac {b-a}{f(b)-f(a)}}}$ (secant method)
    end if
    if (condition 1) s is not between $(3a+b)/4$ and b or
       (condition 2) (mflag is set and |s−b| ≥ |b−c|/2) or
       (condition 3) (mflag is cleared and |s−b| ≥ |c−d|/2) or
       (condition 4) (mflag is set and |b−c| < |δ|) or
       (condition 5) (mflag is cleared and |c−d| < |δ|) then
        ${\textstyle s:={\frac {a+b}{2}}}$ (bisection method)
        set mflag
    else
        clear mflag
    end if
    calculate f(s)
    d := c  (d is assigned for the first time here; it won't be used above on the first iteration because mflag is set)
    c := b
    if f(a)f(s) < 0 then
        b := s 
    else
        a := s 
    end if
    if |f(a)| < |f(b)| then
        swap (a,b) 
    end if
end repeat
output b or s (return the root)

Пример

Предположим, что мы ищем ноль функции, определяемой формулой f ( x ) = ( x + 3) ( x − 1) ² .

Мы берем [ a ₀ , b ₀ ] = [−4, 4/3] в качестве начального интервала.

Имеем f ( a ₀ ) = −25 и f ( b ₀ ) = 0,48148 (все числа в этом разделе округлены), поэтому условия f ( a ₀ ) f ( b ₀ ) < 0 и | ж ( б ₀ )| ≤ | ж ( а ₀ )| удовлетворены.

График f ( x ) = ( x + 3)( x - 1) ²

1. На первой итерации мы используем линейную интерполяцию между ( b ₋₁ , f ( b ₋₁ )) = ( a ₀ , f ( a ₀ )) = (−4, −25) и ( b ₀ , f ( b _{0 )} )) = (1,33333, 0,48148), что дает s = 1,23256. Оно лежит между (3 a ₀ + b ₀ ) / 4 и b ₀ , поэтому это значение принимается. Кроме того, f (1,23256) = 0,22891, поэтому мы положили a ₁ = a ₀ и b ₁ = s = 1,23256.
2. Во второй итерации мы используем обратную квадратичную интерполяцию между ( a ₁ , f ( a ₁ )) = (−4, −25) и ( b ₀ , f ( b ₀ )) = (1,33333, 0,48148) и ( b ₁ , ж ( б ₁ )) = (1,23256, 0,22891). Это дает 1,14205, которое находится между (3 a ₁ + b ₁ )/4 и b ₁ . При этом неравенство |1.14205 − b ₁ | ≤ | б ₀ - б _-1 | / 2 удовлетворено, поэтому это значение принимается. Кроме того, f (1,14205) = 0,083582, поэтому мы полагаем a ₂ = a ₁ и b ₂ = 1,14205.
3. На третьей итерации мы используем обратную квадратичную интерполяцию между ( a ₂ , f ( a ₂ )) = (−4, −25) и ( b ₁ , f ( b ₁ )) = (1,23256, 0,22891) и ( b ₂ , f ( b ₂ )) = (1,14205, 0,083582). Это дает 1,09032, которое находится между (3 a ₂ + b ₂ )/4 и b ₂ . Но здесь вступает в силу дополнительное условие Брента: неравенство |1,09032 − b ₂ | ≤ | б ₁ - б ₀ | / 2 не удовлетворено, поэтому это значение отклоняется. средняя точка m = -1,42897 интервала [ a ₂ , b ₂ Вместо этого вычисляется ]. У нас f ( m ) = 9,26891, поэтому мы положили a ₃ = a ₂ и b ₃ = −1,42897.
4. На четвертой итерации мы используем обратную квадратичную интерполяцию между ( a ₃ , f ( a ₃ )) = (−4, −25) и ( b ₂ , f ( b ₂ )) = (1,14205, 0,083582) и ( b ₃ , f ( b ₃ )) = (−1,42897, 9,26891). Это дает 1,15448, что не находится в интервале между (3 a ₃ + b ₃ )/4 и b ₃ ). Следовательно, она заменяется серединой m = −2,71449. У нас f ( m ) = 3,93934, поэтому мы положили a ₄ = a ₃ и b ₄ = −2,71449.
5. На пятой итерации обратная квадратичная интерполяция дает −3,45500, что лежит в требуемом интервале. Однако предыдущая итерация была шагом пополам, поэтому неравенство |−3,45500 − b ₄ | ≤ | б ₄ - б ₃ | / 2 должны быть удовлетворены. Это неравенство неверно, поэтому мы используем среднюю точку m = −3,35724. У нас есть f ( m ) = −6,78239, поэтому m становится новой контраточкой ( a ₅ = −3,35724 ), а итерация остается той же ( b ₅ = b ₄ ).
6. На шестой итерации мы не можем использовать обратную квадратичную интерполяцию, поскольку b ₅ = b ₄ . Следовательно, мы используем линейную интерполяцию между ( a ₅ , f ( a ₅ )) = (-3,35724, -6,78239) и ( b ₅ , f ( b ₅ )) = (-2,71449, 3,93934). Результат: s = −2,95064, что удовлетворяет всем условиям. Но поскольку итерация не изменилась на предыдущем шаге, мы отвергаем этот результат и возвращаемся к делению пополам. Мы обновляем s = -3,03587 и f ( s ) = -0,58418.
7. На седьмой итерации мы снова можем использовать обратную квадратичную интерполяцию. Результат: s = −3,00219, что удовлетворяет всем условиям. Теперь f ( s ) = −0,03515, поэтому мы устанавливаем a ₇ = b ₆ и b ₇ = −3,00219 ( a ₇ и b ₇ меняются местами так, что условие | f ( b ₇ )| ≤ | f ( a ₇ ) | удовлетворен). ( Правильно : линейная интерполяция ⁠ $s=-2.99436,f(s)=0.089961$⁠ )
8. На восьмой итерации мы не можем использовать обратную квадратичную интерполяцию, поскольку a ₇ = b ₆ . Линейная интерполяция дает s = -2,99994, что принимается. ( Правильно : ⁠ $s=-2.9999,f(s)=0.0016$⁠ )
9. В следующих итерациях корень x = −3 достигается быстро: b ₉ = −3 + 6·10 ⁻⁸ и b ₁₀ = −3 − 3·10 ⁻¹⁵ . ( Правильно : Итер 9: f ( s ) = −1,4 × 10 ⁻⁷ , Итер 10: ж ( с ) = 6,96 × 10 ⁻¹² )

Реализации

• Брент (1973) опубликовал реализацию Алгола 60 .
• Netlib содержит перевод этой реализации на Фортран с небольшими изменениями.
• Метод PARI/ GP solve реализует метод.
• Другие реализации алгоритма (на C++, C и Fortran) можно найти в книгах числовых рецептов .
• Библиотека Apache Commons Math реализует алгоритм на языке Java .
• Модуль оптимизации SciPy реализует алгоритм на Python (языке программирования).
• Стандартная библиотека Modelica реализует алгоритм в Modelica .
• The uniroot функция реализует алгоритм в R (программное обеспечение) .
• The fzero функция реализует алгоритм в MATLAB .
• Boost (библиотеки C++) реализует два алгоритма, основанных на методе Брента на C++ в наборе инструментов Math:
  1. Минимизация функции на сайте minima.hpp с примером поиска функции минимума .
  2. Поиск корня реализует новый TOMS748, более современный и эффективный алгоритм, чем исходный алгоритм Брента, на TOMS748 , а также поиск корня Boost.Math , который использует TOMS748 внутри себя с примерами .
• Пакет Optim.jl реализует алгоритм на Julia (языке программирования).
• Система компьютерной алгебры Emmy (написанная на Clojure (язык программирования) ) реализует вариант алгоритма, предназначенный для минимизации одномерных функций.
• Поиск корня в библиотеке C#, размещенной в проекте кода.

Ссылки

• Брент, Р.П. (1973), «Глава 4: Алгоритм с гарантированной сходимостью для поиска нуля функции», Алгоритмы минимизации без производных , Энглвуд Клиффс, Нью-Джерси: Прентис-Холл, ISBN  0-13-022335-2
• Деккер, Т.Дж. (1969), «Нахождение нуля посредством последовательной линейной интерполяции», в Дежоне, Б.; Хенричи П. (ред.), Конструктивные аспекты фундаментальной теоремы алгебры , Лондон: Wiley-Interscience, ISBN  978-0-471-20300-1

Дальнейшее чтение

• Аткинсон, Кендалл Э. (1989). «Раздел 2.8.». Введение в численный анализ (2-е изд.). Джон Уайли и сыновья. ISBN  0-471-50023-2 .
• Пресс, WH; Теукольский, С.А.; Феттерлинг, WT; Фланнери, BP (2007). «Раздел 9.3. Метод Ван Вейнгаардена–Деккера–Брента» . Численные рецепты: искусство научных вычислений (3-е изд.). Нью-Йорк: Издательство Кембриджского университета. ISBN  978-0-521-88068-8 . Архивировано из оригинала 11 августа 2011 г. Проверено 28 февраля 2012 г.
• Алефельд, GE; Потра, ФА; Ши, Исюнь (сентябрь 1995 г.). «Алгоритм 748: Заключение нулей непрерывных функций» . Транзакции ACM в математическом программном обеспечении . 21 (3): 327–344. дои : 10.1145/210089.210111 . S2CID   207192624 .

Внешние ссылки

• Zeroin.f в Netlib .
• модуль brent на C++ (также C, Fortran, Matlab). Архивировано 5 апреля 2018 г. в Wayback Machine Джоном Буркардтом.
• Реализация GSL .
• Увеличьте реализацию C++.
• Python (Scipy) Реализация