python之使用 scipy.optimize.root 查找根

我正在尝试使用 Python 查找名为 f 的函数的根 y。

这是我的代码:

 def f(y): 
    w,p1,p2,p3,p4,p5,p6,p7 = y[:8]  
    t1 = w - 0.500371726*(p1**0.92894164) - (-0.998515304)*((1-p1)**1.1376649) 
    t2 = w - 8.095873128*(p2**0.92894164) - (-0.998515304)*((1-p2)**1.1376649) 
    t3 = w - 220.2054377*(p3**0.92894164) - (-0.998515304)*((1-p3)**1.1376649) 
    t4 = w - 12.52760758*(p4**0.92894164) - (-0.998515304)*((1-p4)**1.1376649) 
    t5 = w - 8.710859537*(p5**0.92894164) - (-0.998515304)*((1-p5)**1.1376649) 
    t6 = w - 36.66350261*(p6**0.92894164) - (-0.998515304)*((1-p6)**1.1376649) 
    t7 = w - 3.922692207*(p7**0.92894164) - (-0.998515304)*((1-p7)**1.1376649)        
    t8 = p1 + p2 + p3 + p4 + p5 + p6 + p7 - 1 
    return [t1,t2,t3,t4,t5,t6,t7,t8] 
 
 
x0 = np.array([-0.01,0.3,0.1,0.2,0.1,0.1,0.1,0.1]) 
sol = scipy.optimize.root(f, x0, method='lm') 
print sol  
print 'solution', sol.x 
print 'success', sol.success 

无论我在 scipy.optimize.root 中尝试什么方法，Python 都找不到根。

但是有一个，我在Matlab中用函数fsolve找到了它。

它是:

[-0.0622、0.5855、0.087、0.0028、0.0568、0.0811、0.0188、0.1679]。

当我指定 x0 接近根时，python 算法收敛。问题是我不知道根上的先验来指定 x0。实际上，我正在求解许多此类方程式。

我真的很想用 Python。任何人都可以帮助我与 python 融合吗？

请您参考如下方法:

好吧，经过一些闲逛，我们关注好的优化/寻根算法的另一个方面。在上面的评论中，我们反复讨论了使用 scipy.optimize.root() 中的哪个方法。对于近乎无懈可击的“自动”根查找，一个同样重要的问题是确定良好的初始猜测。很多时候，良好的初始猜测实际上根本不接近真实答案。相反，需要对它们进行安排，以便它们自然地引导求解器朝正确的方向发展。

在您的特定情况下，您的猜测实际上是在向奇怪的方向发送算法。

我对你问题的玩具重构是:

import numpy as np 
import scipy as sp 
import scipy.optimize 
def f(y): 
    w,p1,p2,p3,p4,p5,p6,p7 = y[:8] 
    def t(p,w,a):        
        b = -0.998515304 
        e1 = 0.92894164 
        e2 = 1.1376649 
        return(w-a*p**e1 - b*(1-p)**e2) 
    t1 = t(p1,w,1.0) 
    t2 = t(p2,w,4.0)  
    t3 = t(p3,w,16.0) 
    t4 = t(p4,w,64.0) 
    t5 = t(p5,w,256.0) 
    t6 = t(p6,w,512.0) 
    t7 = t(p7,w,1024.0)        
    t8 = p1 + p2 + p3 + p4 + p5 + p6 + p7 - 1.0 
    return(t1,t2,t3,t4,t5,t6,t7,t8) 
guess = 0.0001 
x0 = np.array([-1000.0,guess,guess,guess,guess,guess,guess,guess]) 
sol = sp.optimize.root(f, x0, method='lm') 
print('w=-1000: ', sol.x, sol.success,sol.nfev,np.sum(f(sol.x))) 

请注意，我没有使用您的特定预因子(我想扩大我探索的范围)，尽管我将您的特定指数保留在 p 项上。

真正的 secret 在于最初的猜测，我对所有 p 项进行了相同的猜测。大部分时间都将其设置为 0.1 或更高，因为有些术语希望以一种方式进行，而另一些则进行其他方式。将它降低到 0.01 可以很好地解决这个问题。 (我会注意到 w 项非常稳健——从 -1000 变化到 +1000。对解决方案没有影响)。进一步减少初始猜测对这个特定问题没有影响，但也没有坏处。我会保持它很小。

是的，您知道至少某些项会大得多。但是，您将求解器置于可以清晰直接地朝着真正的解决方案前进的位置。

祝你好运。