輸送方程式_風上差分1次精度近似_2.py

概要

1次元移流方程式を風上差分法（1次精度）で解き、 CFL数を 0.5 に固定したままΔx・Δtを同時に細かくする格子収束の確認を行うコード。前作（_1）がCFL数の違いによる安定性比較だったのに対し、本コードは「同じCFL条件のまま格子を細かくすると数値拡散がどれだけ減少するか」を可視化する。

使用ライブラリ

ライブラリ	用途
numpy	配列演算・初期条件の構築
matplotlib.pyplot	グラフ描画（subplots による横並び比較）

物理背景

対象方程式は1次元線形移流方程式：

∂q/∂t + c · ∂q/∂x = 0　（c = 1 : 移流速度）

解析解は初期ステップ波形が形を変えずに x 方向へ速度 c で平行移動するもの。風上差分1次精度スキームは打ち切り誤差からO(Δx)の数値拡散項を含む。 CFL数を固定したままΔxを小さくすると数値拡散係数 ∝ Δx も小さくなり、格子を細かくするほど解析解に近づく（1次収束）。

数値拡散と格子収束の関係

風上差分法をテイラー展開すると等価な方程式は：

∂q/∂t + c·∂q/∂x = (c·Δx/2)(1 − CFL)·∂²q/∂x²　＋ O(Δx²)

右辺が数値拡散項。CFL を固定した場合、Δx を半分にすると数値拡散係数も半分になる（1次収束）。

計算手法

項目	内容
空間差分	1次精度風上差分（後退差分）：(q[j] − q[j-1]) / Δx
時間積分	前進オイラー法（陽的1次精度）
初期条件	ステップ関数（x < 領域中央で q=1、以降 q=0）
境界条件	両端固定（ループ外で更新なし）
観測時刻	t = 0.1 / 0.2 / 0.3

比較ケース（全てCFL=0.5）

ケース	Δt	Δx	格子点数 jmax	CFL = cΔt/Δx
粗い格子	0.1	0.2	11	0.5
中間格子	0.025	0.05	41	0.5
細かい格子	0.01	0.02	101	0.5

差分式

q[j]ⁿ⁺¹ = q[j]ⁿ − c · Δt/Δx · (q[j]ⁿ − q[j-1]ⁿ)

コードの流れ

グローバルパラメータ c=1 を設定
subplots で3パネルを準備し、zip で (dt, dx, jmax, nmax) の3ケースを一括ループ
各ケースで座標配列 x と初期ステップ関数 q_init を生成（格子数が異なる）
nmax=600 ステップ分時間発展し、t=0.1 / 0.2 / 0.3 時点をプロット
タイトルにΔtとCFL数を表示し、グリッド・凡例を付与して出力

精度・安定性

CFL=0.5 固定のためすべてのケースで安定
格子を細かくするほど（左→右）ステップ波形のなまりが減少し、解析解に近づく
1次精度スキームのため収束は遅く、格子を 10 倍細かくすると誤差が約 1/10 になる
前作（_1）との比較から「CFL 固定での精度向上には格子細分化が必要」という理解が得られる
nmax=600 は固定だが Δt が異なるため物理時間の進み方が異なる点に注意（粗い格子：0〜60、細かい格子：0〜6）

コード全文

import numpy as np
import matplotlib.pyplot as plt

c=1
#dt=0.05

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

for ax, dt, dx, jmax, nmax in zip(axes, [0.1, 0.025, 0.01], [0.2, 0.05, 0.02], [11, 41, 101], [600, 600, 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

    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):
            q[j]=qold[j]-dt*c*(qold[j]-qold[j-1])/(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()

グラフ結果

左：粗い格子（Δx=0.2, jmax=11）　／　中：中間格子（Δx=0.05, jmax=41）　／　右：細かい格子（Δx=0.02, jmax=101）
いずれもCFL=0.5。格子を細かくするにつれてステップ波形のなまりが減少し解析解に収束していく様子が確認できる。

キーワード

移流方程式風上差分法 Upwind Scheme 格子収束 CFL数固定数値拡散数値粘性 1次収束前進オイラー法偏微分方程式