Return to SymPy で理工系の数学B

SymPy でテイラー展開・オイラーの公式

Python で微分積分などの記号計算をする場合は,以下のように必要なパッケージを import します。(弘大 JupyterHub では必要なパッケージはインストール済みです。)

セルをクリックして,上のメニューから「▶︎ Run」をクリックするか,Shift キーを押しながら Enter キー(Return キー)を押すと実行します。

In [1]:
# a Python library for symbolic mathematics
from sympy import *
# 1文字変数の Symbol の定義が省略できる
from sympy.abc import *
# π,ネイピア数,虚数単位
from sympy import pi, E, I

# SymPy Plotting Backends (SPB): グラフを描く際に利用
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

高階導関数

$f(x)$ を $x$ で $n$ 階微分する例。SymPy では diff(f(x), x, n) のように書きます。
例:

$$\frac{d^2}{dx^2} \sin x = \cdots$$

In [2]:
# 手っ取り早く答えを知りたいときは...

diff(sin(x), x, 2)
Out[2]:
$\displaystyle – \sin{\left(x \right)}$
In [3]:
# 何を計算するかを表示させてから実行する場合は...

Eq(Derivative(sin(x), x, 2), 
   diff(sin(x), x, 2))
Out[3]:
$\displaystyle \frac{d^{2}}{d x^{2}} \sin{\left(x \right)} = – \sin{\left(x \right)}$

テイラー展開

テイラー展開は series() 関数を使います。

指数関数のテイラー展開

$f(x) = e^x$ の $x = 0$ のまわりに $x^5$ までテイラー展開する例。5次までの展開なのになぜ 6 かと悩むかもしれないが,これが Python の流儀らしい。

$x^0$ から数えて 6 番目までだから 6 とか,56 未満なので 6 と書くとか,各自納得する方法で覚えておく。

In [4]:
f = exp(x)
series(f, x, 0, 6)
Out[4]:
$\displaystyle 1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + O\left(x^{6}\right)$

$O(x^6)$ は $x^6$ 以上の高次の項という意味。

ちなみに,この式で $x=1$ とおくと
$$ e = \sum_{k=0}^{\infty} \frac{1}{k!}$$ となり,この級数展開を使ってネイピア数 $e$ の近似値を求めることができる。

例として,$\displaystyle \sum_{k=0}^{10} \frac{1}{k!}$ の計算。SymPy では階乗 $k!$ は factorial(k) です。

In [5]:
# 手っ取り早く総和の答えを知りたいときは...

summation(1/factorial(k), (k, 0, 10))
Out[5]:
$\displaystyle \frac{9864101}{3628800}$
In [6]:
# 何を計算するかを表示させてから実行する場合は...

Eq(Sum(1/factorial(k), (k, 0, 10)), 
   Sum(1/factorial(k), (k, 0, 10)).doit())
Out[6]:
$\displaystyle \sum_{k=0}^{10} \frac{1}{k!} = \frac{9864101}{3628800}$
In [7]:
# 浮動小数点数で表示

float(summation(1/factorial(k), (k, 0, 10)))
Out[7]:
2.7182818011463845
In [8]:
# ネイピア数と比較

float(E)
Out[8]:
2.718281828459045

三角関数のテイラー展開

$f(x) = \sin x$ の $x = 0$ のまわりのテイラー展開は
\begin{eqnarray}
\sin x &=& f(0) + f'(0)\, x + \frac{f”(0)}{2!}\,x^2 + \frac{f”'(0)}{3!}\,x^3 +\cdots \\
&=& \sin 0 + \cos 0 \cdot x + \frac{-\sin 0}{2!}\,x^2 + \frac{-\cos 0}{3!}\,x^3 + \cdots \\
&=& x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!} + \cdots\\
&=& \sum_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)!}\,x^{2n+1}
\end{eqnarray}

In [9]:
series(sin(x), x, 0, 8)
Out[9]:
$\displaystyle x – \frac{x^{3}}{6} + \frac{x^{5}}{120} – \frac{x^{7}}{5040} + O\left(x^{8}\right)$
In [10]:
# 階乗の確認

factorial(3), factorial(5), factorial(7)
Out[10]:
(6, 120, 5040)

$f(x) = \cos x$ の $x = 0$ のまわりのテイラー展開は
\begin{eqnarray}
\cos x &=& f(0) + f'(x)\, x + \frac{f”(0)}{2!}\,x^2 + \frac{f”'(0)}{3!}\,x^3 +\cdots \\
&=& \cos 0 – \sin 0 \cdot x + \frac{-\cos 0}{2!}\,x^2 + \frac{\sin 0}{3!}\,x^3 + \cdots \\
&=& 1 – \frac{x^2}{2!} + \frac{x^4}{4!} – \frac{x^6}{6!} + \cdots\\
&=& \sum_{n=0}^{\infty} \frac{(-1)^n}{(2n)!}\,x^{2n}
\end{eqnarray}

In [11]:
series(cos(x), x, 0, 7)
Out[11]:
$\displaystyle 1 – \frac{x^{2}}{2} + \frac{x^{4}}{24} – \frac{x^{6}}{720} + O\left(x^{7}\right)$
In [12]:
# 階乗の確認

factorial(2), factorial(4), factorial(6)
Out[12]:
(2, 24, 720)

参考:テイラー展開した関数のグラフ

$f(x) = e^x$ のテイラー展開

.removeO() で高次の項を取りのぞいて「普通の」関数にしてから plot() します。

In [13]:
expt1 = series(exp(x), x, 0, 2).removeO()
expt2 = series(exp(x), x, 0, 3).removeO()
expt3 = series(exp(x), x, 0, 4).removeO()
expt1, expt2, expt3
Out[13]:
(x + 1, x**2/2 + x + 1, x**3/6 + x**2/2 + x + 1)
In [14]:
plot(exp(x), expt3, expt2, expt1, (x, -2, 2));

$f(x) = \sin x$ のテイラー展開
In [15]:
sint1 = series(sin(x), x, 0, 2).removeO()
sint3 = series(sin(x), x, 0, 4).removeO()
sint5 = series(sin(x), x, 0, 6).removeO()
sint1, sint3, sint5
Out[15]:
(x, -x**3/6 + x, x**5/120 - x**3/6 + x)
In [16]:
plot(sin(x), sint5, sint3, sint1, (x, -pi, pi));

テイラー展開の次数を上げていくと,$|x| < 1$ の範囲では,もとの関数に近づいていくようすがわかります。

人類の至宝:オイラーの公式

オイラーの公式を確認するために,指数関数 $e^{i x}$ を .rewrite(cos) をつけて三角関数で書き直してみます。SymPy では虚数単位 $i$ は I です。

from sympy import pi, E, I
In [17]:
Eq(exp(I*x), 
   exp(I*x).rewrite(cos))
Out[17]:
$\displaystyle e^{i x} = i \sin{\left(x \right)} + \cos{\left(x \right)}$

逆に,三角関数を指数関数を使って表してみます。

In [18]:
cos(x).rewrite(exp)
Out[18]:
$\displaystyle \frac{e^{i x}}{2} + \frac{e^{- i x}}{2}$
In [19]:
sin(x).rewrite(exp)
Out[19]:
$\displaystyle – \frac{i \left(e^{i x} – e^{- i x}\right)}{2}$

以上のことから

\begin{eqnarray}
\cos x &=& \frac{e^{ix} + e^{-ix}}{2} \\
\sin x &=& \frac{e^{ix} – e^{-ix}}{2 i}
\end{eqnarray}

であることが確認できます。

オイラーの等式

$$e^{i\pi} + 1 = 0$$

オイラーの等式のどこがすごいかというと,

  • 可能における単位元ゼロ $0$ と乗法における単位元 $1$ という整数のもっとも基本となる数
  • 無理数の代表選手,ネイピア数 $e$ と円周率 $\pi$
  • そして虚数単位 $i$

という,いずれ名だたる役者達が,加法,乗法および冪乗によって見事に結び付けられているということ!

In [20]:
E**(I * pi) + 1 == 0
Out[20]:
True

三角関数と双曲線関数の関係

三角関数と双曲線関数は,ただまぎらわしいほどに似た表記なだけでなく,アイで結ばれた密接な関係であることを確認する。

In [21]:
cosh(I*x)
Out[21]:
$\displaystyle \cos{\left(x \right)}$
In [22]:
sinh(I*x)
Out[22]:
$\displaystyle i \sin{\left(x \right)}$
In [23]:
cos(I*x)
Out[23]:
$\displaystyle \cosh{\left(x \right)}$
In [24]:
sin(I*x)
Out[24]:
$\displaystyle i \sinh{\left(x \right)}$

… ということで以下のことが確認できた。

\begin{eqnarray}
\cosh i x &=& \cos x \\
\sinh i x &=& i \sin x \\
\cos i x &=& \cosh x \\
\sin i x &=& i \sinh x
\end{eqnarray}

逆三角関数と逆双曲線関数の関係

逆三角関数と逆双曲線関数もまた,ただまぎらわしいほどに似た表記なだけでなく,アイで結ばれた密接な関係であることを確認する。

In [25]:
asin(I*x)
Out[25]:
$\displaystyle i \operatorname{asinh}{\left(x \right)}$
In [26]:
atan(I*x)
Out[26]:
$\displaystyle i \operatorname{atanh}{\left(x \right)}$
In [27]:
acos(I*x)
Out[27]:
$\displaystyle – i \operatorname{asinh}{\left(x \right)} + \frac{\pi}{2}$

… ということで以下のことが確認できた。

\begin{eqnarray}
\sin^{-1} i x &=& i \sinh^{-1} x \\
\tan^{-1} i x &=& i \tanh^{-1} x \\
\cos^{-1} i x &=& – i \sinh^{-1} x + \frac{\pi}{2}
\end{eqnarray}

最後の $\displaystyle \frac{\pi}{2}$ が何かしっくりこないが,これは以下の関係から理解できる。

$$\cos^{-1} x + \sin^{-1} x = \frac{\pi}{2}$$

SymPy でも確認しておく。

In [28]:
Eq(acos(x) + asin(x), 
   acos(x) + asin(x).rewrite(acos))
Out[28]:
$\displaystyle \operatorname{acos}{\left(x \right)} + \operatorname{asin}{\left(x \right)} = \frac{\pi}{2}$