Метод Брента для поиска корней

Введение

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

Пример

Предположим, что мы ищем нуль функции, определенной как f(x) = (x + 3)(x − 1)². Мы берем [a₀, b₀] = [−4, 4/3] в качестве нашего начального интервала. У нас есть f(a₀) = −25 и f(b₀) = 0.48148 (все числа в этом разделе округлены), поэтому условия f(a₀)f(b₀) < 0 и |f(b₀)| ≤ |f(a₀)| выполнены. В первой итерации мы используем линейную интерполяцию между (b₋₁, f(b₋₁)) = (a₀, f(a₀)) = (−4, −25) и (b₀, f(b₀)) = (1.33333, 0.48148), что дает s = 1.23256. Это лежит между (3a₀ + b₀) / 4 и b₀, поэтому это значение принимается. Кроме того, f(1.23256) = 0.22891, поэтому мы устанавливаем a₁ = a₀ и b₁ = s = 1.23256. Во второй итерации мы используем обратную квадратичную интерполяцию между (a₁, f(a₁)) = (−4, −25) и (b₀, f(b₀)) = (1.33333, 0.48148) и (b₁, f(b₁)) = (1.23256, 0.22891). Это дает 1.14205, который лежит между (3a₁ + b₁) / 4 и b₁. Кроме того, неравенство |1.14205 − b₁| ≤ |b₀ − b₋₁| / 2 выполнено, поэтому это значение принято. Кроме того, f(1.14205) = 0.083582, поэтому мы устанавливаем a₂ = a₁ и b₂ = 1.14205. В третьей итерации мы используем обратную квадратичную интерполяцию между (a₂, f(a₂)) = (−4, −25) и (b₁, f(b₁)) = (1.23256, 0.22891) и (b₂, f(b₂)) = (1.14205, 0.083582). Это дает 1.09032, которое лежит между (3a₂ + b₂) / 4 и b₂. Но здесь дополнительное условие Брента вступает в силу: неравенство |1.09032 − b₂| ≤ |b₁ − b₀| / 2 не выполнено, поэтому это значение отклоняется. Вместо этого вычисляется средняя точка m = −1.42897 интервала [a₂, b₂]. У нас есть f(m) = 9.26891, поэтому мы устанавливаем a₃ = a₂ и b₃ = −1.42897. В четвертой итерации мы используем обратную квадратичную интерполяцию между (a₃, f(a₃)) = (−4, −25) и (b₂, f(b₂)) = (1.14205, 0.083582) и (b₃, f(b₃)) = (−1.42897, 9.26891). Это дает 1.15448, которое не находится в интервале между (3a₃ + b₃) / 4 и b₃. Следовательно, оно заменяется средней точкой m = −2.71449. У нас есть f(m) = 3.93934, поэтому мы устанавливаем a₄ = a₃ и b₄ = −2.71449. В пятой итерации обратная квадратичная интерполяция дает −3.45500, которое лежит в требуемом интервале. Однако, предыдущая итерация была шагом бисекции, поэтому неравенство |−3.45500 − b₄| ≤ |b₄ − b₃| / 2 должно быть выполнено. Это неравенство ложно, поэтому мы используем среднюю точку m = −3.35724. У нас есть f(m) = −6.78239, поэтому m становится новой точкой, ограничивающей интервал (a₅ = −3.35724), а итерация остается прежней (b₅ = b₄). В шестой итерации мы не можем использовать обратную квадратичную интерполяцию, потому что 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. В седьмой итерации мы можем снова использовать обратную квадратичную интерполяцию. Результат 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) В восьмой итерации мы не можем использовать обратную квадратичную интерполяцию, потому что a₇ = b₆. Линейная интерполяция дает s = −2.99994, что принимается. (Правильно: s = 2.9999, f(s) = 0.0016) В следующих итерациях корень x = −3 быстро приближается: b₉ = −3 + 6·10⁻⁸ и b₁₀ = −3 − 3·10⁻¹⁵. (Правильно: Iter 9 : f(s) = −1.4 × 10⁻⁷, Iter 10 : f(s) = 6.96 × 10⁻¹²)

Suppose that we are seeking a zero of the function defined by f(x) = (x + 3)(x − 1)2. We take [a0, b0] = [−4, 4/3] as our initial interval. We have f(a0) = −25 and f(b0) = 0.48148 (all numbers in this section are rounded), so the conditions f(a0) f(b0) < 0 and |f(b0)| ≤ |f(a0)| are satisfied. In the first iteration, we use linear interpolation between (b−1, f(b−1)) = (a0, f(a0)) = (−4, −25) and (b0, f(b0)) = (1.33333, 0.48148), which yields s = 1.23256. This lies between (3a0 + b0) / 4 and b0, so this value is accepted. Furthermore, f(1.23256) = 0.22891, so we set a1 = a0 and b1 = s = 1.23256. In the second iteration, we use inverse quadratic interpolation between (a1, f(a1)) = (−4, −25) and (b0, f(b0)) = (1.33333, 0.48148) and (b1, f(b1)) = (1.23256, 0.22891). This yields 1.14205, which lies between (3a1 + b1) / 4 and b1. Furthermore, the inequality |1.14205 − b1| ≤ |b0 − b−1| / 2 is satisfied, so this value is accepted. Furthermore, f(1.14205) = 0.083582, so we set a2 = a1 and b2 = 1.14205. In the third iteration, we use inverse quadratic interpolation between (a2, f(a2)) = (−4, −25) and (b1, f(b1)) = (1.23256, 0.22891) and (b2, f(b2)) = (1.14205, 0.083582). This yields 1.09032, which lies between (3a2 + b2) / 4 and b2. But here Brent's additional condition kicks in: the inequality |1.09032 − b2| ≤ |b1 − b0| / 2 is not satisfied, so this value is rejected. Instead, the midpoint m = −1.42897 of the interval [a2, b2] is computed. We have f(m) = 9.26891, so we set a3 = a2 and b3 = −1.42897. In the fourth iteration, we use inverse quadratic interpolation between (a3, f(a3)) = (−4, −25) and (b2, f(b2)) = (1.14205, 0.083582) and (b3, f(b3)) = (−1.42897, 9.26891). This yields 1.15448, which is not in the interval between (3a3 + b3) / 4 and b3). Hence, it is replaced by the midpoint m = −2.71449. We have f(m) = 3.93934, so we set a4 = a3 and b4 = −2.71449. In the fifth iteration, inverse quadratic interpolation yields −3.45500, which lies in the required interval. However, the previous iteration was a bisection step, so the inequality |−3.45500 − b4| ≤ |b4 − b3| / 2 need to be satisfied. This inequality is false, so we use the midpoint m = −3.35724. We have f(m) = −6.78239, so m becomes the new contrapoint (a5 = −3.35724) and the iterate remains the same (b5 = b4). In the sixth iteration, we cannot use inverse quadratic interpolation because b5 = b4. Hence, we use linear interpolation between (a5, f(a5)) = (−3.35724, −6.78239) and (b5, f(b5)) = (−2.71449, 3.93934). The result is s = −2.95064, which satisfies all the conditions. But since the iterate did not change in the previous step, we reject this result and fall back to bisection. We update s = 3.03587, and f(s) = 0.58418. In the seventh iteration, we can again use inverse quadratic interpolation. The result is s = −3.00219, which satisfies all the conditions. Now, f(s) = −0.03515, so we set a7 = b6 and b7 = −3.00219 (a7 and b7 are exchanged so that the condition |f(b7)| ≤ |f(a7)| is satisfied). (Correct : linear interpolation 1=s = 2.99436, f(s) = 0.089961)
In the eighth iteration, we cannot use inverse quadratic interpolation because a7 = b6. Linear interpolation yields s = −2.99994, which is accepted. (Correct : 1=s = 2.9999, f(s) = 0.0016)
In the following iterations, the root x = −3 is approached rapidly: b9 = −3 + 6·10−8 and b10 = −3 − 3·10−15. (Correct : Iter 9 : f(s) = −1.4 × 10−7, Iter 10 : f(s) = 6.96 × 10−12)