含根式的不等式解值求解与Matplotlib绘图问题求助

问题背景 #

你需要绘制不等式F1(γ)≥0（γ>1.11）的解区域，已通过Sympy构建符号函数F1，但使用scipy.optimize.minimize在特定参数（W=0.003, c=0.2, r=0.11）下无法得到有效解，希望获取该不等式的γ解值并通过Matplotlib实现类似Desmos的绘图效果。

问题分析 #

minimize通过最小化函数绝对值找根的方式，在该参数下因函数特性导致优化器无法收敛。可换用更针对性的单变量根查找工具，或尝试符号解析求解，再结合数值采样绘制解区域。

方法1：使用Scipy根查找工具获取精确根 #

scipy.optimize.root_scalar专门用于单变量函数的根查找，支持二分法等稳定算法，只要找到函数符号变化的区间就能收敛。

import sympy
from sympy import *
import numpy as np
from scipy.optimize import root_scalar
import matplotlib.pyplot as plt

# 1. 构建并简化符号函数
r, γ, W, c = symbols('r γ W c', positive=True, real=True)
D = 2*γ - c - sqrt((2*γ - c)**2 - 4*γ*(1+r)*(1-W))
EU_e = γ - D + D**2/(4*γ)
F = EU_e - (1+r)*W
F1 = F.subs([(W, 0.003), (c, 0.2), (r, 0.11)])

# 2. 转换为数值计算函数
ff = lambdify(γ, F1, modules=['numpy'])

# 3. 确定根的搜索区间（先采样观察函数符号变化）
gamma_test = np.linspace(1.11, 2, 50)
f_vals = ff(gamma_test)
sign_changes = np.where(np.diff(np.sign(f_vals)))[0]

root_gamma = None
if len(sign_changes) > 0:
    # 取第一个符号变化的区间用二分法找根
    a, b = gamma_test[sign_changes[0]], gamma_test[sign_changes[0]+1]
    res = root_scalar(ff, bracket=[a, b], method='bisect')
    root_gamma = res.root
    print(f"找到根：γ = {root_gamma:.6f}")
else:
    # 检查区间内函数是否恒正/恒负
    if np.all(f_vals >= 0):
        print("γ>1.11时，F1(γ)≥0恒成立")
        root_gamma = 1.11
    else:
        print("γ>1.11时，F1(γ)≥0无解")

# 4. 绘制解区域
gamma_range = np.linspace(1.11, 3, 200)
f_vals = ff(gamma_range)

plt.figure(figsize=(8,5))
plt.plot(gamma_range, f_vals, label='$F_1(γ)$')
plt.axhline(y=0, color='r', linestyle='--', label='$F_1(γ)=0$')

if root_gamma is not None:
    plt.scatter(root_gamma, 0, color='g', s=50, zorder=5, label=f'根 γ={root_gamma:.4f}')
    # 填充F1(γ)≥0的区域
    plt.fill_between(gamma_range, f_vals, 0, where=(f_vals >=0), color='lightblue', alpha=0.3, label='$F_1(γ)≥0$')

plt.xlabel('γ')
plt.ylabel('$F_1(γ)$')
plt.xlim(1.11, 3)
plt.legend()
plt.grid(True)
plt.show()

方法2：尝试Sympy符号解析求解 #

若符号函数可化简，Sympy的solve能直接给出解析解，避免数值方法的收敛问题：

# 尝试符号求解F1=0
solutions = solve(F1, γ, dict=True)
print("符号解：", solutions)

# 筛选γ>1.11的实数解
valid_solutions = [sol[γ].evalf() for sol in solutions if sol[γ].is_real and sol[γ].evalf() > 1.11]
if valid_solutions:
    root_gamma = valid_solutions[0]
    print(f"有效符号解：γ={root_gamma:.6f}")
else:
    print("无符合条件的符号解，使用数值方法")

关键说明 #

• 二分法是最稳定的根查找算法之一，只要能定位到函数符号变化的区间，就能保证收敛。
• 先通过数值采样观察函数趋势，可避免盲目设置搜索区间导致的失败。
• Matplotlib的fill_between函数能直观展示F1(γ)≥0的解范围，和Desmos的绘图效果一致。

内容的提问来源于stack exchange，提问作者David s.