sympy.solve对特定系数无返回的问题排查与解决请求

尝试求解一个形式固定、仅指数变化的表达式的临界点，最小复现代码如下：

from sympy import symbols, Mul, diff, solve
from decimal import Decimal


class Histogram:
    def __init__(self, bins, pdf):
        self.bins = bins
        self.pdf = pdf


histogram = Histogram(
    bins=[-1, 0, 1],
    pdf=[Decimal('0.454861111111111104943205418749130330979824066162109375'),
         Decimal('0'),
         Decimal('0.545138888888888839545643349993042647838592529296875')]
)


def computation(histogram):

    A = Decimal('1')
    e = Decimal('1')
    M_factor = Decimal('1.0828567056280801')
    x = symbols('x')
    product_terms = []
    profit = lambda cp: (e * M_factor ** cp - e) * x
    commission = lambda cp: (e * M_factor ** cp + e) * Decimal('0.0005') * x

    for candles_profit, freq in zip(histogram.bins, histogram.pdf):
        profit_term = profit(candles_profit)
        commission_term = commission(candles_profit)
        term = ((A + profit_term - commission_term) / A) ** freq
        product_terms.append(term)

    expression = Mul(*product_terms)
    derivative = diff(expression, x)
    solutions = solve(derivative, x)
    print("Solutions:", solutions)


computation(histogram)

调试中得到的表达式为：

(1.0 - 0.0774785191289289x)^0.454861111111111 * (0.0818152772752661x + 1.0)^0.545138888888889

运行代码时sympy.solve始终无法返回结果。经排查，指数实际是循环小数形式的有理数（131/288、157/288），但SymPy将输入的Decimal转为Float而非Rational，推测这是问题根源。请问该推测是否合理？如何修改代码以正常完成计算？

解答 #

推测合理性 #

你的推测完全合理。SymPy对浮点数的符号计算支持有限，浮点数本质是近似值，会引入精度误差，导致符号求解器无法识别其有理数本质，进而无法完成精确的代数运算，最终导致求解失败。而有理数（Rational类型）是精确的，SymPy能对其进行完整的符号推导和求解。

修改方案 #

核心是将所有Decimal类型的数值转换为SymPy的Rational类型，确保所有运算都基于精确的有理数进行：

1. 替换decimal.Decimal为sympy.Rational；
2. 对于无法直接用分数表示的小数（如M_factor），可以用Rational的字符串构造方法，或通过nsimplify将浮点数转换为最接近的有理数；
3. 移除不必要的Decimal依赖，全程使用SymPy的精确数值类型。

修改后的代码如下：

from sympy import symbols, Mul, diff, solve, Rational, nsimplify


class Histogram:
    def __init__(self, bins, pdf):
        self.bins = bins
        self.pdf = pdf


# 将pdf的Decimal字符串转为Rational
histogram = Histogram(
    bins=[-1, 0, 1],
    pdf=[Rational('0.454861111111111104943205418749130330979824066162109375'),
         Rational('0'),
         Rational('0.545138888888888839545643349993042647838592529296875')]
)


def computation(histogram):
    A = Rational(1)
    e = Rational(1)
    # 使用nsimplify将小数转为接近的有理数，或直接用字符串构造
    M_factor = nsimplify('1.0828567056280801')
    x = symbols('x')
    product_terms = []
    
    # 替换lambda为普通函数，避免闭包问题（可选，但更清晰）
    def profit(cp):
        return (e * M_factor ** cp - e) * x
    
    def commission(cp):
        return (e * M_factor ** cp + e) * Rational('0.0005') * x

    for candles_profit, freq in zip(histogram.bins, histogram.pdf):
        if freq == 0:
            continue  # 指数为0的项等于1，不影响乘积，直接跳过
        profit_term = profit(candles_profit)
        commission_term = commission(candles_profit)
        term = ((A + profit_term - commission_term) / A) ** freq
        product_terms.append(term)

    expression = Mul(*product_terms)
    derivative = diff(expression, x)
    solutions = solve(derivative, x)
    print("Solutions:", solutions)


computation(histogram)

关键说明 #

• Rational('字符串形式的小数')可以精确将小数转换为分数，例如你的pdf值会被自动识别为131/288和157/288；
• nsimplify可以将近似小数转换为最接近的有理数，适合处理M_factor这类无法直接看出分数的数值；
• 跳过指数为0的项（如pdf中的0值），因为任何数的0次方都是1，不会改变乘积结果，减少不必要的计算。

运行修改后的代码，SymPy就能正确识别表达式的有理数结构，完成符号求导并求解出临界点。

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