白虎21のコンピュータ記録: 6月 2022

2022年6月20日月曜日

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

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

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

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

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

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

heyoka

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

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

インストール
実行環境は以下の通りとなります。
- 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」と実行しましょう。

プログラムの実行

適当にプログラムを書きましょう。

以下は公式にあるサンプルです。(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のライブラリでは得られなかった正しい解を得ることができたので、満足です。

2022年6月5日日曜日

VS Code上のjupyter notebookでmatplotlibのアニメーションを再生する方法(How to play matplotlib's animation with jupyter notebook on VS Code)

こんにちは、WhiteTiger-21です。

タイトル通りの答えは(The answer of title is below)

from IPython.display import HTML
from matplotlib import animation, rc
from matplotlib import pyplot as plt

"""
animへアニメーション処理を行う
"""

HTML(anim.to_jshtml())

です。

細かい話を順に追って

matplotlibでアニメーションを表示する際に、普通は

from matplotlib import animation, rc
from matplotlib import pyplot as plt

def init():
	"""
    アニメーションの処理
    """

def animate(i):
	"""
    アニメーションの処理
    """

"""
ゴニョゴニョしたなにか
"""

anim = animation.FuncAnimation(fig, animate, init_func=init,frames=Nt, interval=1000*(t[2]-t[1])*0.8, blit=True)

plt.show()

という感じで書きますよね。

ただ、jupyter notebookでは

from IPython.display import HTML
from matplotlib import animation, rc
from matplotlib import pyplot as plt

def init():
	"""
    アニメーションの処理
    """

def animate(i):
	"""
    アニメーションの処理
    """

"""
ゴニョゴニョしたなにか
"""

anim = animation.FuncAnimation(fig, animate, init_func=init,frames=Nt, interval=1000*(t[2]-t[1])*0.8, blit=True)
HTML(anim.to_html5_video())

という感じで書きますよね。

これでは画面に真っ黒なプレイヤーの画像が出るだけで再生されません。

これの原因は、VS Codeではビデオコーデックが入っていないからだそうです。

Cannot playback video clip in Jupyter within VSCode #7753

その参照で、別のシステム(plotly)を使うことをおすすめされていますが、matplotlibで構築したシステムならば、こちらを使うほうが楽ですよね。

というわけで、

HTML(anim.to_jshtml())

を使うことがVS Code上でjupyter notebookを動かし、アニメーションを再生する方法として楽なものだと思います。