MultiRootFinder

ROOT
1

Hello All,

I am interested in solving (numerically) a bunch of non-linear equations simultaneously using GSL libraries inside ROOT.

I understand there is an example code
root.cern.ch/root/html/tutorials … oot.C.html

that shows how to solve for 2 unknowns.

I would like to solve for 8 unknowns with 8 equations.

If I look here
root.cern.ch/root/html/ROOT__Mat … inder.html

It says it can solve for n unknowns from n equations.
f1(x1,…xn) = 0
f2(x1,…xn) = 0
…
fn(x1,…xn) = 0

I am wondering how can I DEFINE and WRAP my 8 equations so as to pass it to MultiROOTFinder().

Any hint or pointers are very helpful.

Thanks.

2

Here is an example code for 8 equations (adapted from the 2 equation example)

[code]// example of using multiroot finder
// based on GSL algorithm
//
// The MultiRootFinder is based on GSL and it requires the MathMore library
// installed
//
// Usage:
// >.x exampleMultiRoot.C()
// or
// >.x exampleMultiRoot(algoname,printlevel)
//
// where algoname is for an algorithm not using the derivatives:
// hybridS (default) , hybrid, dnewton, broyden
//
#include “RConfigure.h”

#ifdef R__HAS_MATHMORE
#include “Math/MultiRootFinder.h”
#endif
#include “Math/WrappedMultiTF1.h”
//#include “TF2.h”
#include “TError.h”

// example of using multi root finder based on GSL
// need to use an algorithm not requiring the derivative
//like hybrids (default), hybrid, dnewton, broyden

double MyFunc1 (double *x, double *p1 ) {
return ( p1[0] + x[0] + p1[2]*x[1] + p1[3]*x[2] + p1[4]*x[3] + p1[5]*x[4] + p1[6]*x[5] - p1[1] * x[6] * x[7] );
}

double MyFunc2 (double *x, double *p2 ) {
return ( p2[0] + p2[2]*x[0] + p2[3]*x[1] + p2[4]*x[2] + x[3] + p2[5]*x[4] + p2[6]*x[5] + p2[1] * x[6] * x[7] );
}

double MyFunc3 (double *x, double *p3 ) {
return ( p3[0] + p3[2]*x[0] + x[1] + p3[3]*x[2] + p3[4]*x[3] + p3[5]*x[4] + p3[6]*x[5] - p3[1] * x[6] * sqrt(1 - x[7]*x[7] ) );
}

double MyFunc4 (double *x, double *p4 ) {
return ( p4[0] + p4[2]*x[0] + p4[3]*x[1] + p4[4]*x[2] + p4[5]*x[3] + x[4] + p4[6]*x[5] + p4[1] * x[6] * sqrt(1 - x[7]*x[7] ) );
}

double MyFunc5 (double *x, double *p5 ) {
return ( p5[0] + p5[2]*x[0] + p5[3]*x[1] + x[2] + p5[4]*x[3] + p5[5]*x[4] + p5[6]*x[5] - p5[1] * sqrt(1 - x[6]*x[6]) + p5[7]*x[7] );
}

double MyFunc6 (double *x, double *p6 ) {
return ( p6[0] + p6[2]*x[0] + p6[3]*x[1] + p6[4]*x[2] + p6[5]*x[3] + p6[6]*x[4] + x[5] + p6[1] * sqrt(1 - x[6]*x[6]) + p6[7]*x[7] );
}

double MyFunc7 (double *x, double *p7 ) {
return ( sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]) + p7[1]*x[3] + p7[2]*x[4] + p7[3]*x[5] + p7[4]*x[6] + p7[5]*x[7] - p7[0] );
}

double MyFunc8 (double *x, double *p8 ) {
return ( p8[1]*x[0] + p8[2]*x[1] + p8[3]*x[2] + sqrt(x[3]*x[3] + x[4]*x[5] + x[5]*x[5]) + p8[4]*x[6] + p8[5]*x[7] - p8[0] );
}

void exampleMultiRoot(const char * algo = 0, int printlevel = 1) {

#ifndef R__HAS_MATHMORE
Error(“exampleMultiRoot”,“libMathMore is not available - cannot run this tutorial”);
#else

ROOT::Math::MultiRootFinder r(algo);
//defining the function

TF1 * MyFunc1 = new TF1("MyFunc1", MyFunc1, 0, 0, 7);
TF1 * MyFunc2 = new TF1("MyFunc2", MyFunc2, 0, 0, 7);
TF1 * MyFunc3 = new TF1("MyFunc3", MyFunc3, 0, 0, 7);
TF1 * MyFunc4 = new TF1("MyFunc4", MyFunc4, 0, 0, 7);
TF1 * MyFunc5 = new TF1("MyFunc5", MyFunc5, 0, 0, 8);
TF1 * MyFunc6 = new TF1("MyFunc6", MyFunc6, 0, 0, 8);
TF1 * MyFunc7 = new TF1("MyFunc7", MyFunc7, 0, 0, 6);
TF1 * MyFunc8 = new TF1("MyFunc8", MyFunc8, 0, 0, 6);


double mom = 0.267;
MyFunc1->SetParameters( -0.34,  mom, 0., 0., 0., 0., 0.);
MyFunc2->SetParameters( -0.28,  mom, 0., 0., 0., 0., 0.);
MyFunc3->SetParameters(  0.39,  mom, 0., 0., 0., 0., 0.);
MyFunc4->SetParameters(  0.037, mom, 0., 0., 0., 0., 0.);
MyFunc5->SetParameters( -0.74,  mom, 0., 0., 0., 0., 0., 0.);
MyFunc6->SetParameters( -1.02,  mom, 0., 0., 0., 0., 0., 0.);
MyFunc7->SetParameters(  1.11, 0., 0., 0., 0., 0. );
MyFunc8->SetParameters(  0.97, 0., 0., 0., 0., 0. );

// wrap the functions
ROOT::Math::WrappedMultiTF1 g1(*MyFunc1, 8);
ROOT::Math::WrappedMultiTF1 g2(*MyFunc2, 8);
ROOT::Math::WrappedMultiTF1 g3(*MyFunc3, 8);
ROOT::Math::WrappedMultiTF1 g4(*MyFunc4, 8);
ROOT::Math::WrappedMultiTF1 g5(*MyFunc5, 8);
ROOT::Math::WrappedMultiTF1 g6(*MyFunc6, 8);
ROOT::Math::WrappedMultiTF1 g7(*MyFunc7, 8);
ROOT::Math::WrappedMultiTF1 g8(*MyFunc8, 8);

r.AddFunction(g1);
r.AddFunction(g2);
r.AddFunction(g3);
r.AddFunction(g4);
r.AddFunction(g5);
r.AddFunction(g6);
r.AddFunction(g7);
r.AddFunction(g8);

r.SetPrintLevel(printlevel);

// starting point - This starting point is the solution and when used here, program converges.
double x0[8]={ 0.31, -0.62, 0.86, 0.31, 0.2, 0.9, -0.89, 0.12};
/*
root [0] .x exampleMultiRoot.C(“kDNewton”,1)
GSLMultiRootFinder::Solve:dnewton max iterations 100 and tolerance 1e-06
Info in ROOT::Math::GSLMultiRootFinder::Solve: The iteration converged
GSL Algorithm used is : dnewton
Number of iterations = 4
Root values = x[0] = 0.249173 x[1] = -0.557852 x[2] = 0.926722 x[3] = 0.370827 x[4] = 0.130852 x[5] = 0.833278 x[6] = -0.714794 x[7] = 0.475908
Function values = f[0] = -4.58058e-10 f[1] = 4.58043e-10 f[2] = -5.12292e-10 f[3] = 5.12289e-10 f[4] = 2.1217e-10 f[5] = -2.12151e-10 f[6] = 8.80505e-11 f[7] = 1.00293e-10
*/

// With this starting value, program does not converge.

// double x0[8]={ 0.01, -0.02, 0.06, 0.01, 0., 0.01, -0.09, 0.02};
/*
root [0] .x exampleMultiRoot.C(“kDNewton”,1)
GSLMultiRootFinder::Solve:dnewton max iterations 100 and tolerance 1e-06
Info in ROOT::Math::GSLMultiRootFinder::Solve: exceeded max iterations, reached tolerance is not sufficient; absTol = 1e-06
*/
r.Solve(x0);

#endif

}

[/code]

This code runs. See towards the end of the embedded code on how I run it and the results. It is very sensitive to initial value. If I give the solution, it converges but when initial value is guessed to be some numbers, then it does not converge.

I am not sure if I am defining the equations via TF1 function correctly. Can someone please comment?

Thanks