← コード一覧に戻る
輸送方程式_数値流束_流れ方向未知.py
数値流束形式で風上差分（1次・2次）を統一実装し、流速 c の符号に依らず自動で正しい風上方向を選択する
1. 概要
1次元移流方程式を数値流束（conservation form）で解くコード。
流速 c の符号が事前にわからない状況（流れ方向未知）でも対応できるよう、
セル界面フラックスを Rusanov（局所Lax-Friedrichs）型 の形式で構成する。
1次精度（UPWIND1）と2次精度（UPWIND2）の2スキームについて、
c = +1（右向き）と c = −1（左向き）の合計4ケースを計算・可視化する。
2. 使用ライブラリ
ライブラリ 役割
numpy 配列演算・linspace による格子生成 matplotlib.pyplot 時刻歴グラフの描画
3. 物理・数学的背景
支配方程式は1次元線形移流方程式（輸送方程式）：
∂q/∂t + c ∂q/∂x = 0
これを保存型差分形式（conservation form）で離散化する：
q[j]^{n+1} = q[j]^n − (Δt/Δx) ( F_{j+1/2} − F_{j−1/2} )
セル界面フラックス F_{j+1/2} を Rusanov（局所Lax-Friedrichs）型 で定義する：
F_{j+1/2} = 1/2 [ (c·u_R + c·u_L) − |c| (u_R − u_L) ]
ここで u_R・u_L は界面右側・左側の状態量の補間値。
|c| の項が数値粘性（数値拡散）として働き、c の符号によらず 情報伝播方向（風上側）を自動選択する。
c > 0 のとき：|c| = c → F = c·u_L（左側の値を使う ＝ 風上差分）
c < 0 のとき：|c| = −c → F = c·u_R（右側の値を使う ＝ 風上差分）
4. 計算手法
スキーム 補間精度 ステンシル（j+1/2界面） 特徴
UPWIND1
空間1次精度
u_L = q[j]，u_R = q[j+1]
数値拡散大・安定。CFL≤1 で安定
UPWIND2
空間2次精度
u_L = 3/2·q[j] − 1/2·q[j−1]，
数値拡散小・振動が出やすい。より広いステンシルが必要
時間積分は前進オイラー法（FTCS型）。CFL数 = c·Δt/Δx = 1×0.05/0.1 = 0.5 。
5. コードの流れ
1 パラメータ設定 ：dx=0.1, dt=0.05, jmax=21, nmax=10（CFL=0.5）2 init() ：格子点 x と初期値 q を生成（左半分 q1, 右半分 q2 のステップ関数）3 do_computing() ：時間ループを回し、偶数ステップごとにプロット。内部で ff()（数値流束関数）を呼ぶ4 UPWIND1 / UPWIND2 ：界面での左右補間値を計算し、Rusanov 型フラックスを返す5 ケース実行 ：c=+1（q: 1→0, 右向き移流）× 2スキーム、c=−1（q: 0→1, 左向き移流）× 2スキームの計4ケース
6. 解析解との比較・精度
解析解はステップ関数が速度 c で平行移動するだけ（形状維持）
UPWIND1 ：数値拡散でエッジが丸まる（1次精度）。安定かつシンプルUPWIND2 ：拡散が減り波形のシャープさが改善するが、不連続付近で振動が生じやすい（2次精度、非単調）どちらも c の符号を問わず正しく移流する（流れ方向未知問題への対応を確認）
高精度かつ単調性を保つには TVD スキームや MUSCL+リミッタの導入が有効
7. コード全文
import numpy as np
import matplotlib.pyplot as plt
c=1
#dt=0.05
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,order):
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 (order,jmax-order):
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 (f '{ff.__name__} (c={c})' )
plt.legend ()
plt.show ()
def UPWIND1 (q,c,dt,dx,j):
ur=q[j+1 ]
ul=q[j]
fr=c*ur
fl=c*ul
return 1 /2 *((fr+fl)-abs (c)*(ur-ul))
def UPWIND2 (q,c,dt,dx,j):
ur=3 /2 *q[j+1 ]-1 /2 *q[j+2 ]
ul=3 /2 *q[j]-1 /2 *q[j-1 ]
fr=c*ur
fl=c*ul
return 1 /2 *((fr+fl)-abs (c)*(ur-ul))
#***********c=1*****************************
c=1
q1,q2=1 ,0
order=1
x,q=init (q1,q2,dx,jmax)
do_computing (x,q,c,dt,dx,nmax,UPWIND1,order)
order=2
x,q=init (q1,q2,dx,jmax)
do_computing (x,q,c,dt,dx,nmax,UPWIND2,order)
#******************************************
#***********c=-1*****************************
c=-1
q1,q2=0 ,1
order=1
x,q=init (q1,q2,dx,jmax)
do_computing (x,q,c,dt,dx,nmax,UPWIND1,order)
order=2
x,q=init (q1,q2,dx,jmax)
do_computing (x,q,c,dt,dx,nmax,UPWIND2,order)
#******************************************
8. グラフ結果
UPWIND1（1次精度）c = +1（右向き移流）
UPWIND2（2次精度）c = +1（右向き移流）
UPWIND1（1次精度）c = −1（左向き移流）
UPWIND2（2次精度）c = −1（左向き移流）
UPWIND1：数値拡散によりステップエッジが滑らかに丸まる。c の符号によらず正しい方向に移流
UPWIND2：エッジのシャープさが向上するが、不連続付近に振動（アンダー/オーバーシュート）が発生
c = +1 と c = −1 の結果は左右対称となり、フラックス定式化が双方向に対応していることを確認
9. 関連キーワード
移流方程式
数値流束
保存型差分
Rusanov型フラックス
局所Lax-Friedrichs
風上差分
1次精度
2次精度
流れ方向未知
CFL条件
数値拡散
高階関数
