Примеры
Разложение exp(x) с помощью taylor_expand()
Функция taylor_expand принимает в качестве первого аргумента функцию, которую нужно разложить, а в качестве второго — точку, в окрестности которой нужно выполнить разложение. Именованный аргумент order определяет, до какого порядка будет выполняться разложение:
julia> using TaylorSeries
julia> taylor_expand(exp, 0, order=2)
1.0 + 1.0 t + 0.5 t² + 𝒪(t³)Готово! Вычислить простой многочлен Тейлора действительно очень просто. Следующий пример несколько сложнее.
Разложение exp(x) с помощью символьного объекта
Альтернативный способ вычислить однопеременное разложение по Тейлору для функции заключается в определении переменной типа Taylor1 и ее использовании в функции, которую нужно разложить. Аргумент, заданный конструктору Taylor1, — это порядок, до которого нужно выполнить разложение:
julia> using TaylorSeries
julia> x = Taylor1(2)
1.0 t + 𝒪(t³)
julia> exp(x)
1.0 + 1.0 t + 0.5 t² + 𝒪(t³)Давайте также избавимся от выводимой ошибки для следующих нескольких примеров и зададим выводимой независимой переменной значение x:
julia> displayBigO(false)
false
julia> set_taylor1_varname("x")
"x"
julia> exp(x)
1.0 + 1.0 x + 0.5 x²Изменение точки разложения
Переменная, построенная с помощью Taylor1(), автоматически выполняет разложение в окрестности точки x=0. Но что если вы хотите использовать символьный объект для разложения в окрестности точки, отличной от нуля? Поскольку разложение exp(x) в окрестности x=1 точно такое же, как и разложение exp(x+1) в окрестности x=0, просто замените x в выражении на x+1, чтобы выполнить разложение в x=1:
julia> p = exp(x+1)
2.718281828459045 + 2.718281828459045 x + 1.3591409142295225 x²
julia> p(0.01)
2.7456005608350584
julia> exp(1.01)
2.7456010150169163Дополнительные примеры
Вы можете использовать пользовательские функции.
julia> f(a) = 1/(a+1)
f (generic function with 1 method)
julia> f(x)
1.0 - 1.0 x + 1.0 x²Функции могут быть вложенными
julia> sin(f(x))
0.8414709848078965 - 0.5403023058681398 x + 0.11956681346419151 x²и далее усложняться по модульному принципу.
julia> sin(exp(x+2))/(x+2)+cos(x+2)+f(x+2)
0.364113974242596 + 0.41259243488717107 x - 11.843864409375039 x²Тождество четырех квадратов
Первый пример показывает, что тождество четырех квадратов выполняется так:
что первоначально было доказано Эйлером. Код также можно найти в этом тесте пакета.
Сначала мы задаем максимальную степень многочлена равной 4, так как в RHS уравнения априори есть члены четвертого порядка, и определяем 8 независимых переменных.
julia> using TaylorSeries
julia> # Определяем переменные α₁, ..., α₄, β₁, ..., β₄
make_variable(name, index::Int) = string(name, TaylorSeries.subscriptify(index))
make_variable (generic function with 1 method)
julia> variable_names = [make_variable("α", i) for i in 1:4]
4-element Vector{String}:
"α₁"
"α₂"
"α₃"
"α₄"
julia> append!(variable_names, [make_variable("β", i) for i in 1:4])
8-element Vector{String}:
"α₁"
"α₂"
"α₃"
"α₄"
"β₁"
"β₂"
"β₃"
"β₄"
julia> # Создаем переменные TaylorN variables (order=4, numvars=8)
a1, a2, a3, a4, b1, b2, b3, b4 = set_variables(variable_names, order=4)
8-element Vector{TaylorN{Float64}}:
1.0 α₁
1.0 α₂
1.0 α₃
1.0 α₄
1.0 β₁
1.0 β₂
1.0 β₃
1.0 β₄
julia> a1 # переменная a1
1.0 α₁Теперь вычислим каждый член в уравнении. (\ref{eq:Euler})
julia> # левая часть
lhs1 = a1^2 + a2^2 + a3^2 + a4^2 ;
julia> lhs2 = b1^2 + b2^2 + b3^2 + b4^2 ;
julia> lhs = lhs1 * lhs2;
julia> # правая часть
rhs1 = (a1*b1 - a2*b2 - a3*b3 - a4*b4)^2 ;
julia> rhs2 = (a1*b2 + a2*b1 + a3*b4 - a4*b3)^2 ;
julia> rhs3 = (a1*b3 - a2*b4 + a3*b1 + a4*b2)^2 ;
julia> rhs4 = (a1*b4 + a2*b3 - a3*b2 + a4*b1)^2 ;
julia> rhs = rhs1 + rhs2 + rhs3 + rhs4;Теперь сравним две части тождества
julia> lhs == rhs
trueТождество выполнено.
Тест Фейтмана
Ричард Фейтман (Richard J. Fateman) из Беркли в качестве строгой проверки умножения многочленов предложил вычисление , где . Оно реализовано в функции fateman1 ниже. Мы также рассмотрим форму в fateman2, которая предполагает меньше операций (и более справедливо сравнивает то, что делает система Mathematica).
julia> using TaylorSeries
julia> const order = 20
20
julia> const x, y, z, w = set_variables(Int128, "x", numvars=4, order=2order)
4-element Vector{TaylorN{Int128}}:
1 x₁
1 x₂
1 x₃
1 x₄
julia> function fateman1(degree::Int)
T = Int128
s = one(T) + x + y + z + w
s = s^degree
s * ( s + one(T) )
end
fateman1 (generic function with 1 method)(В следующих строках, которые запускаются при сборке документации, по какой-то причине значения времени появляются перед выполнением команды.)
julia> @time fateman1(0);
0.040507 seconds (54.66 k allocations: 28.621 MiB, 16.04% gc time, 57.50% compilation time)
julia> @time f1 = fateman1(20);
1.502543 seconds (1.23 k allocations: 29.089 MiB, 0.45% gc time)Другая реализация того же самого, но с использованием оптимизаций, связанных с ^2, выдает следующее:
julia> function fateman2(degree::Int)
T = Int128
s = one(T) + x + y + z + w
s = s^degree
s^2 + s
end
fateman2 (generic function with 1 method)
julia> fateman2(0);
julia> @time f2 = fateman2(20); # время появляется выше
0.740009 seconds (1.16 k allocations: 27.014 MiB, 0.90% gc time)Заметим, что приведенные выше функции используют разложения в Int128. На самом деле это необходимо, поскольку некоторые коэффициенты больше typemax(Int):
julia> getcoeff(f2, (1,6,7,20)) # коэффициент x y^6 z^7 w^{20}
128358585324486316800
julia> ans > typemax(Int)
true
julia> length(f2)
41
julia> sum(TaylorSeries.size_table)
135751Эти примеры показывают, что fateman2 почти в два раза быстрее, чем fateman1, и что ряд содержит 135751 одночлен в четырех переменных.
Тесты производительности
Описанные выше функции были сопоставлены с системой Mathematica v11.1. Соответствующие файлы, использовавшиеся для тестов производительности, можно найти здесь. При выполнении на MacPro с процессорами Intel-Xeon 2,7 ГГц мы получили, что системе Mathematica требуется в среднем (5 выполнений) 3,075957 секунды на вычисление, а для fateman1 и fateman2 выше мы получили 2,15408 и 1,08337, соответственно.
Затем, в текущей версии TaylorSeries.jl и с использованием Julia версии 0.7.0 наша реализация fateman1 примерно на 30—40 % быстрее. (Исходный тест Фейтмана соответствует fateman1 выше, что позволяет избежать некоторых оптимизаций, связанных с возведением в квадрат. Реализация в системе Mathematica выполнена таким образом, что эти оптимизации отсутствуют.)