2014/10/26(日)rkf45
これは、Runge–Kutta–Fehlberg法という有名なもので、6段5次の公式と5段4次の公式が同時に計算できるように工夫されていて、両者の差を観察することによって(大雑把ではあるが)混入した誤差の評価も出来るという優れものです。
この手法を使ったからと言って精度保証出来るわけでは無いのですが、楽しそうだったのでステップ幅調節らしきものも付けて試しに書いてみて、せっかくなのでkvライブラリ内に混ぜておいたものです。
先日、ちょうどオイラー法やルンゲクッタ法など、初期値問題の解法を講義でやったので、このrkf45を講義内で実演してみたら、なぜか精度が出ない！ ステップ幅調整もしていない素の4段4次ルンゲクッタに負けてる！ 講義中は、「やっぱデモやめた」と誤魔化して先へ進んだのですが、後でチェックしてみたら係数が一箇所間違ってました。精度がでないわけです。
精度保証とは関係ないところですが、ちょっと気分が悪いので、近日中にアップデートしてkv-0.4.11としたいと思います。
2014/10/09(木)kv-0.4.10
webデモ内でlog(対数)が使えなかったバグが見つかったので、修正しました。これはkvライブラリ本体ではなくCGIで書いたpythonからC++へのコンバータのバグなので、本体の問題ではありません。
ライブラリ本体の方はlgammaを追加したり、特異点を持つ関数の数値積分関連を整理したりしてますが、これらはまだ発展途上でドキュメントもありません。更新が落ち着いたらきちんと書きたいと考えています。
2014/10/01(水)山梨精度保証研究会のスライド
yamanashi.pdf
kvライブラリの簡単な説明、使い方 + jsiam2013のスライドからODEの説明を引用、の構成です。他の資料やホワイトボードも使ったのでこのスライドでしゃべった内容を網羅してるわけではありません。(スライドを作っていて途中で力尽きた、とも言う。)
なお、kvライブラリのガンマ関数の部分(gamma, digamma, trigamma)はまだまだ完成度が低く、現在改良中です。今のバージョンはあまり信用しない方がいいかも。
2014/09/15(月)kv-0.4.9
今回は、ポアンカレマップや周期解などの常微分方程式の射撃法による精度保証付き解法の使い方の解説を頑張って書きました。他にも、KrawczykやODEのcallbackの解説を加筆しています。
ライブラリ本体は、ode-affine.hppのステップ幅コントロールのバグfix、gamma関数の改良などです。gamma関数の負の領域の実装が手抜きで極値を含むような大きめの区間を与えたときに非常に大きな区間が返ってきていたのですが、きちんと極値を意識した実装にしました。このためにtrigammaの作成が必要になったのでこれも作成しました。内部ではdigammaのゼロ点をKrawczyk法で精度保証しています。
今日から山梨で研究会です。この研究会についても書けることがあれば書きたいと思います。
2014/08/01(金)kv-0.4.8
newton.hppをちょっと直した、mpfrのラッパーを書いた、ODEの初期値問題の使い方を詳細に書いた、くらいが主な変更です。
ODEの初期値問題はいろいろな解法を試したかったせいで多数のプログラムを作成しており、自分でも思い出しながら使い方を書くのは大変でした。初期値問題のみで、それを使った境界値問題や周期解の問題に対するshooting methodの話はまだです。これらも早く書きたい。
また、MPFRの簡単なラッパーを書き、interval型の両端の数として使えるように整えて、double, dd以外の第3の選択肢を提供しました。以前boost.intervalにmpfrcppを組み合わせて出来ていたことが出来なくなっていたので、早く書きたいと思っていたライブラリの一つです。ラッパーそのものはほぼ一日で書き上がりました。
重めの全解探索問題 (example/test-nishi.cc の5つの解を持つ非線形回路方程式) でベンチマークを取ってみました。
5.96sec (double, -DKV_FASTROUND -DNDEBUG -O3)
29.9sec (dd, -DKV_FASTROUND -DNDEBUG -O3)
326sec (mpfr<53>, -DNDEBUG -O3)
415sec (mpfr<106>, -DNDEBUG -O3)
これを見ると、53bit (doubleと同じ) や106bit (ddと同じ) くらいの精度では遅すぎてMPFRを使う意味はなさそうです。106bitよりも更に高精度が必要な問題に対して仕方なく使うような感じでしょうか。