NumPy と SciPy は、Python で数値計算を行うための強力なライブラリです。NumPy は配列操作や数値演算を効率的に行うためのライブラリで、SciPy は科学技術計算のためのライブラリであり、NumPyの機能を拡張しています。
NumPy を使用するには、まずライブラリをインポートします。
慣習的に np といった省略した名前にします。
import numpy as np
NumPy の中心的なデータ構造は、ndarray
と呼ばれる多次元配列です。
# 1次元配列(ベクトル)
a = np.array([1, 2, 3, 4, 5])
# 2次元配列(行列)
b = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
# 配列の表示
print("1次元配列:", a)
print("2次元配列:", b)
配列の形状、サイズ、要素へのアクセスなどの基本操作を学びます。
# 配列の形状
print("形状:", b.shape)
# 配列のサイズ(要素数)
print("サイズ:", b.size)
# 特定の要素へのアクセス
print("2行3列目の要素:", b[1, 2])
# スライシング
print("1行目:", b[0, :])
print("2列目:", b[:, 1])
NumPy を使って基本的な数学演算を行います。
# 配列の要素ごとの演算
c = np.array([10, 20, 30, 40, 50])
print("要素ごとの足し算:", a + c)
print("要素ごとの掛け算:", a * c)
# 行列の積
d = np.dot(b, b)
print("行列の積:", d)
配列の形状変更や結合、分割などの操作を学びます。
# 配列の形状を変更
e = np.reshape(a, (5, 1))
print("形状変更:", e)
# 配列の結合
f = np.concatenate((a, c))
print("配列の結合:", f)
# 配列の分割
g = np.array_split(f, 5)
print("配列の分割:", g)
SciPy には多くのサブモジュールによって構成されるため目的に合わせて読み込む必要があります。
以下はその例です。
from scipy import linalg # 線形代数用
from scipy import integrate # 数値積分用
from scipy import fft # フーリエ変換等用 (fftpack といったものもあるが古いものなのでこちらを使用する)
SciPyのlinalg
モジュールを使って、線形代数の基本操作を行います。
# 2x2の行列を作成
A = np.array([[1, 2], [3, 4]])
# 行列式の計算
det_A = linalg.det(A)
print("行列式:", det_A)
# 逆行列の計算
inv_A = linalg.inv(A)
print("逆行列:
", inv_A)
# 連立方程式の解
b = np.array([5, 6])
x = linalg.solve(A, b)
print("連立方程式の解:", x)
SciPyのintegrate
モジュールを使って、数値積分を行います。
from scipy import integrate
# 1次元関数の定義
def f(x):
return x**2
# 0から1までの積分
integral, error = integrate.quad(f, 0, 1)
print("積分結果:", integral)
SciPy の fft
モジュールを使って、離散フーリエ変換を行います。 fft
の中には fft が先頭につく関数が多くありますのでscipy.fft
からインポートしています。
以下を実行すると周波数領域で -1 と +1 でピークになるグラフが表示されます(グラフは matplotlib を使用していますので詳しくはそちらを確認してください)。
from scipy.fft import fft, fftfreq, fftshift
import numpy as np
import matplotlib.pyplot as plt
# 刻み幅 0.1、データ数256個、周波数 1 Hz、波形が Sin 波のデータを用意します
timestep = 0.1
data_n = 256
t = timestep*np.arange(data_n)
x = np.sin(2*np.pi*t)
# 離散フーリエ変換を実施します
y = fft(x)
# パワースペクトルを計算します
ps = np.abs(y / len(y))
# 周波数領域の横軸を用意します(データ数と刻み幅から計算)
freq = fftfreq(data_n, d=timestep)
# データの並び方を 0 Hz 中心になるように移動します
freq = fftshift(freq)
ps = fftshift(ps)
# パワースペクトルのグラフを作成します
plt.plot(freq, ps)
plt.show()