最終更新: 2014/9/14

関数オブジェクトによる問題の記述

1 はじめに

本ライブラリで非線形方程式、常微分方程式、数値積分などの問題を解かせるには、解かせたい問題を「関数オブジェクト」として記述する必要がある。

2 関数オブジェクト

関数オブジェクトとは、()演算子をオーバーロードしたクラスのインスタンスである。例えば、

#include <iostream>

struct Func {
    double operator() (double x) {
        return x * x;
    }
};

int main()
{
    Func f;

    std::cout << f(3) << "\n";
}

を実行すると9と表示される。Funcというクラスで()演算子をオーバーロードしているので、そのクラスのインスタンスfに(3)を作用させることが出来る。

これを用いて、簡単なニュートン法を書いた例を示す。

#include <iostream>

struct Func1 {
    double operator() (double x) {
        return x * x - 2.;
    }
};

struct Func2 {
    double operator() (double x) {
        return 2 * x - 1.;
    }
};

template <class F> double newton(F f, double x)
{
    int i;
    double d;

    for (i=0; i<10; i++) {
        d = (f(x + 0.01) - f(x)) / 0.01;
        x = x - f(x) / d;
    }

    return x;
}

int main()
{
    std::cout << newton(Func1(), 5.) << "\n";
    std::cout << newton(Func2(), 5.) << "\n";
}

関数newtonは、解きたい問題の関数オブジェクトを第1引数として受け取っている。解きたい問題毎に型が違うので、テンプレートを用いている。関数ポインタを用いても似たようなことが実現できるが、関数オブジェクトを用いるとポインタを通じた関数呼び出しのオーバーヘッドが無い。この場合、関数newtonはFunc1用とFunc2用の2種類が自動的に生成されることになる。

3 問題の記述

上のnewtonの例のように、例えば

allsol.hpp 内の関数 allsol
kraw-approx.hpp 内の関数 krawczyk_approx
newton.hpp 内の関数 newton
ode-maffine.hpp 内の関数 odelong_maffine
defint.hpp 内の関数 defint_autostep

などの関数は全て解くべき問題を関数オブジェクトで受け取るように設計されている。

上のnewtonの例では、Func1やFunc2はdoubleを受け取ってdoubleを返す関数であった。しかし例えばallsolでは、

#include <kv/allsol.hpp>
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>

namespace ub = boost::numeric::ublas;

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

        y(0) = x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = x(0) - x(1);

        return y;
    }
};

int main()
{
    ub::vector< kv::interval<double> > I;

    std::cout.precision(17);

    I.resize(2);
    I(0) = kv::interval<double>(-10., 10.);
    I(1) = kv::interval<double>(-10., 10.);

    kv::allsol(Func1(), I);
}

のように書く。これは、方程式

x₀² + x₁² - 1 = 0
x₀ - x₁ = 0

の([-10,10],[-10,10])における全ての解を探索する例である。

まず、ベクトル値を受け取ってベクトル値を返すのに、ublasのvectorを用いる仕様になっている。また、そのvectorの内部の型が(double等でなく)Tとtemplate になっている。これは、例えばallsolでは、渡された関数に対して、 (例えば入力ベクトルに渡した値がinterval<double>だったとして)

interval<double>
autodif<interval<double>>

の2通りの型が入力されるので、どちらにも対応できるようにtemplateにする必要がある。

なお、引数を

    ub::vector<T> x

ではなく、

    const ub::vector<T>& x

のようにconst参照にしているのはコピーのオーバーヘッドを無くすためである。

4 doubleで書けない定数の記述

例えば、上の問題をちょっと変えて、

1.18 x₀² + x₁² - 1 = 0
x₀ - x₁ = 0

を解くことを考える。これはどう記述するべきか。

4.1 単純に書くと…

        y(0) = 1.18 * x(0) * x(0) + x(1) * x(1) - 1.;

だと、1.18はdoubleで正確に記述できないので、1.18に最も近いdoubleを係数に持つ式を解いたことになってしまい、厳密で無くなる。

4.2 分数で書いても…

        y(0) = 118./100. * x(0) * x(0) + x(1) * x(1) - 1.;

と書けば、118も100もdoubleで正確に書けるが、除算に誤差が入ってしまうのでやはりまずい。

4.3 区間演算させる

        y(0) = (kv::interval<double>)118./100. * x(0) * x(0) + x(1) * x(1) - 1.;

のように書けば厳密な計算になる。が、この方法だと、内部型をdd型にした場合にエラーになってしまう。コンパイルを通すには、

#include <kv/allsol.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

namespace ub = boost::numeric::ublas;

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

        y(0) = (kv::interval<kv::dd>)118./100. * x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = x(0) - x(1);

        return y;
    }
};

int main()
{
    ub::vector< kv::interval<kv::dd> > I;

    std::cout.precision(34);

    I.resize(2);
    I(0) = kv::interval<kv::dd>(-10., 10.);
    I(1) = kv::interval<kv::dd>(-10., 10.);

    kv::allsol(Func1(), I);
}

のように、型を書き換える必要がある。出来ればFunc1は書き換えずに済ませたい。また、T=doubleの場合 (Func1に対して単に近似計算したい場合) エラーになってしまう。

4.4 文字定数の機能を使う

intervalやdd (とそれを内部に持つautodif等) は文字列で初期化出来るので、

        y(0) = "1.18" * x(0) * x(0) + x(1) * x(1) - 1.;

と書くことが出来る。この方法だと、interval<double>とinterval<dd> で書き分ける必要は無くなる。が、T=doubleでエラーになるのは変わらず。また、文字列からの変換は時間がかかるので、速度が問題となる。

4.5 現状お勧めの書き方

以下は今後変更されるかも知れないが、現状で最も問題の少ない書き方。数学関数中の円周率などの定数を定義するため、constants<T>という構造体を定義していて、constants.hppをincludeすることによって、

#include <kv/constants.hpp>
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main()
{
    std::cout.precision(17);

    std::cout << kv::constants<double>::pi() << "\n";
    std::cout << kv::constants<double>::e() << "\n";
    std::cout << kv::constants<double>::ln2() << "\n";
    std::cout << kv::constants<double>::str("1.18") << "\n";

    std::cout << kv::constants< kv::interval<double> >::pi() << "\n";
    std::cout << kv::constants< kv::interval<double> >::e() << "\n";
    std::cout << kv::constants< kv::interval<double> >::ln2() << "\n";
    std::cout << kv::constants< kv::interval<double> >::str("1.18") << "\n";
    std::cout << kv::constants< kv::interval<double> >::str("1.18", "1.19") << "\n";

    std::cout.precision(34);

    std::cout << kv::constants<kv::dd>::pi() << "\n";
    std::cout << kv::constants<kv::dd>::e() << "\n";
    std::cout << kv::constants<kv::dd>::ln2() << "\n";
    std::cout << kv::constants<kv::dd>::str("1.18") << "\n";

    std::cout << kv::constants< kv::interval<kv::dd> >::pi() << "\n";
    std::cout << kv::constants< kv::interval<kv::dd> >::e() << "\n";
    std::cout << kv::constants< kv::interval<kv::dd> >::ln2() << "\n";
    std::cout << kv::constants< kv::interval<kv::dd> >::str("1.18") << "\n";
    std::cout << kv::constants< kv::interval<kv::dd> >::str("1.18", "1.19") << "\n";
}

のような使い方をすることが出来る (test-constants.cc)。型に応じたpiやeやlog2を得るために作成したものだが、関数strによって、文字列から定数を作成できる。また、interval型については、文字列を2つ指定して区間定数を生成することも出来る。これを利用して、少し冗長になるが、

#include <kv/allsol.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

namespace ub = boost::numeric::ublas;

struct Func1 {
    template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
        ub::vector<T> y(2);
        static T c1 = kv::constants<T>::str("1.18");

        y(0) = c1 * x(0) * x(0) + x(1) * x(1) - 1.;
        y(1) = x(0) - x(1);

        return y;
    }
};

int main()
{
    ub::vector< kv::interval<kv::dd> > I;

    std::cout.precision(34);

    I.resize(2);
    I(0) = kv::interval<kv::dd>(-10., 10.);
    I(1) = kv::interval<kv::dd>(-10., 10.);

    kv::allsol(Func1(), I);
}

のように書くことが出来る。こう書けば、

内部doubleと内部ddで同じ書き方が出来る。
static指定の変数に代入することにより、初回の関数突入時一回だけの評価になるので、文字列の変換の遅さの問題が解消される。
T=doubleでも動く。

と一応今までの問題は全て解消出来る。ただ、定数の定義が「その場に」書けないので、見た目が良くないのが問題。 C++11で導入された「ユーザ定義リテラル」を使うのがベストな解決方法な気がするが、C++11はまだ完全に普及したとは言い難いので、C++11でしかコンパイルできないような書き方に踏み切るにはまだ躊躇している。