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で計算する、というような対策をしました。
2016/03/28(月)kv-0.4.31
具体的には、
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
int main()
{
std::cout.precision(17);
std::cout << kv::interval<double>("1e-1000") << "\n";
std::cout << kv::interval<double>("1e1000") << "\n";
}
のようなプログラムの出力が、[0,0] [inf,inf]となってしまう問題です。これを、正しく
[0,4.9406564584124655e-324] [1.7976931348623157e+308,inf]となるように修正しました。
これは、文字列をdoubleに変換する関数stringtod (kv/conv-double.hpp内にある) で、絶対値が2-1075より小さい値を無条件で0に、絶対値が21024より大きい値を無条件でinfにしてしまっていたせいです。正しくは、丸めモードを見て、絶対値が小さくても上向き丸めなら2-1074を、絶対値が大きくても下向き丸めなら21024-2971(C言語で言うところのDBL_MAX)を返すべきでした。
kv-0.4.31では、kv/conv-double.hppとkv/conv-dd.hppをそのように修正しました。