Примеры

Разложение 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 выполнена таким образом, что эти оптимизации отсутствуют.)