1629-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}

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

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

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

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)

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_3.py", line 9, in <module> exec compile(ur'open("___code___.py","w").write("# -*- coding: utf-8 -*-\n" + _support_.preparse_worksheet_cell(base64.b64decode("cCh4KT1sYWdyYW5nZV9wb2x5KHB0cykKcCh4KQ=="),globals())+"\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single') File "", line 1, in <module> File "/tmp/tmpnv9_yB/___code___.py", line 2, in <module> __tmp__=var("x"); p = symbolic_expression(lagrange_poly(pts)).function(x) NameError: name 'lagrange_poly' is not defined

p(x).expand()

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_9.py", line 9, in <module> exec compile(ur'open("___code___.py","w").write("# -*- coding: utf-8 -*-\n" + _support_.preparse_worksheet_cell(base64.b64decode("cCh4KS5leHBhbmQoKQ=="),globals())+"\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single') File "", line 1, in <module> File "/tmp/tmpxKdKax/___code___.py", line 2, in <module> exec compile(ur'p(x).expand()' + '\n', '', 'single') File "", line 1, in <module> NameError: name 'p' is not defined

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

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{120} x^{4} - \frac{1}{8} x^{3} + \frac{17}{24} x^{2} - \frac{15}{8} x + \frac{137}{60}

Note that Sage automatically evaluates a Polynomial object using Horner's form.  Here we evaluate the polynomial at "$x$".  A polynomial object also uses Horner form for evaluating a polynomial at a number.

L(x)

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{120} \, {\left({\left({\left(x - 15\right)} x + 85\right)} x - 225\right)} x + \frac{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(3))
    html('$f(x)=%s$'%latex(f(x)))
    html('$L(x)=%s=%s$'%(latex(L(x)),latex(L(x).expand())))
    html('|$f(3)-L(3)|= %s$'%(latex(error.n()))) 
    show(dots + p + L_plot+line([(3,f(3)),(3,L(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,-2,2))

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

Sage automatically evaluates a Polynomial object using a technique that is at least 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.