輸送方程式_数値流束.py

1次元移流方程式を数値流束（numerical flux）の形式で定式化し、FTCS・風上差分・Lax法・Lax-Wendroff法の4スキームを比較するコード。

1. 概要

項目	内容
対象方程式	1次元線形移流方程式（スカラー輸送）
比較スキーム	FTCS / 1次風上差分 / Lax法 / Lax-Wendroff法
初期条件	ステップ関数（q=1 / q=0 の2値）
パラメータ	c=1, dx=0.1, dt=0.05, jmax=21, nmax=10（CFL数ν=0.5）
出力	各スキームの時間発展グラフ（n=0, t=0.1, 0.2, …, 0.5）

2. 使用ライブラリ

ライブラリ	用途
numpy	配列演算・等間隔座標生成
matplotlib	時間発展グラフの描画

3. 物理背景

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

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

解析解はステップ関数が速度 c で右に平行移動するもの。数値的には差分スキームによって数値拡散・数値分散・不安定が現れる。

本コードは保存型（フラックス差分型）で離散化する：

q[j]^{n+1} = q[j]^n − (Δt/Δx) × (F_{j+1/2} − F_{j-1/2})

各スキームの違いは数値流束 F の定義にある。

4. 計算手法（数値流束の定義）

スキーム	数値流束 F_{j+1/2}	特徴
FTCS	c/2 × (q[j+1] + q[j])	中心差分。無条件不安定（全CFL数で発散）
UPWIND1	c × q[j]（c>0）	1次風上差分。ν≤1 で安定。数値拡散が大きい
Lax法	c/2 × ((1−1/ν)q[j+1] + (1+1/ν)q[j])	FTCS の q[j]^n を空間平均に置換。ν≤1 で安定。拡散的
Lax-Wendroff法	c/2 × ((1−ν)q[j+1] + (1+ν)q[j])	2次精度（空間・時間とも）。ν≤1 で安定。振動あり

ν = cΔt/Δx（CFL数）。本コードでは ν = 1×0.05/0.1 = 0.5

5. コードの流れ

ステップ	処理
① パラメータ設定	c, dx, dt, jmax, nmax を定数として定義
② `init()`	空間座標 x とステップ関数型初期値 q を生成
③ `do_computing()`	時間ループ内で保存型スキーム q[j] -= dt/dx*(F_{j+1/2} - F_{j-1/2}) を適用。偶数ステップごとに描画
④ 数値流束関数	FTCS / UPWIND1 / LAX / LAXWEN の4関数を定義。ff 引数として do_computing に渡す（関数を引数に渡す高階関数の構造）
⑤ 4回実行	各スキームについて init → do_computing を順番に呼び出し

6. 精度・安定性

スキーム	空間精度	時間精度	安定条件	数値粘性
FTCS	2次	1次	無条件不安定	負（不安定の原因）
UPWIND1	1次	1次	ν ≤ 1	大（拡散的）
Lax	1次	1次	ν ≤ 1	大（拡散的）
Lax-Wendroff	2次	2次	ν ≤ 1	小（振動あり）

7. コード全文

import numpy as np
import matplotlib.pyplot as plt

c=1
dx=0.1
dt=0.05
jmax=21
nmax=10

def init(q1,q2,dx,jmax):
    x=np.linspace(0,dx*(jmax-1),jmax)
    q=np.array([(float(q1)if i< 0.5 *jmax else float(q2))for i in range(jmax)])
    return (x,q)

def do_computing(x,q,c,dt,dx,nmax,ff):
    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):
            ff1=ff(qold,c,dt,dx,j)    ## j+1/2
            ff2=ff(qold,c,dt,dx,j-1)  ## j-1/2
            q[j]=qold[j]-dt/dx*(ff1-ff2)
        if n%2==0:
            plt.plot(x,q,marker='o',lw=2,label=f't={n*dt:.2g}')

    plt.grid(color='black',linestyle='dashed',linewidth=0.5)
    plt.xlabel('x')
    plt.ylabel('q')
    plt.title(ff.__name__)
    plt.legend()
    plt.show()

def FTCS(q,c,dt,dx,j):
    return 0.5*c*(q[j+1]+q[j])

def UPWIND1(q,c,dt,dx,j):
    return c*q[j]

def LAX(q,c,dt,dx,j):
    nu2=1/(c*dt/dx)
    return 0.5*c*((1-nu2)*q[j+1]+(1+nu2)*q[j])

def LAXWEN(q,c,dt,dx,j):
    nu=c*dt/dx
    return 0.5*c*((1-nu)*q[j+1]+(1+nu)*q[j])

c=1
q1,q2=1,0
x,q=init(q1,q2,dx,jmax)
do_computing(x,q,c,dt,dx,nmax,FTCS)
x,q=init(q1,q2,dx,jmax)
do_computing(x,q,c,dt,dx,nmax,UPWIND1)
x,q=init(q1,q2,dx,jmax)
do_computing(x,q,c,dt,dx,nmax,LAX)
x,q=init(q1,q2,dx,jmax)
do_computing(x,q,c,dt,dx,nmax,LAXWEN)

8. グラフ結果

FTCS（無条件不安定）

UPWIND1（1次風上差分・数値拡散大）

LAX法（拡散的・安定）

Lax-Wendroff法（2次精度・振動あり）

9. キーワード

移流方程式数値流束保存型差分 FTCS法風上差分 Lax法 Lax-Wendroff法 CFL数数値拡散数値分散高階関数