最終更新: 2018/11/17

Krawczyk法による非線形方程式の解の精度保証

柏木雅英

1. はじめに

非線形方程式とその近似解が与えられたとき、Krawczyk法でその近似解の近くに真の解があるかどうかテストする。テストする前に軽くニュートン法を行って解を改善することも出来る。

2. ファイル構成

kraw-approx.hpp

下請けとして、

interval.hpp
rdouble.hpp
interval-vector.hpp
autodif.hpp
matrix-inversion.hpp
make-candidate.hpp

を使っている。

3. 使い方

test-kraw-appox.ccを参照。

test/test-kraw-approx.ccは、次のような感じ。非線形方程式

x₀² - x₁ - 1 = 0
(x₀-2)² - x₁ - 1 = 0

に対して近似解(1.01, 0.01)を初期値として5回Newton法を回し、その結果を元にして精度保証を行う例。

#include <kv/kraw-approx.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) {
        ub::vector<T> y(2);
        y(0) = x(0) * x(0) - x(1) - 1.;
        y(1) = (x(0) - 2.) * (x(0) - 2.) - x(1) - 1.;
        return y;
    }
};

int main()
{
    ub::vector<double> x;
    ub::vector<itv> ix;
    bool b;

    std::cout.precision(17);

    x.resize(2);
    ix.resize(2);

    x(0) = 1.01;
    x(1) = 0.01;

    b = kv::krawczyk_approx(Func(), x, ix, 5, 1);
    if (b == false) {
        std::cout << "fail\n";
    }
}

引数の意味は以下の通り。

第1引数: 解きたい問題の関数オブジェクト。書き方は「関数オブジェクトによる問題の記述」を参照。
第2引数: 近似解。型はvector<T>。区間では無いので注意。
第3引数: 精度保証された解が上書きされる。。型はvector<interval<T>>。
第4引数: Newton法の反復回数。default値は2。
第5引数: 1なら画面にいろいろ表示される。0なら表示なし。defaultは1。

戻り値は精度保証成功ならtrue、失敗ならfalse。

4. 注意

単純なプログラムで複雑なことは何もしていない。 Newton法を回し、近似解を元にKrawczyk方が成功しそうな候補者集合を作り、それに対してKrawczyk法を試みるだけ。

なお、候補者集合の作り方は、「その近似解を起点としたNewton法の修正量のノルムの2倍を半径とした超立方体区間」とすることが多いが、ここではそれを少しアレンジした方法を用いている。詳細はこちら(未発表資料, 2008/09/19作成)。

(2018/11/17) Newton法で指定された反復回数未満でも、収束したような様子であれば反復をやめてKrawczykテストに移るようにした。