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

このコードは意図的に不安定なスキームを実装しています

風下差分法（Downwind scheme / 前進差分）は c > 0 の移流問題では無条件不安定であり、格子やΔtの値に関わらず解が発散します。本コードは安定な風上差分法との比較のために「なぜ風下差分を使ってはいけないか」を視覚的に示す教育目的のコードです。

概要

1次元移流方程式を風下差分法（Downwind scheme）で解くコード。空間差分に前進差分（q[j+1]−q[j]）を用いる点が風上差分法と唯一の違いで、それ以外の設定（CFL=0.5 固定、3段階の格子細分化）は風上差分_2 と完全に同じである。格子を細かくしても発散が止まらないことから、スキームの安定性は離散化方向に依存することが直感的に理解できる。

使用ライブラリ

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

物理背景

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

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

移流速度 c > 0 のとき、情報は 上流（左）から下流（右） へ伝わる。差分スキームが安定であるためには、情報が伝わってくる方向（上流側）のデータを使う必要がある。

風上差分（Upwind）：q[j]−q[j-1]（上流 j-1 を参照）→ 安定
風下差分（Downwind）：q[j+1]−q[j]（下流 j+1 を参照）→ 不安定

von Neumann 安定性解析

風下差分のフォン・ノイマン安定性解析では増幅率 G が：

|G|² = 1 + 2·CFL·(1+CFL)·(1−cosθ) ≥ 1　（任意のθ・任意のCFLで）

CFL > 0 かつ θ ≠ 0 のとき常に |G| > 1 となり、CFL の値や格子の細かさに関わらず無条件不安定。

風上差分との差分式比較

スキーム	差分式	安定性
風上差分（Upwind）	q[j]ⁿ⁺¹ = q[j]ⁿ − c·Δt/Δx · (q[j]ⁿ − q[j-1]ⁿ)	CFL ≤ 1 で安定
風下差分（Downwind）	q[j]ⁿ⁺¹ = q[j]ⁿ − c·Δt/Δx · (q[j+1]ⁿ − q[j]ⁿ)	無条件不安定

コード上の違いは 1文字：qold[j-1] → qold[j+1]

計算パラメータ（風上差分_2 と同一設定）

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

コードの流れ

グローバルパラメータ c=1 を設定
subplots で3パネルを準備し、zip で (dt, dx, jmax, nmax) の3ケースを一括ループ
各ケースで座標配列 x と初期ステップ関数 q_init を生成
nmax=600 ステップ分時間発展（前進差分）し、t=0.1 / 0.2 / 0.3 時点をプロット
全ケースで解が発散する様子が観察できる

コード全文

# -*- coding: utf-8 -*-
"""
Created on Sun Mar 29 16:21:27 2026

@author: kishu
"""

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+1]-qold[j])/(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()

グラフ結果

全3ケース（CFL=0.5）で発散。格子を細かくしても不安定性は解消されない。風上差分_2 と同一パラメータだが、差分方向を変えるだけで挙動が正反対になることが確認できる。