Документация Engee

Обзор Roots

Пакет Roots содержит простые процедуры для нахождения нулей непрерывных скалярных функций одной вещественной переменной. Нулем является значение , при котором . Основной интерфейс предоставляет функция find_zero, которая благодаря множественной диспетчеризации может обрабатывать целый ряд различных случаев.

Альтернативный интерфейс доступен в пакете NonlinearSolve.

В дальнейшем мы будем использовать ForwardDiff для вычисления производных.

julia> using Roots, ForwardDiff

Заключение в скобки

Для функции (одномерной, вещественной) скобка — это пара , для которой . То есть значения функции имеют разные знаки в точках и . Если является непрерывной функцией, это гарантирует (Больцано), что в интервале ] будет нуль. Если не является непрерывной, должна существовать точка в ], в которой функция «переходит» через .

Такие значения можно находить вплоть до округления до значения с плавающей запятой. То есть, учитывая f(a) * f(b) < 0, можно найти значение c с a < c < b, где либо f(c) == 0.0, либо f(prevfloat(c)) * f(c) < 0, либо f(c) * f(nextfloat(c)) < 0.

Для примера рассмотрим функцию . С точки зрения тригонометрии мы можем легко увидеть, что ] является скобкой.

Пакет Roots содержит алгоритм поиска делением пополам с помощью find_zero. Мы используем структуру, для которой extrema возвращает (a,b) с a < b, например вектор или кортеж, для задания начального условия и Bisection() для задания алгоритма:

julia> f(x) = cos(x) - x;

julia> x = find_zero(f, (0, pi/2), Bisection())
0.7390851332151607

julia> x, f(x)
(0.7390851332151607, 0.0)

Для этой функции мы видим, что f(x) имеет значение 0.0.

Функции могут быть параметризованы. Далее указана функция, аналогичная вышеприведенной, но в ней является скобкой, когда . Передавая значения p функции find_zero, можно решать различные связанные задачи.

julia> g(x, p=1) = cos(x) - x/p;

julia> x0, M = (0, pi/2), Bisection()
((0, 1.5707963267948966), Bisection())

julia> find_zero(g, x0, M) # как и раньше, решаем cos(x) - x = 0, используя стандартное значение p=1
0.7390851332151607

julia> find_zero(g, x0, M; p=2) # решает cos(x) - x/2 = 0
1.0298665293222589

julia> find_zero(g, x0, M, 2) # позиционный аргумент, полезен при трансляции
1.0298665293222589

Далее рассмотрим . Известным нулем является . Согласно тригонометрии ] будет скобкой. Шаблоном вызова для find_zero является find_zero(f, x0, M; kwargs...), где kwargs может указывать на параметры задачи или допуски для решателя. В этом вызове метод Bisection() не указывается, так как он будет использоваться по умолчанию (поскольку начальное значение не указано, так как число больше Float64 значений:

julia> f(x) = sin(x);

julia> x = find_zero(f, (pi/2, 3pi/2))
3.141592653589793

julia> x, f(x)
(3.141592653589793, 1.2246467991473532e-16)

Это значение x не совсем приводит к нулю, однако оно настолько близко, насколько это возможно:

julia> f(prevfloat(x)) * f(x) < 0.0 || f(x) * f(nextfloat(x)) < 0.0
true

То есть при значении x функция меняет знак.

С точки зрения математики нуль гарантируется для непрерывной функции. Однако компьютерный алгоритм не предполагает непрерывности, он просто ищет изменения знака. Таким образом, алгоритм будет определять разрывы, а не просто нули. Например:

julia> find_zero(x -> 1/x, (-1, 1))
0.0

Значения конечных точек и функции могут быть даже бесконечными для алгоритма по умолчанию Bisection для стандартных типов с плавающей запятой:

julia> find_zero(x -> Inf*sign(x), (-Inf, Inf))  # только для Float64
0.0

Базовый алгоритм, используемый для заключения в скобки в случае, когда значения являются простыми значениями с плавающей запятой, представляет собой модификацию метода деления пополам, где средняя точка берется по битовому представлению a и b.

Для больших значений с плавающей запятой метод делением пополам используется по умолчанию (с ненулевыми допусками), но его применять не рекомендуется. Простой поиск делением пополам по значениям BigFloat может занять много итераций. Для задачи нахождения нуля sin в интервале (big(3), big(4)) поиск делением пополам занимает итераций, тогда как метод A42 —  .

Алгоритмы Алефельда, Потра и Ши, а также хорошо известный алгоритм Брента, начинаются с алгоритма заключения в скобки. Для многих задач эти алгоритмы потребуют выполнения гораздо меньшего количества шагов для достижения сходимости, чем алгоритм метода деления пополам. Их можно вызывать напрямую. Например,

julia> find_zero(sin, (3,4), A42())
3.141592653589793

По умолчанию поиск делением пополам сойдется с машинным допуском. Это может обеспечить большую точность, чем хотелось бы. Можно указать допуск, чтобы завершить работу раньше времени и тем самым использовать меньше ресурсов. Например, для достижения точности в выполняется шага (без указания xatol используется шага):

julia> rt = find_zero(sin, (3.0, 4.0), xatol=1/16)
3.125

julia> rt - pi
-0.016592653589793116

Задачи без заключения в скобки

Методы заключения в скобки имеют гарантированную сходимость, но в общем случае могут потребовать гораздо больше вызовов функции, чем нужно для получения ответа, и не все нули функции заключаются в скобки. Если известно хорошее начальное предположение, функция find_zero предоставляет интерфейс для различных итерационных алгоритмов, которые являются более эффективными. В отличие от метода заключения в скобки, эти алгоритмы могут не сходиться к искомому корню, если начальное предположение выбрано неудачно.

Алгоритм по умолчанию смоделирован на основе алгоритма, используемого в калькуляторах HP-34. Этот алгоритм разработан таким образом, чтобы нестрого оценивать качество начального предположения ценой возможного выполнения большего количества шагов, чем другие алгоритмы, поскольку если алгоритм обнаруживает скобку, будет использован метод заключения в скобки (более эффективный).

Например, на графике видно, что ответ на нашу первоначальную задачу близок к 1. Учитывая это, нуль можно найти следующим образом:

julia> f(x) = cos(x) - x;

julia> x = find_zero(f , 1)
0.7390851332151607

julia> x, f(x)
(0.7390851332151607, 0.0)

Для многочлена разумным кажется начальное значение 2:

julia> f(x) = x^3 - 2x - 5;

julia> x = find_zero(f, 2)
2.0945514815423265

julia> x, f(x), sign(f(prevfloat(x)) * f(nextfloat(x)))
(2.0945514815423265, -8.881784197001252e-16, -1.0)

Для еще большей точности можно использовать числа BigFloat

julia> x = find_zero(sin, big(3))
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

julia> x, sin(x), x - pi
(3.141592653589793238462643383279502884197169399375105820974944592307816406286198, 1.096917440979352076742130626395698021050758236508687951179005716992142688513354e-77, 0.0)

Методы более высокого порядка

Вызов по умолчанию для find_zero использует метод первого порядка, а затем, возможно, заключение в скобки, что потенциально приводит к увеличению количества вызовов функции, чем необходимо. В некоторых случаях требуется более эффективный алгоритм. Здесь лучше всего подойдет метод более высокого порядка. Существуют алгоритмы Order1 (метод секущей), Order2 (Стеффенсена), Order5, Order8 и Order16. Методы 1-го и 2-го порядков обычно достаточно эффективны с точки зрения количества шагов, необходимых для работы со значениями с плавающей запятой. Методы более высоких четных порядков потенциально полезны, когда используется большая точность. Доступ к этим алгоритмам осуществляется путем указания метода после исходной начальной точки:

julia> f(x) = 2x - exp(-x);

julia> x = find_zero(f, 1, Order1())
0.3517337112491958

julia> x, f(x)
(0.3517337112491958, -1.1102230246251565e-16)

Аналогичным образом,

julia> f(x) = (x + 3) * (x - 1)^2;

julia> x = find_zero(f, -2, Order2())
-3.0

julia> x, f(x)
(-3.0, 0.0)
julia> x = find_zero(f, 2, Order8())
1.0000000131073141

julia> x, f(x)
(1.0000000131073141, 6.87206736323862e-16)

Начиная с алгоритм сходится к , показывая, что искомые нули необязательно должны быть простыми нулями. В простом нуле , где . Вообще говоря, ожидается, что непростые нули потребуют гораздо больше вызовов функции, поскольку методы больше не являются суперлинейными. Здесь Order2 использует вызов функции, Order8 —  , а Order0 —  . Метод Roots.Order2B полезен, когда ожидается множественность. Для решения этой задачи требуется вызовов функции.

Чтобы исследовать алгоритм и его сходимость, можно указать аргумент verbose=true. Для хранения промежуточных значений можно использовать объект Roots.Tracks.

Для некоторых функций может потребоваться скорректировать допуски по умолчанию, чтобы обеспечить сходимость. Допуски включают в себя atol и rtol, которые используются для проверки того, что ; xatol и xrtol для проверки того, что ; и maxiters для ограничения количества итераций в алгоритме.

Классические методы

Пакет предоставляет несколько классических методов нахождения корней, таких как Roots.Newton, Roots.Halley и Roots.Schroder. (В настоящее время они не экспортируются, поэтому для их использования необходимо указать в качестве префикса имя пакета.) Мы можем увидеть работу каждого из них на задаче, изученной Ньютоном. Метод Ньютона использует функцию и ее производную:

julia> f(x) = x^3 - 2x - 5;

julia> fp(x) = 3x^2 - 2;

julia> x = Roots.find_zero((f, fp), 2, Roots.Newton())
2.0945514815423265

julia> x, f(x)
(2.0945514815423265, -8.881784197001252e-16)

Функции задаются с помощью кортежа или функции, возвращающей (f(x), f(x)/f'(x)). Последний вариант удобен, когда f' легко вычисляется при наличии f, но в противном случае вычисление может оказаться дорогостоящим.

Метод Галлея имеет кубическую сходимость, в отличие от квадратичной сходимости метода Ньютона. Он также использует вторую производную:

julia> fpp(x) = 6x;

julia> x = Roots.find_zero((f, fp, fpp), 2, Roots.Halley())
2.0945514815423265

julia> x, f(x), sign(f(prevfloat(x)) * f(nextfloat(x)))
(2.0945514815423265, -8.881784197001252e-16, -1.0)

(Метод Галлея выполняется за 3 шага, метод Ньютона — за 4, но в методе Ньютона используется 5 вызовов функции, а в методе Галлея — 10.)

Для многих функций их производные могут быть вычислены автоматически. Пакет ForwardDiff предоставляет соответствующие возможности. Здесь мы определяем оператор D для вычисления производной:

julia> function D(f, n::Int=1)
           n <= 0 && return f
           n == 1 && return x -> ForwardDiff.derivative(f,float(x))
           D(D(f,1),n-1)
       end
D (generic function with 2 methods)

julia> dfᵏs(f,k) = ntuple(i->D(f,i-1), Val(k+1)) # (f, f′, f′′, …)
dfᵏs (generic function with 1 method)
julia> find_zero((f,D(f)), 2, Roots.Newton())
2.0945514815423265

Или для метода Галлея:

julia> find_zero((f, D(f), D(f,2)), 2, Roots.Halley())
2.0945514815423265

В семейство решателей, реализованное в Roots.LithBoonkkampIJzerman(S,D), где S — это количество предыдущих точек, используемых для генерации последующих, а D — количество используемых производных, входит метод секущей (S=2, D=0) и метод Ньютона (S=1, D=1), а также другие. При увеличении объема памяти или добавлении дополнительных производных скорость сходимости увеличивается, но за счет более сложных выражений или большего числа вызовов функции на шаг.

julia> find_zero(dfᵏs(f, 0), 2, Roots.LithBoonkkampIJzerman(3,0)) # аналогично методу секущей
2.0945514815423265

julia> find_zero(dfᵏs(f, 1), 2, Roots.LithBoonkkampIJzerman(2,1)) # аналогично методу Ньютона
2.0945514815423265

julia> find_zero(dfᵏs(f, 2), 2, Roots.LithBoonkkampIJzerman(2,2)) # аналогично методу Галлея
2.0945514815423265

Интерфейс «задача-алгоритм-решение»

Интерфейс «задача-алгоритм-решение» — это популярный шаблон в Julia, реализуемый набором пакетов DifferentialEquations.jl. Шаблон состоит из действий по постановке задачи и ее последующем решении путем указания алгоритма. Это очень похоже на то, что указано в интерфейсе find_zero(f, x0, M), где f и x0 задают задачу, M задает алгоритм, а find_zero вызывает решатель.

Для разбиения этого процесса на шаги Roots содержит методы ZeroProblem и init, solve и solve! из пакета CommonSolve.jl.

Рассмотрим решение с помощью метода Secant, начиная с интервала ].

julia> f(x) = sin(x);

julia> x0 = (3, 4)
(3, 4)

julia> M = Secant()
Secant()

julia> Z = ZeroProblem(f, x0)
ZeroProblem{typeof(f), Tuple{Int64, Int64}}(f, (3, 4))

julia> solve(Z, M)
3.141592653589793

Изменить метод несложно:

julia> solve(Z, Order2())
3.1415926535897944

Интерфейс solve работает и с параметризованными функциями:

julia> g(x, p=1) = cos(x) - x/p;

julia> Z = ZeroProblem(g, (0.0, pi/2))
ZeroProblem{typeof(g), Tuple{Float64, Float64}}(g, (0.0, 1.5707963267948966))

julia> solve(Z, Secant(), 2) # использует p=2 в качестве позиционного аргумента
1.0298665293222589

julia> solve(Z, Bisection(); p=3, xatol=1/16) # p=3; использует ключевые слова для позиций и допусков
1.1959535058744393

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

Примеры

Пересечения

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

julia> k, B, Λ = 3, 1, 1;

julia> f(x) = tan(x); g(x) = x/(B*(Λ*x^2 - 1));

julia> h(x) = f(x) - g(x)
h (generic function with 1 method)

julia> x = find_zero(h, (k*pi, (k + 1/2)*pi)); x, h(x)
(9.530477156207574, 8.326672684688674e-16)

В Julia версии 1.9 существует расширение, позволяющее при загрузке SymPy использовать уравнение для указания левой и правой сторон, как в следующем случае.

using SymPy
@syms x
find_zero(tan(x) ~ x/(B*(Λ*x^2 - 1)), (k*pi, (k + 1/2)*pi))

Обратные функции

С помощью функции find_zero можно определять обратные функции. Предположим, что является монотонной функцией в ]. Тогда обратная функция решает для , если задано . Эта функция выполнит эту задачу и вернет значения в виде функции:

function inverse_function(f, a, b, args...; kwargs...)
    fa, fb = f(a), f(b)
    m, M = fa < fb ? (fa, fb) : (fb, fa)
    y -> begin
        @assert m ≤ y ≤ M
        find_zero(x ->f(x) - y, (a,b), args...; kwargs...)
    end
end
inverse_function (generic function with 1 method)

Тот факт, что монотонна, гарантирует, что скобка ] может быть использована при условии, что находится между и , так что предположение не требуется.

Здесь мы численно находим обратную функцию :

using Plots, Roots;
f(x) = x - sin(x)
a, b = 0, 5pi
plot(inverse_function(f, a, b), f(a), f(b))
inversefunction

Нахождение критических точек

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

julia> f(x) = 1/x^2 + x^3;

julia> find_zero(D(f), 1)
0.9221079114817278

Для более сложных выражений D может потребовать некоторых технических настроек. В этом примере приведена функция, которая моделирует полет стрелы в ветреный день:

julia> function flight(x, theta)
         k = 1/2
         a = 200*cosd(theta)
         b = 32/k
         tand(theta)*x + (b/a)*x - b*log(a/(a-x))
       end
flight (generic function with 1 method)

Общее проделанное расстояние составляет flight(x) == 0.0 для x > 0: Эту задачу можно решить для различных theta с помощью find_zero. Далее мы заметим, что log(a/(a-x)) будет иметь асимптоту при a, поэтому мы начнем поиск с a-5:

julia> function howfar(theta)
         a = 200*cosd(theta)
         find_zero(x -> flight(x, theta), a-5)  # начальная точка имеет тип, определяемый `theta`.
        end
howfar (generic function with 1 method)

Чтобы представить траекторию при выстреле под углом градусов, мы должны иметь следующее.

using Plots;
flight(x, theta) = (k = 1/2; a = 200*cosd(theta); b = 32/k; tand(theta)*x + (b/a)*x - b*log(a/(a-x)))
howfar(theta) = (a = 200*cosd(theta); find_zero(x -> flight(x, theta), a-5))
howfarp(t) = ForwardDiff.derivative(howfar,t)

theta = 45
tstar = find_zero(howfarp, 45)

plot(x -> flight(x,  theta), 0, howfar(theta))
flight

Чтобы максимизировать диапазон, мы решаем задачу об одиночной критической точке howfar в обоснованных начальных точках.

Начиная с Julia версии v"1.9" автоматическое дифференцирование, обеспечиваемое ForwardDiff, будет работать в обход вызова find_zero. До этой версии автоматическое дифференцирование будет работать, если начальная точка имеет правильный тип (в данном случае зависящий от выражения theta). Поскольку мы используем 200*cosd(theta)-5 в качестве отправной точки, это условие выполняется.

julia> (tstar = find_zero(D(howfar), 45)) ≈ 26.2623089
true

На этой схеме будут показаны различия в траекториях:

plot(x -> flight(x, 45), 0, howfar(45))
plot!(x -> flight(x, tstar), 0, howfar(tstar))
flight diff

Чувствительность

В последнем примере, несомненно, важен вопрос о том, как меняется расстояние в зависимости от угла.

В общем случае для функций с параметрами интерес могут представлять производные по переменной (переменным) .

Первой попыткой, как показано выше, может быть попытка автоматического дифференцирования вывода функции find_zero. Например:

f(x, p) = x^2 - p # p является скаляром
p = 2
2
F(p) = find_zero(f, one(p), Order1(), p)
ForwardDiff.derivative(F, p)
Dual{ForwardDiff.Tag{typeof(Main.F), Int64}}(0.3535533905932738,-0.0)

До версии Julia v"1.9" с этим подходом были связаны проблемы, хотя в данном случае он находит правильный ответ, но, как мы увидим, а) он не так производителен, как то, что мы обсудим дальше; б) особое использование one(p) для начальной точки необходимо для гарантии правильного типа для значений ; в) не все алгоритмы будут работать, в частности Bisection не пригоден для этого подхода.

F(p) = find_zero(f, (zero(p), one(p)), Roots.Bisection(), p)
ForwardDiff.derivative(F, 1/2)
Dual{ForwardDiff.Tag{typeof(Main.F), Float64}}(0.7071067811865475,-0.0)

Это будет 0.0, если дифференцирование распространяется через алгоритм. В Julia версии v"1.9" или более поздней производная вычисляется корректно с помощью метода, описанного ниже.

Используя теорему о неявных функциях и следуя этим заметкам, сведениям в этой статье о методе сопряженности или более общих методах, применяемым в пакете ImplicitDifferentiation, мы можем выполнять автодифференцирование без использования этого механизма с помощью функции find_zero.

Решение , выполняемое с помощью функции find_zero, зависит от параметров . Условно это выглядит следующим образом.

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

Поскольку частное в  — скалярная величина, мы можем разделить ее, чтобы решить:

Например, используя ForwardDiff, мы получаем следующее.

xᵅ = find_zero(f, 1, Order1(), p)

fₓ = ForwardDiff.derivative(x -> f(x, p), xᵅ)
fₚ = ForwardDiff.derivative(p -> f(xᵅ, p), p)

- fₚ / fₓ
0.3535533905932738

Конечно, эту задачу можно решить аналитически, чтобы увидеть , так что мы можем легко выполнить сравнение:

ForwardDiff.derivative(sqrt, 2)
0.35355339059327373

Использование с вектором параметров аналогично, только derivative заменяется gradient для переменной p:

f(x, p) = x^2 - p[1]*x + p[2]
p = [3.0, 1.0]
x₀ = 1.0

xᵅ = find_zero(f, x₀, Order1(), p)

fₓ = ForwardDiff.derivative(x -> f(x, p), xᵅ)
fₚ = ForwardDiff.gradient(p -> f(xᵅ, p), p)

- fₚ / fₓ
2-element Vector{Float64}:
 -0.1708203932499369
  0.4472135954999579

Пакет предоставляет расширение пакета для непосредственного использования ForwardDiff для нахождения производных или градиентов, как указано выше, с версией v"1.9" или более поздней Julia, а также реализацию ChainRulesCore.rrule и ChainRulesCore.frule, которая должна разрешить пакетам автоматического дифференцирования на основе ChainRulesCore (например, Zygote), выполнять дифференцирование в p, используя тот же подход. (Спасибо @devmotion за большую помощь по этому вопросу.)

Возможные проблемы

Методы более высоких порядков — это, в сущности, различные версии метода Ньютона без производных (с шагом обновления ). Например, метод Стеффенсена (Order2()) по сути заменяет на . Это разностная аппроксимация вперед для производной, где "$h$" — это , которая, предположительно, уже близка к . Методы с более высоким порядком объединяют этот момент с различными подходами для секущей линии, которые минимизируют количество вызовов функции. Эти методы более высокого порядка могут быть подвержены некоторым обычным проблемам, характерным для метода Ньютона: плохое начальное предположение, малая первая производная или большая вторая производная вблизи нуля.

Когда первая производная близка к , значение следующего шага может быть совсем другим, так как следующий шаг обычно отслеживает точку пересечения с касательной. Мы видим, что начиная с этот поиск становится проблематичным:

julia> try  find_zero(sin, pi/2, Order1()) catch err  "Convergence failed" end
"Convergence failed"

(В то время, начиная с pi/2 + 0.3, где наклон касательной достаточно близок и направлен в сторону , сходимость будет при .)

В качестве классического примера, в котором большая вторая производная является проблемой, у нас есть :

julia> f(x) = cbrt(x);

julia> x = try  find_zero(f, 1, Order2())  catch err  "Convergence failed" end	# все 2, 5, 8 и 16 завершаются сбоем или расходятся к бесконечности
"Convergence failed"

Однако по умолчанию здесь вычисляется корень, так как обнаружена скобка:

julia> x = find_zero(f, 1)
0.0

julia> x,  f(x)
(0.0, 0.0)

Наконец, для многих функций все эти методы требуют хорошего начального предположения. Например, многочленная функция имеет единственный нуль вблизи .

Если мы начнем далеко от нуля, сходимость может произойти, но не гарантированно:

julia> f(x) = x^5 - x - 1;

julia> x0 = 0.1
0.1

julia> try find_zero(f, x0)  catch err  "Convergence failed" end
"Convergence failed"

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

xs = [0.1] # x0
n = 15
for i in 1:(n-1) push!(xs, xs[end] - f(xs[end])/D(f)(xs[end])) end
ys = [zeros(Float64,n)';map(f, xs)'][1:2n]
xs = xs[repeat(collect(1:n), inner=[2], outer=[1])]
plot(f, -1.25, 1.5, linewidth=3, legend=false)
plot!(zero, -1.25, 1.5, linewidth=3)
plot!(xs, ys)
newton

Даже на реальном графике явно различимы лишь несколько шагов, поскольку относительный максимум функции является ловушкой для этого алгоритма Если начать справа от относительного минимума, ближе к нулю, — этой ловушки можно избежать. В стандартном методе используется прием, позволяющий исключить такие неприятные моменты, но он не всегда срабатывает.

Допуски

Математическое нахождение нуля нелинейной функции может оказаться невозможным, поэтому используются численные методы. Однако применение чисел с плавающей запятой для аппроксимации реальных чисел связано с некоторыми нюансами.

Например, рассмотрим многочлен с одним нулем при и эквивалентное ему выражение . Математически это одно и то же, однако при вычислении с плавающей запятой это не так. Здесь мы рассмотрим 21 число с плавающей запятой вблизи :

julia> f(x) = (3x-1)^5;

julia> f1(x) =  -1 + 15*x - 90*x^2 + 270*x^3 - 405*x^4 + 243*x^5;

julia> above = accumulate((x,y) -> nextfloat(x), 1:10, init=1/3);

julia> below = accumulate((x,y) -> prevfloat(x), 1:10, init=1/3);

julia> ns = sort([below...,1/3, above...]); # числа с плавающей запятой вблизи 1/3

julia> maximum(abs.(f.(ns) - f1.(ns))) < 1e-14
true

Экспоненты таковы:

julia> zs .|> abs .|> log10 .|> x -> floor(Int, x)
21-element Vector{Int64}:
 -16
 -16
 -16
 -15
 -15
 -16
 -76
 -77
   ⋮
 -15
 -16
 -75

Мы видим, что значения функции близки для каждой точки, так как максимальная разница составляет . Это примерно соответствует ожиданиям, когда даже одно сложение может внести относительную погрешность со значением , а здесь их несколько.

(Эти значения зависят от особенностей вычисления с плавающей запятой, поэтому могут отличаться в зависимости от базовой архитектуры компьютера.)

Обычно об этой вариативности даже не задумываются, поскольку различия, как правило, незначительны, но когда мы хотим определить, равно ли значение нулю, эти небольшие различия могут иметь значение. Здесь мы посмотрим на знаки значений функции для выполнения, описанного выше.

julia> fs = sign.(f.(ns));

julia> f1s = sign.(f1.(ns));

julia> [ns.-1/3 fs f1s]
21×3 Matrix{Float64}:
 -5.55112e-16  -1.0  -1.0
 -4.996e-16    -1.0  -1.0
 -4.44089e-16  -1.0   1.0
 -3.88578e-16  -1.0   1.0
 -3.33067e-16  -1.0   1.0
 -2.77556e-16  -1.0   1.0
 -2.22045e-16  -1.0  -1.0
 -1.66533e-16  -1.0  -1.0
 -1.11022e-16  -1.0   1.0
 -5.55112e-17  -1.0   1.0
  0.0           0.0  -1.0
  5.55112e-17   0.0   1.0
  1.11022e-16   1.0   1.0
  1.66533e-16   1.0  -1.0
  2.22045e-16   1.0  -1.0
  2.77556e-16   1.0  -1.0
  3.33067e-16   1.0   1.0
  3.88578e-16   1.0   1.0
  4.44089e-16   1.0   1.0
  4.996e-16     1.0  -1.0
  5.55112e-16   1.0   0.0

При анализе возникает ряд неожиданностей. Во-первых, обнаружены два нуля f(x), а не один, как ожидается с математической точки зрения, — значение с плавающей запятой для 1/3 и следующее по величине число с плавающей запятой.

julia> findall(iszero, fs)
2-element Vector{Int64}:
 11
 12

Для f1(x) есть только один нуль, но он не является значением с плавающей запятой для 1/3, а находится на расстоянии 10 чисел с плавающей запятой.

julia> findall(iszero, f1s)
1-element Vector{Int64}:
 21

Кроме того, значения функции f1s несколько раз меняют знак:

julia> findall(!iszero, diff(sign.(f1s)))
9-element Vector{Int64}:
  2
  6
  8
 10
 11
 13
 16
 19
 20

Нет никакой гарантии, что нуль будет присутствовать, но для математической функции, которая меняет знак, будет хотя бы одно изменение знака.

Исходя из этого точным нулем f будет либо то, где iszero(f(x)) истинно, либо то, где функция изменяет знак (либо f(x)*f(prevfloat(x))<0, либо f(x)*f(nextfloat(x)) < 0).

Как уже упоминалось, метод Bisection() по умолчанию для find_zero определяет такие нули для f при условии, что указан начальный заключенный в скобки интервал, когда используются числа Float64. Однако если математическая функция не пересекает ось в нулевой точке, нет никакой гарантии, что значения с плавающей запятой будут удовлетворять одному из этих условий.

Теперь рассмотрим функцию f(x) = exp(x)-x^4. Значение x=8.613169456441398 является нулем в этом смысле, так как происходит изменение знака:

julia> f(x) = exp(x) - x^4;

julia> F(x) = sign(f(x));

julia> x=8.613169456441398
8.613169456441398

julia> F(prevfloat(x)), F(x), F(nextfloat(x))
(-1.0, -1.0, 1.0)

Однако значение f(x) не так мало, как можно было бы изначально ожидать от нуля:

julia> f(x), abs(f(x)/eps(x))
(-2.7284841053187847e-12, 1536.0)

Значение x — это приближение к реальному математическому нулю, назовем его . Существует разница между (математический ответ) и f(x) (ответ с плавающей запятой). Грубо говоря, мы ожидаем, что f(x) будет примерно равен , где  — это разница между x и . Значение будет в интервале abs(x) * eps(), так что в целом мы ожидаем, что ответ будет в диапазоне плюс-минус это значение:

julia> fp(x) = exp(x) - 4x^3; # производная

julia> fp(x) * abs(x) * eps()
5.637565490466956e-12

что примерно соответствует тому, что мы видим.

Метод деления пополам может быть более медленным, чем другие. Для значений Float64 выполнение Bisection() занимает не более 64 шагов, но другие методы могут сходиться к нулю за 4—​5 шагов (при условии, что заданы хорошие начальные значения).

Если требуется меньшее количество вызовов функции, проверка на наличие приближенного нуля может быть предпочтительнее оценки изменения знака, так как в этом случае на каждом шаге потребуется два дополнительных вызова функции. Кроме того, изменение знака не гарантируется для всех нулей. Приближенным нулем будет тот, при котором .

Из вышесказанного следует, что мы должны учитывать соответствующий допуск. Первый пример показывает, что отличия вычислений с плавающей запятой от математических могут вносить ошибки в масштабе eps независимо от размера x. Как видно из второго примера, разница между приближением значений с плавающей запятой к нулю и нулем приводит к потенциальной ошибке, пропорциональной размеру x. Поэтому при допуске могут учитываться оба типа ошибок. Используется как абсолютный допуск, так и относительный, поэтому проверка может выглядеть следующим образом.

abs(f(x)) < max(atol, abs(x) * rtol)

Здесь есть отличие от Julia в isapprox(f(x), 0.0), так как в этом случае в качестве множителя используется abs(f(x)), что делает относительный допуск бесполезным для данного вопроса.

Одна из проблем с относительными допусками заключается в том, что для функций с сублинейным ростом очень большие значения будут считаться нулями. Возвращаясь к предыдущему примеру, в Thukral8 мы имеем неправильно определенный нуль:

julia> find_zero(cbrt, 1, Roots.Thukral8())
1.725042287244107e23

Алгоритм быстро уходит в бесконечность, поэтому относительный допуск , используемый для проверки того, что , велик по сравнению с далеким от нуля значением .

Либо пользователи должны быть проинформированы о такой возможности, либо относительный допуск должен быть установлен равным . В этом случае абсолютный допуск должен быть относительно большим. Консервативный выбор абсолютного допуска может иметь значение sqrt(eps()), или составлять около 1e-8, что, по сути, соответствует выбору в SciPy.

Хотя для очень больших значений этот допуск работать не будет:

julia> find_zero(x -> sqrt(eps()) - eps(x), (0,Inf))
9.981132799999999e7

Это не вариант для Roots. Тот факт, что поиск делением пополам позволяет получить нули настолько точно, насколько это возможно, и тот факт, что ошибка в вычислении функции обычно не входит в диапазон 1e-8, вызывает желание получить большую точность, если она доступна.

В Roots более быстрые алгоритмы используют проверку как размера f(xn), так и размера разницы между двумя последними значениями xn. Проверка f(xn) выполняется с жестким допуском, как и проверка . Если значения функции близки к нулю, объявляется приближенный нуль. Далее, если значения близки друг к другу и значение функции близко к нулю с нестрогим допуском, объявляется приближенный нуль. На практике это работает достаточно хорошо. Нестрогий допуск использует кубический корень абсолютного и относительного допусков.

Поиск всех нулей в интервале

Описанные выше методы используются для нахождения одного из, возможно, нескольких нулей. Функция find_zeros выполняет поиск всех нулей функции в интервале . Она проста в использовании:

julia> f(x) = exp(x) - x^4;

julia> a, b = -10, 10
(-10, 10)

julia> zs = find_zeros(f, a, b)
3-element Vector{Float64}:
 -0.8155534188089606
  1.4296118247255556
  8.613169456441398

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

Алгоритм, используемый для поиска всех нулей в интервале, может работать некорректно из-за следующих нескольких моментов:

  • слишком много нулей в интервале

  • близлежащие нули (значение «близлежащие» зависит от размера , если интервал очень широкий)

Алгоритм адаптивен, поэтому он может завершиться успешно при большом количестве нулей, но может потребоваться увеличить no_pts со значения 12 по умолчанию ценой возможного увеличения времени поиска.

Здесь алгоритм определяет все нули, несмотря на то, что их несколько:

julia> f(x) = cos(x)^2 + cos(x^2); a,b = 0, 10;

julia> rts = find_zeros(f, a, b);

julia> length(rts)
32

Для близлежащих нулей алгоритм работает довольно хорошо, хотя и не идеально.

Здесь мы видим, что для  — с бесконечно большим количеством нулей вблизи  — он находит много экземпляров:

julia> f(x) = iszero(x) ? NaN : sin(1/x);  # избежать ошибки sin(Inf)

julia> rts = find_zeros(f, -1, 1);

julia> length(rts) # найдено 88 нулей
88

Функция выглядит практически аналогично , поэтому find_zeros может завершиться успешно, поскольку два нуля находятся очень близко:

julia> f(x) =  (x-0.5)^3 * (x-0.499)^3;

julia> find_zeros(f, 0, 1)
1-element Vector{Float64}:
 0.5

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

Алгоритм будет выполнен успешно, когда значения показателей степени будут меньше.

julia> f(x) =  (x-0.5)^2 * (x-0.499)^2;

julia> find_zeros(f, 0, 1)
2-element Vector{Float64}:
 0.49899999999999994
 0.5

Он может быть завершен для более близких пар нулей:

julia> f(x) = (x-0.5) * (x - 0.49999);

julia> find_zeros(f, 0, 1)
2-element Vector{Float64}:
 0.49999
 0.5

Сочетания нулей большой (четной) кратности или очень близких нулей могут привести к неправильному нахождению значения.

IntervalRootFinding

Пакет IntervalRootFinding строго определяет изолирующие интервалы для нулей функции. Этот пример, взятый из файла сведений (README) пакета, демонстрирует различия:

julia> using IntervalArithmetic, IntervalRootFinding, Roots

julia> f(x) = sin(x) - 0.1*x^2 + 1
f (generic function with 1 method)

julia> rts = roots(f, -10..10)
4-element Vector{Root{Interval{Float64}}}:
 Root([3.14959, 3.1496], :unique)
 Root([-4.42654, -4.42653], :unique)
 Root([-3.10682, -3.10681], :unique)
 Root([-1.08205, -1.08204], :unique)

julia> find_zeros(f, -10, 10)
4-element Vector{Float64}:
 -4.426534982071949
 -3.1068165552293254
 -1.0820421327607177
  3.1495967624505226

Используя тот факт, что в данном случае интервалы являются заключенными в скобки интервалами f, можно найти нули из вывода roots следующим образом.

julia> [find_zero(f, (interval(u).lo, interval(u).hi)) for u ∈ rts if u.status == :unique]
4-element Vector{Float64}:
  3.1495967624505226
 -4.426534982071949
 -3.1068165552293254
 -1.082042132760718

В версии 1.9 пакета Julia доступно расширение, чтобы при загрузке пакета IntervalRootFinding функция find_zeros вызывала IntervalRootFinding.roots для нахождения изолирующих скобок и find_zero для нахождения корней, когда это возможно, если интервал задан в виде объекта Interval, созданного, скажем, с помощью -1..1.

Добавление решателя

Чтобы добавить решатель, достаточно иметь тип для объявления решателя и метод update_state. В этом примере мы также определяем объект состояния, так как алгоритм использует больше значений, хранящихся в памяти, чем задано по умолчанию.

На странице Википедии, посвященной методу Брента, предлагается современное усовершенствование — метод Чандрапатлы, описанный здесь. Этого описания мы в основном придерживаемся ниже и в реализации пакета Roots.Chandrapatla.

Для реализации алгоритма Чандрапатлы мы сначала определи тип для обозначения метода и объект состояния, который записывает значения , и , необходимые для обратного квадратичного шага.

julia> struct Chandrapatla <: Roots.AbstractBracketingMethod end

julia> struct ChandrapatlaState{T,S} <: Roots.AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    c::T
    fxn1::S
    fxn0::S
    fc::S
end

Метод init_state может использоваться некоторыми методами для добавления дополнительных деталей к базовому объекту состояния. Здесь он запускает старое значение c как a, чтобы определить начальный шаг поиска делением пополам.

julia> function init_state(::Chandrapatla, F, x₀, x₁, fx₀, fx₁)
    a, b, fa, fb = x₁, x₀, fx₁, fx₀
    c, fc = a, fa
    ChandrapatlaState(b, a, c, fb, fa, fc)
end

Основной алгоритм реализуется в методе update_state. Макрос @set! из Setfield.jl используется для изменения объекта состояния, который в противном случае является неизменяемым.

julia> import Roots.Setfield: @set!;

julia> function Roots.update_state(::Chandrapatla, F, o, options, l=NullTracks())

    b, a, c = o.xn1, o.xn0, o.c
    fb, fa, fc = o.fxn1, o.fxn0, o.fc

    # кодирование: a = xₙ, b = xₙ₋₁, c = xₙ₋₂
    ξ = (a - b) / (c - b)
    ϕ = (fa - fb) / (fc - fb)
    ϕ² = ϕ^2
    Δ = (ϕ² < ξ) && (1 - 2ϕ + ϕ² < 1 - ξ) # Неравенство Чандрапатлы для определения следующего шага

    xₜ = Δ ? Roots.inverse_quadratic_step(a, b, c, fa, fb, fc) : a + (b-a)/2

    fₜ = F(xₜ)
    incfn(l)

    if sign(fₜ) * sign(fa) < 0
        a, b, c = xₜ, a, b
        fa, fb, fc = fₜ, fa, fb
    else
        a, c = xₜ, a
        fa, fc = fₜ, fa
    end

    @set! o.xn0 = a
    @set! o.xn1 = b
    @set! o.c = c
    @set! o.fxn0 = fa
    @set! o.fxn1 = fb
    @set! o.fc = fc

    return (o, false)
end

Этот алгоритм выбирает обратный квадратичный шаг или шаг поиска делением пополам в зависимости от соотношения между вычисленными ξ и Φ. Допуски указаны из заданных по умолчанию значений для AbstractBracketingMethod.

Чтобы убедиться в том, что алгоритм работает, нужно выполнить следующее.

julia> find_zero(sin, (3,4), Chandrapatla())
3.1415926535897927

julia> find_zero(x -> exp(x) - x^4, (8,9), Chandrapatla())
8.613169456441398

julia> find_zero(x -> x^5 - x - 1, (1,2), Chandrapatla())
1.1673039782614185