Руководство пользователя
Чтобы ознакомиться с простыми примерами, перейдите в раздел Примеры. Чтобы ознакомиться с подробными инструкциями, читайте дальше.
TaylorSeries.jl — это базовый алгебраический инструмент для работы с многочленами по одной или нескольким переменным. Эти два случая рассматриваются отдельно. Определены три новых типа (Taylor1, HomogeneousPolynomial и TaylorN), которые соответствуют разложениям по одной независимой переменной, однородным многочленам различных переменных и рядам многочленов по множеству независимых переменных, соответственно. Эти типы являются подтипами типа AbstractSeries, который, в свою очередь, является подтипом Number, и определяются параметрически.
Пакет загружается обычным образом:
julia> using TaylorSeriesОдна независимая переменная
Разложения по Тейлору по одной переменной представлены типом Taylor1, который состоит из вектора коэффициентов (имя поля coeffs) и максимального порядка, учитываемого при разложении (имя поля order). Коэффициенты располагаются в порядке возрастания по отношению к степени одночлена, так что coeffs[1] является постоянным членом, coeffs[2] является членом первого порядка (t^1) и т. д. Однако возможно и естественное упорядочение по степени (см. ниже). Это плотное представление многочлена. Порядок многочлена можно не указывать в конструкторе, тогда он фиксируется длиной вектора коэффициентов. Если длина вектора не соответствует order, используется order, что эффективно усекает многочлен до степени order.
julia> Taylor1([1, 2, 3],4) # Многочлен 4-го порядка с коэффициентами 1, 2, 3
1 + 2 x + 3 x²
julia> Taylor1([0.0, 1im]) # Также работает с комплексными числами
( 0.0 + 1.0im ) x
julia> Taylor1(ones(8), 2) # Многочлен, усеченный до 2-го порядка
1.0 + 1.0 x + 1.0 x²
julia> shift_taylor(a) = a + Taylor1(typeof(a),5) ## a + многочлен Тейлора 5-го порядка
shift_taylor (generic function with 1 method)
julia> t = shift_taylor(0.0) # Независимая переменная `t`
1.0 x|
Информация о максимальном рассматриваемом порядке отображается с использованием нотации big-xD835xDCAA. При сочетании различных порядков, а также когда определенные функции используются таким образом, что уменьшают порядок (например, |
В некоторых случаях желательно не отображать нотацию big-xD835__xDCAA. Ее отображение можно контролировать с помощью функции displayBigO.
julia> displayBigO(false) # отключение отображения нотации big O
false
julia> t
1.0 x
julia> displayBigO(true) # включение отображения
true
julia> t
1.0 x + 𝒪(x⁶)Кроме того, управлять некоторыми аспектами формата отображаемых рядов можно с помощью функции use_show_default, use_show_default(true) использует Base.show_default, а use_show_default(false) использует пользовательскую форму отображения (по умолчанию).
julia> use_show_default(true) # использование метода Base.show
true
julia> t
Taylor1{Float64}([0.0, 1.0, 0.0, 0.0, 0.0, 0.0], 5)
julia> use_show_default(false) # использование пользовательского `show`
false
julia> t
1.0 x + 𝒪(x⁶)В определении shift_taylor(a) используется метод Taylor1([::Type{Float64}\], order::Int), который является сокращением для определения независимой переменной разложения по Тейлору заданного типа и порядка (по умолчанию Float64). Предварительное определение независимой переменной является одним из самых простых способов использования пакета.
Обычные арифметические операторы (+, -, *, /, ^, ==) были расширены для работы с типом Taylor1, включая продвижения, предполагающие использование Number. Операции возвращают правильное разложение по Тейлору, соответствующее порядку ряда. Это видно из предпоследнего примера, где первый ненулевой коэффициент находится за пределами порядка разложения, а значит, результат равен нулю.
julia> t*(3t+2.5)
2.5 x + 3.0 x² + 𝒪(x⁶)
julia> 1/(1-t)
1.0 + 1.0 x + 1.0 x² + 1.0 x³ + 1.0 x⁴ + 1.0 x⁵ + 𝒪(x⁶)
julia> t*(t^2-4)/(t+2)
- 2.0 x + 1.0 x² + 𝒪(x⁶)
julia> tI = im*t
( 0.0 + 1.0im ) x + 𝒪(x⁶)
julia> (1-t)^3.2
1.0 - 3.2 x + 3.5200000000000005 x² - 1.4080000000000004 x³ + 0.07040000000000009 x⁴ + 0.011264000000000012 x⁵ + 𝒪(x⁶)
julia> (1+t)^t
1.0 + 1.0 x² - 0.5 x³ + 0.8333333333333333 x⁴ - 0.75 x⁵ + 𝒪(x⁶)
julia> Taylor1(3) + Taylor1(5) == 2Taylor1(3) # применяется соглашение big-_xD835__xDCAA_
true
julia> t^6 # t имеет 5-й порядок
0.0 + 𝒪(x⁶)
julia> t^2 / t # Результат имеет 4-й порядок, а не 5-й
1.0 x + 𝒪(x⁵)Обратите внимание, что последний пример возвращает ряд Taylor1 4-го порядка, а не 5-го. Это соответствует количеству известных коэффициентов возвращаемых рядов, поскольку результат соответствует факторизации t в числителе для отмены того же коэффициента в знаменателе.
Taylor1 также имеет общий порядок, который обеспечивается перегрузкой isless. Упорядочение согласуется с той концепцией, что в алгебре есть бесконечно малые элементы. Подробнее см. в документе M. Berz, Automatic Differentiation as Nonarchimedean Analysis, Computer Arithmetic and Enclosure Methods, (1992), Elsevier, 439-450. Это можно проиллюстрировать на примере.
julia> 1 > t > 2*t^2 > 100*t^3 >= 0
trueЕсли не удается вычислить правильное разложение по Тейлору, возникает ошибка, например когда производная не определена в точке разложения или просто расходится.
julia> 1/t
ERROR: ArgumentError: Division does not define a Taylor1 polynomial;
order k=0 => coeff[0]=Inf.
julia> t^3.2
ERROR: DomainError with 1.0 x + 𝒪(x⁶):
The 0-th order Taylor1 coefficient must be non-zero
to raise the Taylor1 polynomial to a non-integer exponent.
julia> abs(t)
ERROR: DomainError with 1.0 x + 𝒪(x⁶):
The 0th order Taylor1 coefficient must be non-zero
(abs(x) is not differentiable at x=0).Реализовано несколько элементарных функций, коэффициенты которых вычисляются рекурсивно. На момент написания статьи этими функциями являются exp, log, sqrt, тригонометрические функции sin, cos и tan, их инверсии, а также гиперболические функции sinh, cosh и tanh и их инверсии. В будущем будут добавлены и другие функции. Заметим, что этот способ получения коэффициентов Тейлора не является отложенным, в частности для многих независимых переменных. Тем не менее реализация достаточно эффективна, особенно для интегрирования обыкновенных дифференциальных уравнений, что является одним из применений, которые мы имеем в виду (см. TaylorIntegration.jl).
julia> exp(t)
1.0 + 1.0 x + 0.5 x² + 0.16666666666666666 x³ + 0.041666666666666664 x⁴ + 0.008333333333333333 x⁵ + 𝒪(x⁶)
julia> log(1-t)
- 1.0 x - 0.5 x² - 0.3333333333333333 x³ - 0.25 x⁴ - 0.2 x⁵ + 𝒪(x⁶)
julia> sqrt(1 + t)
1.0 + 0.5 x - 0.125 x² + 0.0625 x³ - 0.0390625 x⁴ + 0.02734375 x⁵ + 𝒪(x⁶)
julia> imag(exp(tI)')
- 1.0 x + 0.16666666666666666 x³ - 0.008333333333333333 x⁵ + 𝒪(x⁶)
julia> real(exp(Taylor1([0.0,1im],17))) - cos(Taylor1([0.0,1.0],17)) == 0.0
true
julia> convert(Taylor1{Rational{Int}}, exp(t))
1//1 + 1//1 x + 1//2 x² + 1//6 x³ + 1//24 x⁴ + 1//120 x⁵ + 𝒪(x⁶)Опять же, ошибки возникают, когда это необходимо.
julia> sqrt(t)
ERROR: DomainError with 1.0 x + 𝒪(x⁶):
First non-vanishing Taylor1 coefficient must correspond
to an **even power** in order to expand `sqrt` around 0.
julia> log(t)
ERROR: DomainError with 1.0 x + 𝒪(x⁶):
The 0-th order coefficient must be non-zero in order to expand `log` around 0.Для получения конкретного коэффициента можно использовать getcoeff. Еще одной альтернативой является запрос конкретной степени с использованием векторной нотации, где индекс соответствует степени члена.
julia> expon = exp(t)
1.0 + 1.0 x + 0.5 x² + 0.16666666666666666 x³ + 0.041666666666666664 x⁴ + 0.008333333333333333 x⁵ + 𝒪(x⁶)
julia> getcoeff(expon, 0) == expon[0]
true
julia> rationalize(expon[3])
1//6Обратите внимание, что некоторые арифметические операции или применение некоторых функций могут изменить порядок результата, как упоминалось выше.
julia> t # независимая переменная 5-го порядка
1.0 x + 𝒪(x⁶)
julia> t^2/t # возвращает разложение в ряды 4-го порядка
1.0 x + 𝒪(x⁵)
julia> sqrt(t^2) # возвращает разложение в ряды 2-го порядка
1.0 x + 𝒪(x³)
julia> (t^4)^(1/4) # ряды 1-го порядка
1.0 x + 𝒪(x²)Дифференцирование и интегрирование можно легко выполнить для разложения полиномов по одной переменной с помощью функций differentiate и integrate. (Функция derivative является синонимом differentiate.) Эти функции возвращают соответствующие разложения Taylor1. Порядок производной Taylor1 уменьшается на 1. Для интеграла константа интегрирования может быть задана пользователем (по умолчанию равна нулю). Порядок интегрируемого многочлена для интеграла остается неизменным. Значение -й ( ) производной будет получено с помощью differentiate(n,a), где a — это ряд Тейлора. Аналогичным образом многочлен Taylor1 -й производной будет получен как differentiate(a,n). Итоговый многочлен имеет порядок get_order(a)-n.
julia> differentiate(exp(t)) # exp(t) имеет 5-й порядок; производная имеет 4-й порядок
1.0 + 1.0 x + 0.5 x² + 0.16666666666666666 x³ + 0.041666666666666664 x⁴ + 𝒪(x⁵)
julia> integrate(exp(t)) # полученный ряд Тейлора имеет 5-й порядок
1.0 x + 0.5 x² + 0.16666666666666666 x³ + 0.041666666666666664 x⁴ + 0.008333333333333333 x⁵ + 𝒪(x⁶)
julia> integrate( exp(t), 1.0)
1.0 + 1.0 x + 0.5 x² + 0.16666666666666666 x³ + 0.041666666666666664 x⁴ + 0.008333333333333333 x⁵ + 𝒪(x⁶)
julia> integrate( differentiate( exp(-t)), 1.0 ) == exp(-t)
true
julia> differentiate(1, exp(shift_taylor(1.0))) == exp(1.0)
true
julia> differentiate(5, exp(shift_taylor(1.0))) == exp(1.0) # 5-е дифференцирование от `exp(1+t)`
true
julia> derivative(exp(1+t), 3) # многочлен Taylor1 3-й производной от `exp(1+t)`
2.718281828459045 + 2.718281828459045 x + 1.3591409142295225 x² + 𝒪(x³)Для вычисления ряда Тейлора в заданной точке используется правило Горнера с функцией evaluate(a, dt). Здесь dt — это приращение от точки , в окрестности которой вычисляется разложение по Тейлору для a, то есть ряд вычисляется в точке . Пропуск dt соответствует (см. evaluate).
julia> evaluate(exp(shift_taylor(1.0))) - ℯ # exp(t) в окрестности t0=1 (5-й порядок), вычисляется там (dt=0)
0.0
julia> evaluate(exp(t), 1) - ℯ # exp(t) в окрестности t0=0 (5-й порядок), вычисляется в t=1
-0.0016151617923783057
julia> evaluate(exp( Taylor1(17) ), 1) - ℯ # exp(t) в окрестности t0=0, 17-й порядок
0.0
julia> tBig = Taylor1(BigFloat, 50) # Независимая переменная с BigFloats, 50-й порядок
1.0 x + 𝒪(x⁵¹)
julia> eBig = evaluate( exp(tBig), one(BigFloat) )
2.718281828459045235360287471352662497757247093699959574966967627723419298053556
julia> ℯ - eBig
6.573322999985292556154129119543257102601105719980995128942636339920549561322098e-67Другим способом вычисления значения многочлена Taylor1 p при заданном значении x является вызов p так, как если бы это была функция, то есть p(x):
julia> t = Taylor1(15)
1.0 x + 𝒪(x¹⁶)
julia> p = sin(t)
1.0 x - 0.16666666666666666 x³ + 0.008333333333333333 x⁵ - 0.0001984126984126984 x⁷ + 2.7557319223985893e-6 x⁹ - 2.505210838544172e-8 x¹¹ + 1.6059043836821616e-10 x¹³ - 7.647163731819817e-13 x¹⁵ + 𝒪(x¹⁶)
julia> evaluate(p, pi/2) # значение p при pi/2 с использованием `evaluate`
0.9999999999939766
julia> p(pi/2) # значение p при pi/2 при вычислении p как функции
0.9999999999939766
julia> p(pi/2) == evaluate(p, pi/2)
true
julia> p(0.0)
0.0
julia> p() == p(0.0) # p() — это сокращение для получения коэффициента 0-го порядка для `p`
trueОбратите внимание, что синтаксис p(x) эквивалентен evaluate(p, x), а p() эквивалентен evaluate(p).
Полезными являются функции taylor_expand и update!. Первая возвращает разложение функции в окрестности заданного значения t0, имитируя использование shift_taylor выше. В свою очередь, update! обеспечивает обновление заданного многочлена Тейлора на месте, то есть сдвигает его дальше на заданную величину.
julia> p = taylor_expand( x -> sin(x), pi/2, order=16) # разложение 16-го порядка sin(t) в окрестности pi/2
1.0 + 6.123233995736766e-17 x - 0.5 x² - 1.020538999289461e-17 x³ + 0.041666666666666664 x⁴ + 5.102694996447305e-19 x⁵ - 0.001388888888888889 x⁶ - 1.2149273801065013e-20 x⁷ + 2.48015873015873e-5 x⁸ + 1.6873991390368074e-22 x⁹ - 2.7557319223985894e-7 x¹⁰ - 1.5339992173061888e-24 x¹¹ + 2.08767569878681e-9 x¹² + 9.833328316065313e-27 x¹³ - 1.1470745597729726e-11 x¹⁴ - 4.682537293364434e-29 x¹⁵ + 4.779477332387386e-14 x¹⁶ + 𝒪(x¹⁷)
julia> update!(p, 0.025) # обновляет расширение, заданное p, сдвигая его еще на 0,025
julia> p
0.9996875162757026 - 0.02499739591471227 x - 0.49984375813785126 x² + 0.004166232652452044 x³ + 0.0416536465114876 x⁴ - 0.00020831163262260228 x⁵ - 0.0013884548837162537 x⁶ + 4.959800776728625e-6 x⁷ + 2.4793837209218816e-5 x⁸ - 6.888612189900868e-8 x⁹ - 2.7548708010243125e-7 x¹⁰ + 6.262374718092001e-10 x¹¹ + 2.0870233341100357e-9 x¹² - 4.014342754938811e-12 x¹³ - 1.1467160989730436e-11 x¹⁴ + 1.9117909329549495e-14 x¹⁵ + 4.779477332387386e-14 x¹⁶ + 𝒪(x¹⁷)Множество переменных
Многочлен по переменным может быть представлен (по крайней мере) двумя способами. Как вектор, коэффициенты которого являются однородными многочленами фиксированной степени, или как вектор, коэффициенты которого являются многочленами по переменным . Текущая реализация TaylorSeries.jl соответствует первому варианту, хотя была создана некоторая инфраструктура, позволяющая развивать второй вариант. Элегантная (отложенная) реализация второго представления обсуждалась здесь.
Структура TaylorN строится как вектор параметризованных однородных многочленов, определяемых типом HomogeneousPolynomial, который, в свою очередь, является упорядоченным вектором коэффициентов заданного порядка (степени). Данная реализация обязывает пользователя в начале указать (максимальный) порядок и количество независимых переменных, что удобно сделать с помощью set_variables. Возвращается вектор результирующих переменных Тейлора:
julia> x, y = set_variables("x y")
2-element Vector{TaylorN{Float64}}:
1.0 x + 𝒪(‖x‖⁷)
1.0 y + 𝒪(‖x‖⁷)
julia> typeof(x)
TaylorN{Float64}
julia> x.order
6
julia> x.coeffs
7-element Vector{HomogeneousPolynomial{Float64}}:
0.0
1.0 x
0.0
0.0
0.0
0.0
0.0Как показано, результирующие объекты имеют тип TaylorN{Float64}. В set_variables есть необязательный именованный аргумент order, который используется для указания максимального порядка многочленов TaylorN. Обратите внимание, что переменные можно задавать с помощью вектора символов.
julia> set_variables([:x, :y], order=10)
2-element Vector{TaylorN{Float64}}:
1.0 x + 𝒪(‖x‖¹¹)
1.0 y + 𝒪(‖x‖¹¹)Аналогично, переменные с нижними индексами также доступны при указании имени одной переменной и необязательного именованного аргумента numvars:
julia> set_variables("α", numvars=3)
3-element Vector{TaylorN{Float64}}:
1.0 α₁ + 𝒪(‖x‖¹¹)
1.0 α₂ + 𝒪(‖x‖¹¹)
1.0 α₃ + 𝒪(‖x‖¹¹)В качестве альтернативы set_variables можно использовать get_variables, если вы не хотите изменять внутренние словари. get_variables возвращает вектор независимых переменных TaylorN желаемого порядка (order) (меньше get_order, чтобы не пришлось изменять внутренние компоненты) с длиной и именами переменных, изначально определенными с помощью set_variables.
julia> get_variables(2) # вектор независимых переменных 2-го порядка
3-element Vector{TaylorN{Float64}}:
1.0 α₁ + 𝒪(‖x‖³)
1.0 α₂ + 𝒪(‖x‖³)
1.0 α₃ + 𝒪(‖x‖³)|
Если построение внутренних таблиц не является полностью согласованным, возникает ошибка |
Функция show_params_TaylorN отображает текущие значения параметров в информационном блоке.
julia> show_params_TaylorN()
┌ Info: Parameters for `TaylorN` and `HomogeneousPolynomial`:
│ Maximum order = 10
│ Number of variables = 3
│ Variable names = ["α₁", "α₂", "α₃"]
└ Variable symbols = [:α₁, :α₂, :α₃]Внутреннее изменение параметров (максимального порядка и количества переменных) переопределяет хэш-таблицы, которые преобразуют индекс коэффициентов HomogeneousPolynomial заданного порядка в соответствующие многопеременные одночлены, или наоборот. Фиксация этих значений с самого начала является обязательным условием. Начальными (заданными по умолчанию) значениями являются order = 6 и num_vars=2.
Самый простой способ построения объекта TaylorN заключается в определении независимых переменных. Это можно сделать с помощью функции set_variables, как описано выше, или с помощью метода TaylorN{T<:Number}(::Type{T}, nv::Int) для независимой переменной nv TaylorN{T}. Порядок можно указать с помощью необязательного именованного аргумента order.
julia> x, y = set_variables("x y", numvars=2, order=6);
julia> x
1.0 x + 𝒪(‖x‖⁷)
julia> TaylorN(1, order=4) # переменная 1 4-го порядка
1.0 x + 𝒪(‖x‖⁵)
julia> TaylorN(Int, 2) # переменная 2, тип Int, order=get_order()=6
1 y + 𝒪(‖x‖⁷)Другие способы построения многочленов TaylorN предполагают прямое использование объектов HomogeneousPolynomial, что неудобно.
julia> set_variables(:x, numvars=2); # можно использовать символы
julia> HomogeneousPolynomial([1,-1])
1 x₁ - 1 x₂
julia> TaylorN([HomogeneousPolynomial([1,0]), HomogeneousPolynomial([1,2,3])],4)
1 x₁ + 1 x₁² + 2 x₁ x₂ + 3 x₂² + 𝒪(‖x‖⁵)Разложения по Тейлору реализуются в окрестности 0 для всех переменных. Если требуется разложение в окрестности другого значения, нужно просто преобразовать соответствующую независимую переменную, т. е. .
Как и раньше, функционал обычных арифметических операторов (+, -, *, /, ^, ==) был расширен для работы с объектами TaylorN, включая соответствующие продвижения для работы с обычными числовыми типами. Обратите внимание, что некоторые арифметические операции были расширены для HomogeneousPolynomial, когда результатом является HomogeneousPolynomial. Деление, например, не расширено. При сочетании многочленов TaylorN разного порядка используется то же соглашение, что и для объектов Taylor1. И HomogeneousPolynomial, и TaylorN имеют общий лексикографический порядок, который обеспечивается перегрузкой isless. Упорядочение согласуется с той концепцией, что в алгебре есть бесконечно малые элементы. Подробнее см. в документе M. Berz, Automatic Differentiation as Nonarchimedean Analysis, Computer Arithmetic and Enclosure Methods, (1992), Elsevier, 439-450. По сути, лексикографический порядок работает следующим образом: одночлены меньшего порядка больше, чем одночлены более высокого порядка. Когда порядок одинаков, большие одночлены появляются раньше в хэш-таблицах. Функция show_monomials.
julia> x, y = set_variables("x y", order=10);
julia> show_monomials(2)
1 --> x²
2 --> x y
3 --> y²Тогда имеет место следующее.
julia> 0 < 1e8 * y^2 < x*y < x^2 < y < x/1e8 < 1.0
trueТакже были реализованы элементарные функции путем рекурсивного вычисления их коэффициентов:
julia> exy = exp(x+y)
1.0 + 1.0 x + 1.0 y + 0.5 x² + 1.0 x y + 0.5 y² + 0.16666666666666666 x³ + 0.5 x² y + 0.5 x y² + 0.16666666666666666 y³ + 0.041666666666666664 x⁴ + 0.16666666666666666 x³ y + 0.25 x² y² + 0.16666666666666666 x y³ + 0.041666666666666664 y⁴ + 0.008333333333333333 x⁵ + 0.041666666666666664 x⁴ y + 0.08333333333333333 x³ y² + 0.08333333333333333 x² y³ + 0.041666666666666664 x y⁴ + 0.008333333333333333 y⁵ + 0.001388888888888889 x⁶ + 0.008333333333333333 x⁵ y + 0.020833333333333332 x⁴ y² + 0.027777777777777776 x³ y³ + 0.020833333333333332 x² y⁴ + 0.008333333333333333 x y⁵ + 0.001388888888888889 y⁶ + 0.0001984126984126984 x⁷ + 0.001388888888888889 x⁶ y + 0.004166666666666667 x⁵ y² + 0.006944444444444443 x⁴ y³ + 0.006944444444444443 x³ y⁴ + 0.004166666666666667 x² y⁵ + 0.001388888888888889 x y⁶ + 0.0001984126984126984 y⁷ + 2.48015873015873e-5 x⁸ + 0.0001984126984126984 x⁷ y + 0.0006944444444444445 x⁶ y² + 0.0013888888888888887 x⁵ y³ + 0.0017361111111111108 x⁴ y⁴ + 0.0013888888888888887 x³ y⁵ + 0.0006944444444444445 x² y⁶ + 0.0001984126984126984 x y⁷ + 2.48015873015873e-5 y⁸ + 2.7557319223985893e-6 x⁹ + 2.48015873015873e-5 x⁸ y + 9.920634920634922e-5 x⁷ y² + 0.0002314814814814815 x⁶ y³ + 0.0003472222222222221 x⁵ y⁴ + 0.0003472222222222221 x⁴ y⁵ + 0.0002314814814814815 x³ y⁶ + 9.920634920634922e-5 x² y⁷ + 2.48015873015873e-5 x y⁸ + 2.7557319223985893e-6 y⁹ + 2.7557319223985894e-7 x¹⁰ + 2.755731922398589e-6 x⁹ y + 1.2400793650793652e-5 x⁸ y² + 3.306878306878307e-5 x⁷ y³ + 5.7870370370370366e-5 x⁶ y⁴ + 6.944444444444443e-5 x⁵ y⁵ + 5.7870370370370366e-5 x⁴ y⁶ + 3.306878306878307e-5 x³ y⁷ + 1.2400793650793652e-5 x² y⁸ + 2.755731922398589e-6 x y⁹ + 2.7557319223985894e-7 y¹⁰ + 𝒪(‖x‖¹¹)Функция getcoeff выдает нормализованный коэффициент многочлена, который соответствует одночлену, заданному кортежем или вектором v, содержащим степени. Например, для многочлена exy, приведенного выше, коэффициент одночлена получается с помощью getcoeff(exy, (3,5)) или getcoeff(exy, [3,5]).
julia> getcoeff(exy, (3,5))
0.0013888888888888887
julia> rationalize(ans)
1//720Подобно Taylor1, векторную нотацию можно использовать для запроса определенных коэффициентов объектов HomogeneousPolynomial или TaylorN. Для объектов TaylorN индекс означает степень HomogeneousPolynomial. Для объектов HomogeneousPolynomial индекс означает позицию хэш-таблицы. С помощью функции show_monomials можно получить коэффициент конкретного одночлена с заданной степенью HomogeneousPolynomial.
julia> exy[8] # получение члена 8-го порядка
2.48015873015873e-5 x⁸ + 0.0001984126984126984 x⁷ y + 0.0006944444444444445 x⁶ y² + 0.0013888888888888887 x⁵ y³ + 0.0017361111111111108 x⁴ y⁴ + 0.0013888888888888887 x³ y⁵ + 0.0006944444444444445 x² y⁶ + 0.0001984126984126984 x y⁷ + 2.48015873015873e-5 y⁸
julia> show_monomials(8)
1 --> x⁸
2 --> x⁷ y
3 --> x⁶ y²
4 --> x⁵ y³
5 --> x⁴ y⁴
6 --> x³ y⁵
7 --> x² y⁶
8 --> x y⁷
9 --> y⁸
julia> exy[8][6] # получение 6-го коэффициента члена 8-го порядка
0.0013888888888888887Частичное дифференцирование также реализуется для объектов TaylorN с помощью функции differentiate, указывающей в качестве второго аргумента номер переменной или ее символ.
julia> p = x^3 + 2x^2 * y - 7x + 2
2.0 - 7.0 x + 1.0 x³ + 2.0 x² y + 𝒪(‖x‖¹¹)
julia> q = y - x^4
1.0 y - 1.0 x⁴ + 𝒪(‖x‖¹¹)
julia> differentiate( p, 1 ) # частная производная по первой переменной
- 7.0 + 3.0 x² + 4.0 x y + 𝒪(‖x‖¹¹)
julia> differentiate( q, :y ) # частная производная по :y
1.0 + 𝒪(‖x‖¹¹)Если мы запросим частную производную по неопределенной переменной, будет выведена ошибка.
julia> differentiate( q, 3 ) # ошибка, поскольку мы имеем дело с двумя переменными
ERROR: AssertionError: 1 ≤ r ≤ get_numvars()Для получения более конкретных частных производных у нас есть два специализированных метода, в которых используется кортеж, представляющий собой количество производных по каждой переменной (поэтому длина кортежа должна быть равна фактическому количеству переменных). Эти методы возвращают либо объект TaylorN, либо коэффициент, соответствующий указанному кортежу, нормализованный по факториалам, определенным этим кортежем. Последний по сути является коэффициентом 0-го порядка первого.
julia> differentiate(p, (2,1)) # две производных по :x и одна по :y
4.0 + 𝒪(‖x‖¹¹)
julia> differentiate((2,1), p) # коэффициент 0-го порядка предыдущего выражения
4.0
julia> differentiate(p, (1,1)) # одна производная по :x и одна по :y
4.0 x + 𝒪(‖x‖¹¹)
julia> differentiate((1,1), p) # коэффициент 0-го порядка предыдущего выражения
0.0Интегрирование по r-й переменной для объектов HomogeneousPolynomial и TaylorN выполняется с помощью функции integrate. Обратите внимание, что функция integrate для объектов TaylorN позволяет задать константу интегрирования, которая должна быть независима от интегрируемой переменной. Опять же, эта переменная может быть задана своим символом.
julia> integrate( differentiate( p, 1 ), 1) # интегрирование по первой переменной
- 7.0 x + 1.0 x³ + 2.0 x² y + 𝒪(‖x‖¹¹)
julia> integrate( differentiate( p, 1 ), :x, 2) # интегрирование по :x, константа интегрирования равна 2
2.0 - 7.0 x + 1.0 x³ + 2.0 x² y + 𝒪(‖x‖¹¹)
julia> integrate( differentiate( q, 2 ), :y, -x^4) == q
true
julia> integrate( differentiate( q, 2 ), 2, y)
ERROR: AssertionError: The integration constant ( 1.0 y + 𝒪(‖x‖¹¹)) must be independent of the
y variableevaluate можно использовать и для объектов TaylorN, применяя эту функцию для векторов чисел (Real или Complex). Длина вектора должна совпадать с количеством независимых переменных. Функция evaluate также позволяет указать только одну переменную и значение.
julia> evaluate(exy, [.1,.02]) == exp(0.12)
true
julia> evaluate(exy, :x, 0.0) == exp(y) # вычисление `exy` для :x -> 0
trueАналогично Taylor1, другой способ получить значение многочлена TaylorN p в заданной точке x заключается в его вызове так, как если бы он был функцией: синтаксис p(x) для p::TaylorN эквивалентен evaluate(p,x), а p() эквивалентен evaluate(p).
julia> exy([.1,.02]) == exp(0.12)
true
julia> exy(:x, 0.0)
1.0 + 1.0 y + 0.5 y² + 0.16666666666666666 y³ + 0.041666666666666664 y⁴ + 0.008333333333333333 y⁵ + 0.001388888888888889 y⁶ + 0.0001984126984126984 y⁷ + 2.48015873015873e-5 y⁸ + 2.7557319223985893e-6 y⁹ + 2.7557319223985894e-7 y¹⁰ + 𝒪(‖x‖¹¹)Внутренним образом функция evaluate для TaylorN учитывает отдельно результаты всех объектов HomogeneousPolynomial по порядку (order), которые в итоге складываются после их сортировки на месте (это выполняется по умолчанию) в порядке возрастания по abs2. Это делается для того, чтобы использовать как можно больше значащих цифр всех членов в итоговой сумме, что в итоге должно дать более точный результат. Это значение по умолчанию можно изменить на несортируемое суммирование, которое может быть более производительным или полезным для некоторых подтипов Number, для которых, например, не определено isless. Мотивирующий пример см. в описании этой проблемы. Это можно сделать с помощью ключевого слова sorting в evaluate, которое ожидает Bool, или с помощью логического значения в качестве первого аргумента в функционально-подобном вычислении.
julia> exy([.1,.02]) # по умолчанию используется `sorting=true`
1.1274968515793757
julia> evaluate(exy, [.1,.02]; sorting=false)
1.127496851579376
julia> exy(false, [.1,.02])
1.127496851579376В приведенных примерах первая запись соответствует случаю по умолчанию (sorting=true), который дает тот же результат, что и exp(0.12), а остальные две демонстрируют отключение сортировки членов. Обратите внимание, что результаты не идентичны, поскольку сложение с плавающей запятой не является ассоциативным, что может привести к ошибкам округления.
Функции taylor_expand и update! работают также и для TaylorN.
julia> xysq = x^2 + y^2
1.0 x² + 1.0 y² + 𝒪(‖x‖¹¹)
julia> update!(xysq, [1.0, -2.0]) # разложение в окрестности (1,-2)
julia> xysq
5.0 + 2.0 x - 4.0 y + 1.0 x² + 1.0 y² + 𝒪(‖x‖¹¹)
julia> update!(xysq, [-1.0, 2.0]) # обратный сдвиг
julia> xysq == x^2 + y^2
trueТакже были реализованы функции для вычисления градиента, якобиана и гессиана. Обратите внимание, что эти функции не экспортируются, поэтому для их использования требуется префикс TaylorSeries. Используя многочлены и , определенные выше, можно использовать TaylorSeries.gradient (или ∇). Результаты будут иметь тип Array{TaylorN{T},1}. Чтобы вычислить якобиан и гессиан векторного поля, определенного в точке, мы используем соответственно TaylorSeries.jacobian и TaylorSeries.hessian:
julia> ∇(p)
2-element Vector{TaylorN{Float64}}:
- 7.0 + 3.0 x² + 4.0 x y + 𝒪(‖x‖¹¹)
2.0 x² + 𝒪(‖x‖¹¹)
julia> TaylorSeries.gradient( q )
2-element Vector{TaylorN{Float64}}:
- 4.0 x³ + 𝒪(‖x‖¹¹)
1.0 + 𝒪(‖x‖¹¹)
julia> r = p-q-2*p*q
2.0 - 7.0 x - 5.0 y + 14.0 x y + 1.0 x³ + 2.0 x² y + 5.0 x⁴ - 2.0 x³ y - 4.0 x² y² - 14.0 x⁵ + 2.0 x⁷ + 4.0 x⁶ y + 𝒪(‖x‖¹¹)
julia> TaylorSeries.hessian(ans)
2×2 transpose(::Matrix{Float64}) with eltype Float64:
0.0 14.0
14.0 0.0
julia> TaylorSeries.jacobian([p,q], [2,1])
2×2 transpose(::Matrix{Float64}) with eltype Float64:
13.0 8.0
-32.0 1.0
julia> TaylorSeries.hessian(r, [1.0,1.0])
2×2 transpose(::Matrix{Float64}) with eltype Float64:
-26.0 20.0
20.0 -8.0Другие специфические применения описаны в разделе Examples.
Сочетания
Как было сказано выше, Taylor1{T}, HomogeneousPolynomial{T} и TaylorN{T} являются параметризованными структурами, так что T<:AbstractSeries, последняя является подтипом Number. Тогда мы действительно можем определить разложения по Тейлору по переменным, где одна из переменных (переменная Taylor1) является несколько особенной.
julia> x, y = set_variables("x y", order=3)
2-element Vector{TaylorN{Float64}}:
1.0 x + 𝒪(‖x‖⁴)
1.0 y + 𝒪(‖x‖⁴)
julia> t1N = Taylor1([zero(x), one(x)], 5)
( 1.0 + 𝒪(‖x‖⁴)) x + 𝒪(x⁶)Последняя определяет переменную Taylor1{TaylorN{Float64}}, которая имеет 5-й порядок в t и 3-й порядок в x и y. Затем мы можем вычислить функции с такими многочленами:
julia> cos(2.1+x+t1N)
- 0.5048461045998576 - 0.8632093666488737 x + 0.2524230522999288 x² + 0.14386822777481229 x³ + 𝒪(‖x‖⁴) + ( - 0.8632093666488737 + 0.5048461045998576 x + 0.43160468332443686 x² - 0.0841410174333096 x³ + 𝒪(‖x‖⁴)) x + ( 0.2524230522999288 + 0.43160468332443686 x - 0.1262115261499644 x² - 0.07193411388740614 x³ + 𝒪(‖x‖⁴)) x² + ( 0.14386822777481229 - 0.0841410174333096 x - 0.07193411388740614 x² + 0.0140235029055516 x³ + 𝒪(‖x‖⁴)) x³ + ( - 0.0210352543583274 - 0.03596705694370307 x + 0.0105176271791637 x² + 0.005994509490617178 x³ + 𝒪(‖x‖⁴)) x⁴ + ( - 0.007193411388740615 + 0.00420705087166548 x + 0.0035967056943703073 x² - 0.00070117514527758 x³ + 𝒪(‖x‖⁴)) x⁵ + 𝒪(x⁶)Такого рода разложения представляют интерес при изучении зависимости параметров, например в контексте теории бифуркаций или при рассмотрении зависимости решения дифференциального уравнения от начальных условий в окрестности заданного решения. В этом случае x и y представляют небольшие вариации около заданного значения параметров или определенного начального условия. Такие конструкции используются в пакете TaylorIntegration.jl.