python求函数的根_求函数的根

我曾经遇到过一个类似的多维根查找问题。既然我找到了一个使用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数组中的数据，尽管我从未尝试过。在