6.5. ITP метод#
В общем виде алгоритм с локализацией корня (бисекция, regula falsi, метод Риддерса) можно описать так.
Пусть дана функция \(f\), начальный интервал \((a, b)\), такой что \(f(a) f(b) < 0\) и \(\text{xtol}\). Алгоритм состоит в следующих шагах.
Выбрать некоторое приближение решения \(\tilde{x} \in (a, b)\).
По значению \(f(\tilde{x})\) обновить интервал \((a, b)\) так, чтобы функция имела разные знаки на концах нового интервала.
Если достигнута точность \(|b - a| < \text{xtol}\), завершить, иначе – перейти на шаг 1.
Метод бисекции обладает замечательным свойством – он заранее предсказывает необходимое число итераций \(n_{1/2}\) по заданной точности локализации корня \(\text{xtol}\).
Авторы работы [OT20] показали, что это свойство (называемое minmax optimality) сохраняется во всяком алгоритме с локализацией корня, если новое приближение \(\tilde{x}\) попадает в область \(I = [x_{1/2}-r, x_{1/2}+r]\) вокруг середины интервала \(x_{1/2}\). Ниже приведено выражение для \(r\) на \(i\)-ой итерации алгоритма
Так, можно модифицировать Алгоритм 6.23, добавив после шага 1 проверку: попадает ли приближение \(\tilde{x}\) в область \(I\). Если попадает, то оставить его, как есть, а если не попадает, взять приближение как у бисекции \(\tilde{x} := x_{1/2}\). Таким образом, модифицированный алгоритм гарантированно будет работать не хуже бисекции, а в лучшем случае как алгоритм, взятый за основу (например, метод Риддерса).
6.5.1. Описание метода#
Авторами [OT20] был предложен метод ITP (Interpolate, Truncate and Project).
Interpolate. На этом шаге определяется \(x_f\) приближение по методу regula falsi (пересечение абсциссы секущей).
Truncate. В этом методе, помимо описанной выше модификации присутствует шаг truncate, смещающий найденный при интерполяции корень в сторону \(x_{1/2}\). Делается это в надежде на уменьшение интервала \([a, b]\) более чем в два раза. Результат этого шага точка \(x_t\).
Project. На этом шаге определяется, попадает ли \(x_t\) в область minmax optimality. Если \(x_t\) не попадает, то он перемещается (проецируется) в эту область и получается \(x_{\text{ITP}}\). В конце итерации происходит выбор нового отрезка по точкам \(a_i\), \(b_i\) и \(x_{\text{ITP}}\).
В методе присутствует три параметра. Параметр \(n_0 \ge 0\) определяет число дополнительных итераций, в результате в формуле (6.2) присутствует не \(n_{1/2}\), а \(n_\max = n_{1/2} + n_0\). Параметры \(\kappa_1 > 0\) и \(\kappa_2 \in [1, 1 + \phi)\) (\(\phi \approx 1.618\), золотое сечение) определяют, как происходит truncate шаг. Так, они задают область \([x_{1/2} - \delta, x_{1/2}+\delta]\), где
Если \(x_f\) находится за пределами области \(\delta \le |x_{1/2} - x_f|\), то он смещается в сторону середины \(x_{1/2}\) на \(\delta\). Иначе устанавливается в \(x_t := x_{1/2}\).
Авторами метода показано, что в случае дважды дифференциируемой функции метод сходится к простому корню со скоростью \(\sqrt{\kappa_2}\). Например, в случае \(\kappa_2 = 2\) ожидается сверхлинейная сходимость со скоростью \(\sqrt{2}\).
6.5.2. Реализация#
Метод ITP
"""
Метод ITP поиска корня `f`(x) = 0 c точностью `xtol`.
"""
function itproot(f, x₁, x₂; xtol=eps(), ftol=eps(), κ₁=0.1, κ₂=2, n₀=1)
if x₁ > x₂; x₁, x₂ = x₂, x₁; end
y₁, y₂ = f(x₁), f(x₂)
y₁ * y₂ > 0 && error("Функция должна иметь разные знаки в концах отрезка")
abs(y₁) < ftol && return x₁
abs(y₂) < ftol && return x₂
nbisect = ceil(Int, log2((x₂-x₁)/xtol))
maxiter = nbisect + n₀
brackorig = x₂ - x₁
for i in 1:maxiter
# interpolate
xf = (y₂*x₁ - y₁*x₂)/(y₂ - y₁)
# truncate
xmid = (x₁ + x₂)/2
σ = sign(xmid - xf)
δ = κ₁ * (x₂ - x₁)^κ₂ / brackorig
xt = δ ≤ abs(xmid - xf) ? xf + copysign(δ, σ) : xmid
# project
r = xtol * 2.0^(maxiter - i) - (x₂ - x₁)/2
xnew = abs(xt - xmid) ≤ r ? xt : xmid - copysign(r, σ)
ynew = f(xnew)
if sign(y₂) == sign(ynew)
x₂, y₂ = xnew, ynew
elseif sign(y₁) == sign(ynew)
x₁, y₁ = xnew, ynew
else # ynew == 0
return xnew
end
abs(ynew) < ftol && (x₁ + x₂)/2
abs(x₂ - x₁) < xtol && (x₁ + x₂)/2
end
return (x₁ + x₂)/2
end
Рассмотрим работу метода на функции \(xe^x - 1\).
f = (x) -> x*exp(x) - 1
@show itproot(f, -1, 1; xtol=2e-10, n₀=0)
@show ridders(f, -1, 1; xtol=2e-10, maxiter=4);
itproot(f, -1, 1; xtol = 2.0e-10, n₀ = 0) = 0.5671432904087239
ridders(f, -1, 1; xtol = 2.0e-10, maxiter = 4) = 0.5671432904097838
Если немного поменять функцию itproot, то можно увидеть, что метод справился с поиском корня за 8 итераций, тогда как бисекция требует 34. В то же время метод Риддерса справляется за 4 итерации, но он требует 8 вычислений функции (по два вычисления за итерацию).
На функции \(\ln |x - 10/9|\) ITP выигрывает по числу вычислений функции.
f = (x) -> (log ∘ abs)(x - 10/9)
@show itproot(f, -1, 1; xtol=2e-10, n₀=0)
@show ridders(f, -1, 1; xtol=2e-10, maxiter=5);
itproot(f, -1, 1; xtol = 2.0e-10, n₀ = 0) = 0.11111111111111119
ridders(f, -1, 1; xtol = 2.0e-10, maxiter = 5) = 0.1111111111111111
Здесь ITP потребовалось 8 вычислений, бисекции 34, а методу Риддерса 10.
Сравнение на большем числе функций вы можете посмотреть в оригинальной работе [OT20].