← コード一覧に戻る
輸送方程式_中心差分2次精度近似_2.py
拡張版：dt 3条件・物理時刻で比較
概要
前バージョン（
輸送方程式_中心差分2次精度近似.py）の拡張版。
中心差分（空間2次精度）＋前進オイラー（時間1次精度）のスキームで、
時間刻み Δt を 3 種類（0.1 / 0.05 / 0.01）変えた場合の挙動を
横並び 3 パネルで比較する。
プロットタイミングをステップ番号ではなく
物理時刻（t = 0.1 / 0.2 / 0.3 s）で指定することで、
異なる Δt 同士を同一の物理時刻で公平に比較できる。
- 初期条件：ステップ関数（左半分 q=1、右半分 q=0）
- 格子点数：21点（x = 0 〜 2.0 m）
- 総ステップ数：600（各 dt で同数のループを回す）
- 出力時刻：t = 0.1 / 0.2 / 0.3 s のみ
使用ライブラリ
| numpy | 格子配列・ゼロ配列の生成、copy() による旧値保持、np.isclose による浮動小数点比較 |
|---|
| matplotlib | subplots(1, 3) で 3 パネル横並び表示。各パネルに物理時刻ラベルでプロット |
|---|
物理・数学的背景
支配方程式
∂q/∂t + c ∂q/∂x = 0
解析解は初期分布 q₀ が速度 c で右方向に平行移動するだけ：
q(x, t) = q₀(x − ct)
パラメータ
| c = 1 | 移流速度 [m/s]（全条件共通） |
|---|
| dx = 0.1 | 空間刻み [m]（全条件共通） |
|---|
| jmax = 21 | 格子点数（x = 0 〜 2.0 m） |
|---|
| nmax = 600 | 総時間ステップ数 |
|---|
| dt = 0.1 / 0.05 / 0.01 | 時間刻み（3 パネルで変化させる変数） |
|---|
| plot_times = [0.1, 0.2, 0.3] | プロットする物理時刻 [s] |
|---|
クーラン数とステップ数の比較
| dt | r = c·Δt / Δx | t=0.1s に対応するステップ n |
|---|
| 0.1 | 1.00 | n = 1 |
|---|
| 0.05 | 0.50 | n = 2 |
|---|
| 0.01 | 0.10 | n = 10 |
|---|
計算手法
中心差分（空間2次精度）
∂q/∂x ≈ (q[j+1] − q[j−1]) / (2Δx)
時間更新式（前進オイラー法・陽解法）
q[j]^(n+1) = q[j]^n − c·Δt · (q[j+1]^n − q[j−1]^n) / (2Δx)
物理時刻によるプロット条件
従来の「偶数ステップ（n%2==0）」条件に代わり、物理時刻リストで出力タイミングを管理する。
if any(np.isclose(n*dt, pt) for pt in plot_times)
np.isclose を使う理由：浮動小数点演算では n*dt が厳密に一致しない場合がある
（例：0.1 + 0.1 + 0.1 = 0.30000000000000004）。np.isclose は相対・絶対許容誤差を
考慮して比較するため、意図した時刻を確実に検出できる。
ラベルの有効数字表示
label=f't={n*dt:.2g}'
:.2g は有効数字2桁の汎用フォーマット。末尾のゼロを自動省略するため、
0.10 → 0.1、0.20 → 0.2 のように簡潔に表示される。
コードの流れ
1
パラメータ・格子設定
c, dx, jmax, nmax を定義。np.linspace で x 座標配列を生成
2
初期条件の生成（q_init）
左半分 q=1、右半分 q=0 のステップ関数。各 dt ループで q_init.copy() を使って再利用
3
プロット時刻リストの定義
plot_times = [0.1, 0.2, 0.3]。ここを変更するだけで出力タイミングを自由に指定できる
4
3パネルの figure 作成
plt.subplots(1, 3, figsize=(21, 7))。axes と dt リストを zip で対応付け
5
dt ループ（3回）
各 dt について q = q_init.copy() で初期化し、初期分布（t=0）をプロット
6
時間ループ（n = 1 〜 600）
中心差分で内部格子点を更新。np.isclose で plot_times と一致するステップのみ描画
7
各パネルの装飾
グリッド・軸ラベル・y軸範囲（0〜2.0）・タイトル（dt値）・凡例を設定。plt.tight_layout() で整列
境界の扱い
空間ループは range(1, jmax-1) なので両端（j=0, j=jmax-1）は固定境界として更新されない。
精度と安定性
精度
| 空間差分 | 中心差分 → 2次精度 O(Δx²) |
|---|
| 時間差分 | 前進オイラー（陽解法） → 1次精度 O(Δt) |
|---|
安定性（von Neumann 解析）
中心差分＋前進オイラーの組み合わせは無条件不安定。
クーラン数 r の値によらず数値誤差が増幅し続ける。
本コードでは Δt を変えて同一物理時刻（t=0.1〜0.3 s）での状態を比較することで、
Δt が大きいほど（r が大きいほど）振動の成長が速いことを示す。
3条件の比較まとめ
| dt | クーラン数 r | t=0.3s 時点の挙動 |
|---|
| 0.1 | 1.00 | 振動が最も大きく成長（わずか3ステップで到達） |
|---|
| 0.05 | 0.50 | 中程度の振動（6ステップで到達） |
|---|
| 0.01 | 0.10 | 振動は小さい（30ステップで到達） |
|---|
安定な手法との比較
| 手法 | 空間精度 | 安定条件 |
|---|
| 風上差分 | 1次 | r ≤ 1（条件付き安定） |
|---|
| 中心差分＋前進オイラー | 2次 | 無条件不安定 |
|---|
| Lax-Wendroff法 | 2次 | r ≤ 1（条件付き安定） |
|---|
コード全文
import numpy as np
import matplotlib.pyplot as plt
c = 1 # 移流速度
dx = 0.1 # 空間刻み
jmax = 21 # 格子点数
nmax = 600 # 総時間ステップ数
x = np.linspace(0, dx*(jmax-1), jmax)
# ── 初期条件（ステップ関数） ──────────────────────
q_init = np.zeros(jmax)
for j in range(jmax):
if(j < jmax/2):
q_init[j] = 1
else:
q_init[j] = 0
fig, axes = plt.subplots(1, 3, figsize=(21, 7), dpi=100)
plt.rcParams["font.size"] = 25
plot_times = [0.1, 0.2, 0.3] # プロットする物理時刻 [s]
# ── dt を 3 種類に変えて比較 ──────────────────────
for ax, dt in zip(axes, [0.1, 0.05, 0.01]):
q = q_init.copy()
ax.plot(x, q, marker='o', lw=2, label='t=0.0')
for n in range(1, 1+nmax):
qold = q.copy()
for j in range(1, jmax-1):
q[j] = qold[j] - dt*c*(qold[j+1]-qold[j-1])/(2*dx)
if any(np.isclose(n*dt, pt) for pt in plot_times):
ax.plot(x, q, marker='o', lw=2, label=f't={n*dt:.2g}')
ax.grid(color='black', linestyle='dashed', linewidth=0.5)
ax.set_xlabel('x')
ax.set_ylabel('q')
ax.set_ylim([0, 2.0])
ax.set_title(f'dt = {dt}')
ax.legend()
plt.tight_layout()
plt.show()
グラフ結果
- 3パネルともに t = 0.1 / 0.2 / 0.3 s の同一物理時刻でプロット → 公平な比較が可能
- 左パネル（dt=0.1, r=1.0）：わずか1〜3ステップで振動が大きく成長
- 中央パネル（dt=0.05, r=0.5）：2〜6ステップ、中程度の振動
- 右パネル（dt=0.01, r=0.1）：10〜30ステップ、振動は相対的に小さい
- いずれの Δt でも振動は存在 → 中心差分＋前進オイラーの無条件不安定性を確認
関連キーワード
移流方程式
中心差分
2次精度
前進オイラー法
数値不安定
クーラン数
np.isclose
物理時刻指定
有効数字
subplots
偏微分方程式
numpy
matplotlib