# Examples

## Expanding exp(x) with `taylor_expand()`

The `taylor_expand`

function takes the function to expand as it's first argument, and the point to expand about as the second argument. A keyword argument `order`

determines which order to expand to:

`julia> using TaylorSeries`

`julia> taylor_expand(exp, 0, order=2)`

`1.0 + 1.0 t + 0.5 t² + 𝒪(t³)`

And voìla! It really is that simple to calculate a simple taylor polynomial. The next example is slightly more complicated.

## Expanding exp(x) with a symbolic object

An alternative way to compute the single-variable taylor expansion for a function is by defining a variable of type `Taylor1`

, and using it in the function you wish to expand. The argument given to the `Taylor1`

constructor is the order to expand to:

`julia> using TaylorSeries`

`julia> x = Taylor1(2)`

`1.0 t + 𝒪(t³)`

`julia> exp(x)`

`1.0 + 1.0 t + 0.5 t² + 𝒪(t³)`

Let's also get rid of the printed error for the next few examples, and set the printed independent variable to `x`

:

`julia> displayBigO(false)`

`false`

`julia> set_taylor1_varname("x")`

`"x"`

`julia> exp(x)`

`1.0 + 1.0 x + 0.5 x²`

### Changing point to expand about

A variable constructed with `Taylor1()`

automatically expands about the point `x=0`

. But what if you want to use the symbolic object to expand about a point different from zero? Because expanding `exp(x)`

about `x=1`

is exactly the same as expanding `exp(x+1)`

about `x=0`

, simply replace the `x`

in your expression with `x+1`

to expand about `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`

### More examples

You can even use custum functions

`julia> f(a) = 1/(a+1)`

`f (generic function with 1 method)`

`julia> f(x)`

`1.0 - 1.0 x + 1.0 x²`

Functions can be nested

`julia> sin(f(x))`

`0.8414709848078965 - 0.5403023058681398 x + 0.11956681346419151 x²`

and complicated further in a modular way

`julia> sin(exp(x+2))/(x+2)+cos(x+2)+f(x+2)`

`0.364113974242596 + 0.41259243488717107 x - 11.843864409375039 x²`

## Four-square identity

The first example shows that the four-square identity holds:

\[(a_1^2 + a_2^2 + a_3^2 + a_4^2)\cdot(b_1^2 + b_2^2 + b_3^2 + b_4^2) = \\ \qquad (a_1 b_1 - a_2 b_2 - a_3 b_3 -a_4 b_4)^2 + (a_1 b_2 - a_2 b_1 - a_3 b_4 -a_4 b_3)^2 + \\ \qquad (a_1 b_3 - a_2 b_4 - a_3 b_1 -a_4 b_2)^2 + (a_1 b_4 - a_2 b_3 - a_3 b_2 -a_4 b_1)^2,\]

which was originally proved by Euler. The code can also be found in this test of the package.

First, we reset the maximum degree of the polynomial to 4, since the RHS of the equation has *a priori* terms of fourth order, and define the 8 independent variables.

`julia> using TaylorSeries`

`julia> # Define the variables α₁, ..., α₄, β₁, ..., β₄ 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> # Create the 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 # variable a1`

`1.0 α₁`

Now we compute each term appearing in Eq. (\ref{eq:Euler})

`julia> # left-hand side lhs1 = a1^2 + a2^2 + a3^2 + a4^2 ;`

`julia> lhs2 = b1^2 + b2^2 + b3^2 + b4^2 ;`

`julia> lhs = lhs1 * lhs2;`

`julia> # right-hand side 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;`

We now compare the two sides of the identity,

`julia> lhs == rhs`

`true`

The identity is satisfied. $\square$

## Fateman test

Richard J. Fateman, from Berkeley, proposed as a stringent test of polynomial multiplication the evaluation of $s\cdot(s+1)$, where $s = (1+x+y+z+w)^{20}$. This is implemented in the function `fateman1`

below. We shall also consider the form $s^2+s$ in `fateman2`

, which involves fewer operations (and makes a fairer comparison to what Mathematica does).

`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)`

(In the following lines, which are run when the documentation is built, by some reason the timing appears before the command executed.)

`julia> @time fateman1(0);`

`0.080414 seconds (54.86 k allocations: 41.065 MiB, 16.34% gc time, 50.65% compilation time)`

`julia> @time f1 = fateman1(20);`

`2.692494 seconds (1.86 k allocations: 49.837 MiB, 0.40% gc time)`

Another implementation of the same, but exploiting optimizations related to `^2`

yields:

`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); # the timing appears above`

`1.396940 seconds (1.99 k allocations: 53.986 MiB, 0.38% gc time)`

We note that the above functions use expansions in `Int128`

. This is actually required, since some coefficients are larger than `typemax(Int)`

:

`julia> getcoeff(f2, (1,6,7,20)) # coefficient of x y^6 z^7 w^{20}`

`128358585324486316800`

`julia> ans > typemax(Int)`

`true`

`julia> length(f2)`

`41`

`julia> sum(TaylorSeries.size_table)`

`135751`

These examples show that `fateman2`

is nearly twice as fast as `fateman1`

, and that the series has 135751 monomials in 4 variables.

### Benchmarks

The functions described above have been compared against Mathematica v11.1. The relevant files used for benchmarking can be found here. Running on a MacPro with Intel-Xeon processors 2.7GHz, we obtain that Mathematica requires on average (5 runs) 3.075957 seconds for the computation, while for `fateman1`

and `fateman2`

above we obtain 2.15408 and 1.08337, respectively.

Then, with the current version of `TaylorSeries.jl`

and using Julia v0.7.0, our implementation of `fateman1`

is about 30%-40% faster. (The original test by Fateman corresponds to `fateman1`

above, which avoids some optimizations related to squaring; the implementation in Mathematica is done such that this optimization does not occur.)