デルタ関数とヘヴィサイドの階段関数

ディラックのデルタ関数

デルタ関数の定義

　ディラックのデルタ関数 (Dirac delta function)は $x=0$ で無限大、その他の区間で $0$ の値をもち、全区間での積分が $1$ となるように定義された特殊な関数です。
　
\[\delta(x)=\begin{cases}\infty & (x=0)\\[6pt]
0 & (x\neq 0)\end{cases}\]\[\int_{-\infty}^{\infty}\delta(x)dx=1\]
　デルタ関数を原点を含む区間で積分すると、後述するヘヴィサイドの階段関数が得られます。
　
\[\int_a^{x}\delta(x)dx=H(x)\quad (a\lt 0)\]

sympy.DiracDelta()関数

　SymPyパッケージにはデルタ関数の値を返す sympy.DiracDelta() が用意されています。

import sympy

# 記号xを定義
sympy.var("x")

# y = δ(x)
y = sympy.DiracDelta(x)

dm = y.subs(x, -1)
d0 = y.subs(x, 0)
dp = y.subs(x, 1)

d_list = [dm, d0, dp]

print("[δ(-1), δ(0), δ(1)] =", d_list)

# デルタ関数の不定積分
i = sympy.integrate(y, x)

print("∫δ(x)dx =", i)

[δ(-1), δ(0), δ(1)] = [0, DiracDelta(0), 0]
∫δ(x)dx = Heaviside(x)

ヘヴィサイドの階段関数

階段関数の定義

　$x\lt 0$ で $0$, $0\lt x$ で $1$ となるように定義された不連続関数をヘヴィサイドの階段関数 (Heaviside step function)、あるいはヘヴィサイド関数 (Heaviside function) とよびます。
　
\[H(x)=\begin{cases}0 & (x\lt 0)\\[6pt]
1 & (x\gt 0)\end{cases}\]
　一般に $H(0)$ は定義されませんが、$x=0$ で値をとるように定めることもあり、その場合は適当な実数 $c$ を使って
　
\[H_c(x)=\begin{cases}0 & (x\lt 0)\\[6pt]
c & (x= 0)\\[6pt]
1 & (x\gt 0)\end{cases}\]
のように定義します。このとき、$H_{1/2}(x)$ は符号関数 $\mathrm{sgn}(x)$ と次のような関係で結びつくことになります。
　
\[H_{1/2}(x)=\frac{1+\mathrm{sgn}(x)}{2}\]
　ヘヴィサイド関数を微分するとデルタ関数となります。
　
\[\frac{dH(x)}{dx}=\delta (x)\]
　これはヘヴィサイド関数の原点における接線の傾きが $\infty$ となることを表しています。このように、デルタ関数を用いると不連続関数の微分を求めることができます。
　

numpy.heaviside()

　NumPy にはユニバーサル関数 numpy.heaviside() が用意されています。
　numpy.heaviside() は上記の $H_c(x)$ の定義にしたがうヘヴィサイド関数であり、第 2 引数で $x=0$ でとるべき値を指定します。たとえば、numpy.heaviside([-1, 0, 1], 3) と記述した場合、
　
\[H_3(x)=\begin{cases}0 & (x\lt 0)\\[6pt]
3 & (x= 0)\\[6pt]
1 & (x\gt 0)\end{cases}\]
によって定義されるヘヴィサイド関数の $x=-1,\ 0,\ 1$ における値を返します。

# ヘヴィサイド関数

import numpy as np

x = [-1, 0, 1]

a = np.heaviside(x, 3)

print(a)

[0. 3. 1.]

sympy.Heaviside()関数

　SymPy には sympy.Heaviside()関数が組込まれています。
　sympy.Heaviside() は $x=0$ で値をもたないヘヴィサイド関数です。

# ヘヴィサイド関数

import sympy

# 記号xを定義
sympy.var("x")

# H(-1)
hm = sympy.Heaviside(-1)

# H(0)
h0 = sympy.Heaviside(0)

# H(1)
hp = sympy.Heaviside(1)

h_list = [hm, h0, hp]

print("[H(-1), H(0), H(1)] =", h_list)

# ヘヴィサイド関数をxで微分
dy_dx = sympy.diff(y, x)

print("dH(x)/dx =", dy_dx)

[H(-1), H(0), H(1)] = [0, Heaviside(0), 1]
dH(x)/dx = DiracDelta(x)

　sympy.plottingモジュールを使ってヘヴィサイド関数のグラフを描いてみます。

# ヘヴィサイド関数のグラフを描画

import sympy
from sympy.plotting import plot

# 記号xを定義
sympy.var("x")

# y = H(x)
y = sympy.Heaviside(x)

#[-5, 5]の範囲でヘヴィサイドの階段関数をプロット
p = plot(y, (x, -5, 5), xlabel = "x", ylabel = "y")