Solve system of 6 quadratic equations

Newbie
1

Hi all,

I’m attempting to solve a system of quadratic equations using the GSLMultiRootFinder (ROOT::Math::GSLMultiRootFinder) following this example:
ROOT: tutorials/math/exampleMultiRoot.C File Reference

I can easily extend the example to a set of three quadratic equations

ROOT::Math::MultiRootFinder* r = new ROOT::Math::MultiRootFinder();

   TF3 * f1 = new TF3("f1","[0]*[0] - x*x - y*y");
   TF3 * f2 = new TF3("f2","[0]*[0] - x*x - z*z");
   TF3 * f3 = new TF3("f3","[0]*[0] - y*y - z*z");

   double sigma_01 = 1.23803;  double err_sigma_01 = 0.00668558;
   double sigma_02 = 0.846751; double err_sigma_02 = 0.00586008;
   double sigma_12 = 1.24057;  double err_sigma_12 = 0.00926783;

   f1->SetParameter(0,sigma_01);
   f2->SetParameter(0,sigma_02);
   f3->SetParameter(0,sigma_12);

   // wrap the functions
   ROOT::Math::WrappedMultiTF1 g1(*f1,3);
   ROOT::Math::WrappedMultiTF1 g2(*f2,3);
   ROOT::Math::WrappedMultiTF1 g3(*f3,3);
   r->AddFunction(g1);
   r->AddFunction(g2);
   r->AddFunction(g3);
   r->SetPrintLevel(printlevel);

   // starting point
   double x0[3]={0.5, 1.0, 0.6};
   r->Solve(x0);

   std::cout << "Testing result: " << f1->Eval(x0[0], x0[1], x0[2]) << std::endl;
   std::cout << "Testing result: " << f2->Eval(x0[0], x0[1], x0[2]) << std::endl;
   std::cout << "Testing result: " << f3->Eval(x0[0], x0[1], x0[2]) << std::endl;

However, I don’t know how to treat the functions in case of six because there is no such thing as a TF6.
How can I feed the 6-dimensional functions to the root finder?

    // This does not work -- how to initiate a multi-dimensional function?
    ROOT::Math::IMultiGradFunction* f1 = new ROOT::Math::IMultiGradFunction("f1","[0]*[0] - x[0]*x[0] - (x[1]*x[1] + x[2]*x[2] + x[3]*x[3] + x[4]*x[4] + x[5]*x[5])/5");
    ROOT::Math::IMultiGradFunction* f2 = new ROOT::Math::IMultiGradFunction("f2","[0]*[0] - x[1]*x[1] - (x[0]*x[0] + x[2]*x[2] + x[3]*x[3] + x[4]*x[4] + x[5]*x[5])/5");
    ROOT::Math::IMultiGradFunction* f3 = new ROOT::Math::IMultiGradFunction("f3","[0]*[0] - x[2]*x[2] - (x[0]*x[0] + x[1]*x[1] + x[3]*x[3] + x[4]*x[4] + x[5]*x[5])/5");
    ROOT::Math::IMultiGradFunction* f4 = new ROOT::Math::IMultiGradFunction("f4","[0]*[0] - x[3]*x[3] - (x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[4]*x[4] + x[5]*x[5])/5");
    ROOT::Math::IMultiGradFunction* f5 = new ROOT::Math::IMultiGradFunction("f5","[0]*[0] - x[4]*x[4] - (x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[3]*x[3] + x[5]*x[5])/5");
    ROOT::Math::IMultiGradFunction* f6 = new ROOT::Math::IMultiGradFunction("f6","[0]*[0] - x[5]*x[5] - (x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[3]*x[3] + x[4]*x[4])/5");

    double sigma_0 = 1.5;
    double sigma_1 = 1.4;
    double sigma_2 = 1.5;
    double sigma_3 = 1.5;
    double sigma_4 = 1.6;
    double sigma_5 = 1.5;
    f1->SetParameter(0,sigma_0);
    f2->SetParameter(1,sigma_1);
    f3->SetParameter(2,sigma_2);
    f4->SetParameter(3,sigma_3);
    f5->SetParameter(4,sigma_4);
    f6->SetParameter(5,sigma_5);

    ROOT::Math::WrappedMultiTF1 func1(*f1,6);
    ROOT::Math::WrappedMultiTF1 func2(*f2,6);
    ROOT::Math::WrappedMultiTF1 func3(*f3,6);
    ROOT::Math::WrappedMultiTF1 func4(*f4,6);
    ROOT::Math::WrappedMultiTF1 func5(*f5,6);
    ROOT::Math::WrappedMultiTF1 func6(*f6,6);

    r->AddFunction(func1);
    r->AddFunction(func2);
    r->AddFunction(func3);
    r->AddFunction(func4);
    r->AddFunction(func5);
    r->AddFunction(func6);
    
    // starting point
    double x0[6]={1.5, 1.5, 1.5, 1.5, 1.5, 1.5};
    r->Solve(x0);
2

Hi @jenskroeger,

I am inviting @moneta and @jonas to this thread; I am sure both can help you with your issue.

Cheers,
J.

3

Thanks @jalopezg for the quick reply.

@moneta, @jonas, could you have a quick look at this?
I’m sure it’s a trivial issue and a small technicality only…

Many thanks,
Jens

4

Hi @jenskroeger,

indeed you can solve your issue easily. First of all, you can avoid creating the function wrapper yourself by using the templated MultiRootFinder::AddFunction() version, where you also pass the dimension directly. This will create a (WrappedMultiFunction)[root/WrappedFunction.h at master · root-project/root · GitHub] for you, and you can pass any function that can be called with a const double * x array of input values and returns a double.

To create the functions that you are wrapping, you can’t use the ROOT::Math::IMultiGradFunction class, because it’s an abstract class. Instead, you can create lambda functions, if you want to use some modern C++:

void solveDim6EquationSystem() {

    ROOT::Math::MultiRootFinder* r = new ROOT::Math::MultiRootFinder();

    double sigma_0 = 1.5;
    double sigma_1 = 1.4;
    double sigma_2 = 1.5;
    double sigma_3 = 1.5;
    double sigma_4 = 1.6;
    double sigma_5 = 1.5;

    auto f1 = [sigma_0](const double * x)
        { return sigma_0*sigma_0 - x[0]*x[0] - (x[1]*x[1] + x[2]*x[2]
                 + x[3]*x[3] + x[4]*x[4] + x[5]*x[5])/5; };
    auto f2 = [sigma_1](const double * x)
        { return sigma_1*sigma_1 - x[1]*x[1] - (x[0]*x[0] + x[2]*x[2]
                 + x[3]*x[3] + x[4]*x[4] + x[5]*x[5])/5; };
    auto f3 = [sigma_2](const double * x)
        { return sigma_2*sigma_2 - x[2]*x[2] - (x[0]*x[0] + x[1]*x[1]
                 + x[3]*x[3] + x[4]*x[4] + x[5]*x[5])/5; };
    auto f4 = [sigma_3](const double * x)
        { return sigma_3*sigma_3 - x[3]*x[3] - (x[0]*x[0] + x[1]*x[1]
                 + x[2]*x[2] + x[4]*x[4] + x[5]*x[5])/5; };
    auto f5 = [sigma_4](const double * x)
        { return sigma_4*sigma_4 - x[4]*x[4] - (x[0]*x[0] + x[1]*x[1]
                 + x[2]*x[2] + x[3]*x[3] + x[5]*x[5])/5; };
    auto f6 = [sigma_5](const double * x)
        { return sigma_5*sigma_5 - x[5]*x[5] - (x[0]*x[0] + x[1]*x[1]
                 + x[2]*x[2] + x[3]*x[3] + x[4]*x[4])/5; };

    r->AddFunction(f1,6);
    r->AddFunction(f2,6);
    r->AddFunction(f3,6);
    r->AddFunction(f4,6);
    r->AddFunction(f5,6);
    r->AddFunction(f6,6);

    // starting point
    double x0[6]={1.5, 1.5, 1.5, 1.5, 1.5, 1.5};

    int printlevel = 3;
    r->SetPrintLevel(printlevel);
    
    r->Solve(x0);
}

Is this minimal solution fine for you? If you have further questions and need advice, please feel free to ask!

Cheers,
Jonas

5

Great, thanks a lot! This was exactly what I was looking for :+1: