Murman-Cole.py 解説

ライブラリ	用途
`numpy`	配列演算・格子生成
`matplotlib.pyplot`	時刻歴グラフの描画

対象方程式は無粘性バーガース方程式（非線形スカラー保存則）：

∂q/∂t + ∂f/∂x = 0, f(q) = q²/2

移流速度が q 自身に等しいため、波形は自己変形する。特性線の傾き dX/dt = q であり：

q₁ > q₂（左が大、右が小）：特性線が収束 → 衝撃波を形成
q₁ < q₂（左が小、右が大）：特性線が発散 → 膨張波（希薄波）
q が 0 を跨ぐ場合（q₁=1, q₂=−1 など）：音速点（q=0）近傍でエントロピー条件の扱いが重要

Murman-Coleスキームは、もともと超音速〜亜音速（transonic）流れの混合型方程式に対して提案された手法であり、音速点での差分の切り替え（elliptic↔hyperbolic）を行う。本コードでは数値流束 f(q_j) = q_j²/2 を j+1/2 の流束として用いる形で実装されている。

格子設定

パラメータ	値	意味
`xs`	−1	格子左端
`dx`	0.1	空間刻み幅
`dt`	0.05	時間刻み幅
`jmax`	21	格子点数（x = −1〜1）
`nmax`	10	時間ステップ数
CFL数	dt/dx × max\|q\|	初期値依存（例: q₁=2 → CFL=1.0）

差分スキーム（Murman-Cole型数値流束）

数値流束は以下のように定義される：

f̂_{j+1/2} = f(q_j) = q_j² / 2

時間更新：

q_j^{n+1} = q_j^n − (Δt/Δx) × (f̂_{j+1/2} − f̂_{j−1/2})

これは後退（風上側）の値を用いた1次精度保存型差分であり、正の q の流れに対して風上差分として機能する。

初期条件

格子中央（j = jmax/2）を境にステップ関数を設定する：

q(x, 0) = q₁ (x < 0), q₂ (x ≥ 0)

init(q1, q2, dx, jmax) — 格子配列 x とステップ関数初期値 q を生成
MC(q, c, dt, dx, j) — 数値流束 f̂_{j+1/2} = 0.5 × q[j]² を返す
do_computing(...) — 時間ループで q を更新し、2ステップごとにプロット
6通りの (q₁, q₂) の組み合わせに対して繰り返し実行

ケース	q₁	q₂	物理的意味
1	2	0	強い衝撃波（q₁ > q₂）
2	1	2	膨張波（q₁ < q₂、全域正）
3	1	−1	音速点通過・衝撃波（エントロピー条件の確認が重要）
4	−1	1	音速点通過・膨張波
5	1	0	衝撃波（弱い）
6	−1	0	音速点付近・膨張波

空間・時間とも1次精度
正の q では後退差分（風上差分）として機能し数値安定
負の q では前進差分（風下差分）となり、数値不安定のリスクあり
音速点（q = 0 付近）では流束の符号が変わるため、エントロピー違反（非物理的膨張衝撃波）が生じる場合がある

⚠ ケース3（q₁=1, q₂=−1）では Rankine-Hugoniot 条件を満たす衝撃波が現れるが、エントロピー条件（特性線が衝撃波に流れ込む）の検証が必要。

import numpy as np
import matplotlib.pyplot as plt

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

def init(q1,q2,dx,jmax):
    x=np.linspace(xs,xs+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,q1,q2):
    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__} (q1={q1},q2={q2})')
    plt.legend()
    plt.show()


def MC(q,c,dt,dx,j):
    return 0.5*q[j]**2


c=1
order=1
q1,q2=2,0;  x,q=init(q1,q2,dx,jmax); do_computing(x,q,c,dt,dx,nmax,MC,q1,q2)
q1,q2=1,2;  x,q=init(q1,q2,dx,jmax); do_computing(x,q,c,dt,dx,nmax,MC,q1,q2)
q1,q2=1,-1; x,q=init(q1,q2,dx,jmax); do_computing(x,q,c,dt,dx,nmax,MC,q1,q2)
q1,q2=-1,1; x,q=init(q1,q2,dx,jmax); do_computing(x,q,c,dt,dx,nmax,MC,q1,q2)
q1,q2=1,0;  x,q=init(q1,q2,dx,jmax); do_computing(x,q,c,dt,dx,nmax,MC,q1,q2)
q1,q2=-1,0; x,q=init(q1,q2,dx,jmax); do_computing(x,q,c,dt,dx,nmax,MC,q1,q2)

q1=1, q2=-1 — ケース3: q₁=1, q₂=−1（音速点通過・衝撃波）

非線形移流方程式バーガース方程式 Murman-Coleスキーム衝撃波膨張波音速点エントロピー条件保存型差分数値流束 Rankine-Hugoniot条件偏微分方程式

Murman-Cole.py

使用ライブラリ

物理背景

計算手法

格子設定

差分スキーム（Murman-Cole型数値流束）

初期条件

コードの流れ

精度・安定性

コード全文

グラフ結果

キーワード