Taylor models tutorial

Download the notebook for this tutorial.


The TaylorModels.jl package can be installed with

julia> using Pkg; Pkg.add("TaylorModels")

Once the package is installed, it can be imported. Note that you will need also the IntervalArithmetic.jl package.

using IntervalArithmetic, TaylorModels


Given a function ff which is continuously derivable n+1n+1 times, its nth-order taylor approximation around the point x0x_0 is given by

Pn(x)=i=0nf(n)(x0)i!(xx0)i,P_n(x) = \sum_{i=0}^n\frac{f^{(n)}(x_0)}{i!}(x-x_0)^i,

where Pn(x)P_n(x) is a polynomial of degree nn. It is well known that Taylor expansion gives a good approximation in a small neighbourhood around x0x_0. However it does not give any information on the approximation error.

TaylorModels.jl uses interval arithmetic to give a rigorous polynomial approximation of ff over the interval [a,b][a, b]. More formally, the rigorous approximation is given by the pair (Pn(x),Δ)(P_n(x), \Delta), where

  • Pn(x)P_n(x) is a polynomial of degree nn, whose coefficients are intervals.

  • Δ\Delta is the interval remainder, and it provides an estimate of all numerical errors (rounding, truncation) encoutered during Taylor expansion computation, i.e. it is guaranteed that x[a,b]: P(x)f(x)Δ\forall x\in[a, b]:~|P(x)-f(x)|\in\Delta.

Taylor models for univariate function

Let us now describe how to get the rigorous Taylor approximation for a function ff in an interval aa around the point x0x_0. Let us define our parameters

f(x) = x*(x-1.1)*(x+2)*(x+2.2)*(x+2.5)*(x+3)*sin(1.7*x+0.5)
a = -0.5..1
x0 = mid(a) # expansion around the middle point of the domain

To create a Taylor model of order nn in the interval aa around the point x0x_0 we need the function TaylorModel1 from TaylorModels.jl, (note the 1 stands for 1 dimensional). The signature of this function is Note that x_0 must be given as an interval, even though it is a single point. Observe the following example, where we create a Taylor approximation of order 6 and 7

tm6 = TaylorModel1(6, interval(x0), a)
tm7 = TaylorModel1(7, interval(x0), a)
@show tm6
@show tm7
tm6 =  [0.25, 0.25] + [1, 1] t + [0, 0]
tm7 =  [0.25, 0.25] + [1, 1] t + [0, 0]

Now to compute the Taylor approximation of our functions we can simply give the previously constructed Taylor model to the function as input

ftm6 = f(tm6)
ftm7 = f(tm7)
@show ftm6
@show ftm7
ftm6 =  [-8.36112, -8.36111] + [-47.0612, -47.0611] t + [-38.3545, -38.3544] t² + [90.4376, 90.4377] t³ + [102.183, 102.184] t⁴ + [-20.5814, -20.5813] t⁵ + [-58.0084, -58.0083] t⁶ + [-1.97165, 2.71142]
ftm7 =  [-8.36112, -8.36111] + [-47.0612, -47.0611] t + [-38.3545, -38.3544] t² + [90.4376, 90.4377] t³ + [102.183, 102.184] t⁴ + [-20.5814, -20.5813] t⁵ + [-58.0084, -58.0083] t⁶ + [-14.3341, -14.334] t⁷ + [-0.281401, 1.03908]

The last term of the expression is the interval remainder and the previous ones represent the Taylor polynomial. The result is an object of type TaylorModel1 with the following attributes

  • dom: inverval over which the Taylor approximation is computed (aa in our example)

  • pol: polynomial approximation of the function

  • rem: interval remainder

  • x0: center of Taylor expansion (x0x_0 in our example)

Now we can visualize how well our taylor model works using the plot function.

using Plots
plot(range(inf(a), stop=sup(a), length=1000), f, lw=2, xaxis="x", yaxis="f(x)", label="f(x)")
plot!(ftm6, label="6th order")
plot!(ftm7, label="7th order")

As you can notice, the Taylor model produces a "band" around the function and the polynomial approximation is guaranteed to be somewhere inside that bar. As you may expect, the higher the order of the model, the narrower the band will be, as the following animation shows

orders = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
anim = @animate for n in orders
    tm = TaylorModel1(n, interval(x0), a)
    ftm = f(tm)
    plot(range(inf(a), stop=sup(a), length=1000), f, lw=2, xaxis="x", yaxis="f(x)",
         label="f(x)", ylims=(-30, 10))
    plot!(ftm, title="$(n)th order")