решить систему нелинейных уравнений, используя scipy fsolve (обнаружена ошибка математической области)

Question

Я пытался использовать 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 являются константами в этом случае. Основная цель - найти x, y, которые удовлетворяют одновременно f1(x) = f1(y) и f2(x) = f2(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)

Когда я выполнял код, я всегда сталкивался с ошибкой математического домена. Однако, если бы я попытался решить ее, используя MAPLE fsolve. Я могу получить ответы (0.2542983070, 0.9999999274).

Подключив их обратно к уравнениям, я получу (1.3598772108380786e-09, -1.5558330624053502e-09, -8.434988430355375, -8.435122589529684), что говорит о правильности ответов. Я не знаю, как заставить scipy fsolve работать. Любые предложения будут с благодарностью.

Answer 1

В этом случае вы можете использовать функцию log из numpy.lib.scimath, которая возвращает комплексное число, когда его аргумент отрицательный.

Вместо использования scipy.optimize.fsolve используйте scipy.optimize.root и измените метод на lm, который решает систему нелинейных уравнений в смысле наименьших квадратов, используя модификацию алгоритма Левенберга-Марквардта. Дополнительные методы см. В документации .

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: : array([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])