Python scipy dblquad积分良定义函数报sqrt无效值错误求解

scipy.integrate.dblquad采用自适应求积算法，采样计算点时不会严格限制在传入的积分上下限范围内，会因浮点数舍入误差、自适应步长外推逻辑，采集到少量略超出积分区间的点。
你定义的p积分上限为50/101≈0.4950495，q积分上限为51/101≈0.5049505，当算法采到p比上限大1e-12量级、q比上限大1e-12量级的越界点时，sqrt(0.4950495... - p)、sqrt(0.5049505... - q)的自变量会变成绝对值极小的负数，就会触发"invalid value encountered in sqrt"错误。
你之前逐点验证函数合法性时，选取的都是严格落在定义域内的点，自然无法复现这个采样越界导致的问题。

• 优先选择对被积函数的sqrt输入做钳位处理，用np.clip把所有sqrt的自变量限制在非负范围，越界产生的极小负值截断为0即可，对最终积分精度的影响可以忽略。
• 不推荐直接全局关闭numpy浮点错误警告，这种方式会掩盖真正的非法输入问题，不利于后续调试。
• 原代码中重复导入了两次matplotlib.pyplot，属于冗余代码，可以直接删除。

import sympy as sym
from scipy import integrate
import numpy as np
from numpy import sqrt

beta,h_1,h_2,h_3,h_4,a_1,a_2,m,p,q=sym.symbols('beta h_1 h_2 h_3 h_4 a_1 a_2 m p q', real=True)

f_1=sym.diff((1/beta)*sym.log( (sym.sqrt(p)+4/5*sym.sqrt(50/101-p))**2*(sym.exp(beta*(m+h_1)*a_1)+sym.exp(beta*(m+h_1)*a_2))+
        +(4/5*sym.sqrt(p)+sym.sqrt(50/101-p))**2*(sym.exp(beta*(m+h_2)*a_1)+sym.exp(beta*(m+h_2)*a_2))
        +(sym.sqrt(q)+4/5*sym.sqrt(51/101-q))**2*(sym.exp(beta*(m+h_3)*a_1)+sym.exp(beta*(m+h_3)*a_2))
        +(4/5*sym.sqrt(q)+sym.sqrt(51/101-q))**2*(sym.exp(beta*(m+h_4)*a_1)+sym.exp(beta*(m+h_4)*a_2))),m).subs({beta:10,h_1:0,h_2:1.01/sym.sqrt(2),h_3:-1.01/sym.sqrt(2),h_4:0,a_1:1/sym.sqrt(2),a_2:-1/sym.sqrt(2),m:0})

print(f_1)

P_MAX = 50/101
Q_MAX = 51/101
# 对所有sqrt的输入做非负钳位，避免越界采样触发错误
f=lambda p,q: (
    780.080275764744*sqrt(2)*(0.8*sqrt(np.clip(p, 0, P_MAX)) + sqrt(np.clip(P_MAX - p, 0, P_MAX)))**2 
    - 780.080275764744*sqrt(2)*(sqrt(np.clip(q, 0, Q_MAX)) + 0.8*sqrt(np.clip(Q_MAX - q, 0, Q_MAX)))**2
)/(
    10*(
        156.028873819841*(0.8*sqrt(np.clip(p, 0, P_MAX)) + sqrt(np.clip(P_MAX - p, 0, P_MAX)))**2 
        + 2*(sqrt(np.clip(p, 0, P_MAX)) + 0.8*sqrt(np.clip(P_MAX - p, 0, P_MAX)))**2 
        + 2*(0.8*sqrt(np.clip(q, 0, Q_MAX)) + sqrt(np.clip(Q_MAX - q, 0, Q_MAX)))**2 
        + 156.028873819841*(sqrt(np.clip(q, 0, Q_MAX)) + 0.8*sqrt(np.clip(Q_MAX - q, 0, Q_MAX)))**2
    )
)

I=integrate.dblquad(f, 0, P_MAX, 0, Q_MAX)
print(I)

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