HarmonicBalance.jl: A Julia suite for nonlinear dynamics using harmonic balance

技術要約：「HarmonicBalance.jl：調和バランス法を用いた非線形力学のためのJuliaスイート」

問題提起
非線形常微分方程式（ODE）は、流体力学、ロボティクス、光力学、振動子ネットワークなど、多様な物理システムの時間発展を支配しています。線形システムは通常、一意の定常状態に収束しますが、非線形システムはしばしば複数の定常解を示し、ヒステリシス、自発的対称性の破れ、ノイズ誘起の切り替えといった現象を引き起こします。これらのシステムを解析する上での大きな課題は、標準的な数値ODEソルバー（初期値問題）は、初期条件に依存した単一の定常状態しか見つけられないことです。その結果、解のランドスケープを完全に探索するには、初期条件空間の無限のサンプリングが必要となりますが、これは計算量的に不可能です。さらに、駆動系を回転座標系に変換することで駆動項を静止させることはできますが、得られた結合非線形代数方程式のすべての根を見つけることは、多変数システムにおいては依然として困難です。なぜなら、標準的な根探索アルゴルズム（例：ニュートン・ラフソン法）は、初期推定値の近傍にある解しか特定できないためです。

手法
著者らは、調和駆動される非線形システムのすべての定常解を見つけるために設計されたJuliaパッケージ、HarmonicBalance.jlを導入します。その手法は、以下の3つのコアコンポーネントを統合しています。

1. 調和バランス近似（Harmonic Balance Ansatz）: 本パッケージは、十分な時間が経過した後、システムの応答が、緩やかに変化する振幅を持つ有限個の調和項の和で近似できると仮定します。調和的な時間依存性を持つ$N$個の2次常微分方程式のシステムに対して、解 $x_i(t)$ は以下のように展開されます：
   $$x_i(t) = \sum_{j=1}^{M} [u_{i,j}(T)\cos(\omega_{i,j}t) + v_{i,j}\sin(\omega_{i,j}t)]$$
   ここで、$T$は粗視化されたタイムスケールです。この近似を元のODEに代入し、高速な振動成分を除去する（振幅の2次時間微分を無視する）ことで、問題は調和変数 $u$ および $v$ に関する一連の自律的な1次常微分方程式へと簡約されます。定常状態は、これらの変数の時間微分をゼロと置くことで求められ、結果として結合多項式代数方程式のシステムが得られます。

2. ホモトピー継続法（Homotopy Continuation）: 得られた多項式システムを解くために、本パッケージはHomotopyContinuation.jlライブラリを介したホモトピー継続法を採用しています。反復的な根探索法とは異なり、ホモトピー継続法は、解析的に解けるシステム（既知の根の数を持つ）をターゲットとなるシステムへと変形させます。このアプローチは、すべての複素数解の発見を保証します。パッケージは、計算コストを管理するために2段階のアルゴリズムを利用します。まず「ウォームアップ」フェーズで、一般的なシステムからパスを追跡して特異なパスをフィルタリングし、次にパラメータ・ホモトピー・フェーズで、残りの非特異なパスを物理的な解へと追跡します。

3. 記号および数値積分: 本パッケージは、調和方程式およびヤコビ行列を自動的に導出するために、記号操作を行うSymbolics.jlを使用します。また、時間依存のシミュレーションのためにDifferentialEquations.jlと連携しており、ユーザーは定常状態の結果を検証したり、断熱的なパラメータ掃引（ヒステリシス）を研究したりすることが可能です。

主な貢献

• 包括的な解の探索: 本パッケージの主要な貢献は、初期条件に依存した単一の解ではなく、調和駆動される非線形システムのすべての定常解（実数解および複素数解）を数値的に決定できる能力です。
• 統一されたフレームワーク: 記号代数、調和バランス、およびホモトピー継続を、高性能かつ柔軟な単一のオープンソースJulia環境へと統合しました。
• 安定性とゆらぎの解析: 解をヤコビ行列の固有値に基づいて、実数か複素数か、あるいは安定か不安定かを自動的に分類するツールを提供します。また、安定な固定点に対する線形応答スペクトルの計算もサポートしています。
• 可視化: グラフィカルユーザーインターフェースとプロット機能を備え、1次元および2次元の位相図、解の枝（ブランチ）、および分岐構造を可視化できます。これにより、複雑な解のトポロジーの探索が容易になります。

結果と例
論文では、以下の例を通じてパッケージの能力を実証しています。

• 駆動型ダフィング振動子（Driven Duffing Oscillator）: 単一のダフィング振動子の複数の定常解を特定することに成功しました。これには、非線形性によって生成される高次調和応答（$3\omega_d$）も含まれます。標準的な共鳴曲線や、1つ、2つ、または3つの解が存在する領域を示す位相図を再現しています。
• 結合ダフィング共振器: 2つの結合振動子システムについて、パラメータ空間における解のランドスケープをマッピングし、複雑なトポロジーを明らかにしています。ホモトピーによって見出された「完全な」解の集合と、断熱的な時間掃引によって到達可能な限定的な部分集合を対比させ、ヒステリシスや枝間のジャンプを実証しています。
• パラメトリック駆動振動子: パラメトリック駆動ダフィング振動子（マチュー方程式）を分析し、単一および二重の安定解が存在する領域を特定し、安定状態の数が変化するピッチフォーク分岐を特徴付けています。
• 性能のスケーリング: $N$個の結合ダフィング振動子の鎖に関するベンチマークでは、ベズート境界（理論的な最大解数）は指数関数的（$3^{2N}$）に増大しますが、物理的な実数解の数はそれよりも大幅に少ないことが示されました。本パッケージは、$N=5$（数百のパスを含む）のシステムを、単一のCPUコア上で数分以内に解くことができ、パス追跡戦略の効率性を際立たせています。

意義
著者らは、HarmonicBalance.jlを、理論的解析と実験的観測の間の溝を埋めるためのツールとして位置づけています。すべての解のランドスケープの探索を可能にすることで、本パッケージは、標準的な時間発展シミュレーションでは見落とされがちな、ヒステリシスや対称性の破れといった実験的に観測可能な現象を特定することを可能にします。このソフトウェアは、任意の運動方程式を特別なGUIの制約なしに定義できる「アブイニシオ（原理的）」なアプローチを提供しており、学術コミュニティにとってアクセシブルな設計となっています。著者らは、このツールが構造力学、非線形光学、光力学、アイシングマシンなど、幅広い分野に適用可能であると述べています。現在のバージョンは決定論的な力学と定常状態に焦溶していますが、モジュール設計により、リミットサイクル検出や高次の平均化法などの将来的な拡張をサポートする意図があります。