使用Sympy求多项式Y1在[0,1]区间最大值遇TypeError求助

我在IPython中使用Sympy时碰到了一个无法理解的错误。下述Y1是关于x的多项式，在区间[0,1]内存在最大值（注：图表的y轴是倒置的）。测试简单多项式时一切正常，比如maximum(x*(1-x), x, Interval(0, 1))能正确返回1/4，但对Y1执行maximum(Y1, x, Interval(0, 1))时触发了TypeError，具体信息如下：

In [8]: plot((Y1, interval1), (Y2, interval2), backend=InvertYAxis);

（图表显示Y1在[0,1]区间内有明确的极值）

In [9]: Y1
Out[9]: x**4/24 - 5*x**3/64 + 7*x/192

In [10]: maximum(Y1, x, Interval(0, 1))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 maximum(Y1, x, Interval(0, 1))

File /usr/lib/python3.11/site-packages/sympy/calculus/util.py:792, in maximum(f, symbol, domain)
    789     if domain is S.EmptySet:
    790         raise ValueError("Maximum value not defined for empty domain.")
--> 792     return function_range(f, symbol, domain).sup
    793 else:
    794     raise ValueError("%s is not a valid symbol." % symbol)

File /usr/lib/python3.11/site-packages/sympy/calculus/util.py:208, in function_range(f, symbol, domain)
    203     raise NotImplementedError(
    204             'Infinite number of critical points for {}'.format(f))
    206 critical_points += solution
--> 208 for critical_point in critical_points:
    209     vals += FiniteSet(f.subs(symbol, critical_point))
    211 left_open, right_open = False, False

File /usr/lib/python3.11/site-packages/sympy/sets/sets.py:1561, in Intersection.__iter__(self)
   1559 if not candidates:
   1560     raise TypeError("None of the constituent sets are iterable")
--> 1561 raise TypeError(
   1562     "The computation had not completed because of the "
   1563     "undecidable set membership is found in every candidates.")

TypeError: The computation had not completed because of the undecidable set membership is found in every candidates.

问题原因 #

这个错误是因为Sympy在求解Y1的临界点时，得到的临界点集合和定义域的交集无法被自动判定成员关系，导致function_range函数无法遍历临界点，进而抛出异常。

解决方案 #

换一种思路，手动求导找临界点，再计算临界点和区间端点的函数值，取最大值即可。具体步骤如下：

• 对Y1求导，得到导数表达式
• 求解导数等于0的方程，得到所有临界点
• 筛选出落在区间[0,1]内的临界点
• 计算这些临界点和区间端点（x=0、x=1）处的Y1值
• 比较所有值，取最大的那个

代码实现 #

from sympy import symbols, diff, solve, Interval, N

x = symbols('x')
Y1 = x**4/24 - 5*x**3/64 + 7*x/192

# 1. 求导
Y1_prime = diff(Y1, x)
# 2. 解导数为0的方程，得到临界点
critical_points = solve(Y1_prime, x)
# 3. 筛选区间[0,1]内的有效临界点
valid_cps = [cp for cp in critical_points if Interval(0, 1).contains(cp)]
# 4. 加入区间端点
all_eval_points = valid_cps + [0, 1]
# 5. 计算每个点的函数值
function_values = [Y1.subs(x, point) for point in all_eval_points]
# 6. 取最大值
max_value = max(function_values)

# 输出精确值和近似数值
print("精确最大值:", max_value)
print("近似数值:", N(max_value))

运行这段代码会得到精确结果和近似数值，和图表观察到的极值一致。

如果需要更高效的数值计算，也可以把Sympy表达式转成数值函数后用Scipy的优化工具求解，但手动求导的方法更直接，依赖更少。

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