2016/05/08(日)新しいODE solverの性能評価
x''- μ(1-x2)x'+x = 0この方程式はμ≥0の値によって計算のしやすさが変化し、μが1くらいまでなら計算しやすいnon-stiffな方程式ですが、μが大きくなると極めて計算しづらいstiffな方程式になります。ここでは、μ=1とμ=100でこの方程式の軌道を計算し、ステップ幅の変化を調べてみました。ステップ幅はある方法で局所誤差がmachine epsilon程度になるように計算していますが、stiffなODEの場合そこを変えて大きなステップ幅にしても精度保証の根幹となる不動点定理が成立しなくなってしまい小さなステップ幅への修正を強いられてしまいます。初期値は
x(0) = x'(0) = 1とし、t=200まで計算しました。この軌道を相図(phase diagram)で示すと、次のような感じです。まずμ=1の場合:
周期解(limit cycle)に巻き付いている様子がよく分かります。次にμ=100の場合:
こちらは急激な変化を含んだ周期の長い軌道になります。
さて、これを、
- ode-maffine (従来の標準的なアルゴリズム)
- ode-maffine2 (従来の高速アルゴリズム。初期値に関する微分が出来ない欠点がある)
- ode-maffine3 (新しいアルゴリズム)
- awa (Lohnerによる有名な精度保証プログラム)
横軸は時刻tで、縦軸は上のグラフがx、下のグラフがODE Solverが選んだステップ幅です。non-stiffな方程式なので、どれでもあまり大きな違いがないことが分かります。計算時間は、
| アルゴリズム | 計算時間(sec) |
|---|---|
| maffine | 2.276 |
| maffine2 | 1.240 |
| maffine3 | 1.798 |
| awa | 3.815 |
どの方法もμ=1の場合に比べると刻み幅が小さくなっています。変化の激しさに応じて刻み幅を頑張って適応させていることがよく分かります。刻み幅はawaが一番小さく、次にmaffine、maffine2とmaffine3はほぼ同じで比較的大きい刻み幅で計算できていることが分かります。結果として、計算時間にも次のように大きな差が出ました。
| アルゴリズム | 計算時間(sec) |
|---|---|
| maffine | 53.820 |
| maffine2 | 9.287 |
| maffine3 | 13.478 |
| awa | 176.78 |
- maffine → 消滅 (別名で残すかも)
- maffine2 → 従来通り
- maffine3 → maffineに改名
2016/05/06(金)kv-0.4.35
まず、前回(0.4.34)に混ぜてしまったバグの修正です。1.00e-05みたいな数値が1.00e0-5のように表示されてしまうというバグを直しました。
また、常微分方程式のsolverに新しいアルゴリズム(maffine3)を追加しました。kv-0.4.30とstiffなODEで日記に書いたように、maffineはstiffなODEに弱く、maffine2はstiffに強いが初期値に関する精度保証された微分が得られないという状況でしたが、今回、「stiffに強くて初期値に関する微分もできる」maffine3を追加しました。今年のGWはほぼこの実装に費やしてしまいました。
実装したアルゴリズムのアイデアは極めて単純です。常微分方程式の初期値問題を精度よく解くためには、いわゆるWrapping Effectを防ぐために解の初期値に関する微分を利用することがほぼ必須です。それを計算するのに必要なのが「初期値に関する変分方程式」というものです。すなわち、
x'(t)= f(x(t),t)に対して、x*(t)をこの式の真の解としたy(t)に関する行列微分方程式
x(t0) = x0
y'(t) = fx(x*(t), t)y(t)を初期値に関する変分方程式といいます。ある時刻t1におけるyの値y(t1)は、x(t0)にx(t1)を対応させる写像のヤコビ行列と一致します。
y(t0) = I (単位行列)
このyを計算するのに、従来は上の2つの式を連立させて、
x'(t)= f(x(t),t)のようなxとyの組を未知関数とするn+n×n次元の新しい初期値問題を作成し、これを解くという方法を使っていました。その方が実装が楽だったという事情もあります。
y'(t) = fx(x(t), t)y(t)
x(t0) = x0
y(t0) = I (単位行列)
それを変更して、まずxだけの式を解いて真の解x*を区間関数として求め、それを代入してyの式を解く、素直な2段階法が、新しいアルゴリズムです。stiffでない普通の方程式だとどちらを使っても大差ありませんが、stiffな場合だと後者の方が刻み幅が大きくなり、結果として高速になるようです。
アルゴリズムの詳細や性能評価はまた改めて書く予定です。問題が無ければ、将来のバージョンではこの新しいmaffine3をmaffineとし、旧maffineは抹消してしまうことも考えています。
2016/04/28(木)kv-0.4.34
今回は関根晃太さんの指摘で気づいた問題点の修正です。dd.hppでsqrtやlogに負の値を入れたときに正しくエラーになっていなかった問題を修正しました。これは、interval<dd>では問題は発生せず、ddを単体で使ったとき(つまり精度保証でなく近似計算)の場合のみ発生します。
また、conv-double.hpp, conv-dd.hppで、手抜きしてsprintfを使っていたのを、C++らしく書き直しました。
更に、現在ubuntu 16.04を試用中なので、gcc 5.3とclang 3.8で全てのサンプルをコンパイルしてみました。細かく数値まで精査したわけではありませんが、問題なく動いているようです。
2016/04/15(金)kv-0.4.33
今回は、Visual C++でconv-double.hpp, conv-dd.hppがコンパイル出来ないというご指摘を頂いたので修正しました。この2つのファイルは元々luaというスクリプト言語でプロトタイプを開発し、それをC++に移植したものです。luaのif文の中のand(論理積)を&&に直すのを忘れていた箇所があって、何故か古いVCやgccではこれがコンパイルエラーにならなかったので気づかなかった、というものです。@quartorzさん、@nikqさん、ご指摘ありがとうございました。
その他、自動微分の使い方の説明に少し追記しました。
某S先生がC++11の機能を使っていろいろ作ってるのが楽しそうで、kvもC++11に移行したい気分が高まっています。しかしほぼ全体に及ぶ大手術になるので大変だなあ。
2016/04/02(土)kv-0.4.32
またまたバグフィックスです。intervalとddのsinhとcoshに、引数が負で絶対値が大きい場合に、ゼロ除算が発生してしまうバグがありました。
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
typedef kv::interval<double> itv;
int main()
{
std::cout << sinh(itv(-710.)) << "\n";
}
のようなプログラムを実行すると、terminate called after throwing an instance of 'std::domain_error' what(): interval: division by 0 中止 (コアダンプ)のようになってしまっていました。これを、正しく
[-inf,-8.98846e+307]となるように修正しました。誤った数値を黙って返してしまうのは精度保証的には最悪ですが、停止するのでそれよりマシとは言えみっともないバグではあります。
この原因は、sinh(x)=(exp(x)-exp(-x))/2を、expの計算を一回で済まそうとして(exp(x)-1/exp(x))/2で計算していたため、exp(x)が(値がが小さすぎて)ゼロを含むような場合にゼロ除算が起きてしまうというものでした。xの正負で場合分けして、xが負なら(1/exp(-x)-exp(-x))/2で計算する、というような対策をしました。