続。きれいだな。花火を見るゲーム。

前回のあらすじ。

楕円形の確率分布を持つガウシアン爆発の2次元シミュレーションを行い、5つのタイムステップでその結果をプロットを以下に示します。NumPyを使用して乱数を生成し、Matplotlibでプロットを行います。

パーティクルの飛び散る図。

圧力分布を楕円形のガウス分布に基づいてプロットを以下に示します。圧力は、パーティクルの密度を基に計算され、各タイムステップごとにプロットされます。

パーティクルの飛び散る図のコード。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
num_particles = 1000
num_steps = 5

# 楕円形のガウス分布のパラメータ
sigma_x = 1.0  # 横方向の広がり
sigma_y = 0.5  # 縦方向の広がり

# タイムステップごとにパーティクルの位置を記録するリスト
positions = []

# パーティクルの初期位置を原点に設定
x_initial = np.zeros(num_particles)
y_initial = np.zeros(num_particles)

# パーティクルの速度をガウス分布に基づいて生成
vx = np.random.normal(0, sigma_x, num_particles)
vy = np.random.normal(0, sigma_y, num_particles)

# シミュレーションの実行
for step in range(num_steps):
    # パーティクルの位置を更新
    x = x_initial + vx * step
    y = y_initial + vy * step
    positions.append((x, y))

# プロットの設定
fig, axes = plt.subplots(1, num_steps, figsize=(15, 3))

for i, ax in enumerate(axes):
    x, y = positions[i]
    ax.scatter(x, y, s=1)
    ax.set_xlim(-5, 5)
    ax.set_ylim(-5, 5)
    ax.set_title(f"Step {i+1}")
    ax.set_aspect('equal')

plt.tight_layout()
plt.show()

コードの説明
num_particles: シミュレーションで使用するパーティクルの数です。ここでは1000個のパーティクルを使用しています。
sigma_x, sigma_y: 楕円形のガウス分布の横方向と縦方向の広がりを決めるパラメータです。これにより、横方向には広がりが大きく、縦方向には狭い爆発をシミュレートしています。
positions: 各タイムステップにおけるパーティクルの位置を記録するリストです。
vx, vy: 横方向と縦方向の速度成分を生成しています。これらはガウス分布に基づいています。
プロット: 各タイムステップでのパーティクルの位置をプロットします。5つのサブプロットを横に並べて表示します。
実行結果
このコードを実行すると、5つのタイムステップにおける楕円形のガウス分布に基づく爆発の様子がプロットされます。それぞれのタイムステップで、パーティクルが横方向に広がりやすく、縦方向にはそれほど広がらないことが視覚的に確認できます。

圧力分布を楕円形のガウス分布に基づいてプロットのコード。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# パラメータ設定
num_particles = 1000
num_steps = 5
grid_size = 100
sigma_x = 1.0  # 横方向の広がり
sigma_y = 0.5  # 縦方向の広がり

# グリッドを作成
x_grid = np.linspace(-5, 5, grid_size)
y_grid = np.linspace(-5, 5, grid_size)
x, y = np.meshgrid(x_grid, y_grid)
pos = np.dstack((x, y))

# パーティクルの初期位置を原点に設定
x_initial = np.zeros(num_particles)
y_initial = np.zeros(num_particles)

# パーティクルの速度をガウス分布に基づいて生成
vx = np.random.normal(0, sigma_x, num_particles)
vy = np.random.normal(0, sigma_y, num_particles)

# 圧力分布を保存するリスト
pressure_distributions = []

# シミュレーションの実行
for step in range(num_steps):
    # パーティクルの位置を更新
    x_positions = x_initial + vx * step
    y_positions = y_initial + vy * step

    # 各パーティクルの位置に基づいて圧力分布を計算
    pressure = np.zeros((grid_size, grid_size))
    for i in range(num_particles):
        # 各パーティクルがグリッドに与える影響をガウスカーネルで計算
        particle_pos = np.array([x_positions[i], y_positions[i]])
        rv = multivariate_normal(mean=particle_pos, cov=[[sigma_x**2, 0], [0, sigma_y**2]])
        pressure += rv.pdf(pos)
    
    pressure_distributions.append(pressure)

# プロットの設定
fig, axes = plt.subplots(1, num_steps, figsize=(15, 3))

for i, ax in enumerate(axes):
    ax.imshow(pressure_distributions[i], extent=(-5, 5, -5, 5), origin='lower', cmap='hot')
    ax.set_title(f"Step {i+1}")
    ax.set_aspect('equal')

plt.tight_layout()
plt.show()

コードの説明
grid_size: 圧力分布を計算するグリッドのサイズを設定します。ここでは100x100のグリッドを使用しています。
multivariate_normal: scipy.statsからの多変量正規分布を使用して、各パーティクルの影響をグリッドに広げます。楕円形の分布に基づいて圧力分布を計算します。
pressure_distributions: 各タイムステップにおける圧力分布を記録するリストです。
プロット: 各タイムステップごとに計算された圧力分布をヒートマップとしてプロットします。5つのサブプロットを横に並べて表示します。
実行結果
このコードを実行すると、各タイムステップにおける圧力分布がプロットされます。楕円形のガウス分布に基づくため、横方向に圧力が広がりやすく、縦方向にはそれほど広がらない分布が得られます。ヒートマップで圧力の強さが視覚化され、時間経過とともに圧力分布が変化する様子が確認できます。

参考。

シンプルにグリッド内にいくつのパーテクルが存在するかをカウントして圧力と見なしている簡易的な計算モデルです。