Polynomial Approximation with Derivatives

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from math import factorial

A Brief Intro to sympy

Here we import sympy, a package for symbolic computation with Python.

In [2]:
import sympy as sp
sp.init_printing()

Next, we make a (symbolic) variable $x$ from which we can then build more complicated expressions:

In [3]:
sp.var("x")
x
Out[3]:
$$x$$

Build up an expression with $x$. Assign it to expr. Observe that this expression isn't evaluated--the result of this computation is some Python data that represents the computation:

In [4]:
g = sp.sin(sp.sqrt(x)+2)**2
g
Out[4]:
$$\sin^{2}{\left (\sqrt{x} + 2 \right )}$$

Next, take a derivative, using .diff(x).

In [5]:
g.diff(x)
Out[5]:
$$\frac{1}{\sqrt{x}} \sin{\left (\sqrt{x} + 2 \right )} \cos{\left (\sqrt{x} + 2 \right )}$$
In [6]:
g.diff(x, 4)
Out[6]:
$$\frac{1}{2 x^{2}} \sin^{2}{\left (\sqrt{x} + 2 \right )} - \frac{1}{2 x^{2}} \cos^{2}{\left (\sqrt{x} + 2 \right )} - \frac{15}{8 x^{3}} \sin^{2}{\left (\sqrt{x} + 2 \right )} + \frac{15}{8 x^{3}} \cos^{2}{\left (\sqrt{x} + 2 \right )} + \frac{3}{x^{\frac{5}{2}}} \sin{\left (\sqrt{x} + 2 \right )} \cos{\left (\sqrt{x} + 2 \right )} - \frac{15}{8 x^{\frac{7}{2}}} \sin{\left (\sqrt{x} + 2 \right )} \cos{\left (\sqrt{x} + 2 \right )}$$

Use .subs(x, ...) and .evalf() to evaluate your expression for $x=1$.

In [7]:
g.subs(x,1)
Out[7]:
$$\sin^{2}{\left (3 \right )}$$
In [8]:
g.subs(x, 1).evalf()
Out[8]:
$$0.019914856674817$$

Helper function that takes symbolic functions as argument and plot them

In [9]:
def plot_sympy(my_f, my_pts, **kwargs):
    f_values = np.array([my_f.subs(x, pt) for pt in my_pts])
    plt.plot(pts, f_values, **kwargs)

Polynomial Approximation

In [10]:
f = 1/(20*x-10)
f
Out[10]:
$$\frac{1}{20 x - 10}$$
In [11]:
f.diff(x)
Out[11]:
$$- \frac{20}{\left(20 x - 10\right)^{2}}$$
In [12]:
f.diff(x,2)
Out[12]:
$$\frac{4}{5 \left(2 x - 1\right)^{3}}$$

Write out the degree 2 Taylor polynomial about 0:

In [13]:
taylor2 = (
    f.subs(x, 0)
    + f.diff(x).subs(x, 0) * x
    + f.diff(x, 2).subs(x, 0)/2 * x**2
)
taylor2
Out[13]:
$$- \frac{2 x^{2}}{5} - \frac{x}{5} - \frac{1}{10}$$

Plot the exact function f and the taylor approximation taylor2

In [14]:
pts = np.linspace(-0.4,0.4)

plot_sympy(taylor2, pts, label="taylor2")
plot_sympy(f, pts, label="f")
plt.legend(loc="best")
plt.axis('equal')
plt.grid()
plt.xlabel('$x$')
plt.ylabel('function values')
Out[14]:
Text(0, 0.5, 'function values')

Behavior of the Error

Let's write the taylor approximation for any degree n, and define the error as f - tn

In [15]:
n = 2

tn = 0
for i in range(n+1):
    tn += f.diff(x, i).subs(x, 0)/factorial(i) * x**i
In [16]:
plot_sympy(tn, pts, label="taylor")
plot_sympy(f, pts, label="f")
plt.legend(loc="best")
plt.axis('equal')
plt.grid()
plt.xlabel('$x$')
plt.ylabel('function values')
Out[16]:
Text(0, 0.5, 'function values')
In [17]:
error = f - tn
In [18]:
plot_sympy(error, pts, label="error")
plt.legend(loc="best")
plt.ylim([-1.3, 1.3])
plt.axis('equal')
plt.grid()
plt.xlabel('$x$')
plt.ylabel('error')
Out[18]:
Text(0, 0.5, 'error')

To get a better idea of what happens close to the center, use a log-log plot:

In [19]:
# plot only points close to zero [10^(-3),10^(0.4)]
plt.figure(figsize=(10,6))
pos_pts = 10**np.linspace(-3, 0.4) 
err_values = [abs(error.subs(x, pt)) for pt in pos_pts]
plt.loglog(pos_pts, err_values)
plt.grid()
plt.xlabel("$x$")
plt.ylabel("Error")
Out[19]:
Text(0, 0.5, 'Error')