User guide

For some simple examples, head over to the examples section. For a detailed guide, keep reading.

TaylorSeries.jl is a basic polynomial algebraic manipulator in one or more variables; these two cases are treated separately. Three new types are defined, Taylor1, HomogeneousPolynomial and TaylorN, which correspond to expansions in one independent variable, homogeneous polynomials of various variables, and the polynomial series in many independent variables, respectively. These types are subtypes of AbstractSeries, which in turn is a subtype of Number, and are defined parametrically.

The package is loaded as usual:

julia> using TaylorSeries

One independent variable

Taylor expansions in one variable are represented by the Taylor1 type, which consists of a vector of coefficients (fieldname coeffs) and the maximum order considered for the expansion (fieldname order). The coefficients are arranged in ascending order with respect to the degree of the monomial, so that coeffs[1] is the constant term, coeffs[2] gives the first order term (t^1), etc. Yet, it is possible to have the natural ordering with respect to the degree; see below. This is a dense representation of the polynomial. The order of the polynomial can be omitted in the constructor, which is then fixed by the length of the vector of coefficients. If the length of the vector does not correspond with the order, order is used, which effectively truncates polynomial to degree order.

julia> Taylor1([1, 2, 3],4) # Polynomial of order 4 with coefficients 1, 2, 3 1 + 2 x + 3 x²
julia> Taylor1([0.0, 1im]) # Also works with complex numbers ( 0.0 + 1.0im ) x
julia> Taylor1(ones(8), 2) # Polynomial truncated to order 2 1.0 + 1.0 x + 1.0 x²
julia> shift_taylor(a) = a + Taylor1(typeof(a),5)  ## a + taylor-polynomial of order 5shift_taylor (generic function with 1 method)
julia> t = shift_taylor(0.0) # Independent variable `t` 1.0 x

In some cases, it is desirable to not display the big-𝒪 notation. The function displayBigO allows to control whether it is displayed or not.

julia> displayBigO(false) # turn-off displaying big O notationfalse
julia> t 1.0 x
julia> displayBigO(true) # turn it ontrue
julia> t 1.0 x + 𝒪(x⁶)

Similarly, it is possible to control some aspects of the format of the displayed series through the function use_show_default; use_show_default(true) uses the Base.show_default, while use_show_default(false) uses the custom display form (default).

julia> use_show_default(true) # use Base.show methodtrue
julia> tTaylor1{Float64}([0.0, 1.0, 0.0, 0.0, 0.0, 0.0], 5)
julia> use_show_default(false) # use custom `show`false
julia> t 1.0 x + 𝒪(x⁶)

The definition of shift_taylor(a) uses the method Taylor1([::Type{Float64}], order::Int), which is a shortcut to define the independent variable of a Taylor expansion, of given type and order (the default is Float64). Defining the independent variable in advance is one of the easiest ways to use the package.

The usual arithmetic operators (+, -, *, /, ^, ==) have been extended to work with the Taylor1 type, including promotions that involve Numbers. The operations return a valid Taylor expansion, consistent with the order of the series. This is apparent in the penultimate example below, where the fist non-zero coefficient is beyond the order of the expansion, and hence the result is zero.

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-𝒪 convention appliestrue
julia> t^6  # t is of order 5 0.0 + 𝒪(x⁶)
julia> t^2 / t # The result is of order 4, instead of 5 1.0 x + 𝒪(x⁵)

Note that the last example returns a Taylor1 series of order 4, instead of order 5; this is be consistent with the number of known coefficients of the returned series, since the result corresponds to factorize t in the numerator to cancel the same factor in the denominator.

Taylor1 is also equiped with a total order, provided by overloading isless. The ordering is consistent with the interpretation that there are infinitessimal elements in the algebra; for details see M. Berz, "Automatic Differentiation as Nonarchimedean Analysis", Computer Arithmetic and Enclosure Methods, (1992), Elsevier, 439-450. This is illustrated by:

julia> 1 > t > 2*t^2 > 100*t^3 >= 0true

If no valid Taylor expansion can be computed an error is thrown, for instance, when a derivative is not defined at the expansion point, or it simply diverges.

julia> 1/tERROR: ArgumentError: Division does not define a Taylor1 polynomial;
order k=0 => coeff[0]=Inf.
julia> t^3.2ERROR: 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).

Several elementary functions have been implemented; their coefficients are computed recursively. At the moment of this writing, these functions are exp, log, sqrt, the trigonometric functions sin, cos and tan, their inverses, as well as the hyperbolic functions sinh, cosh and tanh and their inverses; more functions will be added in the future. Note that this way of obtaining the Taylor coefficients is not a lazy way, in particular for many independent variables. Yet, the implementation is efficient enough, especially for the integration of ordinary differential equations, which is among the applications we have in mind (see 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.0true
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⁶)

Again, errors are thrown whenever it is necessary.

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.

To obtain a specific coefficient, getcoeff can be used. Another alternative is to request the specific degree using the vector notation, where the index corresponds to the degree of the term.

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

Note that certain arithmetic operations, or the application of some functions, may change the order of the result, as mentioned above.

julia> t #  order 5 independent variable 1.0 x + 𝒪(x⁶)
julia> t^2/t # returns an order 4 series expansion 1.0 x + 𝒪(x⁵)
julia> sqrt(t^2) # returns an order 2 series expansion 1.0 x + 𝒪(x³)
julia> (t^4)^(1/4) # order 1 series 1.0 x + 𝒪(x²)

Differentiating and integrating is straightforward for polynomial expansions in one variable, using differentiate and integrate. (The function derivative is a synonym of differentiate.) These functions return the corresponding Taylor1 expansions. The order of the derivative of a Taylor1 is reduced by 1. For the integral, an integration constant may be set by the user (the default is zero); the order of the integrated polynomial for the integral is kept unchanged. The value of the $n$-th ($n \ge 0$) derivative is obtained using differentiate(n,a), where a is a Taylor series; likewise, the Taylor1 polynomial of the $n$-th derivative is obtained as differentiate(a,n); the resulting polynomial is of order get_order(a)-n.

julia> differentiate(exp(t)) # exp(t) is of order 5; the derivative is of order 4 1.0 + 1.0 x + 0.5 x² + 0.16666666666666666 x³ + 0.041666666666666664 x⁴ + 𝒪(x⁵)
julia> integrate(exp(t))  # the resulting TaylorSeries is of order 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-th differentiate of `exp(1+t)`true
julia> derivative(exp(1+t), 3)    # Taylor1 polynomial of the 3-rd derivative of `exp(1+t)` 2.718281828459045 + 2.718281828459045 x + 1.3591409142295225 x² + 𝒪(x³)

We also have methods to invert a Taylor1 polynomial, using either inverse, which uses Lagrange inversion, or inverse_map, which uses an algorithm developed by M. Berz; the latter can also be used for TaylorN polynomials; see below.

To evaluate a Taylor series at a given point, Horner's rule is used via the function evaluate(a, dt). Here, dt is the increment from the point $t_0$ around which the Taylor expansion of a is calculated, i.e., the series is evaluated at $t = t_0 + dt$. Omitting dt corresponds to $dt = 0$; see evaluate.

julia> evaluate(exp(shift_taylor(1.0))) - ℯ    # exp(t) around t0=1 (order 5), evaluated there (dt=0)0.0
julia> evaluate(exp(t), 1) - ℯ                 # exp(t) around t0=0 (order 5), evaluated at t=1-0.0016151617923783057
julia> evaluate(exp( Taylor1(17) ), 1) - ℯ     # exp(t) around t0=0, order 170.0
julia> tBig = Taylor1(BigFloat, 50)            # Independent variable with BigFloats, order 50 1.0 x + 𝒪(x⁵¹)
julia> eBig = evaluate( exp(tBig), one(BigFloat) )2.718281828459045235360287471352662497757247093699959574966967627723419298053556
julia> ℯ - eBig6.573322999985292556154129119543257102601105719980995128942636339920549561322098e-67

Another way to evaluate the value of a Taylor1 polynomial p at a given value x, is to call p as if it was a function, i.e., 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) # value of p at pi/2 using `evaluate`0.9999999999939766
julia> p(pi/2) # value of p at pi/2 by evaluating p as a function0.9999999999939766
julia> p(pi/2) == evaluate(p, pi/2)true
julia> p(0.0)0.0
julia> p() == p(0.0) # p() is a shortcut to obtain the 0-th order coefficient of `p`true

Note that the syntax p(x) is equivalent to evaluate(p, x), whereas p() is equivalent to evaluate(p).

Useful shortcuts are taylor_expand and update!. The former returns the expansion of a function around a given value t0, mimicking the use of shift_taylor above. In turn, update! provides an in-place update of a given Taylor polynomial, that is, it shifts it further by the provided amount. Note that the type of the Taylor1 polynomial and the shifted point must be compatible, in the sense that the latter must be convertable into the parametric type of the former.

julia> p = taylor_expand( x -> sin(x), pi/2, order=16) # 16-th order expansion of sin(t) around 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) # updates the expansion given by p, by shifting it further by 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¹⁷)

Many variables

A polynomial in $N>1$ variables can be represented in (at least) two ways: As a vector whose coefficients are homogeneous polynomials of fixed degree, or as a vector whose coefficients are polynomials in $N-1$ variables. The current implementation of TaylorSeries.jl corresponds to the first option, though some infrastructure has been built that permits to develop the second one. An elegant (lazy) implementation of the second representation was discussed here.

The structure TaylorN is constructed as a vector of parameterized homogeneous polynomials defined by the type HomogeneousPolynomial, which in turn is an ordered vector of coefficients of given order (degree). This implementation imposes the user to specify the (maximum) order considered and the number of independent variables at the beginning, which can be conveniently done using set_variables. A vector of the resulting Taylor variables is returned:

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.order6
julia> x.coeffs7-element Vector{HomogeneousPolynomial{Float64}}:
    0.0
  1.0 x
    0.0
    0.0
    0.0
    0.0
    0.0

As shown, the resulting objects are of TaylorN{Float64} type. There is an optional order keyword argument in set_variables, used to specify the maximum order of the TaylorN polynomials. Note that one can specify the variables using a vector of symbols.

julia> set_variables([:x, :y], order=10)2-element Vector{TaylorN{Float64}}:
  1.0 x + 𝒪(‖x‖¹¹)
  1.0 y + 𝒪(‖x‖¹¹)

Similarly, subindexed variables are also available by specifying a single variable name and the optional keyword argument numvars:

julia> set_variables("α", numvars=3)3-element Vector{TaylorN{Float64}}:
  1.0 α₁ + 𝒪(‖x‖¹¹)
  1.0 α₂ + 𝒪(‖x‖¹¹)
  1.0 α₃ + 𝒪(‖x‖¹¹)

Alternatively to set_variables, get_variables can be used if one does not want to change internal dictionaries. get_variables returns a vector of TaylorN independent variables of a desired order (lesser than get_order so the internals doesn't have to change) with the length and variable names defined by set_variables initially.

julia> get_variables(2) # vector of independent variables of order 23-element Vector{TaylorN{Float64}}:
  1.0 α₁ + 𝒪(‖x‖³)
  1.0 α₂ + 𝒪(‖x‖³)
  1.0 α₃ + 𝒪(‖x‖³)

The function show_params_TaylorN displays the current values of the parameters, in an info block.

julia> show_params_TaylorN()┌ Info: Parameters for `TaylorN` and `HomogeneousPolynomial`:
│ Maximum order       = 10
│ Number of variables = 3
│ Variable names      = ["α₁", "α₂", "α₃"]
└ Variable symbols    = [:α₁, :α₂, :α₃]

Internally, changing the parameters (maximum order and number of variables) redefines the hash-tables that translate the index of the coefficients of a HomogeneousPolynomial of given order into the corresponding multi-variable monomials, or the other way around. Fixing these values from the start is imperative; the initial (default) values are order = 6 and num_vars=2.

The easiest way to construct a TaylorN object is by defining the independent variables. This can be done using set_variables as above, or through the method TaylorN{T<:Number}(::Type{T}, nv::Int) for the nv independent TaylorN{T} variable; the order can be also specified using the optional keyword argument order.

julia> x, y = set_variables("x y", numvars=2, order=6);
julia> x 1.0 x + 𝒪(‖x‖⁷)
julia> TaylorN(1, order=4) # variable 1 of order 4 1.0 x + 𝒪(‖x‖⁵)
julia> TaylorN(Int, 2)    # variable 2, type Int, order=get_order()=6 1 y + 𝒪(‖x‖⁷)

Other ways of constructing TaylorN polynomials involve using HomogeneousPolynomial objects directly, which is uncomfortable.

julia> set_variables(:x, numvars=2); # symbols can be used
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‖⁵)

The Taylor expansions are implemented around 0 for all variables; if the expansion is needed around a different value, the trick is a simple translation of the corresponding independent variable, i.e. $x \to x+a$.

As before, the usual arithmetic operators (+, -, *, /, ^, ==) have been extended to work with TaylorN objects, including the appropriate promotions to deal with the usual numeric types. Note that some of the arithmetic operations have been extended for HomogeneousPolynomial, whenever the result is a HomogeneousPolynomial; division, for instance, is not extended. The same convention used for Taylor1 objects is used when combining TaylorN polynomials of different order. Both HomogeneousPolynomial and TaylorN are equiped with a total lexicographical order, provided by overloading isless. The ordering is consistent with the interpretation that there are infinitessimal elements in the algebra; for details see M. Berz, "Automatic Differentiation as Nonarchimedean Analysis", Computer Arithmetic and Enclosure Methods, (1992), Elsevier, 439-450. Essentially, the lexicographic order works as follows: smaller order monomials are larger than higher order monomials; when the order is the same, larger monomials appear before in the hash-tables; the function show_monomials.

julia> x, y = set_variables("x y", order=10);
julia> show_monomials(2) 1  -->  x²
 2  -->  x y
 3  -->  y²

Then, the following then holds:

julia> 0 < 1e8 * y^2 < x*y < x^2  < y < x/1e8 < 1.0true

The elementary functions have also been implemented, again by computing their coefficients recursively:

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‖¹¹)

The function getcoeff gives the normalized coefficient of the polynomial that corresponds to the monomial specified by the tuple or vector v containing the powers. For instance, for the polynomial exy above, the coefficient of the monomial $x^3 y^5$ is obtained using getcoeff(exy, (3,5)) or getcoeff(exy, [3,5]).

julia> getcoeff(exy, (3,5))0.0013888888888888887
julia> rationalize(ans)1//720

Similar to Taylor1, vector notation can be used to request specific coefficients of HomogeneousPolynomial or TaylorN objects. For TaylorN objects, the index refers to the degree of the HomogeneousPolynomial. In the case of HomogeneousPolynomial the index refers to the position of the hash table. The function show_monomials can be used to obtain the coefficient a specific monomial, given the degree of the HomogeneousPolynomial.

julia> exy[8] # get the 8th order term 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] # get the 6th coeff of the 8th order term0.0013888888888888887

Partial differentiation is also implemented for TaylorN objects, through the function differentiate, specifying the number of the variable, or its symbol, as the second argument.

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 )   # partial derivative with respect to 1st variable - 7.0 + 3.0 x² + 4.0 x y + 𝒪(‖x‖¹¹)
julia> differentiate( q, :y )  # partial derivative with respect to :y 1.0 + 𝒪(‖x‖¹¹)

If we ask for the partial derivative with respect to a non-defined variable, an error is thrown.

julia> differentiate( q, 3 )   # error, since we are dealing with 2 variablesERROR: AssertionError: 1 ≤ r ≤ get_numvars()

To obtain more specific partial derivatives we have two specialized methods that involve a tuple, which represents the number of derivatives with respect to each variable (so the tuple's length has to be the same as the actual number of variables). These methods either return the TaylorN object in question, or the coefficient corresponding to the specified tuple, normalized by the factorials defined by the tuple. The latter is in essence the 0-th order coefficient of the former.

julia> differentiate(p, (2,1)) # two derivatives on :x and one on :y 4.0 + 𝒪(‖x‖¹¹)
julia> differentiate((2,1), p) # 0-th order coefficient of the previous expression4.0
julia> differentiate(p, (1,1)) # one derivative on :x and one on :y 4.0 x + 𝒪(‖x‖¹¹)
julia> differentiate((1,1), p) # 0-th order coefficient of the previous expression0.0

Integration with respect to the r-th variable for HomogeneousPolynomials and TaylorN objects is obtained using integrate. Note that integrate for TaylorN objects allows to specify a constant of integration, which must be independent from the integrated variable. Again, the integration variable may be specified by its symbol.

julia> integrate( differentiate( p, 1 ), 1) # integrate with respect to the first variable - 7.0 x + 1.0 x³ + 2.0 x² y + 𝒪(‖x‖¹¹)
julia> integrate( differentiate( p, 1 ), :x, 2) # integration with respect to :x, constant of integration is 2 2.0 - 7.0 x + 1.0 x³ + 2.0 x² y + 𝒪(‖x‖¹¹)
julia> integrate( differentiate( q, 2 ), :y, -x^4) == qtrue
julia> integrate( differentiate( q, 2 ), 2, y)ERROR: AssertionError: The integration constant ( 1.0 y + 𝒪(‖x‖¹¹)) must be independent of the
y variable

evaluate can also be used for TaylorN objects, using it on vectors of numbers (Real or Complex); the length of the vector must coincide with the number of independent variables. evaluate also allows to specify only one variable and a value.

julia> evaluate(exy, [.1,.02]) == exp(0.12)true
julia> evaluate(exy, :x, 0.0) == exp(y)  # evaluate `exy` for :x -> 0true

Analogously to Taylor1, another way to obtain the value of a TaylorN polynomial p at a given point x, is to call it as if it were a function: the syntax p(x) for p::TaylorN is equivalent to evaluate(p,x), and p() is equivalent to 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‖¹¹)

Internally, evaluate for TaylorN considers separately the contributions of all HomogeneousPolynomials by order, which are finally added up after sorting them in place (which is the default) in increasing order by abs2. This is done in order to use as many significant figures as possible of all terms in the final sum, which then should yield a more accurate result. This default can be changed to a non-sorting sum thought, which may be more performant or useful for certain subtypes of Number which, for instance, do not have isless defined. See this issue for a motivating example. This can be done using the keyword sorting in evaluate, which expects a Bool, or using a that boolean as the first argument in the function-like evaluation.

julia> exy([.1,.02]) # default is `sorting=true`1.1274968515793757
julia> evaluate(exy, [.1,.02]; sorting=false)1.127496851579376
julia> exy(false, [.1,.02])1.127496851579376

In the examples shown above, the first entry corresponds to the default case (sorting=true), which yields the same result as exp(0.12), and the remaining two illustrate turning off sorting the terms. Note that the results are not identical, since floating point addition is not associative, which may introduce rounding errors.

The functions taylor_expand and update! work as well for TaylorN.

julia> xysq = x^2 + y^2 1.0 x² + 1.0 y² + 𝒪(‖x‖¹¹)
julia> update!(xysq, [1.0, -2.0]) # expand around (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]) # shift-back
julia> xysq == x^2 + y^2true

Functions to compute the gradient, Jacobian and Hessian have also been implemented; note that these functions are not exported, so its use require the prefix TaylorSeries. Using the polynomials $p = x^3 + 2x^2 y - 7 x + 2$ and $q = y-x^4$ defined above, we may use TaylorSeries.gradient (or ∇); the results are of type Array{TaylorN{T},1}. To compute the Jacobian and Hessian of a vector field evaluated at a point, we use respectively TaylorSeries.jacobian and 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

Other specific applications are described in the Examples.

Mixtures

As mentioned above, Taylor1{T}, HomogeneousPolynomial{T} and TaylorN{T} are parameterized structures such that T<:AbstractSeries, the latter is a subtype of Number. Then, we may actually define Taylor expansions in $N+1$ variables, where one of the variables (the Taylor1 variable) is somewhat special.

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

The last line defines a Taylor1{TaylorN{Float64}} variable, which is of order 5 in t and order 3 in x and y. Then, we can evaluate functions involving such polynomials:

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

This kind of expansions are of interest when studying the dependence of parameters, for instance in the context of bifurcation theory or when considering the dependence of the solution of a differential equation on the initial conditions, around a given solution. In this case, x and y represent small variations around a given value of the parameters, or around some specific initial condition. Such constructions are exploited in the package TaylorIntegration.jl.