プログラム一覧に戻る

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt("kadai1-1.txt")  # データファイルの読み込み
N = (data[0, :].size - 1) // 4  # 粒子数
v = data[:, 2 * N + 1 : 4 * N + 1]  # 各時刻での各粒子の速度

# どの時間帯のデータを使うか
v_plot = v[0, :].flatten()  # 初期時刻での各粒子の速度
# v_plot = v[1:11, :].flatten()  # はじめの方の時間帯の各粒子の速度
# v_plot = v[101:, :].flatten()  # 終わりの方の時間帯の各粒子の速度： 101のところの数字は自分で調整する

# 速度のヒストグラム（規格化したもの）をプロット
plt.hist(v_plot, bins=100, range=(-2, 2), density=True)

T = np.sum(v_plot ** 2) / v_plot.size  # 使用した時間帯の温度を等分配則から計算
v_compare = np.linspace(-2, 2, 100)  # 比較のプロットのための横軸データ
# マクスウェル分布（比較のプロットのための縦軸データ）
f_Maxwell = np.exp(-(v_compare ** 2) / (2 * T)) / np.sqrt(2 * np.pi * T)
plt.plot(v_compare, f_Maxwell, linestyle="dashed", color="black")  # 比較のためのマクスウェル分布を黒破線でプロット

plt.xlabel("velocity")  # x軸の名前を設定
plt.ylabel("distribution")  # y軸の名前を設定
plt.xlim(-2, 2)  # 横軸の範囲を[-2, 2]に設定

plt.show()  # グラフを表示