方程式，解くべし

方程式の解を求める

GSLを用いて, 方程式 f(x)=0 の解を求めましょう. なお, f(x)は連続であるとします. GSLにおいては, 方程式の解は2種類に分類でき, それぞれに適用できる解法が異なります:

• 奇数重解   この場合の特徴は, 解の近傍で f(x-ε)f(x+ε)<0 であることです.
  • bracketing solver: 中間値の定理を使って解を求めていきます. f(x)<0 となるxと, f(x)>0 となるxを与え, その区間を縮小していきます. その区間に存在する, いずれかの解が得られるでしょう. 発散することはまずありません.
    • gsl_root_fsolver_bisection: 2分木探索. 基本的だが遅い
    • gsl_root_fsolver_falsepos:  区間を2等分するのではなく線形補間を使う. 少し速い. 
    • gsl_root_fsolver_brent: Brent-Dekker法で高速.
• 偶数重解   この場合の特徴は, 解の近傍で f(x-ε)f(x+ε)>0 であることです.
  • root polishing solver: f(x)とf'(x)を使って解を求めます. 奇数重解でも効果を持ちます. 解の存在範囲を指定できないので，あなたが想定した解とは異なる解が得られても, 文句を言ってはいけません. この方法で得られるのは, f(x)=0 となる何かのx, ということだけです. 微係数を利用するので, ちょっと関数が妙な形をしていると, すぐに発散するでしょう.
    • gsl_fdfsolver_newton: いわゆるNewton法.
    • gsl_fdfsolver_secant: 簡単化したNewton法.
    • gsl_fsfsolver_steffenson: Aitken加速したNewton法. 最も高速だが, 不安定化するケースがある. 簡単な問題だとNewton法に負ける.

まあ面倒くせえので, gsl_fsfsolver_steffensonを使う例を示します.

#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
void main(){//メインプログラム
	//ソルバーを準備
	auto solver=gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_steffenson);
	gsl_function_fdf FDF;
	//関数のパラメータを準備. structでもOK
	class my_param {
	public:
		double a, b, c;
	};
	my_param params;
	params.a=1.0; params.b=0.0; params.c=-5.0;
	FDF.params = &params;//パラメータを登録
	//関数fだけが計算できるように準備. ラムダ式の例. もちろん関数ポインタでもよい
	FDF.f =  [](double x, void *params){
		auto& p = *(class my_param *) params;
		return ( p.a * x + p.b) * x + p.c;
	};
	//f'を準備
	FDF.df = [](double x, void *params){
		auto& p = *(class my_param *) params;
		return 2.0 * p.a * x + p.b;
	};
	//fとf'を同時に計算する場合を準備
	FDF.fdf = [](double x, void *params,double *y, double *dy){
		auto& p = *(class my_param *) params;
		*y = (p.a * x + p.b) * x + p.c;
		*dy = 2.0 * p.a * x + p.b;
	};
	double x_value=5.0,err_abs=0.0,err_rel=1e-10;//解と精度
	int status;//ステータスコード
	int iter = 0, max_iter = 100;//繰り返しの準備
	gsl_root_fdfsolver_set (solver, &FDF, x_value);
	std::cout << "Using solver [" << gsl_root_fdfsolver_name(solver)
       << "]" << std::endl;
	do {
		iter++;
		status = gsl_root_fdfsolver_iterate (solver);
		double x_prev = x_value;
		x_value = gsl_root_fdfsolver_root (solver);
		status = gsl_root_test_delta (x_value, x_prev, err_abs, err_rel);
		std::cout <<   iter << " " << x_value << " " 
             << x_value-x_prev << std::endl;
	} while (status == GSL_CONTINUE && iter < max_iter);
	//よいこは, 後片付けを忘れません
	gsl_root_fdfsolver_free (solver);
}

Linuxでコンパイルする場合は, ic++ -std=c++11 --gsl ファイル名.cpp です. んまあこんな感じ

[sugimoto@sun2 C++]$ ic++ gsl.cpp --gsl -std=c++11 -o gsl
[sugimoto@sun2 C++]$ ./gsl
Using solver [steffenson]
1 3 -2
2 2.33333 -0.666667
3 2.22222 -0.111111
4 2.23602 0.0138026
5 2.23607 4.31324e-05
6 2.23607 4.16014e-10
7 2.23607 0
[sugimoto@sun2 C++]$

になる.