用scipy求多元函数的根

Question

我在寻找一个函数的根，我知道它有根。在根代码下面的两个细节中找不到scipy.optimize.root. 我的问题涉及到如何找到根（我已经尝试过更改起始点、方法、maxiter和xtol/ftol，但是由于我是python新手，所以我可能没有足够的效率来让程序工作）。在

第一种情况：根在未指定jacobian时找到，但在指定时不存在。雅各布怎么了？ 另外，虽然很明显找不到根，但是算法返回success，为什么以及如何修复它？在

第二种情况：根永远找不到，我该如何修复它？在

from __future__ import division
import numpy as np
import datetime
import pandas 
import scipy
from scipy import optimize

# 1st case
# rap=np.array([.5,12.8,44,10.7,5,12.3])
# 2nd case
rap=np.array([1.2,83,5,41,7,50])

def f2(wc, a1, p0, a1pr, p0pr, atheta, x):
    A = (p0/(1-p0))**(1-a1)
    Apr = (p0pr/(1-p0pr))**(1-a1pr)
    return wc - (((A * x**a1)/(((1-x)**a1) + A * x**a1))*  (1-np.exp(-rap*atheta))/atheta) - (((Apr * (1-x)**a1pr)/((x**a1pr) + Apr * (1-x)**a1pr))*((1-np.exp(atheta))/atheta))

# define the derivative of f2 with respect to x (checked with sympy)
def f2p(wc, a1, p0, a1pr, p0pr, atheta, x):
    A = (p0/(1-p0))**(1-a1)
    Apr = (p0pr/(1-p0pr))**(1-a1pr)
    return - (((A*a1 * x**(a1-1) * (1-x)**(a1-1)) / (((1-x)**a1 + A * x**a1)**2))  \
        *(1-np.exp(-rap*atheta)/atheta)) \
         + (((Apr * a1pr * (1-x)**(a1pr-1) * x**(a1pr-1) ) / ((x**a1pr + Apr * (1-x)**a1pr)**2))  \
         *((1-np.exp(atheta))/atheta))

# f1 defines the function to find the root of
def f1(x, y):
    a1, p0, a1pr, p0pr, atheta = y[:5]
    tosolve = np.full_like(x, x[1:].sum()-1)
    tosolve[:-1] = f2(x[0], a1, p0, a1pr, p0pr, atheta, x[1:])
    return tosolve

# define the jacobian of f1 
def jac(x, y):
    a1, p0, a1pr, p0pr, atheta = y[:5]
    partial = f2p(x[0], a1, p0, a1pr, p0pr, atheta, x[1:])
    matrix = np.diag(partial)
    vector1 = np.full_like(rap,1)
    matrix1 = np.vstack((matrix, vector1))
    vector = np.full_like(x,1)
    vector[-1]=0
    matrix2 = np.column_stack((vector, matrix1))
    return matrix2

# f is finding the roots of f1
def f(y):
    a1, p0, a1pr, p0pr, atheta = y[:5] 
    x0 = [0.000001,0.000001,0.000001,0.000001,0.000001,0.000001,0.000001] 
    print 'x0', x0
    sol = scipy.optimize.root(f1, x0, args=(y), method='lm')
#    sol = scipy.optimize.root(f1, x0, args=(y), method='lm', jac=jac,options={'maxiter': 5000})
    return sol

# 1st case   
# result = f([.8,0.03,0.8,0.3,-0.01])
# 2nd case
result = f([1.99,0.01,1.99,0.989,-0.43])
print result