pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

実践例で学ぶSymPy, NumPy, Matplotlibの使い方:指数関数の解析

数学の具体的な計算にPython(SymPy, NumPy, Matplotlib)を使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使った実践例として指数関数を解析してみます。とくに数式処理ライブラリであるSymPyと数値計算ライブラリのNumPyの両方を使って、シンボリックで解析的な計算と数値的な計算とを組み合わせて対象を解析する方法を具体的にみていきます。


今回の記事は実践例ということで、Pythonの書式についての詳しい説明は省いています。詳細が気になった方は以下の過去記事を参照してみてください:

pianofisica.hatenablog.com

pianofisica.hatenablog.com

pianofisica.hatenablog.com

pianofisica.hatenablog.com

pianofisica.hatenablog.com

pianofisica.hatenablog.com



指数関数

指数関数は自然対数の底(Napir数:ネイピア数 e の冪乗のベキを変数とする関数です:

\quad \displaystyle{\exp(x):\ x\rightarrow e^x}

ここでネイピア数  e は(さまざまな定義がありますが)次の極限

\quad \displaystyle{e=\lim_{n\to\infty}\left(1+\frac{1}{n}\right)^n}

で定まる実数です:Pythonでその数値を求めてみますと

n = 10**8
(1+1/n)**n

から、その近似値が(8桁程度の精度で)

\quad \displaystyle{e\sim 2.7182817983473577}

と求まります。実はNumPyには、あらかじめその値が組み込まれていて

import numpy as np
np.e

から、有効数字の範囲内で

\quad \displaystyle{e= 2.718281828459045}

です。

指数関数は次の指数法則を満たします:

\quad \displaystyle{\exp(x+y)=\exp(x)\exp(y)}


微分公式

指数関数の重要な性質のひとつとして、次の微分公式

\quad \displaystyle{\frac{d}{dx}\exp(x)=\exp(x)}

があげられます。これは以下のようにして確かめることができます。まず、微分の定義、指数法則、ネイピア数の定義式から

\quad \displaystyle{\begin{aligned}\frac{d}{dx}\exp(x)&=\lim_{h\to0}\frac{\exp(x+h)-\exp(x)}{h}\\&=\lim_{h\to0}\frac{\exp(x)\exp(h)-\exp(x)}{h}\\&=\exp(x)\,\lim_{h\to0}\frac{\exp(h)-1}{h}\\&=\exp(x)\lim_{h\to0}\lim_{n\to\infty}\frac{\left(1+\frac{1}{n}\right)^{hn}-1}{h}\end{aligned}}

となります。いま

 \quad \displaystyle{\left(1+\frac{1}{n}\right)^{hn}\sim 1+ hn\times \frac{1}{n}}

に注意します。というのも、ベキ関数  (1+x)^a導関数

 \quad \displaystyle{\frac{d}{dx}(1+x)^a=a(1+x)^{a-1}}

で与えられ、再びベキ関数になり、  x=0 まわりのTaylor展開の公式

\quad \displaystyle{f(x)=f(0)+x\frac{d}{dx}f(0)+\frac{x^2}{2!}\frac{d^2}{dx^2}f(0)+\frac{x^3}{3!}\frac{d^3}{dx^3}f(0)+\dots}

にあてはめると

\quad \displaystyle{(1+x)^a=1+ax+a(a-1)\frac{x^2}{2!}+a(a-1)(a-2)\frac{x^3}{3!}+\dots}

となります。いま考えている状況は  x が十分小さい場合です。したがって、 x のベキが小さい項からの寄与が主要になるため、上の近似でよいことがわかります。実際、たとえば

\quad \displaystyle{x=0.2\qquad a=0.3}

などとしてみます。左辺を計算すると

x = 0.2
a = 0.3
(1+x)**a

により、その値は

\quad \displaystyle{(1+0.2)^{0.3}=1.0562199684392581}

ですが、右辺の最初の数項だけ取り入れた値は

[1+a*x, 1+a*x+a*(a-1)*x**2/2, 1+a*x+a*(a-1)*x**2/2+a*(a-1)*(a-2)*x**3/6]

から

\quad \displaystyle{\begin{aligned}1+ax&=1.06\\ 1+ax+a(a-1)\frac{x^2}{2!}&=1.0558\\ 1+ax+a(a-1)\frac{x^2}{2!}+a(a-1)(a-2)\frac{x^3}{3!}&=1.056276 \end{aligned}}

などとなって、関数値への主要な寄与が  x=0 まわりのTaylor展開のはじめのほうにあることがわかります。あるいは、いま考えている状況により即して、上の式

 \quad \displaystyle{\left(1+\frac{1}{n}\right)^{hn}\sim 1+ h}

\quad \displaystyle{n=20000\qquad h=0.0003}

としてみると

n = 20000
h = 0.0003
(1+1/n)**(h*n) - (1+h)

から

\quad \displaystyle{\left(1+\frac{1}{n}\right)^{hn}-\left( 1+ h\right)=3.750250066048011\times10^{-8}}

と計算され、その差が  h^2\sim h/n\sim 1/n^2 程度であることがわかります。よって上の式の2重極限の値が1と求まります。以上により、公式の成立がわかります。

実際、Pythonから指数関数の導関数を求めると

import sympy as sp

sp.var('x')
sp.diff(sp.exp(x), x, 1)

から公式が確かめられます。



Taylor展開

指数関数の微分公式から、一般の正整数  n に対して

\quad \displaystyle{\frac{d^n}{dx^n}\exp(x)=\exp(x)}

がわかります。また  \exp(0)=1 ですから、  x=0 まわりのTaylor展開の公式

\quad \displaystyle{f(x)=f(0)+x\frac{d}{dx}f(0)+\frac{x^2}{2!}\frac{d^2}{dx^2}f(0)+\frac{x^3}{3!}\frac{d^3}{dx^3}f(0)+\dots}

にあてはめると

\quad \displaystyle{\exp(x)=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\dots=\sum_{n=0}^\infty\frac{x^n}{n!}}

がわかります。実際、

import sympy as sp

sp.var('x')
sp.series(sp.exp(x), x, 0, 5)

から

\quad \displaystyle{\exp(x)=1+x+\frac{x^{2}}{2}+\frac{x^{3}}{6}+\frac{x^{4}}{24}+\mathcal{O}\left(x^{5}\right)}

を得ます。

多項式による指数関数の近似

Taylor級数による指数関数の展開式を  x n 次の項までで打ち切ったとき、その多項式がどれくらい指数関数を近似しているでしょうか?その様子を視覚的に実感するために、グラフを描画してみます。そこでPythonでどのようなことをする必要があるかというと、①指数関数のTaylor級数を各次数ごとに打ち切った多項式をSymPyで求め、②NumPyで各多項式の参照点での値を数値的に求めてリスト化し、③Matplotlibでそのリストをプロットする、という大まかに3つのステップを踏む必要があります。今回の記事のハイライトとして、この一連の作業をなるべくひとまとめに簡潔に書いてみましょう。

import sympy as sp   # ステップ①

sp.var('x')
def c(n,m): return sp.series(sp.exp(x),x,0,n+1).coeff(x,m)   # Taylor展開の係数

Pmax = 8   # Taylor展開の最高次数
Ex = []   # Taylor展開の各次数で打ち切って得られる多項式を要素とするリスト
for n in range(Pmax+1):
   Tayl = []   # Taylor展開の各項を要素とするリスト
   for m in range(n+1):
      Tayl.append(c(n,m)*x**m)   # リスト Tayl に要素を入れる
   def tayl(n): return sum(Tayl)
   Ex.append(tayl(n))   # リスト Ex に要素を入れる
   Tayl = []   # リスト Tayl を初期化する
   
import numpy as np   # ステップ②

N = 1000
xmin = 0
xmax = 3
h = (xmax - xmin)/N
xn = np.linspace( xmin, xmax, N+1)   # xmin から xmax までを N 等分した数(参照点)を要素とするリスト
def F(n,k): return Ex[n].subs(x, xn[k])   # リスト xn の各要素(参照点)での多項式の値を返す関数
G = [[],[],[],[],[],[],[],[],[]]   # 各次数ごとの「多項式の参照点での値のリスト」を要素とするリスト
for p in range(Pmax):
   for k in range(N+1):
      G[p].append(F(p,k))   # p 次までで展開を打ち切った多項式の、各参照点での値をリストに加える
      
import matplotlib.pyplot as plt   # ステップ③

plt.plot(xn, np.exp(xn))
plt.plot(xn, G[1])
plt.plot(xn, G[2])
plt.plot(xn, G[3])
plt.plot(xn, G[5])
plt.show()

以上から得られるのが次のグラフです:

f:id:pianofisica:20200925163305p:plain

青線が指数関数そのもののプロットで、橙、緑、赤、紫はそれぞれ1次、2次、3次、5次までの項で展開を打ち切って得られる多項式のプロットです。




キーワードPython、SymPy、NumPy、Matplotlib