import numpy as np
import matplotlib.pyplot as plt
from math import factorial
sympy
¶Here we import sympy
, a package for symbolic computation with Python.
import sympy as sp
sp.init_printing()
Next, we make a (symbolic) variable $x$ from which we can then build more complicated expressions:
x = sp.Symbol("x")
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:
expr = sp.sin(sp.sqrt(x)+2)**2
expr
Next, take a derivative, using .diff(x)
.
expr.diff(x)
Take 10 derivatives:
expr.diff(x, 10)
Use .subs(x, ...)
and .evalf()
to evaluate your expression for $x=1$.
expr.subs(x, 1).evalf()
Here are some functions to play with:
#f = sp.exp(x)
f = sp.sqrt(1-x**2)
#f = 1/(20*x-10)
Plot f
:
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)
plt.ylim([-1.3, 1.3])
pts = np.linspace(-1, 1, 1000)
plot_sympy(f, pts)
The Taylor polynomial uses derivatives to automatically get close to the function--the closer the higher you increase the degree.
Write out the degree 4 Taylor polynomial about 0, call it taylor4
:
taylor4 = (
f.subs(x, 0)
+ f.diff(x).subs(x, 0) * x
+ f.diff(x, 2).subs(x, 0)/2 * x**2
+ f.diff(x, 3).subs(x, 0)/6 * x**3
+ f.diff(x, 4).subs(x, 0)/24 * x**4
)
taylor4
plot_sympy(taylor4, pts, label="taylor4")
plot_sympy(f, pts, label="f")
plt.legend(loc="best")
Now write code to do this for any order n
. Put the result in taylor
:
n = 20
taylor = 0
for i in range(n):
taylor += f.diff(x, i).subs(x, 0)/factorial(i) * x**i
taylor
plot_sympy(taylor, pts, label="taylor")
plot_sympy(f, pts, label="f")
plt.legend(loc="best")
Let's repeat the setup, for clarity:
f = 1/(20*x-10)
n = 4
taylor = 0
for i in range(n):
taylor += f.diff(x, i).subs(x, 0)/factorial(i) * x**i
error = taylor - f
plot_sympy(error, pts, label="error")
plt.legend(loc="best")
Let's look back at this Taylor expansion:
taylor
next_taylor
.next_taylor = 0
for i in range(n+1):
next_taylor += f.diff(x, i).subs(x, 0)/factorial(i) * x**i
taylor-next_taylor
How does that difference compare to the actual error?
plot_sympy(taylor-next_taylor, pts, label="diff to next")
plot_sympy(error, pts, label="error")
plt.legend(loc="best")
pos_pts = 10**np.linspace(-3, 0.5)
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")
f = sp.sqrt(x-10)
t = sum([f.diff(x, i).subs(x, 12)/factorial(i) * (x-12)**i for i in range(4)])
t
f.subs(x, 12.5)
t.subs(x, 12.5).evalf()
error1 = f.subs(x, 12.5) - t.subs(x, 12.5).evalf()
abs(error1)
Now predict the error at $12.25$:
abs(error1) * (0.25/0.5)**4
error2 = f.subs(x, 12.25) - t.subs(x, 12.25).evalf()
abs(error2)