コード解説：test1.py

初速度 v₀ = 10 m/s で鉛直上向きに投げ上げた物体の高さを、解析解（青線）と前進オイラー法による数値解（黒点）の2通りで求めてグラフ表示するプログラム。数値解法の基礎を学ぶための入門的な実装。

使用ライブラリ

物理・数学的背景

numpy	数値計算用ライブラリ。`np.linspace` で時刻配列を生成し、解析解を配列演算で一括計算する
matplotlib	グラフ描画ライブラリ。解析解を折れ線、数値解を散布点でプロットし、凡例・グリッドを付与する

支配方程式

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

d²h/dt² = −g

これを2回積分すると解析解が得られる：

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

v(t) = −g t + v₀

ここで h₀ は初期高さ（= 0）、v₀ は初速度（= 10 m/s）。

最高点と飛行時間

v(t) = 0 となる時刻が最高点 → t_peak = v₀/g ≈ 1.02 s
h に戻るのは t_land = 2v₀/g ≈ 2.04 s（コードの計算範囲の上限）

計算手法

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

微分方程式 dh/dt = f(t, h) を時間ステップ Δt で近似的に解く最もシンプルな数値解法。

h(t + Δt) ≈ h(t) + f(t, h) · Δt

本コードでは速度 v(t) = −gt + v₀ を使い、各ステップで高さを更新する：

h_new = h + (−g·t + v₀) · Δt

特徴とトレードオフ

実装が簡単で計算コストが低い
時間ステップ Δt が大きいほど誤差が蓄積する（1次精度）
本コードは Δt = 0.5 s と粗く、誤差が顕著に現れる
精度向上には Δt を小さくするか、ルンゲ＝クッタ法（4次精度）などを使う

コードの流れ

パラメータ設定
重力加速度 g=9.8、初速度 v₀=10、初期高さ h₀=0、時間刻み Δt=0.5 を定義

解析解の計算・プロット
np.linspace(0, 2v₀/g, 100) で細かい時刻配列を生成し、h(t) = −½gt² + v₀t を一括計算して青線で描画

数値解のループ計算
t=0, h=h₀ から出発し、while h >= h₀ の間ループ。毎ステップ現在の (t, h) を黒点でプロットし、オイラー法で h と t を更新

グラフ装飾・表示
グリッド・軸ラベル・凡例（Analytical / Numerical）を付けて plt.show() で表示

ループ終了条件の注意点

while h >= h₀ は「高さが初期値（地面）以上の間」継続する条件。Δt が大きいと h が負になった時点でループが止まるため、着地点が実際より手前または奥にずれる。

解析解との比較・精度

Δt による誤差の影響

前進オイラー法の局所打ち切り誤差は O(Δt²)、全体誤差は O(Δt)（1次精度）。本コードの Δt = 0.5 s は飛行時間（≈2 s）の25%にあたるため、誤差が大きく視覚的にはっきり確認できる。

最高点付近で数値解は解析解より低くなる傾向（速度の過小評価）
着地タイミングがずれ、数値解のループが早めに終了する

改善の方向

Δt を小さくする（例：0.01 s）→ 計算点が増えるが誤差は大幅改善
ルンゲ＝クッタ法（RK4）に切り替え → 同じ Δt でも4次精度で高精度
後退オイラー法（陰解法）→ 安定性が高くなる

コード全文

import numpy as np import matplotlib.pyplot as plt g = 9.8 # 重力加速度 [m/s²] v0 = 10 # 初速度 [m/s] h0 = 0 # 初期高さ [m] dt = 0.5 # 時間刻み [s] plt.figure(figsize=(7, 7), dpi=100) plt.rcParams["font.size"] = 25 # ── 解析解 ────────────────────────────── t = np.linspace(0, 2*v0/g, 100) h = -0.5*g*t**2 + v0*t + h0 la, = plt.plot(t, h, color='blue') # ── 数値解（前進オイラー法） ───────────── t = 0 h = h0 while h >= h0: ln = plt.scatter(t, h, marker='o', c='black') h += (-g*t + v0) * dt # オイラー更新 t += dt # ── グラフ装飾 ─────────────────────────── plt.grid(color='black', linestyle='dashed', linewidth=0.5) plt.xlabel('Time') plt.ylabel('Height') plt.legend(handles=[la, ln], labels=['Analytical', 'Numerical']) plt.show()

グラフ結果

青線（Analytical）：解析解 h(t) = −½gt² + v₀t。滑らかな放物線
黒点（Numerical）：前進オイラー法（Δt = 0.5 s）による数値解。点の間隔が粗く、解析解からのずれが視覚的に確認できる
最高点付近で数値解が解析解より低くなり、着地タイミングも早めにずれている（Δt = 0.5 s の誤差蓄積の影響）

test1.py ― 放物運動の解析解と数値解の比較

概要