Computing $\pi$ with Taylor

In [16]:
from math import factorial

import numpy as np
import matplotlib.pyplot as plt

import sympy as sp
sp.init_printing()

x = sp.Symbol("x")
In [17]:
f = sp.sqrt(1-x**2)
In [18]:
pts = np.linspace(0, 1, 100)
f_values = np.array([f.subs(x, pt) for pt in pts])
plt.plot(pts, f_values)
Out[18]:
[<matplotlib.lines.Line2D at 0x7f14ec47f828>]

Get the Taylor coefficients into a list coeffs:

In [43]:
n = 50
In [44]:
coeffs = []
for i in range(n):
    coeffs.append(f.diff(x, i).subs(x, 0)/factorial(i))

coeffs
Out[44]:
$$\left [ 1, \quad 0, \quad - \frac{1}{2}, \quad 0, \quad - \frac{1}{8}, \quad 0, \quad - \frac{1}{16}, \quad 0, \quad - \frac{5}{128}, \quad 0, \quad - \frac{7}{256}, \quad 0, \quad - \frac{21}{1024}, \quad 0, \quad - \frac{33}{2048}, \quad 0, \quad - \frac{429}{32768}, \quad 0, \quad - \frac{715}{65536}, \quad 0, \quad - \frac{2431}{262144}, \quad 0, \quad - \frac{4199}{524288}, \quad 0, \quad - \frac{29393}{4194304}, \quad 0, \quad - \frac{52003}{8388608}, \quad 0, \quad - \frac{185725}{33554432}, \quad 0, \quad - \frac{334305}{67108864}, \quad 0, \quad - \frac{9694845}{2147483648}, \quad 0, \quad - \frac{17678835}{4294967296}, \quad 0, \quad - \frac{64822395}{17179869184}, \quad 0, \quad - \frac{119409675}{34359738368}, \quad 0, \quad - \frac{883631595}{274877906944}, \quad 0, \quad - \frac{1641030105}{549755813888}, \quad 0, \quad - \frac{6116566755}{2199023255552}, \quad 0, \quad - \frac{11435320455}{4398046511104}, \quad 0, \quad - \frac{171529806825}{70368744177664}, \quad 0\right ]$$

For comparison: Here's the full Taylor polynomial:

In [45]:
taylor = 0
for i in range(n):
    taylor += f.diff(x, i).subs(x, 0)/factorial(i) * x**i

taylor
Out[45]:
$$- \frac{171529806825 x^{48}}{70368744177664} - \frac{11435320455 x^{46}}{4398046511104} - \frac{6116566755 x^{44}}{2199023255552} - \frac{1641030105 x^{42}}{549755813888} - \frac{883631595 x^{40}}{274877906944} - \frac{119409675 x^{38}}{34359738368} - \frac{64822395 x^{36}}{17179869184} - \frac{17678835 x^{34}}{4294967296} - \frac{9694845 x^{32}}{2147483648} - \frac{334305 x^{30}}{67108864} - \frac{185725 x^{28}}{33554432} - \frac{52003 x^{26}}{8388608} - \frac{29393 x^{24}}{4194304} - \frac{4199 x^{22}}{524288} - \frac{2431 x^{20}}{262144} - \frac{715 x^{18}}{65536} - \frac{429 x^{16}}{32768} - \frac{33 x^{14}}{2048} - \frac{21 x^{12}}{1024} - \frac{7 x^{10}}{256} - \frac{5 x^{8}}{128} - \frac{x^{6}}{16} - \frac{x^{4}}{8} - \frac{x^{2}}{2} + 1$$

Now compute an approximate integral:

In [46]:
tmp = 0
for i in range(n):
    tmp += (coeffs[i] * 1/(i+1) * (1**(i+1) - 0**(i+1))).evalf()
In [47]:
my_pi = 4*tmp
print(my_pi, np.pi)
3.14468451860141 3.141592653589793
  • What do you think?
  • Did that work well?
  • Did it have a right to work well?
  • One important aspect of studying numerical methods:

    • When do things like this work well/efficiently?
    • When do they not?
    • Why?
In [ ]: