scipy.optimize.fsolve 收敛错误？

Question

这是代码

import scipy as sc
import scipy.optimize as sco

def g(rho):
    return 0.5 * (rho**2 + rho) * sc.exp(-rho)

p = 0.01017036
guess = 1.5879245860401234
sol = sco.fsolve(lambda rho: g(rho) - p, guess, full_output=True)

print g(sol[0]) - p
print sol

输出结果是

[ 0.40970908]
(array([ 1.58792459]), {'qtf': array([-0.40970908]), 'nfev': 4, 'r': array([  9.52007670e+26]), 'fjac': array([[-1.]]), 'fvec': array([ 0.40970908])}, 1, 'The solution converged.')

它说已经收敛了，但显然没有，因为我应该得到的 g(sol[0]) - p 应该离零更近，而不是 0.4

我觉得用来检查收敛的测试有问题，导致出现错误。我知道我可以改变猜测的值来得到正确的结果，但这不是我想要的情况（我需要找到成千上万的根，而我得到了很多错误的根），问题在于这个错误捕捉的算法不可靠。我是不是做错了什么？？

提前谢谢你。

Answer 1

你现在是在最小化一个目标函数，而不是在找一个根，所以你应该使用 optimize.fmin 来做这个：

import scipy as sc
import scipy.optimize as sco

def g(rho):
    return 0.5 * (rho**2 + rho) * sc.exp(-rho)

p = 0.01017036
guess = 1.5879245860401234
sol = sco.fmin(lambda rho: (g(rho)-p)**2, guess)

print sol
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 23
         Function evaluations: 46
[ 8.2240227] 

这可能也跟你最开始的猜测不太好有关：

p = 0.01017036
guess = 1.5879245860401234
print sco.fsolve(f, 10, fprime=gp, args=(p,), full_output=1)
print sco.fsolve(lambda rho: g(rho) - p, 10, full_output=True)

print '\n==================bad_guess====================='
print sco.fsolve(f, guess, fprime=gp, args=(p,), full_output=1)
print sco.fsolve(lambda rho: g(rho) - p, guess, full_output=True)
#(array([ 8.22399478]), {'fvec': array([ -1.90472638e-15]), 'qtf': array([  1.87372265e-10]), 'nfev': 10, 'r': array([ 0.00783117]), 'fjac': array([[-1.]]), 'njev': 1}, 1, 'The solution converged.')
#(array([ 8.22399478]), {'qtf': array([  1.87372918e-10]), 'nfev': 11, 'r': array([ 0.00783117]), 'fjac': array([[-1.]]), 'fvec': array([ -1.90472638e-15])}, 1, 'The solution converged.')
#
#==================bad_guess=====================
#(array([ 1.58792459]), {'fvec': array([ 0.40970908]), 'qtf': array([-0.40970908]), 'nfev': 3, 'r': array([  9.52013062e+26]), 'fjac': array([[-1.]]), 'njev': 1}, 1, 'The solution converged.')
#(array([ 1.58792459]), {'qtf': array([-0.40970908]), 'nfev': 4, 'r': array([  9.52007670e+26]), 'fjac': array([[-1.]]), 'fvec': array([ 0.40970908])}, 1, 'The solution converged.')

可能收敛性是由雅可比矩阵决定的，而你的 guess 可能让雅可比矩阵落在了一个奇怪的地方。注意，rho=8.22399478 和 rho=1.58792459 的雅可比矩阵是一样的，而 1.58792459 正好是你最开始的猜测。看起来 fsolve 由于某种原因从来没有离开过这个初始猜测。