2834-Taylor Polynomials.sagews

Taylor Polynomials

Since Taylor polynomials concentrate their accuracy around a point (making sure the function and polynomial have the same derivatives up to the $n$th derivative at that point), a Taylor polynomial may differ from a function significantly away from the point.

In the example below, we calculate and plot the Taylor polynomial $\hat{f}$ approximating the function $f(x)=1/x$ around the point $x=1$.  We also print the error in the approximation at $x=3$, which is $|f(3)-\hat{f}(3)|$.  Higher order Taylor polynomials usually give better approximations farther from the point, but in this case, as we increase the order of the approximation, the error gets much, much worse.  Check this out by moving the slider to change the order of the approximation.

x0  = 1
f(x)= 1/x
p   = plot(f,-1,5, thickness=2)
dot = point((x0,f(x=x0)),pointsize=80,color='red')
@interact
def _(order=(1..12)):
  taylor = f.taylor(x,x0,order)
  taylor_plot = plot(taylor,-1, 5, color='green', thickness=2)
  error = abs(f(3)-taylor(3))
  html('$f(x)=%s$'%latex(f(x)))
  html('$\hat{f}(x;%s)=%s+\mathcal{O}(x^{%s})$'%(x0,latex(taylor(x)),order+1))
  html('|$f(3)-\hat{f}(3)|= %s$'%(error.n())) 
  show(dot + p + taylor_plot+line([(3,f(3)),(3,taylor(3))],color='purple',linestyle=':'), ymin=-50,ymax=50)

Lagrange Polynomials

Instead of concentrating our accuracy around a specific point, let's make sure that our approximation function is exact at a number of specified points.  The Lagrange polynomial for $n+1$ data points is the unique $n$th degree polynomial that passes through those points.  We define the polynomial by taking linear combinations of functions so that we have the right values at the points we specify.

pts=[(i,1/i) for i in [1..5]]
html.table(pts)

1	1
2	\frac{1}{2}
3	\frac{1}{3}
4	\frac{1}{4}
5	\frac{1}{5}

#pts=[(0,0),(1,1),(2,3)]

def switch_poly(pts,point_to_delete):
    xvalues=[p[0] for p in pts]
    deleted_point=xvalues.pop(point_to_delete)
    numerator=prod((x-pt) for pt in xvalues)
    denominator=numerator(x=deleted_point)
    return numerator/denominator

For each $x$-coordinate, there is a polynomial that is 1 at that $x$-coordinate and zero at all of the others.  Here is the polynomial that is zero at $x_0,x_1,x_3,x_4$ and 1 at $x_2$.

plot(switch_poly(pts,2),(x,0.8,5.2),gridlines='major')

Here is a plot of all of the polynomials, one for each $x_i$.

point_plot=points([(p[0],1) for p in pts],color='red',pointsize=30)
curve_plot=plot([switch_poly(pts,i) for i in range(len(pts))],(x,0.8,5.2),fill=True)
point_plot+curve_plot

The Lagrange polynomial is just a linear combination of each of switching polynomials.  This guarantees that we have the correct values at the points $x_0$, $x_1$, $x_2$, $x_3$, and $x_4$.

def lagrange_poly(pts):
    return sum(switch_poly(pts,i)*pts[i][1] for i in range(len(pts)))

p(x)=lagrange_poly(pts)
p(x)

1/24*(x - 5)*(x - 4)*(x - 3)*(x - 2) - 1/12*(x - 5)*(x - 4)*(x - 3)*(x - 1) + 1/12*(x - 5)*(x - 4)*(x - 2)*(x - 1) - 1/24*(x - 5)*(x - 3)*(x - 2)*(x - 1) + 1/120*(x - 4)*(x - 3)*(x - 2)*(x - 1)

p(x).expand()

1/120*x^4 - 1/8*x^3 + 17/24*x^2 - 15/8*x + 137/60

Sage also has a builtin function to generate Lagrange polynomials as part of its polynomial objects.

R = PolynomialRing(QQ, 'x')
L = R.lagrange_polynomial(pts)
L

1/120*x^4 - 1/8*x^3 + 17/24*x^2 - 15/8*x + 137/60

Sage automatically evaluates a Polynomial object using a technique that is at least as efficient as Horner's form (if there are lots of zero coefficients, then it is more efficient by constructing a graph; type "sage.rings.polynomial.polynomial_compiled.CompiledPolynomialFunction??" to see the details and code).  Here we evaluate the polynomial at "$x$", and you can see that it looks like Horner's form is being used.

L(x)

1/120*(((x - 15)*x + 85)*x - 225)*x + 137/60

Let's plot this approximation to see how good it is.

x0  = 1
f(x)= 1/x
p   = plot(f,-1,5, thickness=2)
R = PolynomialRing(QQ, 'x')
@interact
def _(x_values=input_box([1,2,3,4,5])):
    pts=[(QQ(i),f(i)) for i in x_values]
    dots = point(pts,pointsize=80,color='red')
    L = R.lagrange_polynomial(pts)
    L_plot = plot(L,(x,-1, 5), color='green', thickness=2)
    error = abs(f(3)-L(x=3))
    html('$f(x)=%s$'%latex(f(x)))
    html('$L(x)=%s=%s$'%(latex(L(x=x)),latex(L(x=x).expand())))
    html('|$f(3)-L(x=3)|= %s$'%(latex(error.n()))) 
    show(dots + p + L_plot+line([(3,f(3)),(3,L(x=3))],color='purple',linestyle=':'), ymin=-10,ymax=10)

These approximations can be pretty bad for some functions.  A famous example is the Runge function $$\frac{1}{1+25x^2}$$

plot(1/(1+25*x^2), (x,-1,1))

R = PolynomialRing(QQ, 'x')

f(x)=1/(1+25*x^2)
R.lagrange_polynomial([(i,QQ(f(i))) for i in [0,0.1,0.2,0.3]])

34.6153846153846*x^3 - 15.3846153846154*x^2 - 0.807692307692308*x + 1.00000000000000

x0  = 1
f(x)= 1/(1+25*x^2)
x_bounds=(-1,1)
p   = plot(f,x_bounds, thickness=2)
R = PolynomialRing(QQ, 'x')
@interact
def _(x_values=input_box([-1,..,1,step=1/20])):
    pts=[(i,QQ(f(i))) for i in x_values]
    dots = point(pts,pointsize=40,color='red')
    L = R.lagrange_polynomial(pts)
    L_plot = plot(L(x),x_bounds, color='green', thickness=2)
    error = abs(f(.995)-L(.995))
    html('$f(x)=%s$'%latex(f(x)))
    html('$L(x)=%s$'%(latex(L(x).expand())))
    html('|$f(.995)-L(.995)|= %s$'%(latex(error.n()))) 
    show(dots + p + L_plot, ymin=-10,ymax=10)

It's actually much better to take points that are projections of points evenly spaced around the unit circle (the roots of the "Chebyshev" polynomial).  Try putting in "[cos(i/20*pi.n()) for i in [1,3,..,19]]" for the x_values above.

[cos(i/20*pi.n()) for i in [1,3,..,19]]

[0.987688340595138, 0.891006524188368, 0.707106781186548, 0.453990499739547, 0.156434465040231, -0.156434465040231, -0.453990499739547, -0.707106781186547, -0.891006524188368, -0.987688340595138]

This work by Jason Grout is licensed under a Creative Commons Attribution-Share Alike 3.0 United States License.