コード解説：輸送方程式_中心差分2次精度近似.py

1次元移流方程式（輸送方程式）を中心差分（空間2次精度）＋前進オイラー（時間1次精度）で数値的に解くプログラム。
初期条件としてステップ関数（左半分が1、右半分が0）を与え、時間発展とともに分布がどう変化するかをグラフで確認する。

使用ライブラリ

物理・数学的背景

numpy	格子配列・ゼロ配列の生成、配列演算
matplotlib	各時間ステップの分布をプロット。偶数ステップのみ描画

支配方程式

∂q/∂t + c ∂q/∂x = 0

解析解は初期分布 q₀ が速度 c で右方向に平行移動するだけ：

q(x, t) = q₀(x − ct)

パラメータ

c = 1	移流速度 [m/s]
dt = 0.05	時間刻み [s]
dx = 0.1	空間刻み [m]
jmax = 21	格子点数（x = 0 〜 2.0 m）
nmax = 6	総時間ステップ数

計算手法

中心差分（空間2次精度）

空間微分を左右対称の差分で近似する：

∂q/∂x ≈ (q[j+1] − q[j−1]) / (2Δx)

前後2点を使うため、片側差分（1次精度）より精度が高い（2次精度）。

時間更新式（前進オイラー法・陽解法）

前進オイラー法は陽解法（explicit method）であり、次の時刻の値 q^(n+1) を現在時刻 n の値だけで直接計算できる。連立方程式を解く必要がなく実装が簡単だが、安定条件（CFL条件）を満たす必要がある。

q[j]^(n+1) = q[j]^n − c·Δt · (q[j+1]^n − q[j−1]^n) / (2Δx)

クーラン数

r = c·Δt / Δx = 1 × 0.05 / 0.1 = 0.5

コードの流れ

パラメータ・格子設定
c, dt, dx, jmax, nmax を定義。np.linspace で x 座標配列を生成

初期条件の設定
j < jmax/2（左半分）は q=1、それ以外は q=0 のステップ関数

初期分布をプロット（n=0）

時間ループ（n = 1 〜 nmax）
qold = q.copy() で更新前の値を保存し、中心差分で全格子点を更新。
偶数ステップ（n=2, 4, 6）のみグラフに追加

グラフ表示
各時刻の分布を重ねてプロット。凡例は n=0, 2, 4, 6

境界の扱い

空間ループは range(1, jmax-1) なので、両端 j=0 と j=jmax-1 は更新されない（固定境界）。

qold = q.copy() の役割

更新中に q[j] を書き換えると、隣の計算に新しい値が混入してしまう。copy() で旧値を保持することで、全格子点を同一時刻の値から計算できる。

精度と安定性

精度

空間差分	中心差分 → 2次精度 O(Δx²)
時間差分	前進オイラー（陽解法） → 1次精度 O(Δt)

安定性

中心差分＋前進オイラーの組み合わせは無条件不安定（von Neumann 安定性解析より）。
クーラン数 r がいかなる値でも数値振動が成長し、時間が経つにつれて解が発散する。
グラフではステップの不連続部付近に振動（数値拡散ではなく数値分散）が現れる。

安定な手法との比較

手法	空間精度	安定条件
風上差分	1次	r ≤ 1（安定）
中心差分＋前進オイラー	2次	無条件不安定
Lax-Wendroff法	2次	r ≤ 1（安定）

コード全文

import numpy as np import matplotlib.pyplot as plt c = 1 # 移流速度 dt = 0.05 # 時間刻み dx = 0.1 # 空間刻み jmax = 21 # 格子点数 nmax = 6 # 総時間ステップ数 x = np.linspace(0, dx*(jmax-1), jmax) # ── 初期条件（ステップ関数） ──────────────── q = np.zeros(jmax) for j in range(jmax): if(j < jmax/2): q[j] = 1 else: q[j] = 0 plt.figure(figsize=(7, 7), dpi=100) plt.rcParams["font.size"] = 25 plt.plot(x, q, marker='o', lw=2, label='n=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 n%2==0: plt.plot(x, q, marker='o', lw=2, label=f'n={n}') # ── グラフ装飾 ─────────────────────────────── plt.grid(color='black', linestyle='dashed', linewidth=0.5) plt.xlabel('x') plt.ylabel('q') plt.legend() plt.show()

グラフ結果

n=0：初期ステップ関数（左半分=1、右半分=0）
n=2, 4, 6：時間発展後の分布
不連続部付近に数値振動（リップル）が現れる → 中心差分＋前進オイラーの不安定性
理想的には分布が右に平行移動するだけだが、振動が成長していく様子が確認できる

輸送方程式_中心差分2次精度近似.py ― 1次元移流方程式の中心差分法（2次精度）

概要