如何用Python求解calculate_density(x)方程中的x值？fsolve结果不佳

问题背景 #

你有如下密度计算函数：

a = 5108.30
b = 3967550
c = -281375000
pressure = 100
def calculate_density(x):
    density = (a * x + b * x**2 + c * x**3) * pressure
    return density

需求是已知目标密度值，反向求解方程中的x，但使用fsolve得到的结果并非最优解。

核心分析 #

你的问题本质是求解三次代数方程，先整理为标准形式：
c * pressure * x³ + b * pressure * x² + a * pressure * x - target_density = 0
三次方程最多有3个实根，fsolve这类迭代法依赖初始猜测，容易收敛到局部解，需针对性调整解法。

优化解法 #

方法1：直接求解所有实根后筛选最优 #

利用三次方程的解析特性，先求出所有实根，再根据物理意义筛选符合要求的解：

import numpy as np
from numpy.polynomial import Polynomial

target_density = 123456  # 替换为你的目标密度值

# 构造三次多项式系数：[常数项, x项系数, x²项系数, x³项系数]
coeffs = [
    -target_density,
    a * pressure,
    b * pressure,
    c * pressure
]

# 生成多项式并求解所有根
poly = Polynomial(coeffs)
roots = poly.roots()

# 筛选实根（排除虚部可忽略的复数）
real_roots = [root.real for root in roots if np.isclose(root.imag, 0)]

# 根据业务逻辑筛选有效解（比如x>0，且计算密度与目标误差符合要求）
valid_roots = []
for x in real_roots:
    if x > 0:  # 根据实际场景调整约束
        calc_density = calculate_density(x)
        if np.isclose(calc_density, target_density, rtol=1e-5):
            valid_roots.append(x)

# 获取最优解
x_best = valid_roots[0] if valid_roots else None

方法2：优化fsolve的使用逻辑 #

如果坚持用迭代法，通过调整初始值和求解逻辑提升稳定性：

from scipy.optimize import fsolve, root_scalar

target_density = 123456

# 明确方程：返回计算密度与目标密度的差值
def equation(x, a, b, c, pressure, target_density):
    return (a * x + b * x**2 + c * x**3) * pressure - target_density

# 1. 优化初始猜测值：先通过试值找到密度接近目标的x区间，再选初始点
x_initial_guess = 0.05  # 替换为更接近真实解的初始值
x_solution = fsolve(equation, x_initial_guess, args=(a, b, c, pressure, target_density))

# 验证结果
calc_density = calculate_density(x_solution[0])
print(f"求解x: {x_solution[0]}, 计算密度: {calc_density}, 目标密度: {target_density}")

# 2. 用区间约束求解（更稳定）：已知x范围时用root_scalar的brentq方法
result = root_scalar(equation, args=(a, b, c, pressure, target_density), 
                     method='brentq', bracket=[0, 0.1])  # 替换为实际区间
if result.converged:
    x_best = result.root
else:
    print("求解未收敛")

关键注意事项 #

• 添加物理约束：根据实际场景过滤无意义的根（比如x不能为负）。
• 验证误差：所有求解结果必须代入原函数验证，确保计算密度与目标值的误差在允许范围内。

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