最終更新: 2018/12/10

常微分方程式の初期値問題の精度保証

柏木 雅英

1. はじめに

常微分方程式の初期値問題の精度保証を行う。

ベキ級数演算を利用しており、基本的にアルゴリズムは、

に従っている。

いろいろな手法の優劣を比較するための多数のファイルがある。

2. ファイル構成

多数のファイルがあるが、最も標準的なodelong_maffine関数を構成している ファイルは、

ode.hpp
ode-autodif.hpp
ode-maffine.hpp
ode-param.hpp
ode-callback.hpp

など。下請けとして、

psa.hpp
interval.hpp
interval-vector.hpp
autodif.hpp
affine.hpp
make-candidate.hpp

などが使われている。

以下は実験的アルゴリズムのためのファイル。

ode-affine.hpp
ode-affine-wrapper.hpp
ode-qr.hpp
ode-lohner.hpp
ode-qr-lohner.hpp
ode-maffine2.hpp

3. 基本的な使い方

とりあえず初期値問題を解くには、 KVライブラリのWebデモ の中の、 Verified ODE Solver を使える。それで生成されるSource Codeも参考になるかも。

以下、最も基本的なodelong_maffineの使い方を説明する。 例えば、常微分方程式

dx0/dt = x1
dx1/dt = -x0

を、初期値

x0(0) = 0
x1(0) = 1

としてt=[0,10]で解くには、次のようなファイルを作ればよい (ode-sample1.cc)。 

#include <iostream>
#include <kv/ode-maffine.hpp>

namespace ub = boost::numeric::ublas;
typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
        ub::vector<T> y(2);

        y(0) = x(1);
        y(1) = -x(0);

        return y;
    }
};

int main()
{
    ub::vector<itv> x;
    itv end;
    int r;

    std::cout.precision(17);

    x.resize(2);
    x(0) = 0.;
    x(1) = 1.;

    end = 10.;
    r = kv::odelong_maffine(Func(), x, itv(0.), end);

    if (r == 0) {
        std::cout << "Cannot calculate solution.\n";
    } else if (r == 1) {
        std::cout << "Solution calculated until t = " << end << ".\n";
        std::cout << x << "\n";
    } else {
        std::cout << "Solution calculated.\n";
        std::cout << x << "\n";
    }
}
これを実行すると、 
Solution calculated.
[2]([-0.54402111088937289,-0.54402111088936644],[-0.83907152907645577,-0.83907152907644888])
のようにx0とx1のt=10における値が表示される。 これはsin(10)とcos(10)である。

まず、上の例のFuncのように、解きたい問題の右辺を関数オブジェクトとして 記述する。この問題は自励系なのでtを使っていないが、強制振動系ならば tも使って記述する。

第1引数に解きたい問題の右辺、第2引数に初期値、第3引数に開始時刻、 第4引数に終了時刻を指定する。 計算結果 (終了時の値) は第2引数を上書きし、また計算が途中で終了した 場合はその終了時刻で第4引数を上書きするので、 第2引数及び第4引数は書き込み可能な変数でなければならない。

戻り値(int型)は、次の意味である。

第4引数に無限大を指定すると、出来る限り長い区間で計算する。 このときの戻り値は0または1になる。

第5引数で動作のカスタマイズが出来る。第5引数のデフォルト値は 

ode_param<T>()
というode_param型のインスタンスであり、ODE solverの動作を制御するための パラメータをパックした構造体である。 ode-param.hppで定義されている。 例えば、verbose(饒舌)にして計算中の経過を表示させるには、 次のようにすればよい。 (ode-sample2.cc)。 
#include <iostream>
#include <kv/ode-maffine.hpp>

namespace ub = boost::numeric::ublas;
typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
        ub::vector<T> y(2);

        y(0) = x(1);
        y(1) = -x(0); 
        return y;
    }
};

int main()
{
    ub::vector<itv> x;
    itv end;
    int r;

    std::cout.precision(17);

    x.resize(2);
    x(0) = 0.;
    x(1) = 1.;

    end = 10.;
    r = kv::odelong_maffine(Func(), x, itv(0.), end, kv::ode_param<double>().set_verbose(1));
}
これを実行すると、 
t: [2.0060630097815228,2.0060630097815229]
[2]([0.90675762706802143,0.90675762706802277],[-0.42165223318982914,-0.42165223318982825])
t: [3.8601860547497217,3.8601860547497218]
[2]([-0.65832653095074523,-0.65832653095074278],[-0.75273247481848471,-0.75273247481848204])
t: [5.6744656548386611,5.6744656548386612]
[2]([-0.57181755641836896,-0.57181755641836473],[0.82038081533622209,0.82038081533622587])
t: [7.4500625086106451,7.4500625086106452]
[2]([0.91952774507809498,0.91952774507809965],[0.39302509593101991,0.39302509593102525])
t: [9.2020179791437577,9.2020179791437578]
[2]([0.22092224750854755,0.22092224750855389],[-0.97529142339906616,-0.97529142339905972])
t: [10,10]
[2]([-0.54402111088937289,-0.54402111088936644],[-0.83907152907645577,-0.83907152907644888])
のように途中経過が表示される。 カスタマイズ可能なパラメータは、

パラメータ名 意味 デフォルト値
order order of Taylor expansion 24
epsilon used for deciding error tolerance of each step numeric_limits<T>::epsilon()
iteration maximum number of iterative refinement 2
restart_max 精度保証に失敗した時ステップ幅を半分にしてやり直すが、そのやりなおしの最大回数。 1
ep_reduce affine arithmeticのepsilonの削減アルゴリズムの動作を設定。0なら削減アルゴリズムを使わない。正なら、設定された値がそのままepsilonの数の上限値の設定となる。 0
ep_reduce_limit ep_reduce=0なら意味なし。ep_reduce>0のとき、「epsilonの数がep_reduce_limitを超えたらep_reduceまでepsilonを減らす」という動作になる。ep_reduce_limit ≥ ep_reduceである必要がある。 0
verbose if verbose==1 then become verbose 0

がある。複数のパラメータをセットしたいときは、 

    r = kv::odelong_maffine(Func(), x, itv(0.), end, kv::ode_param<double>().set_verbose(1).set_order(12));
のようにドットを連続して書けばよい。

なお、ode-paramのテンプレートパラメータは、ode_longに渡す値(区間)の 両端に使う型と一致していなければならない。 すなわち、interval<dd>で計算したいなら、ode_param<dd>を 使う必要がある。

odelong_*の名前の関数では、基本的に積分区間を自動分割して計算を行う。 ステップ幅の自動調節の詳細は、 ベキ級数演算について の9章の通り。

4. callback関数を指定する

前節のようにset_verbose(1)を使えば途中経過を表示できる。 しかし、表示せずに更に加工したいとか、あるいはsolverが勝手に刻んだ 分点以外の点での値が欲しいなど、更にカスタマイズしたいことも多いだろう。 そこで、1ステップ進む毎にユーザが指定した関数を呼び出させる機能がある。

例えば、値を「密出力」する(ユーザが指定した時間間隔で結果を表示する) ことができる。

odelong_maffineは更に第6引数を持っており、 そのdefault値はode_callback<T>()である。 ode_callbackはode-callback.hで定義されており、その中身は次のようになっている。 

template <class T> struct ode_callback {
    virtual bool operator()(const interval<T>& start, const interval<T>& end, const ub::vector< interval<T> >& x_s, const ub::vector< interval<T> >& x_e, const ub::vector< psa< interval<T> > >& result) const {
        return true;
    }
};
odelong_maffineでは、1ステップ進む毎に、この第6引数に指定された オブジェクトを という引数で呼び出している。 このdefaultの場合だと何もしない(trueを返すのみ)関数になっているが、 これを継承して別のクラスを作り、それを第6引数に指定することで 1step進む毎に何らかの動作をさせることができる仕組みになっている。

例えば、 

template <class T> struct ode_callback_sample : ode_callback<T> {
    virtual bool operator()(const interval<T>& start, const interval<T>& end, const ub::vector< interval<T> >& x_s, const ub::vector< interval<T> >& x_e, const ub::vector< psa< interval<T> > >& result) const {
        std::cout << "t: " << end << "\n";
        std::cout << x_e << "\n";
        return true;
    }
};
のようにode_callbackを継承してode_callback_sampleを作りoperator()を カスタマイズし、 
    r = kv::odelong_maffine(Func(), ix, itv(0.), end, kv::ode_param<double>(), ode_callback_sample<double>());
のようにカスタマイズした方のインスタンスを渡すと、計算中に 1ステップ進む毎に時刻と値を表示するようになる。

test/test-ode-callback.ccにいくつかカスタマイズのサンプルがある。

注意1 (version 0.4.40で変更された機能) version 0.4.39以前は、このcallback関数の戻り値はvoidだったが、 version 0.4.40でにboolに変更された。odelong_*系の関数はこの戻り値を調べ、もしfalseだったら (指定された終了時刻を無視して)そこまでで計算を中止するようになった。その場合の 戻り値は"3"である。これにより、 version 0.4.39以前のversionのkvを利用してcallback機能を使っているプログラムは コンパイル出来なくなってしまった。 直すには、単に継承して作成したoperator()の戻り値をvoidからboolに変え、 関数本体の最後にreturn true;を補えばよい。 なお、test/test-ode-stop.ccに、この機能を利用して計算を途中で中止するサンプル (ロジスティック方程式で個体数がある上限に達したら中止)がある。

注意2 version 0.4.39以前でtest/test-ode-callback.ccの中にサンプルとして提供されていた

の3つのカスタマイズ例を、version 0.4.40で正式機能に昇格させ、 kv/ode-callback.hppの中に含めるようにした。 これらのサンプルをそのまま同名でコピーして使っていたような場合は、 継承によるクラス定義部を削除すればよい。

5. 初期値に関する微分を得る

射撃法を使う場合など、解の初期値に関する微分が欲しいことがあるが、 それを精度保証付きで得ることが出来る。 前述のodelong_maffineは初期値として区間のベクトルの代わりに 「内部に区間型を持つ自動微分型のベクトル」を指定できるversionがあり、 これを使えば解の初期値に関する微分を得ることが出来る。 例えば次のようにする (ode-sample3.cc)。 
#include <iostream>
#include <kv/ode-maffine.hpp>
#include <kv/autodif.hpp>

namespace ub = boost::numeric::ublas;
typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
        ub::vector<T> y(2);

        y(0) = x(1);
        y(1) = -x(0); 
        return y;
    }
};

int main()
{
    ub::vector<itv> x;
    ub::vector< kv::autodif<itv> > ax;
    itv end;
    int r;
    ub::vector<itv> y;
    ub::matrix<itv> dy;

    std::cout.precision(17);

    x.resize(2);
    x(0) = 0.;
    x(1) = 1.;

    ax = kv::autodif<itv>::init(x);

    end = 10.;
    r = kv::odelong_maffine(Func(), ax, itv(0.), end);

    kv::autodif<itv>::split(ax, y, dy);

    std::cout << y << "\n";
    std::cout << dy << "\n";
}
実行結果は次の通り。 
[2]([-0.54402111088937289,-0.54402111088936644],[-0.83907152907645577,-0.83907152907644888])
[2,2](([-0.83907152907645888,-0.83907152907644544],[-0.54402111088937655,-0.54402111088936322]),([0.54402111088936322,0.54402111088937655],[-0.83907152907645888,-0.83907152907644544]))
dyが、初期値(t=0の値)に対してy(t=10の値)を対応させる関数のヤコビ行列 になっている。

odelong_maffineは高精度化のために自動微分versionでなくても内部で このような量を計算しており、自動微分versionはそれを明示的に 取り出せるようにしただけのことである。

6. 1ステップだけ進む

odelong_maffineは(というかodelong_*は)、1ステップだけ進む関数を 複数回呼び出して繋ぐことにより実現されている。 前述のcallbackシステムを使えば1ステップ毎に何かをさせることは出来るが、 更に凝ったことをしたい場合は1ステップだけ進む関数を単独で使いたい 場合もあろう。

test/test-ode.ccやtest/test-ode-autodif.ccに簡単な使い方が書いてある。

最も簡単な1step進む関数odeは、ode.hで定義されており、 例えば次のように使う (ode-sample4.cc)。 

#include <iostream>
#include <kv/ode.hpp>

namespace ub = boost::numeric::ublas;
typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
        ub::vector<T> y(2);

        y(0) = x(1);
        y(1) = -x(0); 
        return y;
    }
};

int main()
{
    ub::vector<itv> x;
    itv end;
    int r;

    std::cout.precision(17);

    x.resize(2);
    x(0) = 0.;
    x(1) = 1.;

    end = 10.;
    r = kv::ode(Func(), x, itv(0.), end);

    if (r == 0) {
        std::cout << "Cannot calculate solution.\n";
    } else if (r == 1) {
        std::cout << "Solution calculated until t = " << end << ".\n";
        std::cout << x << "\n";
    } else {
        std::cout << "Solution calculated.\n";
        std::cout << x << "\n";
    }
}
実行すると次のようになる。 
Solution calculated until t = [2.0060630097815228,2.0060630097815229].
[2]([0.90675762706802143,0.90675762706802266],[-0.42165223318982914,-0.42165223318982825])
引数の意味は第5引数までodelong_*とほぼ同じ。 第1引数は解きたい問題の右辺、第2引数は初期値、第3引数は開始時刻、 第4引数は終了時刻、第5引数はパラメータ。 計算結果 (終了時の値) は第2引数を上書きし、また指定時刻まで計算出来なかった 場合はその終了時刻で第4引数を上書きするので、 第2引数及び第4引数は書き込み可能な変数でなければならない。

戻り値の意味もodelong_*も同じ。

第4引数に無限大を指定すると、出来る限り長い区間で計算する。 このときの戻り値は0または1になる。

第5引数はパラメータであり、指定の仕方もodelong_*と同様。 使えるパラメータは少し異なり、以下の通り。

パラメータ名 意味 デフォルト値
order order of Taylor expansion 24
epsilon used for deciding error tolerance of each step numeric_limits<T>::epsilon()
iteration maximum number of iterative refinement 2
restart_max 精度保証に失敗した時ステップ幅を半分にしてやり直すが、そのやりなおしの最大回数。 2
autostep trueなら(第4引数に指定された終了時刻を超えない範囲で)step sizeを自動調節し、誤差がepsilon程度になるようにする。falseならepsilonを無視し、指定された終了時刻で精度保証しようとする。 true

なお、第6引数に vector< psa < interval <T> > > 型への 変数へのポインタを渡すと、それがNULLで無かった場合には、 開始時刻から終了時刻までの真の解を含む区間多項式を 受け取ることが出来る。defaultはNULLなので、内部で生成された 区間多項式のデータは捨てられる。 以下に、真の解を含む多項式を取得してそれを表示する例を示す。 ただし、表示が長くなるので次数を5次とした (ode-sample5.cc)。 

#include <iostream>
#include <kv/ode.hpp>

namespace ub = boost::numeric::ublas;
typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
        ub::vector<T> y(2);

        y(0) = x(1);
        y(1) = -x(0); 
        return y;
    }
};

int main()
{
    ub::vector<itv> x;
    itv end;
    int r, i;
    ub::vector< kv::psa<itv> > result_psa;

    std::cout.precision(17);

    x.resize(2);
    x(0) = 0.;
    x(1) = 1.;

    end = 10.;
    r = kv::ode(Func(), x, itv(0.), end, kv::ode_param<double>().set_order(5), &result_psa);

    if (r == 0) {
        std::cout << "Cannot calculate solution.\n";
    } else if (r == 1) {
        std::cout << "Solution calculated until t = " << end << ".\n";
        std::cout << x << "\n";
    } else {
        std::cout << "Solution calculated.\n";
        std::cout << x << "\n";
    }

    for (i=0; i<result_psa.size(); i++) {
        std::cout << result_psa(i) << "\n";
    }
}
実行結果は次の通り。 
Solution calculated until t = [0.0021614881244146895,0.0021614881244146896].
[2]([0.0021614864413251936,0.0021614864413251941],[0.99999766398545342,0.99999766398545354])
[[0,0],[1,1],[-0,-0],[-0.16666666666666669,-0.16666666666666665],[0,0],[0.0083333322518444833,0.008333333333333335]]
[[1,1],[-0,-0],[-0.5,-0.5],[0,0],[0.041666666666666664,0.041666666666666672],[-3.0020668394648476e-06,-0]]
最終項だけ幅の大きい区間となるようなsinとcosのTaylor展開が得られているのが分かる。

なお、関数名は同じodeだが第2引数が自動微分型になっており、初期値に関する 微分を得られる関数が、ode-autodif.hppで定義されている。 使い方は上述のode.hpp内のodeとほぼ同じ。

7. その他の関数(実験的関数を含む)

ode_*系 (1ステップ進む) と、odelong_*系 (多ステップ進む) の関数は、 優劣を比較するために多数の実装がある。 詳細は述べないが、表にまとめておく。

まずはode_* (1ステップ進む) 系。

関数名 header file 説明
(1) ode ode.hpp 最も普通の1step進む関数。
(2) ode ode-autodif.hpp (1)の第2引数をautodifにしたもの。
(3) ode_affine ode-affine.hpp 第2引数がaffine型。affine型を係数に持つpsaを用いた非常に実験的なもの。高性能だが遅いのであまり実用的ではないと思われる。
(4) ode_wrapper ode-affine-wrapper.hpp (3)を少し変形したもの。入力されたaffineベクトルをいったん区間ベクトルに直し、改めてaffine化したもので計算する。affineのダミー変数の数が方程式の次元nより多いとき、(3)より性能が落ちるが速い。
(5) ode_maffine ode-maffine.hpp 第2引数がaffine型。内部では(1)と(2)を呼び出して与式と与式の変分方程式をそれぞれ(1)で解いて平均値形式を構成し、それをaffine型に見せかけている。
(6) ode_maffine2 ode-maffine2.hpp (5)を高速化したもの。初期値に関する微分が出来ないのが欠点。stiffなODEでも比較的大きなstep sizeが取れる。
(7) ode_lohner ode-lohner.hpp Lohner法を実装したもの。
(8) ode_lohner ode-lohner.hpp (7)の第2引数をautodifにしたもの。

以下はodelong_* (多ステップ進む) 系。

関数名 header file 説明
(9) odelong ode.hpp (1)を単に接続したもの。wrapping effectを防ぐ工夫が皆無なので実用性は無い。
(10) odelong ode-autodif.hpp (2)を単に接続したもの。wrapping effectを防ぐ工夫が皆無なので実用性は無い。
(11) odelong_affine ode-affine.hpp (3)を単に接続したもの。第2引数はaffine型。affine型を使っているのでwrapping effectは起きない。
(12) odelong_affine ode-affine.hpp (11)のwrapperで第2引数がただのintervalになっているもの。
(13) odelong_wrapper ode-affine-wrapper.hpp (4)を単に接続したもの。第2引数はaffine型。
(14) odelong_wrapper ode-affine-wrapper.hpp (13)のwrapperで第2引数がただのintervalになっているもの。
(15) odelong_maffine ode-maffine.hpp (5)を単に接続したもの。第2引数はaffine型。affine型を使っているのでwrapping effectは起きない。
(16) odelong_maffine ode-maffine.hpp (15)のwrapperで第2引数がただのintervalになっているもの。
(17) odelong_maffine ode-maffine.hpp (15)のwrapperで第2引数がautodifになっているもの。
(18) odelong_maffine2 ode-maffine2.hpp (6)を単に接続したもの。第2引数はaffine型。affine型を使っているのでwrapping effectは起きない。
(19) odelong_maffine2 ode-maffine2.hpp (18)のwrapperで第2引数がただのintervalになっているもの。
(20) odelong_lohner ode-lohner.hpp (7)を単に接続したもの。wrapping effectを防ぐ仕組みが皆無なので実用性はない。
(21) odelong_lohner ode-lohner.hpp (8)を単に接続したもの。wrapping effectを防ぐ仕組みが皆無なので実用性はない。
(22) odelong_qr_lohner ode-qr-lohner.hpp (7),(8)で平均値形式を構成し、QR分解を用いた方法で接続したもの。第2引数はただのinterval。
(23) odelong_qr_lohner ode-qr-lohner.hpp (22)と同じだが第2引数がautodif。
(24) odelong_qr ode-qr.hpp (1),(2)で平均値形式を構成し、QR分解を用いた方法で接続したもの。第2引数はただのinterval。
(25) odelong_qr ode-qr.hpp (24)と同じだが第2引数がautodif。

たくさんあるが、実際に使用するときにお勧めできるのは上で説明した (1), (16), (17)あたりか。あるいは(16)の高速版としての(19)も勧められる。

8. epsilonの削減アルゴリズムを使う

odelong_maffine (odelong_maffine2, odelong_affineも) では、 affine arithmeticを用いて区間幅の増大を抑制している。 しかし、値を表現するためのダミー変数 epsilon の数が徐々に 増えて計算が遅くなる欠点がある。 具体的には、n変数の常微分方程式では1step毎にepsilonがnずつ増えていく。 Lorenz方程式などchaoticな方程式ではepsilonが増えて計算が遅くなる前に 区間幅が増大して計算が終了してしまうので目立たないが、 chaoticでない普通の方程式では、非常に長い時間に渡って計算すると 徐々に計算が遅くなる。

これを防ぐために、あまり役に立っていないepsilonを削減して速度低下を 抑制する方法が実装されている。 

    r = kv::odelong_maffine(Func(), x, itv(0.), end, kv::ode_param<double>().set_ep_reduce(250).set_ep_reduce_limit(300));
のようにパラメータを設定する。 「epsilonの数が300個を超えたら250個まで減らす。」という意味である。 epsilonの削減作業そのものに時間がかかるため、 削減が発動する頻度を減らすためにこのように少し差を付けている。

例えば、次のようなプログラムで計算速度を比較してみる (ode-sample6.cc)。 

#include <iostream>
#include <kv/ode-maffine.hpp>

namespace ub = boost::numeric::ublas;
typedef kv::interval<double> itv;

struct Func {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
        ub::vector<T> y(2);

        y(0) = x(1);
        y(1) = -x(0); 
        return y;
    }
};

int main()
{
    ub::vector<itv> x;
    itv end;
    int r;

    std::cout.precision(17);

    x.resize(2);
    x(0) = 0.;
    x(1) = 1.;

    end = 10000.;

    // r = kv::odelong_maffine(Func(), x, itv(0.), end);

    r = kv::odelong_maffine(Func(), x, itv(0.), end, kv::ode_param<double>().set_ep_reduce(250).set_ep_reduce_limit(300));

    std::cout << x << "\n";
}
あるPC (core i7 4770K) で、削減無しで16秒かかったのが削減ありだと 5.1秒となった。その分最終結果 (sin(10000)とcos(10000)) の精度は わずかに悪くなる。 250と300は適当に選んだもので、最適パラメータは要調査。

9. 常微分方程式の計算の高速化

ベキ級数演算について の8章にあるように、常微分方程式の近似解を求める際に うまくプログラムを実装するとベキ級数演算の計算量を大きく減らすことが出来る。 特にTaylor展開の次数が高い場合の効果が大きい。 この方法は既に実装されていて、コンパイル時に-DODE_FAST=1とすれば有効、 -DODE_FAST=0とすれば無効になる。デフォルトでは有効になっている。

しかし、実装のまずさもあって安定性に問題があり、実際にODE_FAST絡みの バグが何度も見つかっている。原理的にODE_FAST=0とODE_FAST=1では 全く同じ計算結果になるはずであり、普通に使っていて挙動不審だと感じたら -DODE_FAST=0でコンパイルして計算結果を比較してみて欲しい。 そしてもし異なる結果になった場合は是非作者まで報告して欲しい。 ODE_FAST=0の方は十分安定してると思う。

10. その他

続きはまだ書いてない。