使用 SciPy 求多维标量函数根的最佳方法

Question

假设我有一个函数，其范围是标量，但其域是向量。例如：

def func(x):
  return x[0] + 1 + x[1]**2

找到这个函数的根 的好方法是什么？scipy.optimize.fsolve并scipy.optimize.root期望func返回一个向量（而不是标量），并且scipy.optimize.newton只接受标量参数。我可以重新定义func为

def func(x):
  return [x[0] + 1 + x[1]**2, 0]

然后root和fsolve可以找到根，但是雅可比行列式中的零意味着它并不总是能很好地完成工作。例如：

fsolve(func, array([0,2]))
=> array([-5,  2])

它只会改变第一个参数，但不会改变第二个参数，这意味着它经常会找到很远的零。

---

编辑：看起来下面的 func 重新定义效果更好：

def func(x):
  fx = x[0] + 1 + x[1]**2
  return [fx, fx]

fsolve(func, array([0,5]))
=>array([-16.27342781,   3.90812331])

所以它现在愿意改变这两个参数。但代码仍然有点难看。

Answer 1

因为——对于我的问题——我有一个很好的初步猜测和一个非疯狂的函数，牛顿的方法效果很好。对于标量、多维函数，牛顿法变为：

这是一个粗略的代码示例：

def func(x): #the function to find a root of
  return x[0] + 1 + x[1]**2

def dfunc(x): #the gradient of that function
  return array([1, 2*x[1]])

def newtRoot(x0, func, dfunc):
  x = array(x0)
  for n in xrange(100): # do at most 100 iterations
    f  = func(x)
    df = dfunc(x)

    if abs(f) < 1e-6: # exit function if we're close enough
      break

    x = x - df*f/norm(df)**2 # update guess
  return x

正在使用：

nsolve([0,2],func,dfunc)
=> array([-1.0052546 ,  0.07248865])

func([-1.0052546 ,  0.07248865])
=> 4.3788225025098715e-09

不错！当然，这个函数很粗糙，但是你明白了。它也不适用于“棘手”的函数或您没有良好的起始猜测的情况。我想我会使用类似的方法，但如果牛顿方法不收敛，那么我会回退到fsolveor 。root