Finding roots of a function

Question

I would like to re-write a code from python to cython, and so far I cythonized all the parts which I simplified in this example by not using them. Therefore I could not keep the python shape of this function. However, I need to estimate the root of the function which before the scipy.optimize library has been used. I am wondering what I could substitute to find the roots for this function in cython. Could gsl also provide better tool to find the roots? How should it be done?

def RsMassInsideR(mass, R):
    def f(x):

        xp = R/x

        return (np.log(1+xp) - (xp/(1+xp)))*4*np.pi*delta*rho*x**3 - mass #rho and delta are constant

    try:

        rs = scipy.optimize.brenth(f, 0.01, 10.)

    except ValueError, e:
        print '!!!!!!!!!!!'
        print mass, f(0.01), f(10.)
        raise e

    return rs

score 2 · Accepted Answer · answered Jul 12 '14 at 18:05

I once was facing a similar problem for multidimensional root finding. Since I found a solution using GSL, I'd like to share my code here. The pyx file I used was

cdef extern from "gsl/gsl_errno.h":
    char * gsl_strerror(int gsl_errno)

cdef extern from "gsl/gsl_vector.h":
    ctypedef struct gsl_vector:
        pass
    ctypedef gsl_vector* const_gsl_vector "const gsl_vector*"
    gsl_vector* gsl_vector_alloc(size_t n)
    void gsl_vector_free(gsl_vector* v)
    void gsl_vector_set(gsl_vector* v, size_t i, double x)
    double gsl_vector_get(const_gsl_vector v, size_t i)

cdef extern from "gsl/gsl_multiroots.h":
    # structures
    ctypedef struct gsl_multiroot_function:
        int (*f) (const_gsl_vector x, void* params, gsl_vector* f)
        size_t n
        void* params
    ctypedef struct gsl_multiroot_fsolver_type:
        pass
    ctypedef gsl_multiroot_fsolver_type* const_gsl_multiroot_fsolver_type "const gsl_multiroot_fsolver_type*"
    ctypedef struct gsl_multiroot_fsolver:
        gsl_multiroot_fsolver_type* type
        gsl_multiroot_function* function
        gsl_vector* x
        gsl_vector* f
        gsl_vector* dx
        void* state

    # variables
    gsl_multiroot_fsolver_type* gsl_multiroot_fsolver_hybrids

    # functions
    gsl_multiroot_fsolver* gsl_multiroot_fsolver_alloc(
                                    gsl_multiroot_fsolver_type* T, size_t n)
    void gsl_multiroot_fsolver_free(gsl_multiroot_fsolver* s)
    int gsl_multiroot_fsolver_set(gsl_multiroot_fsolver* s, 
                                  gsl_multiroot_function* f, 
                                  const_gsl_vector x)
    int gsl_multiroot_fsolver_iterate(gsl_multiroot_fsolver* s)
    int gsl_multiroot_test_residual(const_gsl_vector f, double epsabs)


DEF GSL_SUCCESS = 0
DEF GSL_CONTINUE = -2


cdef size_t n = 2

cdef int rosenbrock_f(const_gsl_vector x, void* params, gsl_vector* f):
    cdef double a = 1.
    cdef double b = 10.

    cdef double x0 = gsl_vector_get(x, 0)
    cdef double x1 = gsl_vector_get(x, 1)

    cdef double y0 = a * (1 - x0)
    cdef double y1 = b * (x1 - x0 * x0)

    gsl_vector_set(f, 0, y0)
    gsl_vector_set(f, 1, y1)

    return GSL_SUCCESS


def gsl_find_root():
    #cdef const_gsl_multiroot_fsolver_type T
    cdef gsl_multiroot_fsolver* s

    cdef int status = GSL_CONTINUE
    cdef size_t i, iter = 0

    cdef gsl_multiroot_function func
    func.f = &rosenbrock_f
    func.n = n

    cdef double x_init[2]
    x_init[0] = -10.0
    x_init[1] = -5.0

    cdef gsl_vector* x = gsl_vector_alloc(n)

    gsl_vector_set (x, 0, x_init[0])
    gsl_vector_set (x, 1, x_init[1])

    s = gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_hybrids, n)
    gsl_multiroot_fsolver_set(s, &func, x)

    while iter < 100 and status == GSL_CONTINUE:
        status = gsl_multiroot_fsolver_iterate(s)
        if status != GSL_SUCCESS:
            break

        status = GSL_CONTINUE
        print "%d: %f, %f" % (iter, gsl_vector_get (s.x, 0), gsl_vector_get (s.x, 1))

        status = gsl_multiroot_test_residual(s.f, 1e-7)

        iter += 1

    print("status = %s" % gsl_strerror(status))

    gsl_multiroot_fsolver_free(s)
    gsl_vector_free(x)

It should contain the necessary concepts and it should be easy enough to adapt it to the easier case of one-dimensional root-finding. Note that I here only work with gsl_vector and not numpy arrays. It is easy enough to copy a numpy array into a gsl_vector and it should even be possible to use a gsl_vector working directly on the data in a numpy array, although I've never tried this.

Finding roots of a function

1 Answers1