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

概要

1次元移流方程式（線形移流方程式）を風上差分法（Upwind scheme）で解くコード。 空間差分は1次精度の風上差分、時間差分は前進オイラー法（陽的1次精度）を用いる。 Δt を 0.05 / 0.1 / 0.2 の3種類に変えてクーラン数（Courant数）が解の安定性に与える影響を横並びグラフで比較する。

使用ライブラリ

物理背景

対象方程式は1次元線形移流方程式（スカラー輸送方程式）：

この方程式は流体中のスカラー量（温度・濃度・放射性核種等）が速度 c で上流から下流へ流れる現象を表す。 解析解は初期波形が形を変えずに x 方向へ速度 c で平行移動するものとなる。

クーラン数（Courant–Friedrichs–Lewy 数）

風上差分法の安定条件は CFL 条件で与えられる：

• CFL = 1 のとき：数値拡散ゼロの完全移流（数値的に最も正確）
• 0 < CFL < 1 のとき：安定だが数値拡散（数値粘性）あり
• CFL > 1 のとき：不安定（解が発散）

計算手法

差分式

コードの流れ

1. パラメータ設定：c=1、Δx=0.1、jmax=21、nmax=600
2. 座標配列 x の生成（0 〜 2.0）
3. ステップ関数の初期条件 q_init を設定
4. subplots でΔt=0.05 / 0.1 / 0.2 の3パネルを並べて描画
5. 各パネルで nmax ステップ分時間発展し、t=0.2 / 0.4 / 0.6 時点をプロット
6. タイトルにΔtとCFL数を表示し、グリッド・凡例を付与して出力

精度・安定性

• 空間・時間ともに1次精度のため、数値拡散（数値粘性）が大きく、ステップ波形は時間とともになまる
• 中心差分（2次精度）と異なり 無条件不安定にならない（CFL ≤ 1 で安定）
• CFL = 1.0（Δt=0.1）では数値拡散が最小となり、解析解に最も近い結果が得られる
• CFL = 2.0（Δt=0.2）では安定条件を超えて解が発散する様子が確認できる
• 中心差分との比較から、風上差分の「安定性と精度のトレードオフ」が直感的に理解できる

コード全文

import numpy as np
import matplotlib.pyplot as plt

c=1
#dt=0.05
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.2, 0.4, 0.6]

for ax, dt in zip(axes, [0.05, 0.1, 0.2]):
    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()

グラフ結果

左：CFL=0.5（Δt=0.05）安定・数値拡散あり　／　中：CFL=1.0（Δt=0.1）安定・数値拡散最小　／　右：CFL=2.0（Δt=0.2）不安定・発散

キーワード