ROOTでroot finding

CERNのROOTライブラリを使ってroot findingします。
ほとんどGSLのラッパみたいだけどね。

ほぼドキュメントのRoot Finder AlgorithmsとMultiRootFinderのテストの写経ですが少し自分好みに直しています。
かなり簡単だしなかなか使い勝手がよさそうだった。コードが有ればほぼ説明不要かと思われる。

ハマったところは

• RootFinderAlgorithmsはベクトル値関数用のものがなくMultiRootFinderを使う必要があった
• MultiRootFinder::AddFunctionに第二引数を渡してはいけなかった！

とか。

# include "stdafx.h"                     //プリコンパル済みヘッダ
# include <Math/Functor.h>
# include <Math/RootFinderAlgorithms.h>  //ROOT::Math::Roots::Newton
# include <Math/MultiRootFinder.h>       //ROOT::Math::MultiRootFinder
# include <TMath.h>
# include <boost/format.hpp>
# include <iomanip>
# pragma comment(lib, "libCore.lib")
# pragma comment(lib, "libMathMore.lib")

int _tmain(int argc, _TCHAR* argv[])
{
	//1変数1式
	const auto myfunc0       = [](const double x){return sin(x);};
	const auto myfunc0_deriv = [](const double x){return cos(x);};
	ROOT::Math::GradFunctor1D f0(myfunc0, myfunc0_deriv);

	ROOT::Math::Roots::Newton rf;
	rf.SetFunction(f0, TMath::Pi() - 1);
	rf.Solve();
	const auto a = rf.Root();
	std::cout
		<< boost::format("円周率はおよそ%1%です") % boost::io::group(std::setprecision(3), a)
		<< std::endl;


	//2変数2式
	const auto myfunc1       = [](const double* xs){return 1. - xs[0];};
	const auto myfunc2       = [](const double* xs){return 10. * (xs[1] - xs[0]*xs[0]);};
	const auto myfunc1_deriv = [](const double* xs, const unsigned int icoord)->const double{
		if(icoord == 0){
			return -1.;
		}else{
			return 0.;
		}
	};
	const auto myfunc2_deriv = [](const double* xs, const unsigned int icoord)->const double{
		if(icoord == 0){
			return 10.* -2. * xs[0];
		}else{
			return 10.;
		}
	};

	ROOT::Math::GradFunctor f1(myfunc1, myfunc1_deriv, 2);
	ROOT::Math::GradFunctor f2(myfunc2, myfunc2_deriv, 2);
	ROOT::Math::MultiRootFinder rf_multi(ROOT::Math::MultiRootFinder::kNewton);
	rf_multi.AddFunction(f1);
	rf_multi.AddFunction(f2);
	const double xs0[2] = {-1, -1};
	rf_multi.SetPrintLevel(1);
	rf_multi.Solve(xs0);


	return 0;
}

相変わらずMSVCでしか動かしていない。
とりあえずニュートン法でやった。