test2.py ― 時間刻みの違いによる数値解の比較（Δt = 0.5 vs 0.1）

概要
使用ライブラリ
物理・数学的背景
計算手法
コードの流れ
解析解との比較・精度
コード全文
グラフ結果

概要

test1.py（単一 Δt での比較）を拡張し、Δt = 0.5 s と Δt = 0.1 s の2通りの時間刻みで前進オイラー法を実行し、横並び2グラフで解析解との比較を行うプログラム。
時間刻みが数値解の精度に与える影響を視覚的に確認できる。

使用ライブラリ

numpy	`np.linspace` で解析解用の時刻配列を生成し、解析解を配列演算で一括計算する
matplotlib	`plt.subplots(1, 2)` で2つのグラフ領域を生成。解析解を折れ線、数値解を散布点でプロット

物理・数学的背景

支配方程式

重力加速度 g（= 9.8 m/s²）のみが働く1次元鉛直運動：

d²h/dt² = −g

解析解：

h(t) = −½ g t² + v₀ t + h₀

v(t) = −g t + v₀

初速度 v₀ = 10 m/s、初期高さ h₀ = 0。飛行時間 t_land = 2v₀/g ≈ 2.04 s。

計算手法

前進オイラー法（Explicit Euler Method）

速度 v(t) = −gt + v₀ を用いて高さを逐次更新する：

h(t + Δt) ≈ h(t) + (−g·t + v₀) · Δt

2種類の Δt による比較

Δt	飛行時間に対する割合	計算点数（目安）	精度
0.5 s	約 25%	4〜5 点	誤差が大きく視覚的にずれる
0.1 s	約 5%	20〜21 点	誤差が小さく解析解に近い

コードの流れ

パラメータ設定
g, v₀, h₀ を定義。Δt はループ内で変化させるため、ここでは定義しない

2グラフ領域の生成
plt.subplots(1, 2) で横並び2枚の Axes を取得。axes[0] が左（Δt=0.5）、axes[1] が右（Δt=0.1）

解析解の事前計算
Δt に依存しないため、ループの外で1回だけ計算

ループで各 Δt のグラフを描画
for ax, dt in zip(axes, [0.5, 0.1]) で各グラフ領域と Δt を対応させながら解析解・数値解をプロット

グラフ装飾・表示
各グラフにグリッド・軸ラベル・タイトル（dt = ○○）・凡例を付与。plt.tight_layout() で重なりを防ぐ

zip を使ったペアリング

zip(axes, [0.5, 0.1]) は2つのリストを1対1で対応させる。
1回目：ax = axes[0], dt = 0.5 → 2回目：ax = axes[1], dt = 0.1

解析解との比較・精度

前進オイラー法の全体誤差は O(Δt)（1次精度）。Δt を1/5にすると誤差もおよそ1/5になる。

Δt = 0.5 s：最高点付近で解析解より高さが低くなり、着地も早めにずれる
Δt = 0.1 s：点が密になり解析解の放物線に沿って分布する。誤差が目視では小さくなる

このコードで学べること：
同じアルゴリズムでも Δt の違いが結果に大きく影響することを、グラフの並列比較によって直感的に理解できる。

コード全文

import numpy as np
import matplotlib.pyplot as plt

g  = 9.8   # 重力加速度 [m/s²]
v0 = 10    # 初速度 [m/s]
h0 = 0     # 初期高さ [m]

plt.rcParams["font.size"] = 25
fig, axes = plt.subplots(1, 2, figsize=(14, 7), dpi=100)

# ── 解析解（Δt に依存しないので外で計算） ───────
t_ana = np.linspace(0, 2*v0/g, 100)
h_ana = -0.5*g*t_ana**2 + v0*t_ana + h0

# ── 各 Δt でグラフを描画 ────────────────────────
for ax, dt in zip(axes, [0.5, 0.1]):
    la, = ax.plot(t_ana, h_ana, color='blue')

    t = 0
    h = h0
    while h >= h0:
        ln = ax.scatter(t, h, marker='o', c='black')
        h += (-g*t + v0) * dt   # 前進オイラー更新
        t += dt

    ax.grid(color='black', linestyle='dashed', linewidth=0.5)
    ax.set_xlabel('Time')
    ax.set_ylabel('Height')
    ax.set_xlim([0, 2.5])
    ax.set_ylim([0, 10])
    ax.set_title(f'dt = {dt}')
    ax.legend(handles=[la, ln], labels=['Analytical', 'Numerical'])

plt.tight_layout()
plt.show()

グラフ結果

青線（Analytical）：解析解 h(t) = −½gt² + v₀t。両グラフ共通
黒点（Numerical, Δt = 0.5）：点が粗く、解析解から明らかにずれている
黒点（Numerical, Δt = 0.1）：点が密になり、解析解の放物線に沿っている
Δt を小さくすることで数値解の精度が向上することが視覚的に確認できる