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

当サイトではアフィリエイトプログラムを利用して商品を紹介しています。

ディラックのデルタ関数

ディラックデルタ関数は $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() が用意されています。

# SYMPY_DIRACDELTA

# In[1]

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$ となるように定義された不連続関数をヘヴィサイドの階段関数、あるいはヘヴィサイド関数とよびます。
 \[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$ における値を返します。

# NUMPY_HEAVISIDE

# In[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$ で値をもたないヘヴィサイド関数です。

# SYMPY_HEAVISIDE

# In[1]

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)

# ヘヴィサイド関数を定義
y = sympy.Heaviside(x)

# ヘヴィサイド関数を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モジュールを使ってヘヴィサイド関数のグラフを描いてみます。

# In[2]

from sympy.plotting import plot

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

# ヘヴィサイド関数
y = sympy.Heaviside(x)

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

Python ヘビサイドの階段関数グラフ

コメント

  1. HNaito より:

    HEAVISIDE_01プログラムでは、y = sympy.Heaviside(x) の定義が必要だと思いますので、ご確認ください。

    • あとりえこばと より:

      申し訳ないです。m(_ _)m
      HEAVISIDE_01 に y=sympy.Heaviside(x) を追加しておきました。

  2. HNaito より:

    次の記事「ディオファントス方程式」でも、Colab の SymPy 1.7.1 では diop_solve のインポートエラーが発生し、1.11 にアップグレード後は動作OKでした。
    下記は誤植と思われますので、ご確認ください。
    EQUATION_02 プログラムの実行結果で、(2t_0+t, -3t_0-5) → {(2t_0+t, -3t_0-5)}
    EQUATION_03 プログラムで、最後に print(sol) 追加。
    EQUATION_05 プログラムで、from sympy.solvers.diophantine → from sympy.solvers.diophantine.diophantine

    • あとりえこばと より:

      修正しました。
      ありがとうございます。

      • HNaito より:

        SymPy 1.7.1 でも下記のように修正することで、実行できました。

        EQUATION_01, 03, 04 プログラムで、
        from sympy.solvers.diophantine import diop_solve →
        from sympy.solvers.diophantine.diophantine import diop_solve

        • あとりえこばと より:

          ありがとうございます。SymPyのバージョンによって不具合が生じることをできるだけ避けるために、EQUATION_01, 03, 04 のコードを
          from sympy.solvers.diophantine.diophantine import diop_solve
          に修正しておきました。助かりました。m(_ _)m

          • HNaito より:

            「ディオファントス方程式」の説明文の最後のほうで、コード修正に合わせて文中でも修正が必要だと思います。
            sympy.solvers.diophantine.diop_solve( ) →
            sympy.solvers.diophantine.diophantine.diop_solve( )

  3. HNaito より:

    次の記事「ディオファントス方程式」で、下記は誤植と思われますのでご確認ください。
    EQUATION_01 プログラムの上の文で、(x, y)=(2t+5, −3t+5) → (x, y)=(2t+5, −3t-5)

    • あとりえこばと より:

      細かい所まで丁寧に見てくださって、ありがとうございます。m(_ _)m
      本当に助かります。

  4. HNaito より:

    3x + 2y =5 の媒介変数方程式を SYMPY_PLOTTING_MODULE_02 プログラムを参考にして描いてみました。
    x = 2t + 5
    y = -3t – 5
    いつも 2 変数: x, y の一方を入力として、他方を出力とみなして描いていたグラフが、別のパラメータを使っても描けるというのが何か新鮮でした。一般解が求められる不定方程式の例が他にもあるようでしたら、いくつかご紹介いただけるとありがたいです。

    • あとりえこばと より:

      有名な二次不定方程式として、ピタゴラス方程式 x^2 + y^2 = z^2 があります。いわゆるピタゴラスの定理(三平方の定理)を満たす (x, y, z) を求める方程式です (x, y, z を三角形の辺の長さと考えるなら、直角三角形となるような辺の組を見つけることに相当します)。この方程式の自然数解 (a, b, c) の組をピタゴラス数とよびます。ピタゴラス数を求めることは単位円上の有理点を求める問題に帰着できますが、結果だけ記すと、任意整数 m, n を用いて、(2mn, m^2-n^2, m^2+n^2) または (m^2-n^2, 2mn, m^2+n^2) で表されることが知られています。たとえば、m = 2, n = 1 とすると、
       (2mn, m^2-n^2, m^2+n^2) = (4, 3, 5)
      を得ますが、これは確かにピタゴラスの定理
       4^2 + 3^2 = 5^2
      を満たしています。当然ですが、(4, 3, 5) の 4 と 3 を入れ替えた (3, 4, 5) もピタゴラス数です。ちなみに、1995年にアンドリュー・ワイルズによって証明された「フェルマーの最終定理」により、3 以上の自然数 n について、
       x^n + y^n = z^n
      を満たす自然数 (x, y, z) は存在しないことがわかっています。

      • HNaio より:

        diophantine(x**2 + y**2 – z**2) で一般解が求められました。diophantine( ) で線形、ピタゴラス、平方和不定方程式の一般解は求められると考えてよいのでしょうか。

        • あとりえこばと より:

          ディオファントス方程式
           f(x_1, x_2, …) = 0
          に分類される方程式で、一般解が存在するなら sympy.diophantine() で求められると思います。たとえば、4 変数ピタゴラス方程式
           x^2 + y^2 + z^2 = w^2
          は一般解が存在するので、diop_solve(x**2 + y**2 + z**2 – w**2) は解を返します(ぜひ試してみてください)。しかし、そもそも不定方程式には解が存在しないことがあります。たとえば、ペルの方程式の一種
           x**2 – 3 * y = -7
          は解をもたないので、sympy.diophantine(x**2 – 3 * y + 7) は空の集合を返します。このようなケースがあるため、右辺の -7 を任意自然数 k に置き換えた
           x**2 – 3 * y = k
          は一般解をもたないことになります (k の値によっては解が存在しないので一般解を k を用いて表すことができません)。したがって、sympy.diophantine(x**2 – 3 * y – k) はエラーとなります。ちなみに、diophantine() と diop_solve() は、ほぼ同じ機能をもつ関数ですが、sympy.diophantine() が解のタプルを set に入れて返すのに対し、diop_solve() は (解が存在しない場合を除いて) タプルだけで返す、という違いがあります。細かい話になりますが、実はここで言うタプルは Python の組み込みのタプルとは異なり、sympy.core.containers.Tuple というクラスのオブジェクトです。

          • HNaito より:

            ご回答ありがとうございました。これで、線形不定方程式、2元2次不定式、ピタゴラス方程式に対して、sympy の diop_solve, find_ND, diop_ND, transformation_to_DN のツールの使い方を学べました。

  5. あとりえこばと より:

    【AI連載小説】科学とコードの交差点(42)「ヘヴィサイド関数」
     
    大学の自習室で、開誠と明信、馬締美純(まじめみすず)がヘヴィサイド関数についての議論を始めました。三人は物理学を学ぶ仲間で、数学的なトピックについて熱心に話し合うことが日常的な活動となっていました。
    開誠:「今日はヘヴィサイド関数について話そう。物理学でよく使われるよな」
    美純:「そう、ヘヴィサイド関数はステップ関数とも呼ばれ、時間によって急激な変化があるシステムをモデル化するのに便利だよね」
    明信:「確かに、特に制御理論や信号処理の分野で頻繁に見かけるよね」
    開誠:「ヘヴィサイド関数は数学的には単純だけど、具体的な物理現象においてどのように応用されるか、それが面白いところだな」
    美純:「実際の応用例として、電気回路のスイッチング動作や慣性制御系の動作などが挙げられるね」
    開誠:「ヘヴィサイド関数を数学的にどう表現するかも面白い」
    明信:「そうだね、数学的な表現もいくつかある。単純に階段関数で表す方法もあれば、定義域を考慮して別の方法もある」
    美純:「あと、ディラックのデルタ関数を使っても表現できることがあるよね」
    三人はヘヴィサイド関数の定義や応用について情熱的に議論し、ノートやホワイトボードに数式を書きながら理解を深めていきました。物理学と数学が交差する美しい瞬間が、その自習室に広がっていました。