![](https://www.math-joy-life.com/wp-content/uploads/2023/12/icon1.webp)
![](https://www.math-joy-life.com/wp-content/uploads/2023/12/icon2.webp)
目次
SciPyで特殊関数
Pythonで特殊関数を出力するためにSciPyを使用します。特殊関数とは、数学や物理学においてよく使用される関数で、ガンマ関数、ベッセル関数、ベータ関数、誤差関数などが含まれます。それぞれの特殊関数のPythonコードを分けて、説明を追加します。
![](https://www.math-joy-life.com/wp-content/uploads/2023/12/icon2.webp)
![](https://www.math-joy-life.com/wp-content/uploads/2023/12/icon1.webp)
Pythonプログラム
scipyのinstall
まず、SciPyライブラリをインストールします。以下のコマンドを使用します。
pip install scipy
ガンマ関数
\[ \Gamma(z) = \int_0^\infty t^{z-1} e^{-t} \, dt \]
ガンマ関数は、階乗に関連がある関数です。SciPyのscipy.special
モジュールを使用して計算します。gammaメソッドを呼び出すことで計算ができます。
import scipy.special as sp
# ガンマ関数の計算
gamma_value = sp.gamma(5)
print(f"Gamma(5) = {gamma_value}")
Gamma(5) = 24.0
![](https://www.math-joy-life.com/wp-content/uploads/2023/12/icon2.webp)
ベータ関数
\[ B(x, y) = \int_0^1 t^{x-1} (1-t)^{y-1} \, dt \]
または
\[ B(x, y) = \frac{\Gamma(x) \Gamma(y)}{\Gamma(x + y)} \]
ベータ関数は、ガンマ関数を使用して定義される関数です。betaメソッドを呼び出すことで計算ができます。
import scipy.special as sp
# ベータ関数の計算
beta_value = sp.beta(2, 3)
print(f"Beta(2, 3) = {beta_value}")
Beta(2, 3) = 0.08333333333333333
![](https://www.math-joy-life.com/wp-content/uploads/2023/12/icon2.webp)
エアリー関数
\[ \text{Ai}(x) = \frac{1}{\pi} \int_0^\infty \cos\left(\frac{t^3}{3} + xt\right) \, dt \]
\[ \text{Bi}(x) = \frac{1}{\pi} \int_0^\infty \left[ \exp\left(-\frac{t^3}{3} + xt\right) + \sin\left(\frac{t^3}{3} + xt\right) \right] dt \]
エアリー関数は、線形微分方程式の解として現れる関数です。Ai(x)とBi(x)を計算してグラフにプロットします。airyメソッドを呼ぶことで計算することができます。
import scipy.special as sp
import numpy as np
import matplotlib.pyplot as plt
# エアリー関数の計算とプロット
x = np.linspace(-10, 10, 400)
ai, aip, bi, bip = sp.airy(x)
plt.plot(x, ai, label='Ai(x)')
plt.plot(x, bi, label='Bi(x)')
plt.title("Airy functions")
plt.ylim(-1,1)
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid(True)
plt.show()
ベッセル関数(第二種、第三種)
第二種ベッセル関数
\[ Y_\nu(x) = \frac{J_\nu(x) \cos(\nu \pi) - J_{-\nu}(x)}{\sin(\nu \pi)} \]
第三種ベッセル関数(修正ベッセル関数)
\[ K_\nu(x) = \frac{\pi}{2} \frac{I_{-\nu}(x) - I_\nu(x)}{\sin(\nu \pi)} \]
ベッセル関数は、波動方程式の解として現れる関数です。ここでは第二種と第三種のベッセル関数を計算して、グラフにプロットします。yvとkvメソッドを呼ぶことで計算することができます。
import scipy.special as sp
import numpy as np
import matplotlib.pyplot as plt
# ベッセル関数の計算とプロット
x = np.linspace(0, 10, 100)
y2 = sp.yv(2, x) # 第二種ベッセル関数
y3 = sp.kv(3, x) # 第三種ベッセル関数(修正ベッセル関数)
plt.plot(x, y2, label='Y2(x)')
plt.plot(x, y3, label='K3(x)')
plt.title("Bessel functions (second and third kind)")
plt.xlabel("x")
plt.ylabel("y")
plt.ylim(-1,1)
plt.legend()
plt.grid(True)
plt.show()
ゼータ関数
\[ \zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s} \]
解析接続により次の積分表示もあります。
\[ \zeta(s) = \frac{1}{\Gamma(s)} \int_0^\infty \frac{t^{s-1}}{e^t - 1} \, dt \]
リーマンゼータ関数は、解析数論で重要な役割を果たします。zetaメソッドを呼ぶことで計算することができます。
import scipy.special as sp
# ゼータ関数の計算
zeta_value_2 = sp.zeta(2)
zeta_value_m_2 = sp.zeta(-2)
print(f"Zeta(2) = {zeta_value_2}")
print(f"Zeta(-2) = {zeta_value_m_2}")
Zeta(2) = 1.6449340668482264
Zeta(-2) = 0.0
楕円関数
第一種完全楕円積分
\[ K(k) = \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1 - k^2 \sin^2(\theta)}} \]
第二種完全楕円積分
\[ E(k) = \int_0^{\frac{\pi}{2}} \sqrt{1 - k^2 \sin^2(\theta)} \, d\theta \]
楕円積分は、楕円関数や楕円運動に関連する積分です。ellipe,ellipkメソッドを呼ぶことで計算することができます。
import scipy.special as sp
import numpy as np
import matplotlib.pyplot as plt
# 楕円積分の計算とプロット
x = np.linspace(0, np.pi/2, 100)
ellipk_value = sp.ellipk(x)
ellipe_value = sp.ellipe(x)
plt.plot(x, ellipk_value, label='Elliptic integral K(x)')
plt.plot(x, ellipe_value, label='Elliptic integral E(x)')
plt.title("Elliptic integrals")
plt.xlabel("x")
plt.ylabel("Integral value")
plt.legend()
plt.grid(True)
plt.show()
ルジャンドル関数
\[ P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} \left[(x^2 - 1)^n\right] \]
ルジャンドル多項式は、球面調和関数などの特殊関数に現れます。
![](https://www.math-joy-life.com/wp-content/uploads/2023/12/icon2.webp)
import scipy.special as sp
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-1, 1, 100)
for i in range(1,6):
# ルジャンドル多項式の計算とプロット
legendre_p = sp.legendre(i)
y_legendre = legendre_p(x)
plt.plot(x, y_legendre, label=f'P{i}(x)')
plt.title("Legendre polynomial")
plt.xlabel("x")
plt.ylabel("y_legendre")
plt.legend()
plt.grid(True)
plt.show()