Location : Home > Languages > Perl > Package
Title : Math::Function::Roots
 		Toolbox Logo

名称

　Math::Function::Roots - 関数の根の探索

バージョン

　Version 0.03

概要

　任意の（数学的な意味における）関数の根を探索する（Perl の意味における）関数を集めたものである。関数はサブ参照及び範囲（range）または解の推定値（guess of the answer）をとり、関数の根を返す。

use Math::Function::Roots qw(bisection epsilon max_iter);

epsilon(0); # 望む精度の設定
max_iter(50_000) # 実行回数の上限設定

# (-5,5) の範囲から 2x+1 の根を探索
my $result = bisection( sub {2*shift()+1}, -5, 5);
# $result == -.5

# epsilon 及び max_iter を設定する別の方法
my $result2 = bisection( sub {2*shift()+1}, -5, 5, 
    epsilon=>.00001, 
    max_iter=>20);

説明

　数値解析は正確な解を求めることが困難な問題の解の近似値を求めるためにしばしば繰り返しをすることで求めるアルゴリズムを用いる方法である。根探索アルゴリズム（Root Finding Algorithms）は関数の根を探索するために利用される。これは（あらゆる x に対し１つの f(x) の値が対応する）数学的連続関数を扱う。根とは関数の値が 0 となる値であり、すなわち x 軸との交点である。複数の根を探索するアルゴリズムもあるが、多くは１つだけ探索するものである。
　しかし根探索アルゴリズムがなんであるかを知っているのであればそれで十分である。ここでは以下に述べるアルゴリズムの実装から開始している。基本的であるのはアルゴリズム（関数・推定値）であり、以下の各関数は元のアルゴリズムにちなんで名づけられている。

パラメータ

　全てのアルゴリズムは類似のパラメータを利用しているため、１度しか記述しない。根を探索すべき関数は常に必須で指定しなければならない。

function

　Function はコード参照として渡される。"\&Function" または "sub{...}" といった形状で記述される。単純な関数であれば sub{} という形でインラインで表示し、複雑な場合は参照で記述することが推奨される。

# f(x) = 2x - 4
# sub{ 2*shift() - 4 used as
my $root = bisection( sub{ 2*shift() - 4 }, -10, 10 );

　複数の変数を持つ関数を利用することが多いだろう。これはラッパ関数として以下のように利用される。

sub foo{
   my ($x1,$x2) = @_;
   return $x1**2 + $x2**2;
}

my $wrapper = sub {
   my $x2 = shift;
   return foo( 5, $x2 );
}

my $result = bisection( $wrapper, -10, 10 );

　サブルーチンが渡されれば、１引数で呼び出され、１つの結果を返すことが期待される。説明にフィットしない関数はラッパを必要とする。
　ここでは f(x) = 5**2 + x**2 の根を求めるものとする。関数のあるタイプにはよいアルゴリズムとなるであろう。

guess, min/max

　たいていのアルゴリズムは初期範囲または予測値を必要とする。予測値が複数指定されても解は [guess1,guess2] の間にある必要はないが、範囲または最小値・最大値が指定されれば、解は [min, max] の間になければならない。

epsilon, max_iter

　Epsilon（e）は要求する精度を設定するために用いられる。精度を大きくすると繰り返し回数も少なく計算も速く終了する。一般に e は真値と得られた解との最大乖離を示す。もし小数点第6位までの精度がほしければ e = .000_000_1 と設定する。この値がデフォルトである。計算丸め誤差を防ぐために、要求される精度よりも細かい一般的に推奨される精度で計算する。epsilon はアルゴリズム実行時に e に設定されたパラメータの名称で、常に以下のように設定すべき値である。

my $result = bisection( sub{...}, -10, 10, epsilon => .01 );

　epsilon() 関数は e をグローバルに定義することになるので注意が必要である。
　epsilon と同様に、アルゴリズムの繰り返し回数の最大値を max_iter で設定する。グローバルに設定するには max_iter(i) と定義する。この最大値はエラーを把握する際に用いられる。すなわち所与の関数が収束しないか、バグがある場合（あぁ）である。これを実行時の制御に用いてはならない。解を速く得たい場合には epsilon を大きくすればよい。この値を変更する唯一の理由はゆっくり収束させることで、そのことによりよい値を得ることでアルゴリズムが許す最大限度まで計算を行う。デフォルト値は 5,000 である。

パフォーマンスの確認

last_iter( )

　最終結果を得るための繰り返し計算回数を返す。これはデータに対しどの程度アルゴリズムが有効に機能したかを示す。

アルゴリズム

　以下に可能なアルゴリズムを列挙する。関数には、特に根の周辺部分で多くの制約がある。

bisection( function, min, max )

　２分法を用いる。平均的なパフォーマンスを示すが、信頼性が高い。根の存在する範囲を min 及び max で指定する。f(min) と f(max) の符号は逆でなければならない（両者の間に少なくとも１箇所 0 を横切るため）。所与の範囲に複数の根がある場合にはたいていの場合うまく動かない。このメソッドは関数の形にかかわらず信頼性がある。関数の形状から判断するアルゴリズムに比べて幾分遅いかも知れない。

fixed_point( fixed point function, guess )

　定点法アルゴリズム（The Fixed-Point Iteration algorithm）は高速で頑強なメソッドであるが、不幸にも、適用できる問題の領域が限られており、いくつかの幾何を要請する。長所は高速収束であり、根の範囲は必ずしも既知でなくてよく、最終的に収束すればよい。
　不動点とは g(x) = x である点である。メソッドは f(x) が根を持つ不動点を持つ関数 g(x) を探索する。g(x) = x - f(x) とすれば求めることができる。一般的なケースとして、f(x) から導出された x = ff(x)と同じで、 x を因数に持つ g(x) = x = ff(x) により解く。
　上述したように g(x) の選択に関しては制限があり、g(x) を微分した関数の絶対値は１より小さくなければならない； |g'(x)| < 1。g'(x) が 0 のに近いほど収束は速い。
　不動点を含む範囲 [a,b] で |g'(x)| < 1 を満たす関数を考える。これは無限の範囲または関数の一部かも知れない。初期推定値がこの範囲にある限りはアルゴリズムは収束する。
　guess は解の推定値である。範囲 [a,b] 内にさえあれば、guess と真値との関係にかかわらず収束する。
　関数によっては容易に不動点関数に変換できる。また初期値にかかわらず微分値が 0 近傍であれば収束は速い。

secant( function, guess1, guess2 )

　割線法（The secant method）は、関数の根のよりよい推定値に関数の微分を使うニュートン法の単純化である。割線法は傾き（関数上の２点を結んだ線分）を既知または関数の微分から計算された値を代入して求めるものである。
　通常、関数が与えられ、２つの推定値が与えられる。２分法と異なり、これには括弧は必要ない。傾きが 0 に近い、極小または極大値近辺ではこのアルゴリズムはうまく機能しない。しかし２つの推計値が真値に近ければこのアルゴリズムは高速で収束する。

false_position( function, min, max )

　挟み撃ち法（False Position）は割線法と似ており、関数の傾きをよりよい推計値として用いる。違いは２分法の括弧付けと割線法の速度を併せ持った点である。
　括弧付けは望ましい属性であり、アルゴリズムの信頼性を高めている。括弧付けによりその範囲内に解があることを示している。これは高次関数で厳密に近い解を求めたいときに有用である。
　制限は、 [min,max] の範囲内で微分値が 0 とならないことであることと、範囲内に根は１つしかないことである。また f(min) と f(max) は符号が逆でなければならない。

将来の改善

　最優先事項はアルゴリズムの追加である。それから与えられたデータに対して最適なアルゴリズムを選択する機構を組み込むことに興味を持って言えう。最後に Inline::C または XS を使って C におけるアルゴリズムを書き換え、望ましいパフォーマンスを得ることである。理想的には C コンパイラが利用可能であればこれを稼動させ、 C バージョンをコンパイルして利用し、そうでなければ Perl バージョンが利用されるようなものを考えている。このような例を見たことはあるが、どのように実現するのか今のところ理解できていないのでほうっておいたままである。
　それからテストカバー率を上げることである。

著者

　Spencer Ogden, <spencer@spencerogden.com>

バグ

　バグ報告や機能の要望は bug-algorithm-bisection@rt.cpan.org にメールを送るか、http://rt.cpan.org を通じて行ってほしい。通知してもらえればバグに関する変更を自動的に報告する。

著作権とライセンス

　Copyright 2005 Spencer Ogden, All Rights Reserved.

　本プログラムはフリーソフトであり、 Perl 本体と同等の条件で修正／再配布してもよい。

Toolbox Logo
 Updated : 2007/04/23