使用scipy fsolve（遇到数学域错误）求解非线性方程组

Question

我试图用Scipy的fsolve找到两个非线性方程组的答案。这两个方程是：

f1 = math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2
f2 = math.log(1 - x) - (m - 1)*x + chi*m*x**2

在这种情况下，m和chi是常量。基本目标是找到同时满足f1（x）=f1（y）和f2（x）=f2（y）的x，y。我知道x，y的初始猜测是0.3和0.99。下面是我的代码。

from scipy.optimize import fsolve
import math

# some global variables
m = 46.663
chi = 1.1500799949128826

def binodal_fsolve():
    def equations(p):
        x, y = p
        out = []
        out.append(math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2 - (math.log(y) + 1. - ((1. + (m - 1)*y) / m) + chi * (1 - y)**2))
        out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))

        return out

    initial_guess = [0.3, 0.99]
    ans = fsolve(equations, initial_guess)

    return ans

def test_answers(phiL, phiR):
    def functions(x):
        return math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2, math.log(1 - x) - (m - 1)*x + chi*m*x**2

    return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]

print (test_answers(0.2542983070, 0.9999999274))
# (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
res = binodal_fsolve()
print (res)

当我执行代码时，我总是遇到数学域错误。但是，如果我尝试使用MAPLEfsolve解决它。我可以得到答案（0.2542983070，0.999999 274）。

通过将这些插入方程，我得到（1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)，这表明答案是正确的。我不知道如何让西皮解决问题。任何建议将不胜感激。

Answer 1

在这种情况下，您可以使用numpy中的log功能。自由党。scimath当其参数为负数时返回复数。

代替使用scipy.optimize.fsolve，使用scipy.optimize.root并将该方法更改为lm，该方法使用Levenberg-Maraqudt的修改在最小二乘意义上求解非线性方程组算法。有关更多方法，请参阅留档。

from scipy.optimize import root

import numpy.lib.scimath as math

# some global variables
m = 46.663
chi = 1.1500799949128826

def binodal_fsolve():
    def equations(p):
        x, y = p
        out = []
        out.append(math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2 - (math.log(y) + 1. - ((1. + (m - 1)*y) / m) + chi * (1 - y)**2))
        out.append(math.log(1 - x) - (m - 1)*x + chi*m*x**2 - (math.log(1 - y) - (m - 1)*y + chi*m*y**2))

        return out

    initial_guess = [0.3, 0.99]
    #ans = fsolve(equations, initial_guess)
    ans = root(equations, initial_guess, method='lm')

    return ans

def test_answers(phiL, phiR):
    def functions(x):
        return math.log(x) + 1. - ((1. + (m - 1)*x) / m) + chi * (1 - x)**2, math.log(1 - x) - (m - 1)*x + chi*m*x**2

    return functions(phiL)[0], functions(phiR)[0], functions(phiL)[1], functions(phiR)[1]

print (test_answers(0.2542983070, 0.9999999274))
# (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
res = binodal_fsolve()
print (res)

其中给出了以下根x和y：：数组（[0.25429812, 0.99999993]）。

完整输出：

(1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684)
/home/user/.local/lib/python3.6/site-packages/scipy/optimize/minpack.py:401: ComplexWarning: Casting complex values to real discards the imaginary part
  gtol, maxfev, epsfcn, factor, diag)
   cov_x: array([[6.49303571e-01, 8.37627537e-07],
       [8.37627537e-07, 1.08484856e-12]])
    fjac: array([[ 1.52933340e+07, -1.00000000e+00],
       [-1.97290115e+01, -1.24101235e+00]])
     fun: array([-2.22945317e-07, -7.20367503e-04])
    ipvt: array([2, 1], dtype=int32)
 message: 'The relative error between two consecutive iterates is at most 0.000000'
    nfev: 84
     qtf: array([-0.00338589,  0.00022828])
  status: 2
 success: True
       x: array([0.25429812, 0.99999993])