Нахождение корней функции

Question

Я хотел бы переписать код с python на cython, и до сих пор я цитировал все части, которые я упростил в этом примере, не используя их. Поэтому я не мог сохранить форму питона этой функции. Тем не менее, мне нужно оценить корень функции, которая перед scipy.optimize библиотека была использована. Мне интересно, что я мог бы заменить, чтобы найти корни этой функции в Cython. Мог gsl также предоставить лучший инструмент, чтобы найти корни? Как это должно быть сделано?

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

Answer 1

Однажды я столкнулся с подобной проблемой для многомерного поиска корней. Поскольку я нашел решение с использованием GSL, я хотел бы поделиться своим кодом здесь. pyx файл, который я использовал, был

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)

Он должен содержать необходимые понятия и достаточно легко адаптировать его к более простому случаю одномерного поиска корней. Обратите внимание, что я здесь только работаю с gsl_vector а не numy массивы. Достаточно просто скопировать пустой массив в gsl_vector и даже должно быть возможно использовать gsl_vector работая непосредственно с данными в массиве Numpy, хотя я никогда не пробовал это.