数値微分

この記事では、微分係数の定義と数値微分の考え方、Pythonにおける中心差分公式の実装例、SciPyやSymPyを用いて関数を微分する方法などについて解説します。

【Python】微分係数と中心差分公式

関数 $y=f(x)$ について、点 $a$ における微分係数は
 \[f'(a)=\lim_{h\rightarrow 0}\frac{f(a+h)-f(a)}{h}\tag{1}\]
によって与えられますが、数値計算においては厳密に $h\rightarrow 0$ の極限をとることはできません。そこで $h$ を十分に小さな有限数として、
 \[\frac{f(a+h)-f(a)}{h}\tag{2}\]
を計算すれば微分係数の近似値を得られます。(2) は前進差分公式とよばれますが、一般的には
 \[f'(a)=\frac{f(a+h)-f(a-h)}{2h}\tag{3}\]
で表される中心差分公式のほうが精度が高く、数値計算で多用されています。実装は簡単です。

# PYTHON_DERIVATIVE_01

# In[1]

# 微分係数を返す関数を定義
def derivative(func, a, h):
    return (func(a + h) - func(a - h)) / (2 * h)

試しに、二次関数 $f(x)=x^2+1$ の $x=1$ における微分係数を求めてみます。$f'(x)=2x$ なので真値は $f'(1)=2$ です。

# In[2]

# h=0.00001
h = 1e-5

# 二次関数y=x**2+1のx=1における微分係数(真値は2)
val = derivative(lambda x : x**2 + 1, 1, h)

print(val)
# 2.000000000002

若干の誤差がありますが、まずまずの精度で微分係数を計算できました。次は math モジュールをインポートして、$f(x)=\sin x$ の $x=\pi/3$ における微分係数を求めてみましょう。

# In[3]
import math

# y=sinxのx=π/3における微分係数(真値は0.5)
val = derivative(math.sin, math.pi/3, 1e-5)

print(val)
# 0.4999999999921733

真値は $0.5$ なので、やはり高い精度で微分係数を求めることができました。NumPy をインポートすることで、derivative() の引数 a に連続データを渡すことができます。このとき、derivative() は func の近似的な導関数を返すといえます。$f(x)=\sin x$ の導関数を得て、Matplotlib でプロットしてみましょう。

# In[4]

import numpy as np
import matplotlib.pyplot as plt

# 描画領域の作成と設定
fig = plt.figure()
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("Derivative of sinx")
ax.set_xlabel("x", fontsize = 16)
ax.set_ylabel("y", fontsize = 16)

# 連続データの作成
x = np.linspace(-2*np.pi, 2*np.pi, 65)

# derivativeに連続データxを渡すことで、
# 各点の微分係数、つまり導関数を得る(真関数はcosx)
dfunc = derivative(np.sin, x, h)

# 導関数をプロット
ax.plot(x, dfunc, color="blue")

SciPy には任意関数 func を受け取って微分係数を返す derivative() 関数が用意されています (SciPy v1.10.0 で実行したところ、この関数は非推奨となっており、SciPy v1.12.0 から完全に削除される予定です)。この関数は引数 n で微分の階数を指定できます (n 次微分係数を取得できます)。一例として、
 \[f(x)=\frac{1}{1+x^2}\tag{4}\]
の一次導関数と二次導関数のデータを得て、グラフにプロットしてみます。

# In[5]

from scipy.misc import derivative

# 描画領域の作成と設定
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("f(x)=1/(1+x**2), f'(x), f''(x)",
             fontsize=15, pad=10)
ax.set_xlabel("x", fontsize = 16)
ax.set_ylabel("y", fontsize = 16)
ax.set_aspect('equal')
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)

# 連続データの作成
x = np.linspace(-2*np.pi, 2*np.pi, 129)

def func(x):
    return 1 / (1 + x**2)

# funcの一次導関数
dfunc_1 = derivative(func, x, h)

# funcの二次導関数
dfunc_2 = derivative(func, x, h, n=2)

# 原始関数と導関数をプロット
ax.plot(x, func(x), color="red", label="f(x)")
ax.plot(x, dfunc_1, color="blue", label="f'(x)")
ax.plot(x, dfunc_2, color="green", label="f''(x)")

ax.legend(fontsize=12, shadow=True)

SymPy を用いると解析的な微分を実行できます。私たちがノートで計算をするのと同じように微分を扱えるということです。最初に SymPy をインポートし、数式で使いそうな文字を定め、f と g を x の抽象的な関数 (Functionオブジェクト) として定義しておきます。

# PYTHON_DERIVATIVE_02

# In[1]

import sympy as sym

# 実行結果をTeXで表示する
sym.init_printing()

# 数式で使用するシンボルを定義
sym.var('x y', real=True)

# fとgを関数記号として定義
f = sym.Function('f')
g = sym.Function('g')

SymPy の数式や Function オブジェクトは導関数を得る diff() メソッドを備えています。メソッドの第 1 引数には微分する変数を指定します。積の微分 $(fg)’$ を確認してみましょう。

# In[2]

# 積の微分公式(fg)'
(f(x) * g(x)).diff(x)

第 2 引数には導関数の次数 (何回微分するか) を指定します。$fg$ の二次導関数を得てみます。

# In[3]

# 積の2階微分
(f(x) * g(x)).diff(x, 2)

次は商の微分公式 (g/f)’ を確認してみます。

# In[4]

# 商の微分公式
expr = (g(x) / f(x)).diff(x)

# 式を整理する
sym.simplify(expr)

次は合成関数 $f(g(x))$ の微分公式です。

# In[5]

# 合成関数の微分
f(g(x)).diff(x)

次は具体的な計算例です。$y=x\sin x$ を微分して、もとの関数と導関数のグラフを同時プロットしてみます。

# In[6]

# y=xsinx
y = x * sym.sin(x)

# xsinxの導関数
dy_dx = y.diff(x)

# 導関数を表示
display(dy_dx)

# 関数と導関数をプロット
ax = sym.plot(y, dy_dx, (x,-5,5), legend=True, show=False)
ax[0].line_color = "navy"
ax[1].line_color = "orangered"

ax.show()

$f(x,y)=0$ の形で $x$ と $y$ の関係が定められた関数を陰関数とよびます。たとえば、よく知られた単位円の方程式 $x^2+y^2=1$ も陰関数です。このような関数について導関数 $dy/dx$ を求めたいときは、idiff() メソッドを用います：

# In[7]

# x,yを再定義
sym.var('x y', real=True)

# 陰関数x**2+y**2=1のdy/dxを得る
sym.idiff(y ** 2 + x ** 2 - 1, y, x)

多変数関数の偏微分も簡単に実行できます。
例として次のような二変数関数を定義します。
 \[f(x,y)=\sqrt{x^{3} + \cos{\left(x y \right)}}\tag{5}\]

# In[8]

# 2変数関数を定義
func = sym.sqrt(x**3 + sym.cos(x*y))

func.diff(x) と記述すると、変数 x で偏微分します。

# In[9]

# funcを変数xで偏微分
func.diff(x)

同様に func.diff(y) と記述すると、変数 y で偏微分を実行します。

# In[9]

# funcを変数yで偏微分
func.diff(y)

func.diff(x,y) は x と y で 1 回ずつ偏微分した、二次偏導関数を返します。

# In[10]

# funcを変数xとyで1回ずつ偏微分
func.diff(x, y)

scipy.misc.derivative()

scipy.misc.derivative() は中心差分公式を使って、受け取った関数の指定点における微分係数を計算します。

scipy.misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3)

func には関数オブジェクト、x0 には微分係数を求める点 (中心差分公式のa) を渡します。dx は中心差分公式における h の値であり、デフォルトでは 1 に設定されていますが、より小さな値に設定しないと意味のある精度の値が得られません (具体的な値はサンプルコードを参照してください)。n は微分の階数であり、デフォルトでは 1 に設定されています。
 
3次関数 $f(x)=x^3$ の 1 階微分係数を求めてみます。

# SCIPY_MISC_DERIVATIVE

# In[1]

from scipy.misc import derivative

# 3次関数を定義
f = lambda x : x**3

# x=1における1階微分係数
df_1 = derivative(f, 1, dx=1e-5)

print(df_1)
# 3.000000000097369

n = 2 を指定すると 2 階微分係数を計算します。

# In[2]

# x=1における2階微分係数
df_2 = derivative(f, 1, dx=1e-5, n=2)

print(df_2)
# 6.000002716888274

simpy.diff()

SymPy には任意関数の導関数を返す diff() 関数が用意されています。引数に渡した関数 func を 変数 var で n階微分するときは次の構文を記述します。

diff(func, var, n)

引数 n を省略すると 1 階微分します。

# SYMPY_DIFF

# In[1]

import sympy

# 記号a, b, c, x, yを定義
sympy.var('a b c x y')

# 関数f(x)を定義
f = a*x**2 + b*x + c

# f(x)をxで微分
dfdx1 = sympy.diff(f, x)

# f(x)をxで2階微分
dfdx2 = sympy.diff(f, x, 2)

# f(x)をxで3階微分
dfdx3 = sympy.diff(f, x, 3)

print("dfdx1 = {}".format(dfdx1))
print("dfdx2 = {}".format(dfdx2))
print("dfdx3 = {}".format(dfdx3))

# dfdx1 = 2*a*x + b
# dfdx2 = 2*a
# dfdx3 = 0

引数 var に複数の記号 x1, x2, … を渡して導関数を得る構文もあります。

diff(func, x1[, x2, ...])

次のサンプルコードでは、$f(x)=\cos(xy)$ を $x$ と $y$ で 1 回ずつ偏微分しています。

# In[2]

# f(x)を定義
f = sympy.cos(x*y)

# f(x)をxとyで1回ずつ偏微分
df_dx_dy = sympy.diff(f, x, y)

print("df_dx_dy = {}".format(df_dx_dy))
# df_dx_dy = -(x*y*cos(x*y) + sin(x*y))

diff()関数の代わりに diff()メソッドを使って微分することもできます。

# In[3]

# f(x) = xexp(x)
f = x * sympy.exp(x)

# f(x)をxで微分
df_dx = f.diff(x, 1)

print(df_dx)
# x*exp(x) + exp(x)