已知积分下限与积分值，求解9阶函数积分上限的实数值方法

我有一组x-y数据，拟合得到一个9阶多项式：25342𝑥^9 −155900𝑥^8+409218𝑥^7 − 599317𝑥^6 + 537190𝑥^5 − 303116𝑥^4 + ... + 274。用SymPy求解积分上限时，得到的10个解全是复数，但已知当下限为0.0925、上限约0.945时积分值约为100，求能获取该积分上限实数解的方法。

相关代码 #

p_optimal = estimate_function_from_data_points()

from sympy import integrate, solve
from sympy.abc import x, u
f = 25342.695344882944*x**9 - 155900.56387247072*x**8 + 409218.9290579793*x**7 - 599317.5264117827*x**6 + 537190.6784517929*x**5 - 303116.0648042093*x**4 + 105493.81468203208*x**3 - 20422.11374996729*x**2 + 1263.9293528900394*x + 274.55521542679185
lower = 0.0925 # Initial
upper = u
eq = integrate(f, (x, lower, upper))
eq, solve(eq - 100, u) # 注：原代码写eq+100，实际应为eq-100=0，对应积分值等于100

Out: [-0.114399781774514 - 0.112912224139529*I,
 -0.114399781774514 + 0.112912224139529*I,
 0.145632609024802 - 0.532284754794354*I,
 0.145632609024802 + 0.532284754794354*I,
 0.646926125977188 - 0.679233975801008*I,
 0.646926125977188 + 0.679233975801008*I,
 1.20499184950745 - 0.534200757552949*I,
 1.20499184950745 + 0.534200757552949*I,
 1.53445822404458 - 0.201823934360761*I,
 1.53445822404458 + 0.201823934360761*I]

补充验证代码 #

lower = 0.0925 # Initial 
upper = 0.945
integrate(f, (x, lower, upper)).evalf()

Out: 100.016292426307

可行的解决方法 #

1. 数值求根法（推荐） #

高次多项式的符号解通常会包含大量复数，直接用数值求根工具针对实数区间求解更高效。利用已知的近似值范围（0.9附近），可以用SciPy的root_scalar：

from scipy.optimize import root_scalar
import sympy as sp

x = sp.symbols('x')
f = 25342.695344882944*x**9 - 155900.56387247072*x**8 + 409218.9290579793*x**7 - 599317.5264117827*x**6 + 537190.6784517929*x**5 - 303116.0648042093*x**4 + 105493.81468203208*x**3 - 20422.11374996729*x**2 + 1263.9293528900394*x + 274.55521542679185
lower = 0.0925

# 定义目标函数：积分结果与100的差值为0
def objective(u):
    return sp.integrate(f, (x, lower, u)).evalf() - 100

# 用二分法，指定解所在的区间[0.8, 1.0]
result = root_scalar(objective, method='bisect', bracket=[0.8, 1.0])
print("精确实数解：", result.root)

运行后会得到精确的实数解，这种方法收敛速度快，且只关注我们需要的实数区间。

2. SymPy解筛选实数解 #

如果一定要用SymPy的符号解，可以筛选出虚部极小的解（因浮点数精度，真实实数解可能带极小虚部）：

from sympy import integrate, solve, N
from sympy.abc import x, u

f = 25342.695344882944*x**9 - 155900.56387247072*x**8 + 409218.9290579793*x**7 - 599317.5264117827*x**6 + 537190.6784517929*x**5 - 303116.0648042093*x**4 + 105493.81468203208*x**3 - 20422.11374996729*x**2 + 1263.9293528900394*x + 274.55521542679185
lower = 0.0925
eq = integrate(f, (x, lower, u)) - 100
solutions = solve(eq, u)

# 设置阈值筛选虚部接近0的解
real_solutions = [N(sol) for sol in solutions if abs(sol.imag) < 1e-6]
print("筛选出的实数解：", real_solutions)

这种方法适合需要保留符号计算流程的场景，但效率不如纯数值方法。

3. 手动计算不定积分后求根 #

先手动推导9阶多项式的不定积分（得到10阶多项式），代入下限得到常数项，转化为求P(u) = 100的实数根，再用数值求根工具求解。本质和方法1一致，但可以避免SymPy符号积分的开销，适合对性能要求较高的场景。

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