Визуализация шага для различных алгоритмов нахождения нулей
Метод Ньютона
Помимо Roots, мы используем пакеты Plots и ForwardDiff:
using Roots
using Plots, ForwardDiff
Base.adjoint(f::Function) = x -> ForwardDiff.derivative(f, float(x)) # f' вычисляет производнуюАлгоритм нахождения нулей решает уравнение или, возможно, для значения . Здесь мы обсуждаем итерационные алгоритмы, которые используют результат одного или нескольких предыдущих шагов для выполнения следующего шага (то есть для некоторой функции , представляющей алгоритм).
Метод Ньютона — это итерационный алгоритм нахождения нулей, с которым знакомятся на уроках математики после ввода понятия касательной.
Значение описывается как точка пересечения оси с касательной в точке . Чтобы представить это в явном виде, мы подставим в уравнение касательной :
Решение дает нам новую формулу:
Сходимость не гарантируется для всех начальных предположений , но для простого нуля непрерывно дифференцируемой функции существует некоторый интервал в окрестности нуля, в котором обеспечивается квадратичная сходимость.
Чтобы представить метод Ньютона геометрически, можно построить график касательной.
Функция не имеет простого аналитического решения для единственного вещественного нуля, так как это многочлен пятой степени. Однако на графике видно, что она имеет один нуль в интервале от до , ближе к . Начиная с мы получаем следующее визуальное представление :
f(x) = x^5 - x - 1
x0 = 1.4
α = find_zero((f, f'), x0, Roots.Newton())
tl(x) = f(x0) + f'(x0)*(x-x0)
x1 = x0 - f(x0)/f'(x0)
p = plot(f, 1.1, 1.5; legend=false, linewidth=3)
plot!(zero)
plot!(tl; color="red", linewidth=3)
scatter!([x0, x1], [0, 0]; markercolor="blue")
annotate!([(x0,0,"x0", :bottom), (x1, 0, "x1", :bottom)])
scatter!([x0], [f(x0)]; markercolor=:blue)
scatter!([α], [0]; markercolor=:blue)
annotate!([(α, 0, "α", :top)])
pДля нахождения нуля мы использовали Roots.Newton().
Метод секущей
Метод секущей гораздо старше метода Ньютона, хотя и похож на него тем, что на следующем шаге алгоритма используется пересечение линии с осью . Угловой коэффициент секущей традиционно вычисляется проще, чем угловой коэффициент касательной, для которого требуется понятие производной. Для метода секущей требуются две начальные точки, и , и используется секущая вместо касательной. Секущая имеет угловой коэффициент . В результате получается следующий алгоритм:
Метод секущей достаточно легко визуализировать. Допустим, мы начинаем с точек и :
x0, x1 = 1.4, 1.3
x2 = x1 - (x1-x0)/(f(x1)-f(x0)) * f(x1)
sl(x) = f(x1) + (f(x1)-f(x0))/(x1-x0) * (x-x1)
p = plot(f, 1.1, 1.5; legend=false, linewidth=3)
plot!(zero)
plot!(sl, color=:red, linewidth=3)
scatter!([x0, x1, x2], [0,0,0]; markercolor=:blue)
annotate!([(x0,0,"x0", :bottom), (x1, 0, "x1", :bottom), (x2,0,"x2", :bottom)])
scatter!([x0, x1], [f(x0), f(x1)]; markercolor=:blue)
scatter!([α],[0]; markercolor=:blue)
annotate!([(α, 0, "α", :top)])
pМетод секущей реализован в Secant().
Метод Стеффенсена (Root.Steffensen()) похож на метод секущей, но в качестве точек берутся не и , а и . При стремлении к стремится к , поэтому этот метод сходится с асимптотической скоростью, как и метод Ньютона (но ценой того, что он требует двух новых вычислений функции за шаг, в то время как метод секущей — только одного).
Обратное представление
Как показано выше, понятие секущей вполне естественное, но на него можно взглянуть немного иначе. Возьмем две точки: и . Две несовпадающие точки определяют прямую. В данном случае мы поменяли местами значения и , то есть обратили координаты прямой. Для нахождения или прямой, заданной в какой-либо иной форме, необходимо решить два уравнения с двумя неизвестными. В каждое уравнение подставляется известная точка:
Эту систему линейных уравнений можно решить, например, с помощью такого кода SymPy:
using SymPy @syms x0, y0, x1, y1, m ,b u = solve([x0 ~ y0 * m + b, x1 ~ y1 * m + b], (m,b))
Результат
Dict{Any, Any} with 2 entries:
b => (-x0*y1 + x1*y0)/(y0 - y1)
m => (x0 - x1)/(y0 - y1)
Значение m — это величина, обратная угловому коэффициенту, так как мы обратили вид. Значение b — это место пересечения обратной прямой с осью , такое же, как при использовании метода секущей:
sm = x1 - y1 * (x1-x0)/(y1-y0) simplify(sm - u[b])
что дает:
0
Обратные квадратичный и кубический методы
Метод Брента (Roots.Brent()) — это ограничивающий метод, предусматривающий использование обратного квадратичного шага для ускорения сходимости по сравнению с методом секущей. На обратном квадратичном шаге используется тот факт, что три (неколлинеарные) точки определяют квадратный многочлен. Как и ранее, берутся точки, обратные точкам , и . Используя представленный выше метод, можно показать, что при :
x0, x1, x2 = xs = 1.4, 1.3, 1.2
fx0, fx1, fx2 = ys = f.(xs)
@syms x[0:2], y[0:2], a, b, c
u = solve([xᵢ ~ a*yᵢ^2 + b * yᵢ + c for (xᵢ, yᵢ) ∈ zip(x, y)], (a, b, c))
x3 = u[c]
for (k, v) ∈ u
for (xᵢ, yᵢ, x,y) ∈ zip(x, y, xs, ys)
v = v(xᵢ => x, yᵢ => y)
end
u[k] = v
end
u[a], u[b], u[c]
В результате возвращается:
(-0.00930682152998560, 0.104752944517765, 1.17057129242798)
Последнее значение, c, также вычисляется в методе (3,0) алгоритма LithBoonkkampIJzerman, который реализует данный метод:
x0, x1, x2 = 1.4, 1.3, 1.2
xs = [x0, x1,x2]
x3 = Roots.lmm(Roots.LithBoonkkampIJzerman{3,0}(), xs, f.(xs))1.1705712924279839Визуализировать это можно следующим образом:
a, b, c = (-0.00930682152998560, 0.104752944517765, 1.17057129242798)
iq(y) = a * y^2 + b * y + c
p = plot(f, 1.1, 1.5; legend=false, linewidth=3)
plot!(zero)
ys′ = range(f(1.1), f(1.5), length=100)
plot!(iq.(ys′), ys′; color=:red, linewidth=3)
scatter!(xs, f.(xs); markercolor=:blue)
scatter!([x3], [f(x3)]; markercolor=:blue)
annotate!([(x0,0,"x0", :bottom), (x1, 0, "x1", :bottom),
(x2,0,"x2", :bottom), (x3,0,"x3", :bottom)])
scatter!(xs, zero.(xs); markercolor=:blue)
scatter!([α],[0]; markercolor=:blue)
annotate!([(α, 0, "α", :top)])
pОбратный кубический метод похож на предыдущий, но для определения следующей точки берутся предыдущие. В данном примере мы решаем полученную систему линейных управлений численно (метод Roots.LithBoonkkampIJzerman{4,0}() реализует алгоритм):
xs = [1.4, 1.35, 1.3, 1.25]
ys = f.(xs)
A = zeros(Float64, 4, 4)
for i ∈ reverse(0:3)
A[:,4-i] .= ys.^i
end
a, b, c, d = A \ xs
ic(y) = a * y^3 + b * y^2 + c * y + d
p = plot(f, 1.1, 1.5; legend=false, linewidth=3)
plot!(zero)
plot!(ic.(ys′), ys′; color=:red, linewidth=3)
x4 = d
scatter!(vcat(xs, x4), zeros(5); markercolor=:blue)
for (i,x) ∈ enumerate(xs)
annotate!([(x, 0, "x$(i-1)", :bottom)])
end
annotate!([(x4, 0, "x4", :bottom)])
scatter!(xs, f.(xs); markercolor=:blue)
scatter!([α],[0]; markercolor=:blue)
annotate!([(α, 0, "α", :top)])
pНа графике видно, что для этой функции и выбранных значений приближения обратным квадратичным и обратным кубическим методами очень близки к фактическому нулю, что свидетельствует о быстрой сходимости. Roots.Brent(), Roots.Chandrapatlu() и Roots.AlefeldPotraShi() — это различные ограничивающие методы, предусматривающие использование обратного квадратичного шага, если это целесообразно, или других методов оценки в противном случае. Аналогичным образом, обратный кубический шаг по возможности используется алгоритмом Roots.A42(). Методы LithBoonkkampIJzerman{S,0} используют предыдущих точек ( ) и соответствующий обратный полиномиальный шаг. Так как они не являются ограничивающими, сходимость алгоритмов гарантируется только для близлежащих начальных предположений.
Разновидности метода Ньютона с производными более высокого порядка
Метод Ньютона можно визуализировать не только как пересечение оси с определенной касательной. Его можно также представить как пересечение оси с прямой с двумя ограничениями:
-
Точка лежит на прямой.
-
Угловой коэффициент прямой совпадает с угловым коэффициентом касательной в этой точке (условие касания).
В совокупности это означает, что если ограниченная условиями прямая имеет угловой коэффициент и проходит через точку , то же самое верно и в отношении касательной. Отсюда вытекает метод Ньютона.
При выборе квадратичной аппроксимации графика в точке возможны и другие алгоритмы.
Геометрические построения итерационных функций для решения нелинейных уравнений, предложенные Аматом (Amat), Бускьер (Busquier) и Гутьерресом (Gutiérrez), предполагают систематический подход, которому мы и будем следовать.
Метод Эйлера
Рассмотрим квадратичное выражение . При том условии, что выражение проходит через точку и выполняются условия касания и , мы получаем в качестве решения многочлен Тейлора второго порядка .
Допустим
Тогда пересечение кривой с осью можно представить так:
Этот метод известен как метод Эйлера или иррациональный метод Галлея и реализован в Roots.IrrationalHalley():
L_f(x) = f(x) * f''(x) / (f'(x))^2
x0 = 1.4
x1 = x0 - 2 / (1 + sqrt(1 - 2L_f(x0))) * f(x0)/f'(x0)
t2(x) = f(x0) + f'(x0)*(x-x0) + f''(x0)/2 * (x - x0)^2
a, b = 1.1, 1.5
p = plot(f, a, b; legend=false, linewidth=3)
plot!(zero)
plot!(t2; color=:red, linewidth=3)
scatter!([x0, x1], [0,0]; markercolor=:blue)
annotate!([(x0,0,"x0", :bottom), (x1, 0, "x1", :bottom)])
scatter!([x0], [f(x0)]; markercolor=:blue)
scatter!([α],[0]; markercolor=:blue)
annotate!([(α, 0, "α", :top)])
pМетод Галлея
Общая форма квадратного уравнения уточняется выше путем задания и с последующим наложения того условия, что точка является решением, а также условий касания. Известный метод Галлея можно представить как специализацию гиперболы: . В результате получается кривая
и итерационный алгоритм:
Визуализировать это можно следующим образом, используя контурный график для представления гиперболы:
x1 = x0 - 2 / (2 - L_f(x0)) * f(x0)/f'(x0)
F(x,y) = y - f(x0) - f'(x0)*(x-x0) - f''(x0)/(2f'(x0)) * (x-x0) * (y-f(x0))
a, b = 1.1, 1.5
p = plot(f, a, b; legend=false, linewidth=3)
plot!(zero)
xs, ys = range(a, b, length=50), range(f(a), f(b), length=50);
zs = [F(x,y) for y ∈ ys, x ∈ xs];
contour!(xs, ys, zs; levels = [0], color=:red, linewidth=3)
scatter!([x0, x1], [0,0]; markercolor=:blue)
annotate!([(x0,0,"x0", :bottom), (x1, 0, "x1", :bottom)])
scatter!([x0], [f(x0)]; markercolor=:blue)
scatter!([α],[0]; markercolor=:blue)
annotate!([(α, 0, "α", :top)])
pРеализация содержится в Roots.Halley().
Метод Чебышева
Метод Чебышева использует обратное квадратичное приближение, специализированное уравнением , для вычисления следующей итерации. Выразить это можно так:
И получаем следующий алгоритм:
Визуализация похожа на предыдущий пример:
x1 = x0 - (1 + 1/2 * L_f(x0)) * f(x0) / f'(x0)
F(x, y) = -f''(x0)/(2f'(x0)^2) * (y-f(x0))^2 + y - f(x0) - f'(x0) * (x- x0)
a, b = 1.1, 1.5
p = plot(f, a, b; legend=false, linewidth=3)
plot!(zero)
xs, ys = range(a, b, length=50), range(f(a), f(b), length=50);
zs = [F(x,y) for y ∈ ys, x ∈ xs];
contour!(xs, ys, zs; levels = [0], color=:red, linewidth=3)
scatter!([x0, x1], [0,0]; markermarkercolor=:blue)
annotate!([(x0,0,"x0", :bottom), (x1, 0, "x1", :bottom)])
scatter!([x0], [f(x0)]; markermarkercolor=:blue)
scatter!([α],[0]; markermarkercolor=:blue)
annotate!([(α, 0, "α", :top)])
pРеализация содержится в Roots.InverseQuadratic(), а более быстрая версия — в Roots.ChebyshevLike().
Амат, Бускьер и Гутьеррес также рассматривают следующую гиперболу:
Так как здесь неизвестные и только ограничения, решение зависит от параметра, который они называют , что дает следующую систему:
В результате получается следующий алгоритм:
При мы приходим к методу Ньютона, при — к методу Чебышева, а при — к методу Галлея.
Суперметод Галлея имеет место при . Визуализировать это можно так:
cn = -f'(x0)
bn = -f''(x0)/f'(x0)
an = -f''(x0)/(2f'(x0)^2) - bn/f'(x0)
x1 = x0 - (1 + 1/2 * (L_f(x0) / (1 + bn * f(x0) / f'(x0)))) * f(x0)/f'(x0)
F(x, y) = x - x0 + (y-f(x0)) * (1 + an * (y - f(x0))) / (bn * (y - f(x0)) + cn)
a, b= 1.1, 1.5
p = plot(f, a, b; legend=false, linewidth=3)
plot!(zero)
xs, ys = range(a, b, length=50), range(f(a), f(b), length=50);
zs = [F(x,y) for y ∈ ys, x ∈ xs];
contour!(xs, ys, zs; levels = [0], color=:red, linewidth=3)
scatter!([x0, x1], [0,0]; markercolor=:blue)
annotate!([(x0,0,"x0", :bottom), (x1, 0, "x1", :bottom)])
scatter!([x0], [f(x0)]; markercolor=:blue)
scatter!([α],[0]; markercolor=:blue)
annotate!([(α, 0, "α", :top)])
pРеализация содержится в Roots.SuperHalley().
Авторы обсуждают использование различных соприкасающихся кривых, например кубического уравнения. Как они отмечают, все их методы принимают следующую форму (в нотации «О большое»):
Алгоритмы, реализованные в методах Roots.LithBoonkkampIJzerman{S,D}(), предполагают подход на основе дифференциальных уравнений к вычислению функции, обратной , в точке . Все эти методы находят обратные полиномиальные аппроксимации . Для эти методы представляют обратную секущую при , обратную квадратичную аппроксимацию при и обратную кубическую аппроксимацию при . При и метод Эйлера обращается в метод Ньютона, а при и применяются обратный квадратичный метод или метод Чебышева.