微分・偏微分

1. Python 数値微分
- 1.1. scipy.misc.derivative()
- 1.2. simpy.diff()

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")

python derivative() 関数で求めた sinx の導関数データ
　
　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)

f(x)=1/(1+x^2) の1階導関数と2階導関数
　
　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)

\[f{\left(x \right)} \frac{d}{d x} g{\left(x \right)} + g{\left(x \right)} \frac{d}{d x} f{\left(x \right)}\]

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

# In[3]

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

\[f{\left(x \right)} \frac{d^{2}}{d x^{2}} g{\left(x \right)} + g{\left(x \right)} \frac{d^{2}}{d x^{2}} f{\left(x \right)} + 2 \frac{d}{d x} f{\left(x \right)} \frac{d}{d x} g{\left(x \right)}\]

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

# In[4]

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

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

\[\frac{f{\left(x \right)} \frac{d}{d x} g{\left(x \right)} - g{\left(x \right)} \frac{d}{d x} f{\left(x \right)}}{f^{2}{\left(x \right)}}\]

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

# In[5]

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

\[\frac{d}{d g{\left(x \right)}} f{\left(g{\left(x \right)} \right)} \frac{d}{d x} g{\left(x \right)}\]

　次は具体的な計算例です。$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()

\[x \cos{\left(x \right)} + \sin{\left(x \right)}\]

sympy diiメソッドによるxsinxの微分
　
　$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)

\[- \frac{x}{y}\]

　多変数関数の偏微分も簡単に実行できます。
　例として次のような二変数関数を定義します。

\[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)

\[\frac{\frac{3 x^{2}}{2} - \frac{y \sin{\left(x y \right)}}{2}}{\sqrt{x^{3} + \cos{\left(x y \right)}}}\]

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

# In[9]

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

\[- \frac{x \sin{\left(x y \right)}}{2 \sqrt{x^{3} + \cos{\left(x y \right)}}}\]

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

# In[10]

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

\[\frac{- 2 x y \cos{\left(x y \right)} + \frac{x \left(3 x^{2} - y \sin{\left(x y \right)}\right) \sin{\left(x y \right)}}{x^{3} + \cos{\left(x y \right)}} - 2 \sin{\left(x y \right)}}{4 \sqrt{x^{3} + \cos{\left(x y \right)}}}\]

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)

コメント（5件）

HNaito より:

2023年3月4日 9:34 PM

DERIVATIVE_01 In[2] プログラムでは、 lambda 関数の使い方を一つ覚えさせていただきました。
“`
def f(x):
　return x**2 +1
“`
と外部で定義するよりも、実行時間もだいぶ短くなるのでしょうか。
　
DERIVATIVE_01 In[5] プログラムの実行結果では、scipy.misc.derivative( ) は非推奨のメッセージが出ました。
DERIVATIVE_02 In[7] プログラムでは、シンボル x, y の再定義( 初期化 ) が必要になりました。

返信
- あとりえこばとより:
  
  2023年3月5日 1:27 PM
  
  def 構文とlambda 式で処理速度に大きな違いはないと思います。念のために、Google Colab で計測した結果を載せておきます。
  ――――――――――――――――――――――――――
  def derivative(func, a, h):
  return (func(a + h) – func(a – h)) / (2 * h)
  　
  def f(x):
  return x**2 + 1
  　
  h = 1e-5
  　
  %timeit derivative(lambda x : x**2 + 1, 1, h)
  # 1.95 µs ± 740 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
  　
  %timeit derivative(f, 1, h)
  # 1.77 µs ± 327 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
  ――――――――――――――――――――――――――
  　lambda 式は記述を簡素にすると同時に、使い捨て関数なので関数オブジェクトがより速くメモリから解放されるという利点があります。ですから、一回限りしかその関数を使わない場合は lambda 式を使ったほうがいいです。
  　
  　おっしゃるように、scipy.misc.derivative() は非推奨となったようです。SciPy v1.12.0 から完全に削除されるようですね。実行結果の警告文には、外部ライブラリの findiff を使うことを勧めています。本文中に非推奨であることを付け加えておきます。また後で findiff の基本的な使い方を確認したうえで加筆したいと思います。
  
  返信
通りすがりより:

2023年5月2日 1:04 PM

重箱の角で恐縮です。
(3)式の分子の
f(a+h)−(a−h)
は
f(a+h)−f(a−h)
（右項もf()）
となりますでしょうか。

返信
- あとりえこばとより:
  
  2023年5月2日 1:50 PM
  
  ありがとうございます。
  ご指摘いただいた箇所を修正しておきました。
  助かりました。m(_ _)m
  
  返信
あとりえこばとより:

2023年8月17日 9:15 PM

【GPTプログラミング】任意関数の微分係数

　
[Q] Pythonで任意関数の微分係数を計算する関数を作ってください。ただし、ライブラリを使わないでください。
　
[A] 以下に、Pythonを使用して任意の関数の微分係数を計算する関数を示します。この関数は、中心差分法を使用して数値微分を行います。
```
def differentiate(func, x, h=1e-5):
    """
    Calculate the derivative of a given function at a specific point using the central difference method.
    
    Parameters:
    func (function): The function for which the derivative needs to be calculated.
    x (float): The point at which the derivative is to be computed.
    h (float): The step size for numerical differentiation. Default is 1e-5.
    
    Returns:
    float: The derivative of the function at the given point.
    """
    derivative = (func(x + h) - func(x - h)) / (2 * h)
    return derivative

# Example usage:
def quadratic(x):
    return x ** 2

derivative_at_2 = differentiate(quadratic, 2)
print("Derivative of x^2 at x = 2:", derivative_at_2)
```
　
この例では、differentiate関数を使用して、与えられた関数の微分係数を計算しています。関数は中心差分法を使用して数値微分を行い、指定された点における微分係数を返します。ご自身の関数に対してこの関数を使用する際は、quadratic関数のように、微分したい関数を定義し、それをdifferentiate関数に渡してください。また、必要に応じてhの値を調整して、求める精度を調整できます。

返信