In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from sympy import *
from functools import reduce
from IPython.display import display, Math
init_printing()
x = Symbol("x")

$x=\frac{\pi}2$ で $1$ になるよう係数 $a$ を調整する。

$\sin\left(\frac{\pi}6\right)=\frac 1 2$ と値を比較する。(緑の縦線)

In [2]:
def plot(x1, x2):
    r = range(x1, x2 + 1)
    t = reduce(lambda t, i: t + (latex(x) if i == 0 else "(" + latex(x - i * pi) + ")"), r, "")
    y = reduce(lambda y, i: y * (x - i * pi), r, 1)
    a = 1 / y.subs([(x, pi / 2)])
    ay = a * y
    
    def f(c, rx):
        c.plot(rx, np.sin(rx))
        c.plot(rx, [ay.evalf(subs={x: i}) for i in rx])
        
    _, cols = plt.subplots(ncols=3, figsize=(15, 4))
    f(cols[0], np.linspace(np.pi * (x1 - 1), np.pi * (x2 + 1), 100))
    f(cols[1], np.linspace(0, np.pi, 100))
    f(cols[2], np.linspace(np.pi / 8, np.pi / 4, 100))
    cols[2].plot([np.pi / 6] * 2, [0.6, 0.4])
    
    display(Math("f(x)=%s%s" % (latex(a), t)))
    v = ay.evalf(subs={x: pi / 6})
    a2 = "\\left(\\frac{\\pi}6\\right)"
    display(Math("f%s=%f,\\quad f%s-\\sin%s=%f" % (a2, v, a2, a2, v - sin(pi / 6))))

奇数の次数だと中心がずれる。

In [3]:
plot(-1, 1)
$$f(x)=- \frac{8}{3 \pi^{3}}(x + \pi)x(x - \pi)$$
$$f\left(\frac{\pi}6\right)=0.432099,\quad f\left(\frac{\pi}6\right)-\sin\left(\frac{\pi}6\right)=-0.067901$$

偶数の次数でどんどん増やして確認する。

In [4]:
plot(0, 1)
$$f(x)=- \frac{4}{\pi^{2}}x(x - \pi)$$
$$f\left(\frac{\pi}6\right)=0.555556,\quad f\left(\frac{\pi}6\right)-\sin\left(\frac{\pi}6\right)=0.055556$$
In [5]:
plot(-1, 2)
$$f(x)=\frac{16}{9 \pi^{4}}(x + \pi)x(x - \pi)(x - 2 \pi)$$
$$f\left(\frac{\pi}6\right)=0.528121,\quad f\left(\frac{\pi}6\right)-\sin\left(\frac{\pi}6\right)=0.028121$$
In [6]:
plot(-2, 3)
$$f(x)=- \frac{64}{225 \pi^{6}}(x + 2 \pi)(x + \pi)x(x - \pi)(x - 2 \pi)(x - 3 \pi)$$
$$f\left(\frac{\pi}6\right)=0.518732,\quad f\left(\frac{\pi}6\right)-\sin\left(\frac{\pi}6\right)=0.018732$$

20次でも誤差はまだかなり残る。

In [7]:
plot(-9, 10)
$$f(x)=\frac{1048576}{428670161650355625 \pi^{20}}(x + 9 \pi)(x + 8 \pi)(x + 7 \pi)(x + 6 \pi)(x + 5 \pi)(x + 4 \pi)(x + 3 \pi)(x + 2 \pi)(x + \pi)x(x - \pi)(x - 2 \pi)(x - 3 \pi)(x - 4 \pi)(x - 5 \pi)(x - 6 \pi)(x - 7 \pi)(x - 8 \pi)(x - 9 \pi)(x - 10 \pi)$$
$$f\left(\frac{\pi}6\right)=0.505583,\quad f\left(\frac{\pi}6\right)-\sin\left(\frac{\pi}6\right)=0.005583$$
In [ ]: