求函数的根

Question

我想把一段代码从python改写成cython，到目前为止，我已经把所有我在这个例子中简化的部分都用cython处理了，所以我没法保持这个函数的python样子。不过，我需要估算这个函数的根，而之前是用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，而不是numpy数组。把numpy数组复制到gsl_vector里是很简单的，甚至可以直接用gsl_vector在numpy数组的数据上操作，尽管我自己从来没有尝试过。