Pythonでマクローリン展開を実装してわかりやすく説明してみた
高校までは奥手だった自分も、大学に入って素敵な恋の展開を迎えることができる
かと思いきや、
なんてことがありがちな理系大学生。
今回の記事ではPythonを使って感覚的にテイラー展開とマクローリン展開を理解していきます。
テイラー展開とは
テイラー展開とは、複雑な数式を多項式で近似する方法のことです。
$f(x)$の$a$を中心としたテイラー展開は、
\begin{eqnarray}
f(x) &=& \sum_{k=0}^{\infty} \frac{f^{(k)}(a)}{k!}(x-a)^k \\
&=& f(a) + \frac{f'(a)}{1!} (x - a) + \frac{f''(a)}{2!} (x - a)^2 + … + \frac{f^{(k)}(a)}{k!} (x - a)^k + …
\end{eqnarray}
と表すことができます。
マクローリン展開とは
テイラー展開の$a=0$、つまり$f(x)$の$x=0$、ようするに原点周りのテイラー展開は、「マクローリン展開」と呼ばれ、以下のように表すことができます。
\begin{eqnarray}
f(x) &=& \sum_{k=0}^{\infty} \frac{f^{(k)}(0)}{k!}(x)^k \\
&=& f(0) + \frac{f'(0)}{1!} (x) + \frac{f''(0)}{2!} (x)^2 + … + \frac{f^{(k)}(0)}{k!} (x)^k + …
\end{eqnarray}
マクローリン展開の何が嬉しいのか
マクローリン展開が複雑な式を原点周りで近似できるのはわかった、で、なにが嬉しいの？
と思った人がいるかと思います。
マクローリン展開の嬉しいところを紹介していきます。
計算が簡単になる
例えば$sin(x)$のグラフを正確に書こうと思ったとき、$sin(0.1)$とか、$sin(0.01)$などを求めることになりますが、それを計算するのはかなり大変です。
そこで、マクローリン展開を使って計算します。計算機でもマクローリン展開を使っているらしいです。
実は、高校物理でも使われています。$x << 1$（$x$が1よりも十分に小さい）の時、
$$ e^{x} = 1 $$
$$ \sin{x} = x $$
$$ \cos{x} = 1 $$
などが成り立つと教えられた記憶があると思います。なので$sin(0.01) = 0.01$としても高校物理ではOKというわけです。1よりも十分に小さい時に上の式でも十分に近似することはこのグラフを見ればわかります。
$x = 1$付近では$0.2$くらい差がありそうです。
もうひとつ項数を増やしてみます。
項数を増やすと、$x=1.0$あたりまでは近似していることがわかります。このように、項数を増やせば増やすほどに近似していきます。
実は上記三つの式はマクローリン展開した時の第一項目を表したものです。
もう一度表記すると、マクローリン展開は、
\begin{eqnarray}
f(x) &=& \sum_{k=0}^{\infty} \frac{f^{(k)}(0)}{k!}(x)^k \\
&=& f(0) + \frac{f'(0)}{1!} (x) + \frac{f''(0)}{2!} (x)^2 + … + \frac{f^{(k)}(0)}{k!} (x)^k + …
\end{eqnarray}
と表すことができます。
マクローリン展開を使うと、$e^x$や$\sin{x}$や$\cos{x}$は以下のように近似できます。
$$ e^x = \sum_{k=0}^{\infty} \frac{x^{k}}{k!} = 1 + \frac{x}{1!} + \frac{x^2}{2!} + … $$
$$ \sin{x} = \sum_{k=0}^{\infty} (-1)^k \frac{x^{2k+1}}{(2k+1)!} = \frac{x}{1!} - \frac{x^3}{3!} + \frac{x^5}{5!} +… $$
$$ \cos{x} = \sum_{k=0}^{\infty} (-1)^k \frac{x^{2k}}{(2k)!} = 1 - \frac{x^2}{2!} + \frac{x^4}{4!} +… $$
マクローリン展開との組み合わせでシンプルな結果を得ることができる。
たとえば積率母関数とマクローリン展開を組み合わせるとこんなことができます。
標準正規分布に従う確率変数Xのモーメントの分布を調べるとき、
直接計算するには、
$$ E(X^k) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} x^k e^{\frac{-x^2}{2}} dx $$
を計算することになりますがこれは非常に面倒です。そこで、
Xの積率母関数（モーメント母関数）$M_X(t) = E(e^{tX})$をマクローリン展開すると、
$$ E(X^k) = \begin{cases} \frac{(2n)!}{2^n n!} (偶数(k=2n)の時) \\ 0 (奇数の時) \end{cases} $$
というシンプルな結果を得ることができます。これは元の積分の式よりも明らかに使いやすいですよね。
※ 詳しくは、ひとつ前の記事に書いてあります。
こんな風にマクローリンはいろんなところで役立つので、解析学の最初にやるわけですね。
マクローリン展開の導出
マクローリン展開は、まず複雑な式を多項式で近似できると仮定、つまり
$$ f(x) = A_0 + A_1 x + A_2 x^2 + … $$
が成り立つとした時、原点周りでは近似すると仮定しているので、$x=0$を代入すると$A_0$が明らかになります。さらに$f(x)を$微分したものも原点周りで近似するはずなので、$f(x)$を微分して$x=0$を代入すると$A_1$が明らかになり、さらにその式を微分して$x=0$を代入して、、、と逐次的に係数が明らかにできることから導出することができます。
すると、一般にマクローリン展開は
\begin{eqnarray}
f(x) &=& \sum_{k=0}^{\infty} \frac{f^{(k)}(0)}{k!}(x)^k \\
&=& f(0) + \frac{f'(0)}{1!} (x) + \frac{f''(0)}{2!} (x)^2 + … + \frac{f^{(k)}(0)}{k!} (x)^k + …
\end{eqnarray}
と表すことができるのです。
$\sin{x}の場合を考えてみます。$、
$$\sin(x) = A_0 + A_1 x + A_2 x^2 + A_3 x^3 + A_4 x^4 + …$$
とした時、$x = 0$を代入すれば、$sin0 = A_0$より、$A_0 = 0$となり、これを微分すると
$$ \sin'(x) = \cos(x) = A_1 + 2A_2 x + 3A_3 x^2 + 4A_4^3 + …$$
に$x = 0$を代入すれば、$\cos0 = A_1$ となるため、$A_1 = 1$となります。さらにこれを微分すると
$$\cos'(x) = - \sin(x) = 2A_2 + 6A_3 x + 12A_4^2 + …$$
となり、これに $x = 0$を代入すれば$- \sin0 = 2A_2$となるため、$A_2 = \frac{0}{2} = 0$。この手順で進めていけば、
$$ \sin{x} = \sum_{k=0}^{\infty} (-1)^k \frac{x^{2k+1}}{(2k+1)!} = \frac{x}{1!} - \frac{x^3}{3!} + \frac{x^5}{5!} +… $$
となることがわかります。
Pythonでマクローリン展開を実装してみた
Pythonの記号計算を行うためのライブラリであるsympyを使って、マクローリン展開を実装してみました。
ライブラリのインポート、各種設定を行います。
import numpy as np
import sympy
from sympy.plotting import plot
sympy.init_printing(use_unicode=True)
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()$$ f(x) = \sum_{k=0}^{\infty} \frac{f^{(k)}(0)}{k!}(x)^k $$
でした。
$k$は大きくすればするほど、元の式に近いものになりますが、ここでは上限を柔軟に変えるために関数を呼び出す時にk_maxで指定するようにします。
def mcLaughlin_expansion(formula,k_max):
f = formula
approximation = np.sum([(sympy.diff(f, x, k).subs(x, 0) / (sympy.factorial(k))) * ((x)**k) for k in range(k_max)])
return approximation$e^x$をマクローリン展開
$e^x$をマクローリン展開してできた多項式を見ていきます。
x = sympy.symbols("x")
formula = sympy.E**x
print("formula: ",formula)
print("mcLaughlin: ",mcLaughlin_expansion(formula,10))formula: 2.71828182845905**x mcLaughlin: 2.75573192239859e-6*x**9 + 2.48015873015873e-5*x**8 + 0.000198412698412698*x**7 + 0.00138888888888889*x**6 + 0.00833333333333333*x**5 + 0.0416666666666667*x**4 + 0.166666666666667*x**3 + 0.5*x**2 + 1.0*x + 1
$\sum_{k=0}^{\infty} \frac{x^{k}}{k!} = 1 + \frac{x}{1!} + \frac{x^2}{2!} + …$になってますね。
それでは、$e^x$と$e^x$をマクローリン展開してできた多項式を可視化して比較していきます。
x = sympy.symbols("x")
formula = sympy.E**x
plt.rcParams['figure.figsize'] = 10, 8
f1 = mcLaughlin_expansion(formula, 3)
f2 = mcLaughlin_expansion(formula, 6)
f3 = mcLaughlin_expansion(formula, 9)
f4 = mcLaughlin_expansion(formula, 12)
f5 = mcLaughlin_expansion(formula, 15)
ax = sympy.plot(f1,f2,f3,f4,f5,formula,
legend=True, show=False)
colors = ["orange", "red", "purple", "blue", "green"]
for i in range(6):
if i == 5:
ax[i].line_color="black"
ax[i].label='$e^x$'
else:
ax[i].line_color=colors[i]
label_i = i+1
ax[i].label='$f%d$'%label_i
ax.show()
このように、k_maxが大きければ大きいほど、つまり項が多ければ多いほど、黒線である元の式$e^x$に近似することがわかります。
$sin(x)$のマクローリン展開
x = sympy.symbols("x")
formula = sympy.sin(x)
print("formula: ",formula)
print("mcLaughlin: ",mcLaughlin_expansion(formula,10))formula: sin(x) mcLaughlin: x**9/362880 - x**7/5040 + x**5/120 - x**3/6 + x
plt.rcParams['figure.figsize'] = 10, 8
f1 = mcLaughlin_expansion(formula, 5)
f2 = mcLaughlin_expansion(formula, 10)
f3 = mcLaughlin_expansion(formula, 15)
f4 = mcLaughlin_expansion(formula, 20)
f5 = mcLaughlin_expansion(formula, 25)
ax = sympy.plot(f1,f2,f3,f4,f5,formula,
ylim=(-5,5), legend=True, show=False)
colors = ["orange", "red", "purple", "blue", "green"]
for i in range(6):
if i == 5:
ax[i].line_color="black"
ax[i].label='$sin(x)$'
else:
ax[i].line_color=colors[i]
label_i = i+1
ax[i].label='$f%d$'%label_i
ax.show()
$cos(x)$のマクローリン展開
x = sympy.symbols("x")
formula = sympy.cos(x)
print("formula: ",formula)
print("mcLaughlin: ",mcLaughlin_expansion(formula,10))formula: cos(x) mcLaughlin: x**8/40320 - x**6/720 + x**4/24 - x**2/2 + 1
plt.rcParams['figure.figsize'] = 10, 8
f1 = mcLaughlin_expansion(formula, 5)
f2 = mcLaughlin_expansion(formula, 10)
f3 = mcLaughlin_expansion(formula, 15)
f4 = mcLaughlin_expansion(formula, 20)
f5 = mcLaughlin_expansion(formula, 25)
ax = sympy.plot(f1,f2,f3,f4,f5,formula,
ylim=(-5,5), legend=True, show=False)
colors = ["orange", "red", "purple", "blue", "green"]
for i in range(6):
if i == 5:
ax[i].line_color="black"
ax[i].label='$sin(x)$'
else:
ax[i].line_color=colors[i]
label_i = i+1
ax[i].label='$f%d$'%label_i
ax.show()
まとめ
テイラー展開、マクローリン展開の有用性を解説したり、Pythonで実装して可視化してみたりしてみました。
可視化すると幾何学的な模様ができて、ワクワクしました。
sympyでハイパボリックタンジェントやシグモイド、ReLUを実装して可視化した記事も書いてます。
