Python中使用fsolve求解含对数二元非线性方程报错问题

现有两个包含变量a、b的非线性方程，方程中对数项通过sympy.log()生成，使用fsolve对该方程组求解、设置初始迭代猜测值为(1,1)时触发报错：提示无法将SymPy表达式转换为浮点数，函数调用返回结果不符合浮点数数组格式要求。需确认该方程组是否具备可解性，以及对应报错的解决方法。

方程1 #

-a*(25570.2642 - 48973.016*a) - a*(48973.016*a - 23402.7518) + b*(25570.2642 - 48973.016*b) + b*(48973.016*b - 23402.7518) + (1 - a)*(25570.2642 - 48973.016*a) + (1 - a)*(48973.016*a - 23402.7518) - (1 - b)*(25570.2642 - 48973.016*b) - (1 - b)*(48973.016*b - 23402.7518) + 4731.31823640719 - (b/(1 - b))**(-9600*b)*((b/(1 - b))**(9600*b)*(1 - b)**9600*(9600*(1 - b)*(b/(1 - b)**2 + 1/(1 - b)) + 9600*log(b/(1 - b))) - 9600*(b/(1 - b))**(9600*b)*(1 - b)**9599)/(1 - b)**9600 + (a/(1 - a))**(-9600*a)*((a/(1 - a))**(9600*a)*(1 - a)**9600*(9600*(1 - a)*(a/(1 - a)**2 + 1/(1 - a)) + 9600*log(a/(1 - a))) - 9600*(a/(1 - a))**(9600*a)*(1 - a)**9599)/(1 - a)**9600

方程2 #

-a*(25570.2642 - 48973.016*a) - a*(48973.016*a - 23402.7518) + (1 - a)*(25570.2642 - 48973.016*a) + (1 - a)*(48973.016*a - 23402.7518) + 10382.138477723 - (a*(1 - a)*(25570.2642 - 48973.016*a) + a*(1 - a)*(48973.016*a - 23402.7518) + 10382.138477723*a - b*(1 - b)*(25570.2642 - 48973.016*b) - b*(1 - b)*(48973.016*b - 23402.7518) - 5650.82024131583*b + log((a/(1 - a))**(9600*a)*(1 - a)**9600) - log((b/(1 - b))**(9600*b)*(1 - b)**9600) + 320.132207782401)/(a - b) + (a/(1 - a))**(-9600*a)*((a/(1 - a))**(9600*a)*(1 - a)**9600*(9600*(1 - a)*(a/(1 - a)**2 + 1/(1 - a)) + 9600*log(a/(1 - a))) - 9600*(a/(1 - a))**(9600*a)*(1 - a)**9599)/(1 - a)**9600

def equations(p):
    x, y = p
    return (eq1, eq2)

x, y =  fsolve(equations, (1, 1))

TypeError                                 Traceback (most recent call last)
/usr/local/lib/python3.7/dist-packages/sympy/core/expr.py in __float__(self)
    348         if result.is_number and result.as_real_imag()[1]:
    349             raise TypeError("can't convert complex to float")
--> 350         raise TypeError("can't convert expression to float")
    351 
    352     def __complex__(self):

TypeError: can't convert expression to float
---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-148-1d80da78acbb> in <module>()
     91     return (eq1, eq2)
     92 
---> 93 x, y =  fsolve(equations, (1, 1))

1 frames
/usr/local/lib/python3.7/dist-packages/scipy/optimize/minpack.py in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
    223             maxfev = 200 * (n + 1)
    224         retval = _minpack._hybrd(func, x0, args, 1, xtol, maxfev,
--> 225                                  ml, mu, epsfcn, factor, diag)
    226     else:
    227         _check_func('fsolve', 'fprime', Dfun, x0, args, n, (n, n))

error: Result from function call is not a proper array of floats.

• 核心问题：fsolve是纯数值求解函数，要求传入的目标函数返回浮点数类型的计算结果。原代码中eq1、eq2是SymPy符号对象，没有将迭代过程中传入的x、y值代入表达式做数值计算，直接返回符号表达式，自然无法被转换为浮点数。
• 初始值不合法：方程中存在1-a、1-b作为分母，且包含log(a/(1-a))、log(b/(1-b))项，有效定义域为0 < a < 1、0 < b < 1，原代码设置的初始值(1,1)刚好在定义域边界，代入后会出现除零、对数真数无意义的计算错误，即便解决符号转换问题，该初始值也无法正常迭代。
• 表达式冗余度高：方程中存在大量可约去的高次幂项（最高次达9600次），直接计算极易出现数值溢出，导致迭代失败。

该方程组为定义域内的连续非线性方程组，只要(0,1)区间内存在满足两个方程残差为0的点，就可以通过数值迭代方法求解。但受高次幂项带来的数值刚性影响，迭代收敛性高度依赖初始值选择，建议先化简表达式降低计算难度，再选择合理初始值求解。

1. 先对SymPy表达式做化简，约去冗余的高次幂项，降低数值溢出风险，可直接调用sympy.simplify()完成化简。
2. 使用sympy.lambdify()将符号表达式转换为兼容NumPy的数值计算函数，替换原代码中直接返回符号对象的逻辑。
3. 调整初始迭代值到(0,1)开区间内（比如(0.5, 0.5)），同时在目标函数中增加定义域判断，避免迭代过程中出现无意义的计算。

修正后的参考代码：

import numpy as np
from scipy.optimize import fsolve
import sympy as sp

# 定义符号变量
a, b = sp.symbols('a b', real=True)

# 代入原方程表达式后先做化简
eq1 = sp.simplify(
    -a*(25570.2642 - 48973.016*a) - a*(48973.016*a - 23402.7518) + b*(25570.2642 - 48973.016*b) + b*(48973.016*b - 23402.7518) + (1 - a)*(25570.2642 - 48973.016*a) + (1 - a)*(48973.016*a - 23402.7518) - (1 - b)*(25570.2642 - 48973.016*b) - (1 - b)*(48973.016*b - 23402.7518) + 4731.31823640719 - (b/(1 - b))**(-9600*b)*((b/(1 - b))**(9600*b)*(1 - b)**9600*(9600*(1 - b)*(b/(1 - b)**2 + 1/(1 - b)) + 9600*sp.log(b/(1 - b))) - 9600*(b/(1 - b))**(9600*b)*(1 - b)**9599)/(1 - b)**9600 + (a/(1 - a))**(-9600*a)*((a/(1 - a))**(9600*a)*(1 - a)**9600*(9600*(1 - a)*(a/(1 - a)**2 + 1/(1 - a)) + 9600*sp.log(a/(1 - a))) - 9600*(a/(1 - a))**(9600*a)*(1 - a)**9599)/(1 - a)**9600
)
eq2 = sp.simplify(
    -a*(25570.2642 - 48973.016*a) - a*(48973.016*a - 23402.7518) + (1 - a)*(25570.2642 - 48973.016*a) + (1 - a)*(48973.016*a - 23402.7518) + 10382.138477723 - (a*(1 - a)*(25570.2642 - 48973.016*a) + a*(1 - a)*(48973.016*a - 23402.7518) + 10382.138477723*a - b*(1 - b)*(25570.2642 - 48973.016*b) - b*(1 - b)*(48973.016*b - 23402.7518) - 5650.82024131583*b + sp.log((a/(1 - a))**(9600*a)*(1 - a)**9600) - sp.log((b/(1 - b))**(9600*b)*(1 - b)**9600) + 320.132207782401)/(a - b) + (a/(1 - a))**(-9600*a)*((a/(1 - a))**(9600*a)*(1 - a)**9600*(9600*(1 - a)*(a/(1 - a)**2 + 1/(1 - a)) + 9600*sp.log(a/(1 - a))) - 9600*(a/(1 - a))**(9600*a)*(1 - a)**9599)/(1 - a)**9600
)

# 转换为数值计算函数
calc_eq1 = sp.lambdify((a, b), eq1, 'numpy')
calc_eq2 = sp.lambdify((a, b), eq2, 'numpy')

def equations(p):
    x, y = p
    # 定义域校验，超出范围返回大残差引导迭代
    if x <= 1e-6 or x >= 1-1e-6 or y <= 1e-6 or y >= 1-1e-6 or abs(x-y) <= 1e-6:
        return [1e12, 1e12]
    return [calc_eq1(x, y), calc_eq2(x, y)]

# 选择定义域内的初始值
sol = fsolve(equations, (0.3, 0.7))
print(f"求解结果：a={sol[0]:.6f}, b={sol[1]:.6f}")
print(f"残差：eq1={calc_eq1(sol[0], sol[1]):.6e}, eq2={calc_eq2(sol[0], sol[1]):.6e}")

内容的提问来源于stack exchange，提问作者Sonic Sharma