グラフ結果
CFL = ±0.5 での計算結果。
- 左図（c = +1）：波が右方向に移流。初期ステップ（x=1 で q が 0→1 に立ち上がり）が右へ移動しながら数値拡散で広がる。
- 右図（c = −1）：波が左方向に移流。同じ初期条件が左へ移動。フラックス分割により c の符号が反転しても安定かつ対称的な拡散挙動が得られる。
- 両ケースとも数値拡散が生じるが、振動・発散はなく CFL ≤ 1 での安定性を確認できる。
1次元移流方程式に対し、移流速度の符号に応じて差分方向を自動的に切り替える「風上下差分（フラックス分割）」を実装し、c = +1（右向き移流）と c = −1（左向き移流）の挙動を横並びで比較するコード。
| 項目 | 内容 |
|---|---|
| 対象方程式 | 1次元移流方程式 ∂q/∂t + c ∂q/∂x = 0 |
| 空間刻み | dx = 0.10（格子点数 jmax = 21） |
| 時間刻み | dt = 0.05（ステップ数 nmax = 600） |
| CFL数 | |ν| = |c|·dt/dx = 0.5 |
| 初期条件 | ステップ関数（x < 1 で q=0、x ≥ 1 で q=1） |
| 出力時刻 | t = 0.1, 0.2, 0.3 |
| 比較ケース | 左図：c = +1（右向き移流） / 右図：c = −1（左向き移流） |
| ライブラリ | 用途 |
|---|---|
numpy | 配列演算・初期条件の設定 |
matplotlib.pyplot | 1×2 サブプロットへの描画 |
1次元移流方程式の解析解は波形が速度 c で平行移動するものとなる：
c > 0 のとき波は右へ、c < 0 のとき波は左へ進む。 差分近似では、波が来る方向（風上側）の情報を使うことが安定性の鍵となる。 符号が固定されていない場合、c の符号に応じて差分方向を自動選択する必要がある。
フラックス分割（風上下差分）の考え方：
移流速度 c を正成分と負成分に分割する：
それぞれに対し、情報が来る方向（風上側）の差分を適用することで、どちらの向きの流れでも安定なスキームを構成できる。
これを展開すると：
精度は時間・空間ともに1次精度。CFL条件 |ν| ≤ 1 のもとで条件付き安定。 ステップ波形は数値拡散により緩やかになる（風上差分の特性）。
for ax, c in zip(axes, [1,-1]) ループで c=+1 と c=−1 の2ケースを処理q_init を作成（x < 1 で q=0、x ≥ 1 で q=1）c2・c3 をその場で計算し、フラックス分割型差分式で q を更新plot_times に一致するステップで描画。plt.tight_layout() でレイアウト整形して出力境界条件：j = 1 ～ jmax-2 のみ更新（両端 j=0, j=jmax-1 は固定値のまま）。
| 項目 | 内容 |
|---|---|
| 空間精度 | 1次精度（後退または前進差分） |
| 時間精度 | 1次精度（前進オイラー法） |
| 安定条件 | |ν| = |c|·dt/dx ≤ 1 （本コードは |ν| = 0.5 で安定） |
| 主な誤差特性 | 数値拡散（差分打ち切り誤差の主成分が拡散項）。ステップ波形が時間とともになだらかになる |
| 改善方向 | Lax-Wendroff 法（2次精度）や MUSCL スキームで数値拡散を低減できる |
c = +1 と c = −1 の両ケースで対称的な拡散挙動が得られる点がフラックス分割の特徴。
import numpy as np
import matplotlib.pyplot as plt
c=1
#dt=0.05
fig, axes = plt.subplots(1, 2, figsize=(14, 7), dpi=100)
plt.rcParams["font.size"]=25
dt=0.05
dx=0.10
jmax=21
nmax=600
for ax, c in zip(axes, [1,-1]):
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]=0
else:
q_init[j]=1
plot_times = [0.1, 0.2, 0.3]
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):
c2=(c+abs(c))/2
c3=(c-abs(c))/2
q[j]=qold[j]-dt*c2*(qold[j]-qold[j-1])/(dx)-dt*c3*(qold[j+1]-qold[j])/(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,1.5])
ax.set_title(f'dt = {dt},Courant={dt*c/dx:.2g}')
ax.legend()
plt.tight_layout()
plt.show()CFL = ±0.5 での計算結果。