C++

C++でRKF45 Methodを実装した

  • このエントリーをはてなブックマークに追加
  • Pocket

RKF45 Method書いたったたたた

https://github.com/LoliGothick/Cranberries/blob/dev/cranberries/rkf45.hpp

Runge-Kutta-Fehlberg 4(5) method というのはルンゲクッタ法の一種で, 陽的埋め込み型ルンゲクッタ法と呼ばれるものの一つである.

ルンゲクッタ法を使うと数値解析において微分方程式の初期値問題に対して近似解を与えることができる.

ありていに言えば, 厳密に解を計算することが困難な数式の解をコンピュータを使って近似的(数値的)に求めることができる.

なぜ作ったのか

Boost.Odientにそれっぽいのあるんだけど、コンパイルエラーで動かんくてな…(´・_・`)

Boostのソースコード読むの嫌すぎて作ったってわけよ.

ルンゲクッタ法

まず、最も単純なルンゲクッタ法を考えよう.

次のように初期問題を考える.

    \[y'=f(t,y), y(t_0)=y_0\]

y(t) が近似的に求めたい未知関数であり, tにおける勾配はf(t,y)によってty(t)の関数として与えられている.

時刻 t_0における初期値はy_0で与えられている.

 

いま、初期値から十分に小さい刻み幅hで勾配

    \[f(t,x)\]

を計算することで

    \[y(t+h)=y(t) + h f(t)\]

を得る.

これは最も単純なルンゲクッタ法で(前進)オイラー法と呼ばれている.

 

4次のルンゲクッタ法(RK4)では厳密解とRK4のテイラー展開を4次の項まで一致させ、1ステップの推定誤差は

    \[\displaystyle{O(h^{5})}\]

のオーダーになる. 目的の時刻の y を求めるのに必要なステップ数は

    \[\displaystyle{O(h^{-1})}\]

になるので, 全体の推定誤差は

    \[\displaystyle{O(h^{4})}\]

になる.

 

埋め込みルンゲクッタ法(Embedded Runge-Kutta method)

またの名を適応型ルンゲクッタ法(adaptive Runge–Kutta method)とも言うらしい.

ルンゲクッタでは各ステップでの局所誤差を精確に計算することは困難であるが, 誤差を許容範囲内にコントロールすることが望ましい.

埋め込みルンゲクッタは4次のルンゲクッタと5次のルンゲクッタの2つを計算することでその差(誤差)が許容範囲に収まるように各ステップで刻み幅hを調整することを可能にする手法である.

y(t)が急激に変化するようなtの範囲ではhを小さくとり誤差を抑え, 変化量が一定なtの範囲ではhを大きくとるというような制御を可能にする.

シキノートさんが詳しいので気になるひとはこちらで(丸投げ

ルンゲクッタ法の説明と刻み幅制御

Cnraberries.RKF45

Cranberries Libの新しいライブラリである.

Boost.Odientとクソ似ている().

多変数の微分方程式を解くことができる.

使い方

まず, 微分方程式をどう書き表すかだ.

例としてローレンツ方程式を用いる.

カオス的振る舞いをする非線形方程式として有名で, カオス研究の発端となった方程式である(というのを聞いたことがあるような無いような.

dx/dt = p(y-x)

dy/dt = -xz + rx – y

dz/dt = xy – bz

というx, y, zの3つの変数についての方程式になっていてシステムの振る舞いはp,r,bという3つの定数によって決定する.

この方程式をコードで表すと

となる.

3次元ベクトルをstd::array<T,3> で表す.

0,1,2のインデックスがそれぞれx,y,zに対応していてる.

x,y,zに関する方程式をそれぞれそのまま記述すれば良い.

t, arrayを引数にとり, arrayを返す関数としてシステムを記述する.

また, それぞれの変数の初期値

    \[x_0=10, y_0=1, z_0=1\]

を与える必要があるので(値は適当),

これもarrayとして

と記述する.

ここまでで下準備は終了, ソルバーを用意すして方程式を解く.

まず, make_rkf45に解きたい方程式を渡してRKF45クラスを作る.

次に各種設定を行う.

set_integrate_rangeではf(t,x,y,z)を求めるtの範囲を与える, 上記の場合0<t<25の範囲でf(t,x,y,z)を計算する設定となる.

set_toleranceでは許容誤差の限界を与える, 上記の場合8.0 \times 10^{-3}が各ステップでの許容される最小の局所誤差となる.

set_step_size_rangeでは刻み幅hの最小値と最大値を与える, 上記の場合は1.0 \times 10^{-6} < h < 0.05である.

許容誤差がある程度小さい場合は刻み幅を小さくしなければならないため下限には注意が必要である(上限はいい感じに設定されるので大きめにとっても良い). 詳しくは後述する.

ここまでで設定は終わりである.

最後に解いてもらうのに呼び出す関数がintegrateである(微分方程式を解くということは積分するということなので).

integrateの第1引数は初期値のarrayである. 第2引数にはオブザーバを渡す, いくら解けても結果が見えなければ意味がないので各ステップでのt, state(x,y,z)を受け取ってゴニョゴニョする関数オブジェクトを渡す.

上記の場合は標準出力に出力しているだけだが, どう考えても数字の羅列を見ても何もわからないため, 当然グラフにしたいのでバイナリファイルに書き出すなどした方がいいだろう.

ということで, 使い方は以上である.

C++でやるのはここまで!

ここからは隣のPython部屋に移動してグラフを見てみよう.

グラフを描画してみる

まず, テスト条件は…書くの面倒なので実際のコードどうぞ(例と変わってないし…)

バイナリで出力する.

プロットはC++でやるとか考えられないのでPythonで.

スクリプトは誰が書いても似たようになるだろう.

そして, 出力されたのが…これ

Wikipediaのローレンツ方程式のページのグラフとそっくりですね.

問題なく解けているようです!

失敗する場合がある

integrateが途中で失敗することがある.

誤差を許容範囲内に抑えるためにコントロールした刻み幅hが設定したh_{min}より小さくなってしまった場合だ、この場合は処理が中断される.

integrateの戻り値は成功すればtrueで失敗すればfalseとなっている.

まとめ

Pythonはいいぞ!

跋文

以上, Crranberries Lib の新しいライブラリの紹介でした.

何かに必要で書き始めたんですが, 何に必要だったのかも忘れたこの折に完成いたしました(´・_・`)

間違いの指摘や質問等はツイッターアカウントをお持ちであれば

Twitter|いなむのみたま

までお願いします.

Twitterがだめならコメント欄もありますのでそちらへ(反応は鈍いと思います

  • このエントリーをはてなブックマークに追加
  • Pocket

コメントを残す

*

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください