md_packets_FuncAnimation.py 技術ドキュメント

プログラムの動作

md_packets_FuncAnimation.py は、1次元の結合振動子系における波動パケットの衝突および散乱現象をシミュレーションし、リアルタイムで可視化するPythonプログラムです。線形および四次の非線形バネモデルを採用しており、多様な散乱体の設定や初期条件に対応しています。

主な機能:

1次元格子モデル: 質点とバネで構成される1次元格子モデルをシミュレートします。
非線形バネ: 質点間の結合に線形（二次）と非線形（四次）の両方のバネ定数を考慮します。
時間積分: Euler法またはVelocity Verlet法による時間積分が選択可能です。
波動パケット生成: ガウス型波束を初期条件として設定し、単一または二つのパケットの衝突をシミュレートできます。
散乱体: 質量欠陥（重い/軽い質量）、結合欠陥（軟らかい/硬い/非線形バネ）、バリアなど、様々な種類の散乱体を格子中に配置できます。
リアルタイム可視化: matplotlib.animation.FuncAnimation を使用して、質点の変位の動的な変化をリアルタイムでグラフ表示します。
GIF保存: シミュレーションアニメーションをGIFファイルとして保存する機能があります。
解析機能:
- フーリエ変換（FFT）により、波動パケットのモードパワースペクトルを可視化します。
- 系の全エネルギーの時間変化をプロットし、シミュレーションの安定性とエネルギー保存則のチェックを行います。

解決する課題:

このプログラムは、1次元の線形および非線形格子における波動の伝播、散乱、およびエネルギー交換といった複雑な物理現象を、インタラクティブかつ視覚的に理解するためのツールとして機能します。特に、非線形性が波動伝播に与える影響や、異なる種類の散乱体が波動パケットに与える影響の学習・研究に役立ちます。

原理

このシミュレーションは、1次元に配置された $N$ 個の質点の運動を記述します。

1. 質点と結合モデル:

各質点 $i$ は質量 $m_{\text{site}}[i]$ を持ちます。
質点 $i$ と質点 $i+1$ の間には、線形バネ定数 $k_{2\text{_bond}}[i]$ と非線形バネ定数 $k_{4\text{_bond}}[i]$ を持つバネが接続されています。

2. ポテンシャルエネルギー: 質点 $i$ の平衡位置からの変位を $u[i]$ とすると、質点 $i$ と $i+1$ 間のバネによるポテンシャルエネルギー $V_i$ は以下の形で表されます。 $$V_i = \frac{1}{2} k_{2\text{_bond}}[i] (u[i+1]-u[i])^2 + \frac{1}{4} k_{4\text{_bond}}[i] (u[i+1]-u[i])^4$$ 系の全ポテンシャルエネルギー $V$ はこれらの和となります。 $$V = \sum_i \left[ \frac{1}{2} k_{2\text{_bond}}[i] (u[i+1]-u[i])^2 + \frac{1}{4} k_{4\text{_bond}}[i] (u[i+1]-u[i])^4 \right]$$

3. 力の計算: 質点 $j$ に作用する力 $F[j]$ は、隣接するバネからの寄与によって計算されます。ここで、$du[i] = u[i+1] - u[i]$ とし、結合 $i$ と $i+1$ 間の「バネの力」に相当する量 $g[i]$ を以下のように定義します。 $$g[i] = k_{2\text{_bond}}[i] du[i] + k_{4\text{_bond}}[i] du[i]^3$$ すると、内部の質点 $j$ に作用する力 $F[j]$ は、右側のバネから受ける力と左側のバネから受ける力の差として与えられます。 $$F[j] = g[j] - g[j-1]$$ 境界条件は「fixed」であり、両端の質点 $0$ と $N-1$ は固定され、それらの変位 $u[0], u[N-1]$ および速度 $v[0], v[N-1]$ は常に $0$ に設定されます。

4. 運動方程式: 各質点 $j$ の運動は、ニュートンの運動方程式によって記述されます。 $$m_{\text{site}}[j] \frac{d^2u[j]}{dt^2} = F[j]$$

5. 時間積分アルゴリズム: 運動方程式を数値的に解くために、以下のいずれかの時間積分法が使用されます。

Euler法 (step_euler): $v_{new} = v + \Delta t \cdot F/m$ $u_{new} = u + \Delta t \cdot v_{new}$ この実装は修正Euler法に似ていますが、Velocity Verlet法とは異なります。
Velocity Verlet法 (step_verlet): より安定でエネルギー保存性に優れています。 $v_{half} = v + 0.5 \Delta t \cdot F(u)/m$ $u_{new} = u + \Delta t \cdot v_{half}$ $F_{new} = F(u_{new})$ $v_{new} = v_{half} + 0.5 \Delta t \cdot F_{new}/m$

6. 波動パケットの初期条件: ガウス型波束が初期変位 $u(x,0)$ および初期速度 $v(x,0)$ として設定されます。変位: $$u(x,0) = A \exp\left(-\frac{(x-x_0)^2}{2 \sigma^2}\right) \cos(k(x-x_0))$$ 速度: $$v(x,0) = A \exp\left(-\frac{(x-x_0)^2}{2 \sigma^2}\right) \left( \frac{v_g (x-x_0)}{\sigma^2} \cos(k(x-x_0)) + \omega \sin(k(x-x_0)) \right)$$ ここで、$A$ は振幅、$\sigma$ は幅、$k$ は中心波数、$x_0$ はパケットの中心位置です。$\omega$ は線形分散関係、 $v_g = d\omega/dk$ は群速度です。線形近似における分散関係 $\omega(k)$ は、背景の線形バネ定数 $k_{2_0}$ と質量 $m_0$ を用いて以下のように表されます。 $$\omega(k) = 2 \sqrt{k_{2\_0}/m_0} |\sin(ka/2)|$$ 群速度 $v_g(k)$ は以下の通りです。 $$v_g(k) = \sqrt{k_{2\_0}/m_0} a \cos(ka/2) \text{ sign}(k)$$

7. 全エネルギー: 系の全エネルギー $E$ は、運動エネルギーとポテンシャルエネルギーの和として計算されます。 $$E = \sum_j \frac{1}{2} m_{\text{site}}[j] v[j]^2 + \sum_i \left[ \frac{1}{2} k_{2\text{_bond}}[i] (u[i+1]-u[i])^2 + \frac{1}{4} k_{4\text{_bond}}[i] (u[i+1]-u[i])^4 \right]$$ Velocity Verlet法は、エネルギー保存則を比較的良好に満たします。

8. FFT解析: シミュレーション中に一定間隔で質点の変位スナップショット $u(t, x)$ を記録します。これらの空間データに対して離散フーリエ変換を行うことで、時間発展するモードパワースペクトル $|u_k(t)|^2$ を可視化し、特定の波数が時間とともにどのように変化するかを分析できます。

必要な非標準ライブラリとインストール方法

このプログラムの実行には、以下のPythonライブラリが必要です。

numpy: 数値計算
matplotlib: グラフ描画とアニメーション
pillow: (オプション) アニメーションをGIF形式で保存する場合に必要

これらのライブラリは pip コマンドでインストールできます。

pip install numpy matplotlib pillow

必要な入力ファイル

md_packets_FuncAnimation.py は、外部からの入力ファイルを必要としません。すべてのシミュレーションパラメータは、コマンドライン引数として直接指定されるか、プログラム内で設定されたデフォルト値が使用されます。

生成される出力ファイル

--save 1 オプションが指定された場合、シミュレーションのアニメーションがGIFファイルとして保存されます。

ファイル名形式:

生成されるGIFファイルの名前は、設定されたシミュレーションパラメータに基づいて自動的に決定されます。 md_packets_FuncAnimation_<method>_<initial_mode>_<scatterer_mode>_<k2_value>_<k4_value>.gif

<method>: 使用された時間積分法 (euler または verlet)
<initial_mode>: 設定された初期波動パケットのモード (left_only または two_packets)
<scatterer_mode>: 配置された散乱体の種類（heavy_mass, soft_bond, none など）
<k2_value>: 背景の線形バネ定数（例: 5p0 は 5.0 を表す）
<k4_value>: 背景の非線形バネ定数（例: 500p0 は 500.0 を表す）

例:

md_packets_FuncAnimation_verlet_left_only_heavy_mass_5p0_500p0.gif

ファイルの内容:

GIFファイルには、シミュレーション中の質点の変位の時間発展がアニメーションとして記録されます。これはリアルタイムで表示されるアニメーションと同じ内容です。

GIFファイル以外に、シミュレーションの実行中に標準出力にステップごとのエネルギー保存則のチェック結果が出力され、シミュレーション終了後にFFTによるモードパワースペクトルとエネルギーの時間変化を示すMatplotlibグラフが新しいウィンドウで表示されます。これらのグラフはファイルとしては保存されません。

コマンドラインでの使用例 (Usage)

プログラムは、以下の形式でコマンドラインから実行できます。

python md_packets_FuncAnimation.py [OPTIONS]

利用可能な主なオプションは以下の通りです。完全なリストは --help オプションで確認できます。

-h, --help: ヘルプメッセージを表示し、終了します。
--n INT: 質点の数（デフォルト: 200）。
--a FLOAT: 格子間隔（デフォルト: 1.0）。
--m0 FLOAT: 背景の質点質量（デフォルト: 1.0）。
--k2 FLOAT: 背景の線形バネ定数（デフォルト: 5.0）。
--k4 FLOAT: 背景の非線形バネ定数（デフォルト: 500.0）。
--method {euler,verlet}: 時間積分法 (euler または verlet、デフォルト: verlet)。
--dt FLOAT: 時間ステップ（デフォルト: 0.05）。
--nstep INT: シミュレーション/アニメーションステップ数（デフォルト: 1000）。
--interval INT: アニメーションのフレーム間隔 [ms]（デフォルト: 1）。
--save_interval INT: FFTおよびエネルギーのスナップショットを保存する間隔（デフォルト: 10）。
--boundary {fixed}: 境界条件（現在 fixed のみサポート、デフォルト: fixed)。
--A FLOAT: パケットの振幅（デフォルト: 0.20）。
--sigma FLOAT: パケットの幅（デフォルト: 2.0）。
--k0 FLOAT: 中心波数（デフォルト: 0.1）。
--initial_mode {left_only,two_packets}: 初期パケット設定 (left_only または two_packets、デフォルト: left_only)。
--x0_left FLOAT: 左パケットの中心。省略すると n*a*0.25。
--x0_right FLOAT: 右パケットの中心。省略すると n*a*0.75。
--scatterer_mode {none,heavy_mass,light_mass,soft_bond,hard_bond,nonlinear_bond,barrier}: 散乱体の種類（デフォルト: heavy_mass)。
--scatterer_center INT: 散乱体の中心インデックス。省略すると n//2。
--mass_width INT: 質量欠陥の幅（デフォルト: 3）。
--heavy_mass_factor FLOAT: 重い質量欠陥の質量倍率（デフォルト: 10.0）。
--light_mass_factor FLOAT: 軽い質量欠陥の質量倍率（デフォルト: 0.2）。
--bond_width INT: 結合欠陥の幅（デフォルト: 4）。
--soft_bond_k2_factor FLOAT: 軟らかい結合のk2倍率（デフォルト: 0.2）。
--hard_bond_k2_factor FLOAT: 硬い結合のk2倍率（デフォルト: 5.0）。
--nonlinear_bond_k4_factor FLOAT: 非線形結合のk4倍率（デフォルト: 10.0）。
--barrier_half_width INT: バリアの半幅（デフォルト: 3）。
--barrier_mass FLOAT: バリアの質量値（デフォルト: 8.0）。
--barrier_k2 FLOAT: バリアのk2値（デフォルト: 20.0）。
--barrier_k4 FLOAT: バリアのk4値。省略すると背景の k4 が使用されます。
--save {0,1}: アニメーション終了後にGIFを保存するかどうか (0: しない, 1: する、デフォルト: 0)。
--gif_fps INT: GIF保存時のフレームレート [fps]（デフォルト: 30）。
--gif_dpi INT: GIF保存時のDPI（デフォルト: 100）。

コマンドラインでの具体的な使用例

例1: デフォルト設定での実行

最も基本的な実行方法です。デフォルトでは、左から右へ進む一つの波動パケットが、中央に配置された重い質量散乱体と衝突する様子をシミュレーションします。

コマンド:

python md_packets_FuncAnimation.py

実行結果の説明:

アニメーションウィンドウ: matplotlib のウィンドウが開き、1次元格子上の質点の変位が時間とともにアニメーション表示されます。左から来た波動パケットが中心の散乱体で一部反射し、一部透過していく様子が確認できます。
コンソール出力: シミュレーション設定の概要、各質点・バネ定数の最小/最大値、波動パケットの設定、および初期エネルギーが表示されます。シミュレーションの進行中、約50ステップごとに現在のステップ数、時間、および全エネルギーの相対ドリフト（初期エネルギーからのずれ）が報告されます。
終了後のグラフ: アニメーションウィンドウを閉じると、FFTによるモードパワースペクトルと、全エネルギーの相対ドリフトの時間変化を示す2つのグラフが追加のウィンドウで表示されます。

例2: 2つの波動パケットの衝突とGIF保存

--two_packets オプションで2つの波動パケットを設定し、--scatterer_mode none で散乱体なし、--save 1 でGIFファイルを保存します。時間積分法は verlet を使用し、ステップ数を増やしてより長いシミュレーションを行います。

コマンド:

python md_packets_FuncAnimation.py --method verlet --initial_mode two_packets --scatterer_mode none --nstep 2000 --save 1 --gif_fps 20 --gif_dpi 150

実行結果の説明:

アニメーションウィンドウ: matplotlib のウィンドウが開き、左から右へ進むパケットと右から左へ進むパケットが中央で衝突し、互いをすり抜けていく様子がアニメーション表示されます（散乱体がないため）。
コンソール出力: 例1と同様のシミュレーション情報と進行状況が出力されます。「animation_phase」が "display" と "save" の2段階で進むことが示されます。
GIFファイルの生成: アニメーションウィンドウを閉じると、指定された設定に基づいた名前（例: md_packets_FuncAnimation_verlet_two_packets_none_5p0_500p0.gif）でGIFファイルがプログラムと同じディレクトリに保存されます。このGIFファイルには、シミュレーションの全過程がアニメーションとして記録されています。
終了後のグラフ: FFTによるモードパワースペクトルと、全エネルギーの相対ドリフトの時間変化を示すグラフが表示されます。verlet 法を使用しているため、エネルギーのドリフトは比較的低い値に抑えられていることが確認できます。

例3: 非線形性の強いバリアによる散乱

特定の散乱体 (barrier) を設定し、その特性（質量、k2、k4）をデフォルト値から変更して、非線形性の効果を強調します。背景の非線形性も調整します。

コマンド:

python md_packets_FuncAnimation.py --scatterer_mode barrier --barrier_mass 20.0 --barrier_k2 100.0 --barrier_k4 2000.0 --k4 100.0 --nstep 1500

実行結果の説明:

アニメーションウィンドウ: 左から進む波動パケットが、質量の重い、結合の硬い、そして強い非線形性を持つバリア領域に進入し、反射と透過が起こる様子がアニメーション表示されます。強い非線形性の影響により、パケットの形状が変化したり、追加の周波数成分が発生したりする可能性があります。
コンソール出力: 設定されたバリアのパラメータが反映された質点質量とバネ定数の情報が表示され、シミュレーションの進行状況が報告されます。
終了後のグラフ: FFTとエネルギーのグラフが表示されます。非線形性が強いため、FFTスペクトルにはより広範囲の波数成分が見られる可能性があり、Verlet法を使用していてもエネルギー保存則からのわずかなずれが見られることがあります。