User guide

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 t + 3 t² + 𝒪(t⁵)

julia> Taylor1([0.0, 1im]) # Also works with complex numbers
 ( 1.0 im ) t + 𝒪(t²)

julia> Taylor1(ones(8), 2) # Polynomial truncated to order 2
 1.0 + 1.0 t + 1.0 t² + 𝒪(t³)

julia> shift_taylor(a) = a + Taylor1(typeof(a),5)  ## a + taylor-polynomial of order 5
shift_taylor (generic function with 1 method)

julia> t = shift_taylor(0.0) # Independent variable `t`
 1.0 t + 𝒪(t⁶)

Note that the information about the maximum order considered is displayed using a big-𝒪 notation. 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 notation
false

julia> t
 1.0 t

julia> displayBigO(true) # turn it on
true

julia> t
 1.0 t + 𝒪(t⁶)

The definition of shift_taylor(a) uses the method Taylor1([::Type{Float64}], [order::Int64=1]), which is a shortcut to define the independent variable of a Taylor expansion, of given type and order (defaults are Float64 and order=1). This is one of the easiest ways to work with 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 of maximum order. This is apparent in the last example below, where the answer is beyond the order of the expansion.

julia> t*(3t+2.5)
 2.5 t + 3.0 t² + 𝒪(t⁶)

julia> 1/(1-t)
 1.0 + 1.0 t + 1.0 t² + 1.0 t³ + 1.0 t⁴ + 1.0 t⁵ + 𝒪(t⁶)

julia> t*(t^2-4)/(t+2)
 - 2.0 t + 1.0 t² + 𝒪(t⁶)

julia> tI = im*t
 ( 1.0 im ) t + 𝒪(t⁶)

julia> (1-t)^3.2
 1.0 - 3.2 t + 3.5200000000000005 t² - 1.4080000000000004 t³ + 0.07040000000000009 t⁴ + 0.011264000000000012 t⁵ + 𝒪(t⁶)

julia> (1+t)^t
 1.0 + 1.0 t² - 0.5 t³ + 0.8333333333333333 t⁴ - 0.75 t⁵ + 𝒪(t⁶)

julia> t^6  # t is of order 5
 0.0 + 𝒪(t⁶)

If no valid Taylor expansion can be computed, an error is thrown, for instance when a derivative is not defined (or simply diverges):

julia> 1/t
ERROR: ArgumentError: Division does not define a Taylor1 polynomial;
order k=0 => coeff[0]=Inf.

julia> t^3.2
ERROR: ArgumentError: The 0th order Taylor1 coefficient must be non-zero
to raise the Taylor1 polynomial to a non-integer exponent.

julia> abs(t)
ERROR: ArgumentError: 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 will be added in the future. Note that this way of obtaining the Taylor coefficients is not the laziest way, in particular for many independent variables. Yet, it is quite efficient, especially for the integration of ordinary differential equations, which is among the applications we have in mind (see also TaylorIntegration.jl).

julia> exp(t)
 1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)

julia> log(1-t)
 - 1.0 t - 0.5 t² - 0.3333333333333333 t³ - 0.25 t⁴ - 0.2 t⁵ + 𝒪(t⁶)

julia> sqrt(1 + t)
 1.0 + 0.5 t - 0.125 t² + 0.0625 t³ - 0.0390625 t⁴ + 0.02734375 t⁵ + 𝒪(t⁶)

julia> imag(exp(tI)')
 - 1.0 t + 0.16666666666666666 t³ - 0.008333333333333333 t⁵ + 𝒪(t⁶)

julia> real(exp(Taylor1([0.0,1im],17))) - cos(Taylor1([0.0,1.0],17)) == 0.0
true

julia> convert(Taylor1{Rational{Int64}}, exp(t))
 1//1 + 1//1 t + 1//2 t² + 1//6 t³ + 1//24 t⁴ + 1//120 t⁵ + 𝒪(t⁶)

Again, errors are thrown whenever it is necessary.

julia> sqrt(t)
ERROR: ArgumentError: First non-vanishing Taylor1 coefficient must correspond
to an **even power** in order to expand `sqrt` around 0.

julia> log(t)
ERROR: ArgumentError: The 0-th order `TaylorN` 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 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)

julia> getcoeff(expon, 0) == expon[0]
true

julia> rationalize(expon[3])
1//6

Differentiating and integrating is straightforward for polynomial expansions in one variable, using derivative and integrate. These functions return the corresponding Taylor1 expansions. The last coefficient of a derivative is set to zero to keep the same order as the original polynomial; for the integral, an integration constant may be set by the user (the default is zero). The order of the resulting polynomial is not changed. The value of the $n$-th ($n \ge 0$) derivative is obtained using derivative(n,a), where a is a Taylor series.

julia> derivative(exp(t))
 1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 𝒪(t⁶)

julia> integrate(exp(t))
 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)

julia> integrate( exp(t), 1.0)
 1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)

julia> integrate( derivative( exp(-t)), 1.0 ) == exp(-t)
true

julia> derivative(1, exp(shift_taylor(1.0))) == exp(1.0)
true

julia> derivative(5, exp(shift_taylor(1.0))) == exp(1.0) # 5-th derivative of `exp(1+t)`
true

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))) - e # exp(t) around t0=1 (order 5), evaluated there (dt=0)
0.0

julia> evaluate(exp(t), 1) - e # exp(t) around t0=0 (order 5), evaluated at t=1
-0.0016151617923783057

julia> evaluate(exp( Taylor1(17) ), 1) - e # exp(t) around t0=0, order 17
0.0

julia> tBig = Taylor1(BigFloat, 50) # Independent variable with BigFloats, order 50
 1.000000000000000000000000000000000000000000000000000000000000000000000000000000 t + 𝒪(t⁵¹)

julia> eBig = evaluate( exp(tBig), one(BigFloat) )
2.718281828459045235360287471352662497757247093699959574966967627723419298053556

julia> e - eBig
6.573322999985292556154129119543257102601105719980995128942636339920549561322098e-67

Another way to obtain the value of a Taylor1 polynomial p at a given value x, is to call p as if it were a function, i.e., p(x):

julia> t = Taylor1(15)
 1.0 t + 𝒪(t¹⁶)

julia> p = sin(t)
 1.0 t - 0.16666666666666666 t³ + 0.008333333333333333 t⁵ - 0.0001984126984126984 t⁷ + 2.7557319223985893e-6 t⁹ - 2.505210838544172e-8 t¹¹ + 1.6059043836821616e-10 t¹³ - 7.647163731819817e-13 t¹⁵ + 𝒪(t¹⁶)

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 function
0.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). For more details about function-like behavior for a given type in Julia, see the Function-like objects section of the Julia manual.

Useful shortcuts are taylor_expand are update!. The former returns the expansion of a function around a given value t0. In turn, update! provides an in-place update of a given Taylor polynomial, that is, it shifts it further by the provided amount.

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 t - 0.5 t² - 1.020538999289461e-17 t³ + 0.041666666666666664 t⁴ + 5.102694996447305e-19 t⁵ - 0.001388888888888889 t⁶ - 1.2149273801065013e-20 t⁷ + 2.48015873015873e-5 t⁸ + 1.6873991390368074e-22 t⁹ - 2.7557319223985894e-7 t¹⁰ - 1.5339992173061888e-24 t¹¹ + 2.08767569878681e-9 t¹² + 9.833328316065313e-27 t¹³ - 1.1470745597729726e-11 t¹⁴ - 4.682537293364434e-29 t¹⁵ + 4.779477332387386e-14 t¹⁶ + 𝒪(t¹⁷)

julia> update!(p, 0.025) # updates the expansion given by p, by shifting it further by 0.025

julia> p
 0.9996875162757026 - 0.02499739591471227 t - 0.49984375813785126 t² + 0.004166232652452044 t³ + 0.0416536465114876 t⁴ - 0.00020831163262260228 t⁵ - 0.0013884548837162537 t⁶ + 4.959800776728625e-6 t⁷ + 2.4793837209218816e-5 t⁸ - 6.888612189900868e-8 t⁹ - 2.7548708010243125e-7 t¹⁰ + 6.262374718092001e-10 t¹¹ + 2.0870233341100357e-9 t¹² - 4.014342754938811e-12 t¹³ - 1.1467160989730436e-11 t¹⁴ + 1.9117909329549495e-14 t¹⁵ + 4.779477332387386e-14 t¹⁶ + 𝒪(t¹⁷)

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 a 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 Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 x + 𝒪(‖x‖⁷)
  1.0 y + 𝒪(‖x‖⁷)

julia> typeof(x)
TaylorSeries.TaylorN{Float64}

julia> x.order
6

julia> x.coeffs
7-element Array{TaylorSeries.HomogeneousPolynomial{Float64},1}:
    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 Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 x + 𝒪(‖x‖¹¹)
  1.0 y + 𝒪(‖x‖¹¹)

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

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

Alternatively to set_variables, get_variables can be used if one doesn't 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(order=2) #vector of independent variables of order 2
3-element Array{TaylorSeries.TaylorN{Float64},1}:
  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      = String["α₁", "α₂", "α₃"]
Variable symbols    = Symbol[:α₁, :α₂, :α₃]

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 numbers. (Some of the arithmetic operations have been extended for HomogeneousPolynomial, whenever the result is a HomogeneousPolynomial; division, for instance, is not extended.)

Also, the elementary functions have been implemented, again by computing their coefficients recursively:

julia> x, y = set_variables("x y", order=10);

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.0013888888888888887 x⁶ + 0.008333333333333331 x⁵ y + 0.020833333333333332 x⁴ y² + 0.027777777777777776 x³ y³ + 0.020833333333333332 x² y⁴ + 0.008333333333333331 x y⁵ + 0.0013888888888888887 y⁶ + 0.00019841269841269839 x⁷ + 0.0013888888888888885 x⁶ y + 0.004166666666666666 x⁵ y² + 0.006944444444444443 x⁴ y³ + 0.006944444444444443 x³ y⁴ + 0.004166666666666666 x² y⁵ + 0.0013888888888888885 x y⁶ + 0.00019841269841269839 y⁷ + 2.4801587301587298e-5 x⁸ + 0.00019841269841269836 x⁷ y + 0.0006944444444444443 x⁶ y² + 0.0013888888888888887 x⁵ y³ + 0.0017361111111111108 x⁴ y⁴ + 0.0013888888888888887 x³ y⁵ + 0.0006944444444444443 x² y⁶ + 0.00019841269841269836 x y⁷ + 2.4801587301587298e-5 y⁸ + 2.7557319223985884e-6 x⁹ + 2.4801587301587295e-5 x⁸ y + 9.920634920634918e-5 x⁷ y² + 0.0002314814814814814 x⁶ y³ + 0.0003472222222222221 x⁵ y⁴ + 0.0003472222222222221 x⁴ y⁵ + 0.0002314814814814814 x³ y⁶ + 9.920634920634918e-5 x² y⁷ + 2.4801587301587295e-5 x y⁸ + 2.7557319223985884e-6 y⁹ + 2.7557319223985883e-7 x¹⁰ + 2.7557319223985884e-6 x⁹ y + 1.2400793650793647e-5 x⁸ y² + 3.306878306878306e-5 x⁷ y³ + 5.787037037037036e-5 x⁶ y⁴ + 6.944444444444443e-5 x⁵ y⁵ + 5.787037037037036e-5 x⁴ y⁶ + 3.306878306878306e-5 x³ y⁷ + 1.2400793650793647e-5 x² y⁸ + 2.7557319223985884e-6 x y⁹ + 2.7557319223985883e-7 y¹⁰ + 𝒪(‖x‖¹¹)

The function getcoeff gives the normalized coefficient of the polynomial that corresponds to the monomial specified by a vector v containing the powers. For instance, for the polynomial exy above, the coefficient of the monomial $x^3 y^5$ is

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.4801587301587298e-5 x⁸ + 0.00019841269841269836 x⁷ y + 0.0006944444444444443 x⁶ y² + 0.0013888888888888887 x⁵ y³ + 0.0017361111111111108 x⁴ y⁴ + 0.0013888888888888887 x³ y⁵ + 0.0006944444444444443 x² y⁶ + 0.00019841269841269836 x y⁷ + 2.4801587301587298e-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 term
0.0013888888888888887

Partial differentiation is also implemented for TaylorN objects, through the function derivative, 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> derivative( p, 1 )   # partial derivative with respect to 1st variable
 - 7.0 + 3.0 x² + 4.0 x y + 𝒪(‖x‖¹¹)

julia> derivative( 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> derivative( q, 3 )   # error, since we are dealing with 2 variables
ERROR: AssertionError: 1 ≤ r ≤ get_numvars()

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( derivative( p, 1 ), 1) # integrate with respect to the first variable
 - 7.0 x + 1.0 x³ + 2.0 x² y + 𝒪(‖x‖¹¹)

julia> integrate( derivative( 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( derivative( q, 2 ), :y, -x^4) == q
true

julia> integrate( derivative( 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 -> 0
true

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.0013888888888888887 y⁶ + 0.00019841269841269839 y⁷ + 2.4801587301587298e-5 y⁸ + 2.7557319223985884e-6 y⁹ + 2.7557319223985883e-7 y¹⁰ + 𝒪(‖x‖¹¹)

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^2
true

Functions to compute the gradient, Jacobian and Hessian have also been implemented. Using the polynomials $p = x^3 + 2x^2 y - 7 x + 2$ and $q = y-x^4$ defined above, we may use 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 jacobian and hessian:

julia> ∇(p)
2-element Array{TaylorSeries.TaylorN{Float64},1}:
  - 7.0 + 3.0 x² + 4.0 x y + 𝒪(‖x‖¹¹)
                    2.0 x² + 𝒪(‖x‖¹¹)

julia> gradient( q )
2-element Array{TaylorSeries.TaylorN{Float64},1}:
  - 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> hessian(ans)
2×2 Array{Float64,2}:
  0.0  14.0
 14.0   0.0

julia> jacobian([p,q], [2,1])
2×2 Array{Float64,2}:
  13.0  8.0
 -32.0  1.0

julia> hessian(r, [1.0,1.0])
2×2 Array{Float64,2}:
 -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 Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 x + 𝒪(‖x‖⁴)
  1.0 y + 𝒪(‖x‖⁴)

julia> t1N = Taylor1([zero(x), one(x)], 5)
 ( 1.0 + 𝒪(‖x‖⁴)) t + 𝒪(t⁶)

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‖⁴)) t + ( 0.2524230522999288 + 0.43160468332443686 x - 0.1262115261499644 x² - 0.07193411388740614 x³ + 𝒪(‖x‖⁴)) t² + ( 0.14386822777481229 - 0.0841410174333096 x - 0.07193411388740614 x² + 0.0140235029055516 x³ + 𝒪(‖x‖⁴)) t³ + ( - 0.0210352543583274 - 0.03596705694370307 x + 0.0105176271791637 x² + 0.005994509490617178 x³ + 𝒪(‖x‖⁴)) t⁴ + ( - 0.007193411388740615 + 0.00420705087166548 x + 0.0035967056943703073 x² - 0.00070117514527758 x³ + 𝒪(‖x‖⁴)) t⁵ + 𝒪(t⁶)

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.