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

Справка по API

Пакет Roots предоставляет несколько разных алгоритмов для решения уравнения f(x)=0.

Roots

Пакет для решения f(x) = 0 для одномерных скалярных функций.

Доступны следующие основные методы:

  • find_zero для использования одного из нескольких методов для нахождения нуля

  • ZeroProblem для нахождения нуля с помощью интерфейса CommonSolve

  • find_zeros для эвристического нахождения всех нулей в заданном интервале

Расширенная справка

Функции нахождения корней в Julia

docs stable blue docs dev blue badge badge

Этот пакет содержит простые процедуры для нахождения корней (или нулей) скалярных функций одной вещественной переменной с помощью математический операций с плавающей запятой. Функция find_zero предоставляет основной интерфейс. Основной вызов имеет следующий вид: find_zero(f, x0, [M], [p]; kws...), где, как правило, f — это функция, x0 — это начальная точка или заключаемый в скобки интервал, M используется для настройки алгоритмов, применяемых по умолчанию, а p можно использовать для передачи параметров.

В число различных алгоритмов входят следующие.

  • Алгоритмы в стиле метода деления пополам. Для функций, в которых известен заключаемый в скобки интервал (когда f(a) и f(b) имеют противоположные знаки), можно задать метод заключения в скобки, например Bisection. По умолчанию для большинства типов чисел с плавающей запятой используется Bisection в соответствии с соглашениями о хранении значений с плавающей запятой. Для других числовых типов (например, BigFloat) по умолчанию используется алгоритм Алефельда, Потра и Ши. Эти заданные по умолчанию методы гарантированно сходятся. Также доступны другие методы заключения в скобки.

  • Несколько алгоритмов без производных. Они задаются с помощью методов Order0, Order1 (метод секущей), Order2 (метод Стеффенсена), Order5, Order8 и Order16. Число приблизительно указывает порядок сходимости. Метод Order0 используется по умолчанию и является самым надежным, но для сходимости может потребоваться больше вызовов функции, так как при возможности применяет метод заключения в скобки. От использования методов более высокого порядка ожидается более высокая сходимость, однако они не всегда выдают результаты с меньшим количеством вызовов функций по сравнению с Order1 или Order2. Методы Roots.Order1B и Roots.Order2B являются суперлинейными и квадратично сходящимися методами, не зависящими от кратности нуля.

  • Существуют исторические алгоритмы, которые требуют указания одной или двух производных: Roots.Newton и Roots.Halley. Roots.Schroder предоставляет квадратичный метод, подобный методу Ньютона, который не зависит от кратности нуля. Он обобщается с помощью Roots.ThukralXB (где X — это 2, 3, 4 или 5).

  • Существует несколько неэкспортируемых алгоритмов, таких как Roots.Brent(), Roots.LithBoonkkampIJzermanBracket и Roots.LithBoonkkampIJzerman.

Дополнительные сведения можно найти в документации по каждому методу.

Некоторые примеры:

julia> using Roots

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

julia> α₀, α₁, α₂ = -0.8155534188089607, 1.4296118247255556, 8.6131694564414;

julia> find_zero(f, (8,9), Bisection()) ≈ α₂ # метод деления пополам с указанной скобкой
true

julia> find_zero(f, (-10, 0)) ≈ α₀ # Поиск делением пополам используется по умолчанию, если x в `find_zero(f, x)` не является скаляром
true


julia> find_zero(f, (-10, 0), Roots.A42()) ≈ α₀ # меньшее количество вычислений функции, чем при поиске делением пополам
true

Для методов без заключения в скобки начальная позиция передается в виде скалярного значения или, возможно, для методов секущих в виде итерируемой переменной типа (x_0, x_1):

julia> find_zero(f, 3) ≈ α₁  # find_zero(f, x0::Number) будет использовать Order0()
true

julia> find_zero(f, 3, Order1()) ≈ α₁ # тот же ответ, другой метод (метод секущей)
true

julia> find_zero(f, (3, 2), Order1()) ≈ α₁ # запуск метода секущей с (3, f(3), (2, f(2))
true


julia> find_zero(sin, BigFloat(3.0), Order16()) ≈ π # 2 итерации до 6 с использованием Order1()
true

Функцию find_zero можно использовать с вызываемыми объектами:

julia> using Polynomials;

julia> x = variable();

julia> find_zero(x^5 - x - 1, 1.0) ≈ 1.1673039782614187
true

Функция должна учитывать единицы пакета Unitful:

julia> using Unitful

julia> s, m  = u"s", u"m";

julia> g, v₀, y₀ = 9.8*m/s^2, 10m/s, 16m;


julia> y(t) = -g*t^2 + v₀*t + y₀
y (generic function with 1 method)

julia> find_zero(y, 1s)  ≈ 1.886053370668014s
true

Метод Ньютона можно использовать без взятия производных вручную. В следующем примере используется пакет ForwardDiff:

julia> using ForwardDiff

julia> D(f) = x -> ForwardDiff.derivative(f,float(x))
D (generic function with 1 method)

Теперь мы имеем:

julia> f(x) = x^3 - 2x - 5
f (generic function with 1 method)

julia> x0 = 2
2

julia> find_zero((f, D(f)), x0, Roots.Newton()) ≈ 2.0945514815423265
true

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

julia> using Statistics: mean, median

julia> as = rand(5);

julia> M(x) = sum((x-a)^2 for a in as)
M (generic function with 1 method)

julia> find_zero(D(M), .5) ≈ mean(as)
true

julia> med(x) = sum(abs(x-a) for a in as)
med (generic function with 1 method)

julia> find_zero(D(med), (0, 1)) ≈ median(as)
true

Интерфейс CommonSolve

Интерфейс DifferentialEquations для постановки задачи, инициализации задачи и ее решения также реализуется с помощью типов ZeroProblem и методов init, solve! и solve (из CommonSolve).

Например, мы можем решить задачу множеством различных методов, как показано ниже.

julia> f(x) = exp(-x) - x^3
f (generic function with 1 method)

julia> x0 = 2.0
2.0

julia> fx = ZeroProblem(f, x0)
ZeroProblem{typeof(f), Float64}(f, 2.0)

julia> solve(fx) ≈ 0.7728829591492101
true

При отсутствии значений по умолчанию и указании одной начальной точки используется метод Order1 по умолчанию. Метод solve позволяет передавать другие методы нахождения корней, а также другие параметры. Например, чтобы использовать метод Order2 с критерием сходимости (см. ниже), согласно которому |xₙ - xₙ₋₁| ≤ δ, можно выполнить следующий вызов:

julia> solve(fx, Order2(); atol=0.0, rtol=0.0) ≈ 0.7728829591492101
true

В отличие от функции find_zero, которая при несходимости выдает ошибку, метод solve при несходимости возвращает NaN.

В следующем примере нулевое значение находится в позиции 0.0, но для большинства начальных значений будет уходить в сторону ±∞, что иногда приводит к тому, что относительный допуск возвращает неправильное значение. Здесь мы можем увидеть различия:

julia> f(x) = cbrt(x) * exp(-x^2)
f (generic function with 1 method)

julia> x0 = 0.1147
0.1147

julia> find_zero(f, x0, Roots.Order1()) ≈ 5.075844588445206 # прекращается, поскольку |f(xₙ)| ≤ |xₙ|ϵ
true

julia> find_zero(f, x0, Roots.Order1(), atol=0.0, rtol=0.0) # ошибка, поскольку нет проверки в `|f(xn)|`
ERROR: Roots.ConvergenceFailed("Algorithm failed to converge")
[...]

julia> fx = ZeroProblem(f, x0);

julia> solve(fx, Roots.Order1(), atol=0.0, rtol=0.0) # NaN, не ошибка
NaN

julia> fx = ZeroProblem((f, D(f)), x0); # методы более высокого порядка позволяют найти нуль этой функции

julia> solve(fx, Roots.LithBoonkkampIJzerman(2,1), atol=0.0, rtol=0.0)
0.0

Функции могут быть параметризованы, как показано ниже:

julia> f(x, p=2) = cos(x) - x/p
f (generic function with 2 methods)

julia> Z = ZeroProblem(f, pi/4)
ZeroProblem{typeof(f), Float64}(f, 0.7853981633974483)

julia> solve(Z, Order1()) ≈ 1.0298665293222586     # использовать p=2 (по умолчанию)
true

julia> solve(Z, Order1(), p=3) ≈ 1.170120950002626 # использовать p=3
true

julia> solve(Z, Order1(), 4) ≈ 1.2523532340025887  # по позиции, использует p=4
true

Множество нулей

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

julia> f(x) = exp(x) - x^4
f (generic function with 2 methods)

julia> find_zeros(f, -10,10) ≈ [α₀, α₁, α₂] # из вышеприведенного
true

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

julia> using IntervalSets

julia> find_zeros(f, -10..10) ≈ [α₀, α₁, α₂]
true

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

Сходимость

Для большинства алгоритмов сходимость определяется, когда

  • значение |f(x_n)| <= tol при tol = max(atol, abs(x_n)*rtol), или

  • значения x_n ≈ x_{n-1} с допусками xatol и xrtol и f(x_n) ≈ 0 с нестрогим допуском на основе atol и rtol.

Алгоритм find_zero останавливается, если

  • обнаруживает NaN или Inf, или

  • количество итераций превышает maxevals

Если алгоритм останавливается и выполняется условие нестрогой сходимости, возвращается предполагаемый нуль. В противном случае возникает ошибка, указывающая на отсутствие сходимости. Для настройки допусков find_zero принимает именованные аргументы atol, rtol, xatol и xrtol, как показано в примерах выше.

Методы Bisection и Roots.A42 гарантированно сходятся, даже если заданы нулевые допуски, поэтому они используются по умолчанию. Для xatol и xrtol можно указать ненулевые значения, чтобы сократить количество вызовов функции, когда требуется более низкая точность.

julia> fx = ZeroProblem(sin, (3,4));

julia> solve(fx, Bisection(); xatol=1/16)
3.125

Альтернативный интерфейс

Эту функциональность обеспечивает функция fzero, знакомая пользователям MATLAB. Roots предоставляет этот альтернативный интерфейс:

  • fzero(f, x0::Real; order=0) вызывает метод без производных с порядком, указывающим один из Order0, Order1 и т. д.

  • fzero(f, a::Real, b::Real) вызывает алгоритм find_zero с методом Bisection.

  • fzeros(f, a::Real, b::Real) вызовет find_zeros.

Примеры использования

julia> f(x) = exp(x) - x^4
f (generic function with 2 methods)

julia> fzero(f, 8, 9) ≈ α₂   # заключение в скобки
true

julia> fzero(f, -10, 0) ≈ α₀
true

julia> fzeros(f, -10, 10) ≈ [α₀, α₁, α₂]
true

julia> fzero(f, 3) ≈ α₁      # по умолчанию используется Order0()
true

julia> fzero(sin, big(3), order=16)  ≈ π # использует метод более высокого порядка
true

Функции find_zero и find_zeros

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

find_zero(f, x0, M, [N::AbstractBracketingMethod], [p=nothing]; kwargs...)

Интерфейс одного из нескольких методов нахождения нулей одномерной функции, например решение .

Аргументы

Позиционные аргументы

  • f: функция (одномерная или f(x,p) с p, содержащим параметры)

  • x0: начальное условие (значение, начальные значения или заключаемый в скобки интервал)

  • M: некий метод AbstractUnivariateZeroMethod, указывающий решатель

  • N: некий метод заключения в скобки, при указании которого создается гибридный метод

  • p: для указания параметра функции f. Также может быть ключевым словом, но позиционный аргумент полезен при трансляции.

Именованные аргументы

  • xatol, xrtol: абсолютный и относительный допуски, чтобы решить, xₙ₊₁ ≈ xₙ

  • atol, rtol: абсолютный и относительный допуски, чтобы решить, f(xₙ) ≈ 0

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

  • verbose::Bool: указывает, следует ли показывать подробную информацию об алгоритме.

  • tracks: позволяет указать спецификацию объектов Tracks

Расширенная справка

Начальное значение

Для большинства методов x0 — это скалярное значение, указывающее начальное значение в итеративной процедуре. (Методы секущих могут иметь кортеж, задающий их начальные значения.) Значения должны быть подтипом Number и иметь методы, определенные для float, real и oneunit.

Для заключаемых в скобки интервалов x0 задается с помощью кортежа, вектора или любого итерируемого объекта с определенными экстремумами (extrema). Заключаемый в скобки интервал ] представляет собой интервал, где и имеют разные знаки.

Возвращаемое значение

Если алгоритм завершается успешно, возвращается найденный приблизительный корень. Если алгоритм завершается сбоем, возникает ошибка ConvergenceFailed . В случае сбоя альтернативная форма solve(ZeroProblem(f,x0), M) возвращает NaN .

Указание метода

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

  • Существуют методы, в которых указывается скобка: Bisection, A42, AlefeldPotraShi, Roots.Brent, помимо прочих. Метод деления пополам используется по умолчанию для базовых типов с плавающей запятой, но для A42 обычно требуется гораздо меньше итераций.

  • Существует несколько методов без производных: ср. Order0, Order1 (также Roots.Secant), Order2 (также Steffensen), Order5, Order8 и Order16, где число указывает порядок сходимости.

  • Есть несколько классических методов, в которых требуется спецификация производных: Roots.Newton, Roots.Halley, Roots.Schroder, помимо прочих.

  • Методы, предназначенные для решения задач с кратностями, включают Roots.Order1B, Roots.Order2B и Roots.ThukralXB для разных значений X.

  • Семейство Roots.LithBoonkkampIJzerman{S,D} для различных S и D использует линейный многошаговый метод поиска корней. Метод (2,0) представляет собой метод деления пополам, (1,1) является методом Ньютона.

Дополнительные сведения см. на странице справки по каждому методу (например, ?Order1). Неэкспортируемые методы должны быть снабжены именем модуля, как в ?Roots.Schroder.

Если метод не указан, используемый по умолчанию метод зависит от x0.

  • Если x0 является скаляром, по умолчанию используется более надежный метод Order0.

  • Если x0 является кортежем, вектором или итерируемым объектом с определенным экстремумом (extrema), , указывающим заключаемый в скобки интервал, то метод Bisection используется для типов Float64, Float32 или Float16. В противном случае используется метод A42.

По умолчанию выбираются надежные методы. Они могут быть не столь эффективны, как некоторые другие.

Указание функции

Функции передаются в качестве первых аргументов.

Для тех немногих методов, которые используют одну или несколько производных (Newton, Halley, Schroder, LithBoonkkampIJzerman(S,D) и т. д.), применяется кортеж функций. Для классических алгоритмов можно использовать функцию, возвращающую (f(x), f(x)/f'(x), [f'(x)/f''(x)]).

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

  • xatol — абсолютный допуск для значений x.

  • xrtol — относительный допуск для значений x.

  • atol — абсолютный допуск для значений f(x).

  • rtol — относительный допуск для значений f(x).

  • maxiters — ограничение на максимальное число итераций.

  • strict — если false (по умолчанию), то при остановке алгоритма возможные нули проверяются с нестрогим допуском.

  • verbose — если true, то после успешного завершения будет показана трассировка алгоритма. Сведения о сохранении этой трассировки см. в справке по внутреннему объекту Roots.Tracks.

См. раздел справки по Roots.assess_convergence со сведениями о сходимости См. страницу справки по Roots.default_tolerances(method) со сведениями о заданных по умолчанию допусках.

В целом при работе с числами с плавающей запятой сходимость не следует понимать как безусловное утверждение. Даже если математически α является ответом, а xstar — реализацией плавающей запятой, может оказаться, что f(xstar) - f(α) ≈ xstar ⋅ f'(α) ⋅ eps(α), поэтому необходимо принимать во внимание допуски, а иногда и указывать их.

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

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

Примечание. На порядок выполнения метода указывает схема именования. Схема имеет порядок r, если при eᵢ = xᵢ - α, eᵢ₊₁ = C⋅eᵢʳ. Если погрешность eᵢ достаточно мала, то, по сути, на каждом шаге она будет набирать в r раз больше начальных нулей. Однако если погрешность велика, этого не произойдет. Без хороших начальных предположений метод высокого порядка может сходиться медленно, если вообще будет сходиться. Методы OrderN имеют ряд эвристических правил, используемых для получения более широкого диапазона сходимости ценой не совсем точной реализации метода, хотя они доступны через неэкспортированные методы.

Примеры:

Методы по умолчанию.

julia> using Roots

julia> find_zero(sin, 3)  # использует Order0()
3.141592653589793

julia> find_zero(sin, (3,4)) # использует Bisection()
3.141592653589793

При указании метода

julia> find_zero(sin, (3,4), Order1())            # можно задать две начальные точки для метода секущей
3.141592653589793

julia> find_zero(sin, 3.0, Order2())              # Использует метод Стеффенсена
3.1415926535897936

julia> find_zero(sin, big(3.0), Order16())        # быстрая сходимость
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

julia> find_zero(sin, (3, 4), A42())              # меньше вызовов функции, чем при использовании метода Bisection(), в данном случае
3.141592653589793

julia> find_zero(sin, (3, 4), FalsePosition(8))   # 1 из 12 возможных алгоритмов ложного положения
3.141592653589793

julia> find_zero((sin,cos), 3.0, Roots.Newton())  # использует метод Ньютона
3.141592653589793

julia> find_zero((sin, cos, x->-sin(x)), 3.0, Roots.Halley())  # использует метод Галлея
3.141592653589793

Изменение допусков.

julia> fn = x -> (2x*cos(x) + x^2 - 3)^10/(x^2 + 1);

julia> x0, xstar = 3.0,  2.9947567209477;

julia> fn(find_zero(fn, x0, Order2())) <= 1e-14  # f(xₙ) ≈ 0, но Δxₙ может быть довольно большим
true

julia> find_zero(fn, x0, Order2(), atol=0.0, rtol=0.0) # ошибка: x_n ≉ x_{n-1}; только f(x_n) ≈ 0
ERROR: Roots.ConvergenceFailed("Algorithm failed to converge")
[...]

julia> fn = x -> (sin(x)*cos(x) - x^3 + 1)^9;

julia> x0, xstar = 1.0,  1.112243913023029;

julia> isapprox(find_zero(fn, x0, Order2()), xstar; atol=1e-4)
true

julia> find_zero(fn, x0, Order2(), maxiters=3)    # необходимо больше шагов для сходимости
ERROR: Roots.ConvergenceFailed("Algorithm failed to converge")
[...]

Трассировка

При передаче verbose=true будет показана подробная информация о шагах алгоритма. Алгоритм tracks позволяет передать объект Roots.Tracks для записи значений x и f(x), используемых в алгоритме.

Сведения об альтернативном интерфейсе см. в описании solve! и ZeroProblem.

 
find_zeros(f, a, [b]; [no_pts=12, k=8, naive=false, xatol, xrtol, atol, rtol])

Выполняет поиск нулей f в интервале [a,b] с помощью эвристического алгоритма.

  • f: функция или вызываемый объект

  • a, b: если указан b, используется интервал ]. Если указан только a, он передается в extrema, чтобы определить интервал, в котором следует выполнить поиск. Предполагается, что ни одна из конечных точек не равна нулю.

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

Расширенная справка

Примеры

julia> using Roots

julia> find_zeros(x -> exp(x) - x^4, -5, 20)        # несколько хорошо расположенных нулей
3-element Vector{Float64}:
 -0.8155534188089606
  1.4296118247255556
  8.613169456441398

julia> find_zeros(x -> sin(x^2) + cos(x)^2, 0, 2pi)  # много нулей
12-element Vector{Float64}:
 1.78518032659534
 2.391345462376604
 3.2852368649448853
 3.3625557095737544
 4.016412952618305
 4.325091924521049
 4.68952781386834
 5.00494459113514
 5.35145266881871
 5.552319796014526
 5.974560835055425
 6.039177477770888

julia> find_zeros(x -> cos(x) + cos(2x), (0, 4pi))    # сочетание простых и непростых нулей
6-element Vector{Float64}:
  1.0471975511965976
  3.141592653589943
  5.235987755982988
  7.330382858376184
  9.424777960769228
 11.519173063162574

julia> f(x) = (x-0.5) * (x-0.5001) * (x-1)          # близлежащие нули
f (generic function with 1 method)

julia> find_zeros(f, 0, 2)
3-element Vector{Float64}:
 0.5
 0.5001
 1.0

julia> f(x) = (x-0.5) * (x-0.5001) * (x-4) * (x-4.001) * (x-4.2)
f (generic function with 1 method)

julia> find_zeros(f, 0, 10)
3-element Vector{Float64}:
 0.5
 0.5001
 4.2

julia> f(x) = (x-0.5)^2 * (x-0.5001)^3 * (x-4) * (x-4.001) * (x-4.2)^2  # трудно найти
f (generic function with 1 method)

julia> find_zeros(f, 0, 10, no_pts=21)                # слишком сложно для заданных по умолчанию
5-element Vector{Float64}:
 0.49999999999999994
 0.5001
 4.0
 4.001
 4.200000000000001

В некоторых случаях количество нулей может быть занижено:

  • если начальный интервал (a,b) слишком широк

  • если нули расположены очень близко друг к другу

  • функция является плоской, например x→0.

Основной алгоритм проверяет наличие нулей в конечных точках, затем разбивает интервал (a,b) на подынтервалы no_pts-1, а затем переходит к поиску нулей с помощью метода деления пополам или метода без производной. Поскольку проверка на наличие заключаемого в скобки интервала относительно дешева, а метод деления пополам гарантированно сходится, в каждом интервале есть k пар промежуточных точек, проверяемых на наличие скобки.

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

Допуски используются для уменьшения интервалов, но не для поиска нулей в диапазоне поиска. При поисках метод деления пополам гарантированно сходится без указанного допуска. Для поиска без производной используется модификация метода Order0, который в худшем случае сравнивает |f(x)| <= 8*eps(x) для нахождения нуля. Алгоритм может определить несколько значений для нуля из-за аппроксимаций с плавающей запятой. Если потенциальная пара нулей удовлетворяет isapprox(a,b,atol=sqrt(xatol), rtol=sqrt(xrtol)), они объединяются.

Алгоритм может выполнять множество вызовов функции. При обнаружении нулей в интервале примитивный поиск выполняется в каждом подынтервале. Чтобы сократить количество вызовов функции, но с повышенной вероятностью пропуска нескольких нулей, адаптивный характер можно отменить, задав аргумент naive=true или уменьшив количество точек.

Алгоритм заимствован из PR у @djsegal..

Пакет IntervalRootFinding предоставляет серьезную альтернативу этой эвристике. Этот пакет использует интервальную арифметику, поэтому может вычислять границы размера изображения интервала для функции f. Если это изображение содержит 0, он может искать нуль. Метод деления пополам, с другой стороны, будет искать нуль только в том случае, если две конечные точки имеют разные знаки, что является гораздо более жестким условием для потенциального нуля.

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

Например, эта функция (из-за @truculentmath) особенно сложна, поскольку она положительна при каждом числе с плавающей запятой, но имеет два нуля (вертикальная асимптота при 15//11 отрицательна только в смежных значениях с плавающей запятой):

julia> using IntervalArithmetic, IntervalRootFinding, Roots

julia> g(x) = x^2 + 1 +log(abs( 11*x-15 ))/99
g (generic function with 1 method)

julia> find_zeros(g, -3, 3)
Float64[]

julia> IntervalRootFinding.roots(g, -3..3, IntervalRootFinding.Bisection)
1-element Vector{Root{Interval{Float64}}}:
 Root([1.36363, 1.36364], :unknown)

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

julia> g(x) = exp(x) - x^5
g (generic function with 1 method)

julia> rts = IntervalRootFinding.roots(g, -20..20)
2-element Vector{Root{Interval{Float64}}}:
 Root([12.7132, 12.7133], :unique)
 Root([1.29585, 1.29586], :unique)

julia> find_zeros(g, -20, 20)
2-element Vector{Float64}:
  1.2958555090953687
 12.713206788867632

Интерфейс CommonSolve

Интерфейс «задача-алгоритм-решение» — это популярный шаблон в Julia, реализуемый набором пакетов DifferentialEquations.jl. Его можно использовать в качестве альтернативы функции find_zero. В отличие от find_zero, solve при несходимости возвращает NaN.

solve!(P::ZeroProblemIterator)
solve(fx::ZeroProblem, [M], [N]; p=nothing, kwargs...)
init(fx::ZeroProblem, [M], [N];
     p=nothing,
     verbose=false, tracks=NullTracks(), kwargs...)

Solve for the zero of a scalar-valued univariate function specified through ZeroProblem or ZeroProblemIterator using the CommonSolve interface.

The defaults for M and N depend on the ZeroProblem: if x0 is a number, then M=Secant() and N=AlefeldPotraShi() is used (Order0); if x0 has 2 (or more values) then it is assumed to be a bracketing interval and M=AlefeldPotraShi() is used.

The methods involved with this interface are:

  • ZeroProblem: used to specify a problem with a function (or functions) and an initial guess

  • solve: to solve for a zero in a ZeroProblem

The latter calls the following, which can be useful independently:

  • init: to initialize an iterator with a method for solution, any adjustments to the default tolerances, and a specification to log the steps or not.

  • solve! to iterate to convergence.

Returns NaN, not an error like find_zero, when the problem can not be solved. Tested for zero allocations.

Examples:

julia> using Roots

julia> fx = ZeroProblem(sin, 3)
ZeroProblem{typeof(sin), Int64}(sin, 3)

julia> solve(fx)
3.141592653589793

Or, if the iterable is required

julia> problem = init(fx);

julia> solve!(problem)
3.141592653589793

keyword arguments can be used to adjust the default tolerances.

julia> solve(fx, Order5(); atol=1/100)
3.1425464815525403

The above is equivalent to:

julia> problem = init(fx, Order5(), atol=1/100);

julia> solve!(problem)
3.1425464815525403

The keyword argument p may be used if the function(s) to be solved depend on a parameter in their second positional argument (e.g., f(x, p)). For example

julia> f(x,p) = exp(-x) - p # to solve p = exp(-x)
f (generic function with 1 method)

julia> fx = ZeroProblem(f, 1)
ZeroProblem{typeof(f), Int64}(f, 1)

julia> solve(fx; p=1/2)  # log(2)
0.6931471805599453

This would be recommended, as there is no recompilation due to the function changing. For use with broadcasting, p may also be the last positional argument.

The argument verbose=true for init instructs that steps to be logged;

The iterator interface allows for the creation of hybrid solutions, such as is used when two methods are passed to solve. For example, this is essentially how the hybrid default is constructed:

julia> function order0(f, x)
           fx = ZeroProblem(f, x)
           p = init(fx, Roots.Secant())
           xᵢ,st = ϕ = iterate(p)
           while ϕ !== nothing
               xᵢ, st = ϕ
               state, ctr = st
               fᵢ₋₁, fᵢ = state.fxn0, state.fxn1
               if sign(fᵢ₋₁)*sign(fᵢ) < 0 # check for bracket
                   x0 = (state.xn0, state.xn1)
                   fx′ = ZeroProblem(f, x0)
                   p = init(fx′, Bisection())
                   xᵢ = solve!(p)
                   break
               end
               ϕ = iterate(p, st)
           end
           xᵢ
       end
order0 (generic function with 1 method)

julia> order0(sin, 3)
3.141592653589793
solve(fx::ZeroProblem, args...; verbose=false, kwargs...)

Выполняет диспетчеризацию в solve!(init(fx, args...; kwargs...)). Дополнительные сведения см. в описании solve!.

ZeroProblem{F,X}

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

Классические методы на основе производных

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

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

Метод Ньютона описывается легко:

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

Для понимания методов, доступных в Roots, полезно знать ряд фактов:

  • Для метода Ньютона существует формула погрешности: при условии, что , где  — нуль, где  — некоторое значение между и .

  • Остаточный член, заданный в виде \|\epsilon_{n+1}\| \leq C\cdot\|\epsilon_n\|^2, можно использовать для определения интервала вокруг , для которого гарантируется сходимость. Такая сходимость называется квадратичной (порядок 2). Для решений с плавающей запятой квадратичная сходимость и правильно выбранная начальная точка могут приводить к сходимости за 4 или 5 итераций. В общем случае сходимость называется имеющей порядок , когда

  • Член указывает на наличие возможных проблем, когда производная слишком велика в окрестности или производная слишком мала в окрестности . В частности, если , квадратичной сходимости может не быть и для достижения сходимости может потребоваться много итераций. Нуль, для которого при , называется простым при и непростым при . Метод Ньютона является квадратичным в окрестности простых нулей, но необязательно является квадратичным в окрестности непростых нулей. Аналогичным образом, если производная слишком велика в окрестности или производная слишком мала в окрестности либо если находится слишком далеко от (то есть ), погрешность может увеличиваться и сходимость не гарантируется.

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

  • Количество вычислений функции за один шаг для метода Ньютона составляет 2.

Roots.Newton()

Реализует метод Ньютона: xᵢ₊₁ = xᵢ - f(xᵢ)/f'(xᵢ). Это квадратично сходящийся метод, требующий одну производную и два вызова функции для каждого шага.

Примеры

julia> using Roots

julia> find_zero((sin,cos), 3.0, Roots.Newton()) ≈ π
true

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

julia> find_zero(x -> (sin(x), sin(x)/cos(x)), 3.0, Roots.Newton()) ≈ π
true

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

Погрешность eᵢ = xᵢ - α может быть выражена в виде eᵢ₊₁ = f[xᵢ,xᵢ,α]/(2f[xᵢ,xᵢ])eᵢ² (методы Сиди, унифицированной обработки ложного положения, Ньютона-Рафсона, секущей и Стеффенсена для нелинейных уравнений).

Roots.Halley()

Реализует метод Галлея `xᵢ₊₁ = xᵢ

  • (f/f')(xᵢ) * (1 - (f/f')(xᵢ) * (f''/f')(xᵢ) * 1/2){caret}(-1)` Этот метод

сходится кубически, но требует вызова функции на шаг. Метод Галлея находит xₙ₊₁ как нуль гиперболы в точке (xₙ, f(xₙ)), совпадающей с первой и второй производными f.

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

Примеры

julia> using Roots

julia> find_zero((sin, cos, x->-sin(x)), 3.0, Roots.Halley()) ≈ π
true

julia> function f(x)
       s,c = sincos(x)
       (s, s/c, -c/s)
       end;

julia> find_zero(f, 3.0, Roots.Halley()) ≈ π
true

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

Погрешность eᵢ = xᵢ - α удовлетворяет условию eᵢ₊₁ ≈ -(2f'⋅f''' -3⋅(f'')²)/(12⋅(f'')²) ⋅ eᵢ³ (все вычисляется по α).

Roots.QuadraticInverse()

Реализует квадратичный обратный метод, также известный как метод Чебышева, xᵢ₊₁ = xᵢ - (f/f')(xᵢ) * (1 + (f/f')(xᵢ) * (f''/f')(xᵢ) * 1/2). Этот метод сходится кубически и требует вызова функции на каждый шаг.

Пример

julia> using Roots

julia> find_zero((sin, cos, x->-sin(x)), 3.0, Roots.QuadraticInverse()) ≈ π
true

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

julia> find_zero(x -> (sin(x), sin(x)/cos(x), -cos(x)/sin(x)), 3.0, Roots.QuadraticInverse()) ≈ π
true

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

Погрешность eᵢ = xᵢ - α удовлетворяет условию eᵢ₊₁ ≈ (1/2⋅(f''/f')² - 1/6⋅f'''/f')) ⋅ eᵢ³ (все вычисляется по α).

CHEBYSHEV-LIKE METHODS AND QUADRATIC EQUATIONS (J. A. EZQUERRO, J. M. GUTIÉRREZ, M. A. HERNÁNDEZ and M. A. SALANOVA)

An acceleration of Newton’s method: Super-Halley method (J.M. Gutierrez, M.A. Hernandez)

Методы Ньютона и Галлея входят в это семейство методов:

LithBoonkkampIJzerman{S,D} <: AbstractNewtonLikeMethod
LithBoonkkampIJzerman(S,D)

Семейство различных методов, включающее метод секущей и метод Ньютона.

Определяет линейный многошаговый решатель с S шагами и D производными, как описано в статье Лит (Lith), Бункампа (Boonkkamp) и Айзермана (IJzerman) (https://doi.org/10.1016/j.amc.2017.09.003).

Расширенная справка

Примеры

julia> using Roots

julia> find_zero(sin, 3, Roots.LithBoonkkampIJzerman(2,0)) ≈ π # метод секущей
true

julia> find_zero((sin,cos), 3, Roots.LithBoonkkampIJzerman(1,1)) ≈ π # метод Ньютона
true

julia> find_zero((sin,cos), 3, Roots.LithBoonkkampIJzerman(3,1)) ≈ π # Быстрая скорость сходимости
true

julia> find_zero((sin,cos, x->-sin(x)), 3, Roots.LithBoonkkampIJzerman(1,2)) ≈ π # метод в стиле Галлея
true

Метод может быть более устойчивым к начальным условиям. Это пример из статьи (стр. 13). Метод Ньютона (случай S=1, D=1) завершается сбоем, если |x₀| ≥ 1.089, но методы с большим объемом памяти работают успешно.

julia> fx =  ZeroProblem((tanh,x->sech(x)^2), 1.239); # нуль в 0.0

julia> solve(fx, Roots.LithBoonkkampIJzerman(1,1)) |> isnan# Ньютон, NaN
true

julia> solve(fx, Roots.LithBoonkkampIJzerman(2,1)) |> abs |> <(eps())
true

julia> solve(fx, Roots.LithBoonkkampIJzerman(3,1)) |> abs |> <(eps())
true

Множество производных можно создать автоматически с помощью автоматического дифференцирования. Например,

julia> using ForwardDiff

julia> function δ(f, n::Int=1)
           n <= 0 && return f
           n == 1 && return x -> ForwardDiff.derivative(f,float(x))
           δ(δ(f,1),n-1)
       end;

julia> fs(f,n) = ntuple(i -> δ(f,i-1), Val(n+1));

julia> f(x) = cbrt(x) * exp(-x^2); # ср. табл. 6 в документе, α = 0

julia> fx = ZeroProblem(fs(f,1), 0.1147);

julia> opts = (xatol=2eps(), xrtol=0.0, atol=0.0, rtol=0.0); # сходимость при |xₙ - xₙ₋₁| <= 2ϵ

julia> solve(fx, Roots.LithBoonkkampIJzerman(1, 1); opts...) |> isnan # NaN — без сходимости
true

julia> solve(fx, Roots.LithBoonkkampIJzerman(2, 1); opts...) |> abs |> <(eps()) # сходится
true

julia> fx = ZeroProblem(fs(f,2), 0.06);                       # требуется лучшая начальная точка

julia> solve(fx, Roots.LithBoonkkampIJzerman(2, 2); opts...) |> abs |> <(eps()) # сходится
true

Для случая D=1 метод заключения в скобки, основанный на этом подходе, реализован в LithBoonkkampIJzermanBracket.

Справка

В статье Лита, Бункампа и Айзермана (https://doi.org/10.1016/j.amc.2017.09.003) приведен анализ скорости сходимости при использовании линейных многошаговых методов для решения 0=f(x) с помощью f⁻¹(0) = x, когда f является достаточно сглаженной линейной функцией. Переформулировка, приписываемая Грау-Санчесу, находит дифференциальное уравнение для f⁻¹: dx/dy = [f⁻¹]′(y) = 1/f′(x) = F в виде x(0) = x₀ + ∫⁰_y₀ F(x(y)) dy.

Для численного решения этого уравнения используется линейный многошаговый метод. Пусть S — количество шагов памяти (S= 1,2,…​), а D — количество используемых производных, тогда при F(x) = dx/dy x_{n+S} = ∑_{k=0}^{S-1} aₖ x_{n+k} +∑d=1^D ∑_{k=1}^{S-1} aᵈ_{n+k}F⁽ᵈ⁾(x_{n+k}). Значения aₖ и aᵈₖ вычисляются на каждом шаге.

Эта таблица взята из таблиц 1 и 3 статьи, и в ней показана скорость сходимости для простых корней, указанных в ней:

s: number of steps remembered
d: number of derivatives uses
s/d  0    1    2    3    4
1    .    2    3    4    5
2    1.62 2.73 3.79 4.82 5.85
3    1.84 2.91 3.95 4.97 5.98
4    1.92 2.97 3.99 4.99 5.996
5    1.97 .    .    .    .

То есть больший объем памяти приводит к более высокой скорости сходимости. Большее количество производных приводит к более высокой скорости сходимости. Однако интервал около α, нуль, где скорость сходимости гарантирована, может уменьшиться.

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

Методы без вычисления производных

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

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

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

Secant()
Order1()
Orderφ()

Метод Order1() является псевдонимом для Secant. Он задает метод секущей. Этот метод хранит в своем состоянии два значения — xₙ и xₙ₋₁. Обновленная точка является точкой пересечения оси с секущей, сформированной из двух точек. Метод секущей использует вычисление функции на каждом шаге и имеет порядок φ≈ (1+sqrt(5))/2.

Погрешность eᵢ = xᵢ - α удовлетворяет условию eᵢ₊₂ = f[xᵢ₊₁,xᵢ,α] / f[xᵢ₊₁,xᵢ] * (xᵢ₊₁-α) * (xᵢ - α).

Secant()
Order1()
Orderφ()

Метод Order1() является псевдонимом для Secant. Он задает метод секущей. Этот метод хранит в своем состоянии два значения — xₙ и xₙ₋₁. Обновленная точка является точкой пересечения оси с секущей, сформированной из двух точек. Метод секущей использует вычисление функции на каждом шаге и имеет порядок φ≈ (1+sqrt(5))/2.

Погрешность eᵢ = xᵢ - α удовлетворяет условию eᵢ₊₂ = f[xᵢ₊₁,xᵢ,α] / f[xᵢ₊₁,xᵢ] * (xᵢ₊₁-α) * (xᵢ - α).

Steffensen()

Квадратично сходящийся метод Стеффенсена используется для алгоритма Order2() без производных. В отличие от квадратично сходящегося метода Ньютона, производные не требуются, но, как и в методе Ньютона, на каждый шаг приходится два вызова функции. Алгоритм Стеффенсена более чувствителен, чем метод Ньютона, к плохим начальным предположениям, когда f(x) велика из-за того, как f'(x) аппроксимируется. Метод Order2 заменяет шаг Стеффенсена шагом секущей, когда значение f(x) велико.

Погрешность eᵢ - α удовлетворяет условию eᵢ₊₁ = f[xᵢ, xᵢ+fᵢ, α] / f[xᵢ,xᵢ+fᵢ] ⋅ (1 - f[xᵢ,α] ⋅ eᵢ²

Order2

Steffensen с защитой на шаге секущей.

Order5()
KumarSinghAkanksha()

Реализует алгоритм 5-го порядка из статьи A New Fifth Order Derivative Free Newton-Type Method for Solving Nonlinear Equations от Manoj Kumar, Akhilesh Kumar Singh и Akanksha, Appl. Math. Inf. Sci. 9, No. 3, 1507-1513 (2015), DOI: 10.12785/amis/090346. Для каждого шага требуется пять вызовов функции. Метод Order5 заменяет шаг Стеффенсена шагом секущей, когда значение f(x) велико.

Погрешность eᵢ = xᵢ - α удовлетворяет условию eᵢ₊₁ = K₁ ⋅ K₅ ⋅ M ⋅ eᵢ⁵ + O(eᵢ⁶)

Order8()
Thukral8()

Реализует алгоритм 8-го порядка из статьи New Eighth-Order Derivative-Free Methods for Solving Nonlinear Equations от Rajinder Thukral, International Journal of Mathematics and Mathematical Sciences Volume 2012 (2012), Article ID 493456, 12 pages DOI: 10.1155/2012/493456. Для каждого шага требуется четыре вызова функции. Метод Order8 заменяет шаг Стеффенсена шагом секущей, когда значение f(x) велико.

Погрешность eᵢ = xᵢ - α выражается как eᵢ₊₁ = K ⋅ eᵢ⁸ в (2.25) статьи для явного K.

Order16()
Thukral16()

Реализует алгоритм 16-го порядка из статьи New Sixteenth-Order Derivative-Free Methods for Solving Nonlinear Equations от R. Thukral, American Journal of Computational and Applied Mathematics p-ISSN: 2165-8935; e-ISSN: 2165-8943; 2012; 2(3): 112-118 DOI: 10.5923/j.ajcam.20120203.08.

Для каждого шага требуется пять вызовов функции. Хотя этот метод быстро сходится, в целом он не быстрее (меньше вызовов функций/шагов) других методов при использовании значений Float64, но может быть полезен для решения задач с BigFloat. Метод Order16 заменяет шаг Стеффенсена шагом секущей, когда значение f(x) велико.

Погрешность eᵢ = xᵢ - α выражается как eᵢ₊₁ = K⋅eᵢ¹⁶ для явного K в уравнении (50) из статьи.

Ограничивающие методы

Методом деления пополам находится нуль непрерывной функции между и , когда и имеют разные знаки. (Интервал ] называется ограничивающим, когда .) Базовый алгоритм очень простой: интервал ] делится по . Либо , либо один из интервалов ] или ] является ограничивающим. Он обозначается ]. Из этого описания видно, что ] имеет длину в раз большую, чем длина ], поэтому в конечном итоге будет найден нуль или либо интервалы сойдутся к нулю. Данная сходимость медленная (эффективность составляет всего ), но гарантированная. Для 16-, 32- и 64-битных значений с плавающей запятой переинтерпретация способа нахождения средней точки ( ) обеспечивает сходимость не более чем за итераций, в отличие от описанного выше способа, который в некоторых случаях может требовать гораздо большего числа шагов для сходимости.

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

При описанном выше способе выбора средней точки по умолчанию учитывается только знак функции . Алгоритмы, учитывающие форму функции, могут быть значительно эффективнее. Например, ограничивающий метод Roots.AlefeldPotraShi по Алефельду, Потра и Ши имеет эффективность . Этот метод также применяется по умолчанию в функции find_zero, когда дана одна начальная точка и определен ограничивающий интервал.

Bisection()

Если возможно, используйте метод деления пополам для значений Float64. Метод заключения в скобки начинается с заключаемого в скобки интервала [a,b] и разбивает его на два интервала [a,c] и [c,b]. Если c не является нулем, то один из этих двух интервалов будет заключен в скобки и процесс продолжится. c вычисляется с помощью _middle, которая по-иному интерпретирует значения с плавающей запятой в целые числа без знака и выполняет разбиение. Его автором является Джейсон Меррилл (Jason Merrill). Этот метод позволяет избежать проблем с плавающей запятой и при нулевом значении допусков (по умолчанию) гарантирует «лучшее» решение (такое, в котором найден нуль или интервал скобок имеет тип [a, nextfloat(a)]).

Если заданы допуски, алгоритм завершается, когда длина интервала меньше или равна допуску max(δₐ, 2abs(u)δᵣ), причем u в {a,b} выбирается меньшим из |f(a)| и |f(b)|, или значение функции меньше max(tol, min(abs(a), abs(b)) * rtol). Последний вариант используется только в том случае, если настроены допуски по умолчанию (atol или rtol).

При решении для с помощью Bisection нельзя брать производную непосредственно через автоматическое дифференцирование, так как алгоритм не является дифференцируемым. Альтернативные варианты см. в разделе Чувствительность документации.

Roots.A42()

Метод заключения в скобки, который позволяет найти корень непрерывной функции в заданном заключаемом в скобки интервале [a, b], не требуя производных. Он основан на алгоритме 4.2, описанном в документе: G. E. Alefeld, F. A. Potra, and Y. Shi, "Algorithm 748: enclosing zeros of continuous functions," ACM Trans. Math. Softw. 21, 327—​344 (1995), DOI: 10.1145/210089.210111. Индекс асимптотической эффективности, , равен .

Автор: Джон Треверс (John Travers).

В статье, приведенной выше, показано, что для непрерывно дифференцируемой в ] с простым корнем алгоритм завершается на нуле, или асимптотически шаги имеют обратный кубический вид (Лемма 5.1) Это доказано в предположении, что четырежды непрерывно дифференцируема.

 
Roots.AlefeldPotraShi()

Следует алгоритму 4.1 в статье ON ENCLOSING SIMPLE ROOTS OF NONLINEAR EQUATIONS от Alefeld, Potra, Shi; DOI: 10.1090/S0025-5718-1993-1192965-2.

Порядок сходимости — 2 + √5. Асимптотически на каждом шаге выполняется 3 вычисления функции. Индекс асимптотической эффективности равен . Менее эффективен, но может работать быстрее, чем связанный метод A42.

Автор: Джон Треверс (John Travers).

Roots.Brent()

Реализация метода Брента (или Брента-Деккера). Этот метод использует обратную квадратичную интерполяцию или шаг секущей, при необходимости возвращаясь к шагу поиска делением пополам.

Roots.Chandrapatla()

Использует алгоритм Чандрапатлы (ср. метод Шеррера) для решения .

Алгоритм Чандрапатлы выбирает обратный квадратичный шаг или шаг поиска делением пополам на основе вычисленного неравенства.

Roots.Ridders()

Реализует метод Риддерса. Этот метод заключения в скобки находит среднюю точку x₁, затем интерполирует экспоненту, а затем использует ложное положение с интерполированным значением для нахождения c. Используется, если c и x₁ образуют скобку, в противном случае применяется подынтервал [a,c] или [c,b].

Пример:

julia> using Roots

julia> find_zero(x -> exp(x) - x^4, (5, 15), Roots.Ridders()) ≈ 8.61316945644
true

julia> find_zero(x -> x*exp(x) - 10, (-100, 100), Roots.Ridders()) ≈ 1.74552800274
true

julia> find_zero(x -> tan(x)^tan(x) - 1e3, (0, 1.5), Roots.Ridders()) ≈ 1.3547104419
true

Риддерс (Ridders) показал, что погрешность eₙ₊₁ ≈ 1/2 eₙeₙ₋₁eₙ₋₂ ⋅ (g^2-2fh)/f удовлетворяет условию для f=F', g=F''/2, h=F'''/6, что предполагает сходимость со скоростью ≈ 1.839.... Использует четыре вычисления функции на каждом шаге, поэтому порядок сходимости составляет ≈ 1.225....

Roots.ITP(;[κ₁-0.2, κ₂=2, n₀=1])

Использует метод заключения в скобки ITP . Этот метод утверждает, что он «является первым алгоритмом поиска корней, который достигает суперлинейной сходимости метода секущей, сохраняя при этом оптимальную производительность метода деления пополам в наихудшем случае».

Значения κ1, κ₂ и n₀ являются параметрами настройки.

  1. Предлагаемое значение для κ₁ составляет 0.2/(b-a), но по умолчанию здесь используется 0.2. Значение κ₂ равно 2, а значением по умолчанию n₀ является 1.

Примечание.

Предложено на форуме discourse пользователем @TheLateKronos, который предоставил исходную версию кода.

FalsePosition([galadino_factor])

Использует метод ложного положения для нахождения нуля для функции f в заключаемом в скобки интервале [a,b].

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

Чтобы ускорить сходимость для вогнутых функций, этот алгоритм реализует понижающих коэффициентов Галдино (A family of regula falsi root-finding methods). Они задаются числом, как в FalsePosition(2), или одним из трех имен FalsePosition(:pegasus), FalsePosition(:illinois) или FalsePosition(:anderson_bjork) (по умолчанию). Вариант по умолчанию, как правило, имеет лучшую производительность, чем остальные, хотя бывают и исключения.

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

Примеры

find_zero(x -> x^5 - x - 1, (-2, 2), FalsePosition())
LithBoonkkampIJzermanBracket()

Метод заключения в скобки, являющийся модификацией метода Брента, согласно статье Лита, Бункампа и Айзермана (IJzerman). Наилучшая возможная скорость сходимости составляет 2,91.

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

Состояние включает 3 точки — скобку [a,b] (в b=xₙ f(b) ближе всего к 0) и c=xₙ₋₁ — и соответствующие значения функции и ее производной в этих трех точках.

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

BracketedHalley

Для скобки [a,b] использует метод Roots.Halley, начиная со значения x, для которого fa или fb ближе всего к 0. Использует платформу Roots.AbstractAlefeldPotraShi для применения метода заключения в скобки, каждый раз выполняя дополнительный шаг двойной секущей.

BracketedChebyshev

Для скобки [a,b] использует метод Roots.QuadraticInverse, начиная со значения x, для которого fa или fb ближе всего к 0. Использует платформу Roots.AbstractAlefeldPotraShi для применения метода заключения в скобки, каждый раз выполняя дополнительный шаг двойной секущей.

BracketedSchroder

Для скобки [a,b] использует метод Roots.Schroder, начиная со значения x, для которого fa или fb ближе всего к 0. Использует платформу Roots.AbstractAlefeldPotraShi для применения метода заключения в скобки, каждый раз выполняя дополнительный шаг двойной секущей.

Непростые нули

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

Для методов без вычисления производных для непростых нулей реализованы следующие типы:

Roots.King()

Суперлинейная (порядка 1.6...) модификация метода секущей для нескольких корней. Представлено в статье A SECANT METHOD FOR MULTIPLE ROOTS, автор RICHARD F. KING, BIT 17 (1977), 321-328

Основная идея аналогична методу Шредера: применение метода секущей к f/f'. Однако здесь используется f' ~ fp = (fx - f(x-fx))/fx (шаг Стеффенсена). В этой реализации Order1B, когда значение fx слишком велико, используется шаг одной секущей f .

Асимптотическая погрешность eᵢ = xᵢ - α задается следующим образом: eᵢ₊₂ = 1/2⋅G''/G'⋅ eᵢ⋅eᵢ₊₁ + (1/6⋅G'''/G' - (1/2⋅G''/G'))^2⋅eᵢ⋅eᵢ₊₁⋅(eᵢ+eᵢ₊₁).

Roots.Order1B()

Метод King с защищенным шагом секущей.

Roots.Esser()

Метод Эссера. Это квадратично сходящийся метод, который, как и метод Шредера, не зависит от кратности нуля. Метод Шердера имеет обновленный шаг x - r2/(r2-r1) * r1, где ri = fⁱ⁻¹/fⁱ. Esser approximates Эссер считает, что f' ~ f[x-h, x+h], f'' ~ f[x-h,x,x+h], where , где h = fx, как и в методе Стеффенсена, требуется 3 вызова функции на каждый шаг. Реализация Order2B использует шаг секущей, когда значение |fx| считается слишком большим.

Esser, H. Computing (1975) 14: 367. DOI: 10.1007/BF02253547 Eine stets quadratisch konvergente Modification des Steffensen-Verfahrens

Примеры

f(x) = cos(x) - x
g(x) = f(x)^2
x0 = pi/4
find_zero(f, x0, Order2(), verbose=true)        #  3 шага / 7 вызовов функции
find_zero(f, x0, Roots.Order2B(), verbose=true) #  4 / 9
find_zero(g, x0, Order2(), verbose=true)        #  22 / 45
find_zero(g, x0, Roots.Order2B(), verbose=true) #  4 / 10
Roots.Order2B()

Метод Esser с защищенным шагом секущей.

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

Roots.Schroder()

Метод Шредера, как и метод Галлея, использует f, f' и f''. В отличие от метода Галлея, он сходится квадратично, но не зависит от кратности нуля (ср. Schröder, E. Über unendlich viele Algorithmen zur Auflösung der Gleichungen. Math. Ann. 2, 317-365, 1870; mathworld).

Метод Шредера применяет метод Ньютона к f/f' — функции со всеми простыми нулями.

Пример

m = 2
f(x) = (cos(x)-x)^m
fp(x) = (-x + cos(x))*(-2*sin(x) - 2)
fpp(x) = 2*((x - cos(x))*cos(x) + (sin(x) + 1)^2)
find_zero((f, fp, fpp), pi/4, Roots.Halley())     # 14 шагов
find_zero((f, fp, fpp), 1.0, Roots.Schroder())    # 3 шага

(В то время как при m=1 в методе Галлея выполняется 2 шага, а в методе Шредера — 3.)

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

find_zero(x -> (sin(x), sin(x)/cos(x), -cos(x)/sin(x)), 3.0, Roots.Schroder())

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

Погрешность eᵢ = xᵢ - α такая же, как в методе Newton, с f замененной f/f'.

Семейство методов для непростых нулей, которые требуют производных порядка и сводятся к методу Шредера при , реализовано в следующем типе:

AbstractThukralBMethod

Абстрактный тип для методов ThukralXB для X, являющихся 2, 3, 4 или 5.

Это семейство методов, которые

  • эффективны (порядок X) для непростых корней (например, Thukral2B — это метод Schroder)

  • выполняют X+1 вызовов функции на каждый шаг

  • требуют X производных. Они могут быть переданы как кортеж функций (f, f', f'', …) или как

функция, возвращающая коэффициенты: x -> (f(x), f(x)/f'(x), f'(x)/f''(x), …).

Примеры

using ForwardDiff
Base.adjoint(f::Function)  = x  -> ForwardDiff.derivative(f, float(x))
f(x) = (exp(x) + x - 2)^6
x0 = 1/4
find_zero((f, f', f''), x0, Roots.Halley())               # 14 итераций; ≈ 48 вычислений функции
find_zero((f, f', f''), big(x0), Roots.Thukral2B())       #  3 итерации; ≈ 9 вычислений функции
find_zero((f, f', f'', f'''), big(x0), Roots.Thukral3B()) #  2 итерации; ≈ 8 вычислений функции

Справка

Introduction to a family of Thukral -order method for finding multiple zeros of nonlinear equations, R. Thukral, JOURNAL OF ADVANCES IN MATHEMATICS 13(3):7230-7237, DOI: 10.24297/jam.v13i3.6146.

Гибридные методы

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

Order0()

Метод Order0 разработан как более надежная, хотя, возможно, и более медленная альтернатива другим методам нахождения корней без производных. Реализация примерно следует алгоритму, описанному в статье Personal Calculator Has Key to Solve Any Equation , кнопка SOLVE из HP-34C. Основная идея заключается в использовании шага секущей. Если по пути будет найдена скобка, следует переключиться на алгоритм заключения в скобки, используя AlefeldPotraShi. Если шаг секущей не уменьшает значение функции, используется шаг квадратичного метода до раз.

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

Скорость сходимости

Порядок метода —  , где e_{i+1} \approx e_i^q. Известно, что метод Ньютона является квадратичным для простых корней; метод секущей порядка . Однако для метода Ньютона требуется вызовов, а для метода секущей — только . Асимптотическая эффективность равна ценой увеличения числа вызовов функций. Есть и другие методы порядка , требующие вызовов функции за шаг, например метод Галлея. Другие же требуют меньше шагов, как показано ниже. Многие методы используют обратные квадратичные шаги, а некоторые — обратные кубические. Они имеют порядок при решении уравнения ( для квадратичного случая). В случае с надежными методами для достижения скорости сходимости, как правило, требуется дополнительный вызов функции, что хорошо видно на примере метода Schroder.

Тип Метод Порядок Вычислений F Асимптотическая эффективность

Гибридный

Order0

Без вычисления производных

Secant

Без вычисления производных

Steffensen

Без вычисления производных

Order5

Без вычисления производных

Order8

Без вычисления производных

Order16

Классический

Newton

Классический

Halley

Классический

QuadraticInverse

Классический

ChebyshevLike

Классический

SuperHalley

Многошаговый

LithBoonkkampIJzerman{S,D}

варьируется, максимум

Ограничивающий

BisectionExact

Ограничивающий

A42

Ограничивающий

AlefeldPotraShi

Ограничивающий

Brent

Ограничивающий

ITP

Ограничивающий

Ridders

Ограничивающий

FalsePosition

Ограничивающий

LithBoonkkampIJzermanBracket

Надежный

King

Надежный

Esser

Надежный

Schroder

Надежный

Thukral3

Надежный

Thukral4

Надежный

Thukral5

Сходимость

Для определения того, сходится ли алгоритм или расходится, требуется указание допусков и условий сходимости.

В случае с точным делением пополам сходимость гарантируется математически. Для чисел с плавающей запятой либо находится точный нуль, либо ограничивающий интервал может быть разделен на ], где и  — соседние значения с плавающей запятой. То есть в случае с числами с плавающей запятой имеет минимальный возможный размер. Это можно рассматривать как условие остановки в . Для преждевременного завершения (обеспечивающего меньшую точность, но и меньшее число вызовов функции) можно задать допуск, чтобы при достаточно небольшом алгоритм успешно завершал выполнение. В случае со значениями с плавающей запятой для оценки равенства требуются два допуска: относительный допуск, так как минимальная разница значений с плавающей запятой зависит от размера и , и абсолютный допуск для значений в окрестности . Для определения близости в функцию Base.isapprox передаются значения xrtol и xatol.

Не во всех задачах можно полагаться на близость двух значений . В некоторых примерах разность может быть весьма малой и даже равна , хотя значение не близко к нулю . В связи с этим для неограничивающих методов также применяется проверка размера . При поиске аппроксимаций с плавающей запятой к , то есть нуля, необходимо учитывать значения, малые при малом . При использовании аппроксимации Тейлора можно ожидать, что они будут в районе . Таким образом, малая величина зависит от размера и производной в точке . Первая зависимость учитывается посредством относительной и абсолютной погрешностей (rtol и atol). Размер зависит от задачи, и для его учета можно увеличивать относительную или абсолютную погрешность.

Когда алгоритм возвращает значение NaN, его выполнение завершается. Это может происходить при приближении к сходимости либо указывать на наличие проблем. При досрочном завершении сходимость проверяется в пределах размера с увеличенной погрешностью, если задан аргумент strict=false (по умолчанию).

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

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

Условия сходимости зависят от метода и оцениваются методами Roots.assess_convergence.

Roots.assess_convergence(method, state, options)

Определяет сходимость алгоритма.

Возвращает флаг сходимости и логическое значение, указывающее, завершился ли алгоритм (со схождением или без схождения).

Если алгоритм не сошелся, возвращается (:not_converged, false).

Если алгоритм остановился или сошелся, возвращается флаг и true. Существуют следующие флаги:

  • :x_converged, если xn1 ≈ xn, обычно с заданными ненулевыми допусками.

  • :f_converged, если |f(xn1)| < max(atol, |xn1|*rtol)

  • :nan или :inf, если fxn1 является NaN или бесконечностью.

  • :not_converged, если алгоритм должен продолжить выполнение.

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

В decide_convergence значения остановки (и :x_converged при strict=false) проверяются на сходимость с нестрогим допуском.

Погрешности по умолчанию задаются посредством методов Roots.default_tolerances.

default_tolerances(M::AbstractUnivariateZeroMethod, [T], [S])

По умолчанию для большинства методов используются допуски xatol=eps(T), xrtol=eps(T), atol=4eps(S) и rtol=4eps(S) с соответствующими единицами (абсолютные допуски имеют единицы x и f(x); относительные допуски единицы не используют). Для значений Complex{T} используется T.

Количество итераций ограничивается maxiters=40.

default_tolerances(M::AbstractBisectionMethod, [T], [S])

Для метода Bisection, когда значения x имеют тип Float64, Float32 или Float16, допуски по умолчанию равны нулю, а количество итераций не ограничено. В этом случае алгоритм гарантированно сходится к точному нулю или точке, в которой функция изменяет знак одного из соседних значений ответа с плавающей запятой.

Для других типов по умолчанию заданы ненулевые допуски для xatol и xrtol.

Упрощенные версии

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

secant_method(f, xs; [atol=0.0, rtol=8eps(), maxevals=1000])

Выполняет метод секущей для решения f(x) = 0.

Метод секущей — это итеративный метод с шагом обновления, заданным с помощью b - fb/m, где m — наклон секущей линии между (a,fa) и (b,fb).

Начальные значения могут быть указаны как пара двух значений, как в (x₀, x₁) или [x₀, x₁], или как одно значение x₁, в этом случае выбирается значение x₀.

Алгоритм возвращает m, когда abs(fm) <= max(atol, abs(m) * rtol). Если это не произойдет до шагов maxevals, или алгоритм столкнется с проблемой, будет возвращено значение NaN. Если выполнено слишком много шагов, текущее значение проверяется на предмет изменения знака для соседних значений с плавающей запятой.

Метод Order1 для find_zero также реализует метод секущей. Этот вариант должен быть немного быстрее, поскольку требуется меньше затрат на настройку.

Примеры:

Roots.secant_method(sin, (3,4))
Roots.secant_method(x -> x^5 -x - 1, 1.1)

Эта функция будет специализировать функцию f, так что начальный вызов может занять больше времени, чем вызов метода Order1(), хотя последующие вызовы будут происходить гораздо быстрее. Использование FunctionWrappers.jl может гарантировать, что начальный вызов будет таким же быстрым, как и последующие.

 
bisection(f, a, b; [xatol, xrtol])

Выполняет метод деления пополам для нахождения нуля непрерывной функции.

Предполагается, что (a,b) является скобкой, то есть функция имеет разные знаки в точках a и b. Интервал (a,b) преобразуется в значение с плавающей запятой и уменьшается, когда a или b являются бесконечными. В типичном случае функция f может быть бесконечной. Если f не является непрерывной, алгоритм может найти не только нули, но и точки перехода через ось x.

Если заданы нетривиальные допуски, процесс завершится, когда скобка (a,b) удовлетворяет условию isapprox(a, b, atol=xatol, rtol=xrtol) Для нулевых допусков, по умолчанию, для значений. For zero tolerances, the default, for Float64, Float32 или значений Float16 процесс завершится при значении x, при котором f(x)=0 или f(x)*f(prevfloat(x)) < 0 или f(x) * f(nextfloat(x)) < 0. For other number types, the Для других числовых типов используется метод Roots.A42.

muller(f, xᵢ; xatol=nothing, xrtol=nothing, maxevals=100)
muller(f, xᵢ₋₂, xᵢ₋₁, xᵢ; xatol=nothing, xrtol=nothing, maxevals=100)

Метод Мюллера обобщает метод секущей, но использует квадратичную интерполяцию между тремя точками вместо линейной интерполяции между двумя. Решение для поиска нулей квадратичного уравнения позволяет находить сложные пары корней. Учитывая три предыдущих предположения о корне xᵢ₋₂, xᵢ₋₁, xᵢ и значения многочлена f в этих точках, выводится следующая аппроксимация xᵢ₊₁.

Отрывок и алгоритм взяты из статьи

W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery Numerical Recipes in C, Cambridge University Press (2002), p. 371

Сходимость здесь определяется тем, что xᵢ₊₁ ≈ xᵢ с использованием указанных допусков, которые оба по умолчанию имеют значение eps(one(typeof(abs(xᵢ))))^4/5 в соответствующих единицах. Каждая итерация выполняет три вычисления f. Первый метод случайным образом выбирает две оставшиеся точки в относительной близости от xᵢ.

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

Примеры:

muller(x->x^3-1, 0.5, 0.5im, -0.5) # → -0.500 + 0.866…im
muller(x->x^2+2, 0.0, 0.5, 1.0) # → ≈ 0.00 - 1.41…im
muller(x->(x-5)*x*(x+5), rand(3)...) # → ≈ 0.00
muller(x->x^3-1, 1.5, 1.0, 2.0) # → 2.0, без сходимости
Roots.newton(f, fp, x0; kwargs...)

Реализация метода Ньютона: xᵢ₊₁ = xᵢ - f(xᵢ)/f'(xᵢ).

Аргументы:

  • f::Function — функция, нуль которой требуется найти.

  • fp::Function — производная f.

  • x0::Number — начальное предположение. Для метода Ньютона это может быть сложно сделать.

Пакет ForwardDiff позволяет вычислять производные автоматически. Например, при определении D(f) = x -> ForwardDiff.derivative(f, float(x)) можно использовать D(f) для первой производной.

Именованные аргументы передаются функции find_zero с помощью метода Roots.Newton().

Более простые реализации см. также в описании Roots.newton((f,fp), x0) и Roots.newton(fΔf, x0).

dfree(f, xs)

Более надежная реализация метода секущей

Решение для f(x) = 0 с использованием алгоритма из статьи Personal Calculator Has Key to Solve Any Equation f(x) = 0, кнопка SOLVE из HP-34C.

Также реализовано как метод Order0 для find_zero.

Начальные значения могут быть указаны как пара двух значений, как в (a,b) или [a,b], или как одно значение, в этом случае вычисляется значение b, возможно, из fb. Основная идея заключается в том, чтобы следовать методу секущих, за исключением случаев, когда

  • найдена скобка, тогда используется AlefeldPotraShi;

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

Сходимость наступает, когда f(m) == 0, происходит смена знака между m и соседним значением с плавающей запятой или f(m) <= 2^3*eps(m).

Значение NaN возвращается, если алгоритм делает слишком много шагов перед нахождением нуля.

Примеры

Roots.dfree(x -> x^5 - x - 1, 1.0)

Интерфейс MATLAB

Согласно исходной схеме именования функция называлась fzero, а не fzeros, то есть так же, как аналогичная функция в MATLAB — fzero. Этот интерфейс не рекомендуется к использованию, но все еще поддерживается.

fzero(f, x0; order=0; kwargs...)
fzero(f, x0, M; kwargs...)
fzero(f, x0, M, N; kwargs...)
fzero(f, x0; kwargs...)
fzero(f, a::Number, b::Number; kwargs...)
fzero(f, a::Number, b::Number; order=?, kwargs...)
fzero(f, fp, a::Number; kwargs...)

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

  • f: скалярная функция или вызываемый объект.

  • x0: начальное предположение, скалярное значение или кортеж из двух значений.

  • order: целое число, символ или строка, указывающая алгоритм, который следует использовать для find_zero. Order0 по умолчанию можно указать напрямую с помощью order=0, order=:0 или order="0"; Order1() с помощью order=1, order=:1, order="1" или order=:secant; Order1B() с помощью order="1B" и т. д.

  • M: конкретный метод, как был бы передан find_zero без использования ключевого слова order.

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

  • a, b: если при передаче двух значений не указано значение order, Bisection будет использоваться в заключаемом в скобки интервале (a,b). Если указано значение order, значение x0 будет установлено равным (a,b), и будет использоваться заданный метод.

  • fp: если задано значение fp (предполагается, что вычисляется производная f), будет использоваться метод Ньютона.

  • kwargs...: спецификацию допусков и других именованных аргументов см. в описании find_zero.

Примеры:

fzero(sin, 3)                  # использование метода Order0(), по умолчанию
fzero(sin, 3, order=:secant)   # использование метода секущей (также только `order=1`)
fzero(sin, 3, Roots.Order1B()) # использование варианта метода секущей для кратных корней.
fzero(sin, 3, 4)               # использование метода деления пополам в (3,4)
fzero(sin, 3, 4, xatol=1e-6)   # использование метода деления пополам до тех пор, пока |x_n - x_{n-1}| <= 1e-6
fzero(sin, 3, 3.1, order=1)    # использование метода секущей с x_0=3.0, x_1 = 3.1
fzero(sin, (3, 3.1), order=2)  # использование метода Стеффенсена с x_0=3.0, x_1 = 3.1
fzero(sin, cos, 3)             # использование метода Ньютона

В отличие от find_zero, fzero не специализирует тип аргумента функции. Преимущество этого способа заключается в том, что первое использование функции f происходит быстрее, а последующие — медленнее.

 
fzeros(f, a, b; kwargs...)
fzeros(f, ab; kwargs...)

Выполняет поиск всех нулей f в интервале (a,b). Предполагается, что ни a, ни b не равны нулю.

Интерфейс совместимости для find_zeros.

Отслеживание итераций

При вызове функции find_zero можно добавить именованный аргумент verbose=true для получения подробных сведений о решении и данных по каждой итерации. Для сохранения этих данных можно передать объект Tracks в аргументе tracks.

Roots.Tracks{T,S}

Экземпляр Tracks используется для записи хода выполнения алгоритма. T — это тип входных данных функции, а S — тип выходных данных функции. Оба по умолчанию имеют тип Float64. Обратите внимание, что поскольку этот тип не экспортируется, вам придется написать Roots.Tracks(), чтобы создать объект Tracks.

По умолчанию при нахождении корня отслеживание не ведется. Для изменения этого поведения создайте объект Tracks и передайте его именованному аргументу tracks. При этом изменится объект Tracks, хранящий значения входных данных и функции при каждой итерации, а также дополнительную информацию о процессе поиска корней.

Объекты Tracks отображаются в удобном для чтения формате. На внутреннем уровне хранится либо кортеж пар (x,f(x)), либо пары (aₙ, bₙ), последние — для методов заключения в скобки. (Эти особенности реализации могут быть изменены без уведомления.) Интерес могут представлять методы empty! для сброса объекта Tracks, get для получения отслеживаний, last для получения значения преобразования.

Если вы хотите только вывести только информацию, но в дальнейшем она вам не понадобится, можно передать verbose=true функции нахождения корней. Это не повлияет на возвращаемое значение, которое по-прежнему будет корнем функции.

Примеры

julia> using Roots

julia> f(x) = x^2-2
f (generic function with 1 method)

julia> tracker = Roots.Tracks()
Algorithm has not been run

julia> find_zero(f, (0, 2), Roots.Secant(), tracks=tracker) ≈ √2
true

julia> tracker
Results of univariate zero finding:

* Converged to: 1.4142135623730947
* Algorithm: Secant()
* iterations: 7
* function evaluations ≈ 9
* stopped as |f(x_n)| ≤ max(δ, |x|⋅ϵ) using δ = atol, ϵ = rtol

Trace:
x₁ = 0,	 fx₁ = -2
x₂ = 2,	 fx₂ = 2
x₃ = 1,	 fx₃ = -1
x₄ = 1.3333333333333333,	 fx₄ = -0.22222222222222232
x₅ = 1.4285714285714286,	 fx₅ = 0.04081632653061229
x₆ = 1.4137931034482758,	 fx₆ = -0.0011890606420930094
x₇ = 1.4142114384748701,	 fx₇ = -6.0072868388605372e-06
x₈ = 1.4142135626888697,	 fx₈ = 8.9314555751229818e-10
x₉ = 1.4142135623730947,	 fx₉ = -8.8817841970012523e-16

julia> empty!(tracker)  # сбрасывает

julia> find_zero(sin, (3, 4), Roots.A42(), tracks=tracker) ≈ π
true

julia> get(tracker)
4-element Vector{NamedTuple{names, Tuple{Float64, Float64}} where names}:
 (a = 3.0, b = 4.0)
 (a = 3.0, b = 3.157162792479947)
 (a = 3.141592614491745, b = 3.1415926926910007)
 (a = 3.141592653589793, b = 3.141592653589794)

julia> last(tracker)
3.141592653589793

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