ITP метод

6.5. ITP метод#

В общем виде алгоритм с локализацией корня (бисекция, regula falsi, метод Риддерса) можно описать так.

Алгоритм 6.23 (Алгоритм с локализацией корня)

Пусть дана функция \(f\), начальный интервал \((a, b)\), такой что \(f(a) f(b) < 0\) и \(\text{xtol}\). Алгоритм состоит в следующих шагах.

  1. Выбрать некоторое приближение решения \(\tilde{x} \in (a, b)\).

  2. По значению \(f(\tilde{x})\) обновить интервал \((a, b)\) так, чтобы функция имела разные знаки на концах нового интервала.

  3. Если достигнута точность \(|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.2)#\[r_i = \text{xtol } 2^{n_{1/2} - i} - \frac{b_i - a_i}{2}, \quad i = 1, 2, ..., n_{1/2}.\]

Так, можно модифицировать Алгоритм 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]\), где

\[\delta = \kappa_1 (b_i - a_i)^{\kappa_2} / (b_0 - a_0).\]

Если \(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. Реализация#

Функция 6.24 (itproot)

Метод 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
Демонстрация 6.25 (Метод ITP)

Рассмотрим работу метода на функции \(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].

Содержание