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のライブラリでは得られなかった正しい解を得ることができたので、満足です。
   

「Ray Tracing in One Weekend」をやってみて

最近レイトレやc++、OpenGLなどを触り始めました。
今回はPeter Shirley氏のRay Tracing in One Weekendについて書いていこうかと思います。
日本語訳はここ。
githubはここ。
pdfをココ最近配布されたのですが、url忘れました。
*2019/3/10追記
pdfは ここ。

技術的な解説は他の方がたくさんされていらっしゃって、非常にわかりやすいので、
レイトレ自体の解説は省き、環境構築、注意事項についてのお話です。

 

環境構築について

私はDebianで開発しています。
理由は以下のとおりです。
1.Linuxは開発環境としてぶっ壊しと再生が楽
2.VisualStudioが重すぎてヤダ
3.他の開発を考えるとクソ楽
の3つです。

1と3は似たような感じなのでまとめますが、私はEmacsをプログラミングをする際に使っています。また、freeglutやglfwを触ったりもしています。そして、gitを頻繁に使っています。
Windows環境下ですとソフトウェアのダウンロードがめんどくさいからです。
今回では特別なライブラリが必要ではないのですが、最低必要なエディタにコンパイラ、gitをインストールするはめになります。

2ですが、私の環境におけるEmacs、g++、gdbなどなどその他諸々を制覇した素晴らしいIDEです。ですが、その分重い。後、Windowsでしか起動できないはず。
叩いて3sくらいで開かないとやる気なくしてしまうところがあるからです。
*2019/3/10追記
ウィンドウにツリーとか表示して邪魔。
自動で{}を縮めてくれるのはいいんですけど。
コード一気にみたい。

注意事項

私自身あるあるでググったら出てくるオチなのですが、 Chapter8でsphere.hのclassにmaterial *mat_ptrを入れなければなりません。

In One Weekend

後、MAXFLOATでダメならば、FLT_MAXにすると走ります。

感想

多分ほぼ初めてc++とレイトレを触ったのに近い感じの私ですが、レイトレはどのように実装されて、c++でどのように書かれているのか少しわかりました。
この一冊では球体しか出すことができませんが、レイトレの基礎を触りだけ簡単に理解できる気がします。
CGの基礎をきちんと実装も含めて書いてあるので気になるなぁというお方はエディタとコンパイラさえあればできるのでされてみるのはいかがかなぁと思います。
macOSやLinuxをお使いの方は何もせずできると言っても過言ではないくらい依存ライブラリがありません。
今現在、私は床井教授が出されているOpenGLのをしているので先になると思うのですが、この次のものもチャレンジしてみようかと思うところです。