Trying to mimic the fsolve function from Python in a cpp scritp for simple functions

  algebra, c++, math, python, symbolic-math

First of all, I’m new to C++. Tbh I find it very hard to "get used to" but the last week I have been trying to "translate" a script from Python to C++ due to computational time requirements.

One of the problem I have is root-finding of simple 1D functions:

First the simple MWE from Python:
Trying to solve ax^2+bx+c = 0 :

def f(x,a,b,c):
   return a*x**2+b*x+c
optimize.fsolve(f,-1.0e2,args=(1,1,-1))

which will return either -1.618 or 0.618 depending on the starting’s guess sign

Ok, in C++ from what I have searched online is a bit more complicated and used the GSL Root finding library.

For simple non-polyonimal functions, it works like a charm, but when second order come in,
the addition of endpoints seem to be a problem when you simpled search for a "quick" solution:

#include <stdio.h>
#include <math.h>
#include <functional>
#include <stdlib.h>
#include <iostream>
#include <vector>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>


struct my_f_params { double a; double b; double c; };

double
my_f (double x, void * p)
  {
    struct my_f_params * params = (struct my_f_params *)p;
    double a = (params->a);
    double b = (params->b);
    double c = (params->c);

    return  (a * x + b) * x + c;
  }

double root (struct my_f_params prms, double r)
{
  int status;
  int iter = 0, max_iter = 50;
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  double x_lo= -5e0, x_hi=11e0;
  gsl_function F;
   

  F.function = &my_f;
  F.params = &prms;

  T = gsl_root_fsolver_falsepos;
  s = gsl_root_fsolver_alloc (T);
  gsl_root_fsolver_set (s, &F, x_lo, x_hi);




  do
    {
      iter++;
      gsl_set_error_handler_off();
      status = gsl_root_fsolver_iterate (s);
      r = gsl_root_fsolver_root (s);
      x_lo = gsl_root_fsolver_x_lower (s);
      x_hi = gsl_root_fsolver_x_upper (s);
      status = gsl_root_test_interval (x_lo, x_hi,
                                       0, 0.001);
    printf("%f  %fn",x_lo,x_hi);
    }
  while (status == GSL_CONTINUE && iter < max_iter);


  return r;
}
int main(int argc, char const *argv[])
{
   struct my_f_params params = {1, 1, -1};
   printf("root of x2+x1-1=0 %fn",  root(params,1.25));
  return 0;
}

now, if the starting x_lo , x_hi cover the 2 solution integrals, it does not proceed to finding the closest one, giving the error

gsl: falsepos.c:74: ERROR: endpoints do not straddle y=0
Default GSL error handler invoked.
Aborted (core dumped)

Have tried a lot of things from google before I posted here.

Really really thank you for your time, anything is appreciated!

Source: Windows Questions C++

LEAVE A COMMENT