C++で常微分方程式を解く方法

私はラグランジュ方程式による二重振り子の計算(たぶん)という常微分方程式をPythonのscipyのsolve_ivpやodeintで解いていましたが、もう少し速度がほしいと考えました。

 

そのため、コンパイル型の言語で常微分方程式を解くことができるライブラリを探していました。

 

そこでC++でBoostのOdeintを用いました。

しかしながら、BoostのOdeintではダメでした。

そこで、ODE-Solverを探したところ、「heyoka」に出会いました。

heyoka

 テイラー展開を手法としては用いているようで、一般的に聞く、オイラーやルンゲクッタ法ではないようです。

このライブラリを使うのに個人的には苦労したので、インストールから実際にプログラムを書くまでを記しておこうかと思います。

1. インストール

   実行環境は以下の通りとなります。
   • Linux:Pop!_OS 22.04LTS(Debian系Ubuntu派生ディストリビューション)
   • gcc:11.2.0
   • Boost:1.74
   • {fmt}:8.1.1
   • LLVM:12~14
   • spdlog:1.9.2
   • TBB:2021.5.0-7
   
   
   これらのインストールは
   「apt install libfmt-dev libspdlog-dev libtbb-dev libboost-serialization-dev clang-13 lldb-13 lld-13」
   としておけば、問題ないかと思われます。
   その後、「git clone https://github.com/bluescarni/heyoka.git」 を実行し、githubからリポジトリをクローンしましょう。
   そして、そのリポジトリのディレクトリに移動し、「cmake .」 を実行し、cmakeの準備をしましょう。
   そして、ライブラリとしてインストールするならば、「cmake --build . --target install」と実行しましょう。
   
   
   

2. プログラムの実行

   適当にプログラムを書きましょう。
   以下は公式にあるサンプルです。(As a simple example, consider the ODE system corresponding to the pendulum)
   single_pendulum.cppと命名。
   
   

   #include <iostream>
   
   #include <heyoka/heyoka.hpp>
   
   using namespace heyoka;
   
   int main()
   {
       // Create the symbolic variables x and v.
       auto [x, v] = make_vars("x", "v");
   
       // Create the integrator object
       // in double precision.
       auto ta = taylor_adaptive <double> {// Definition of the ODE system:
                                         // x' = v
                                         // v' = -9.8 * sin(x)
                                         {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                         // Initial conditions
                                         // for x and v.
                                         {0.05, 0.025}};
   
       // Integrate for 10 time units.
       ta.propagate_for(10.);
   
       // Print the state vector.
       std::cout << "x(10) = " << ta.get_state()[0] << '\n';
       std::cout << "v(10) = " << ta.get_state()[1] << '\n';
   }
     

   
   
   これを
   

   g++ simple_pendulum.cpp -o single_pendulum.out -lheyoka

   と実行。
   こうすると実行ファイルが生成されます。
    
   
   

   まとめ

   今日触ってみた限りでは、あまり高速だとは感じませんでしたし、Pythonとくらべても0.07sしか差がありませんでした。

   しかしながら、バッチ化やもっと正しい時間の進め方があるようなので、もっと高速になるかと思います。

   また、Boostのライブラリでは得られなかった正しい解を得ることができたので、満足です。