如何执行数值积分？量子电容积分计算报错求助

问题分析 #

第一段代码报错原因 #

1. 常量未定义：直接使用e、kb、T、mu但未提前赋值
2. 积分式语法/逻辑错误：np.cosh(((x/(2*kb*T))**2) (x-mu))缺少运算符，公式不符合量子电容的正确形式；同时直接传入整个y_array（DOS数组）与积分变量相乘，而非根据当前能量值插值得到对应DOS
3. 积分区间错误：仅积分-1到1，未覆盖数据的能量范围（-18.68到18.73）

第二段代码警告原因 #

1. 公式错误：误用hbar，量子电容的正确公式为 $C_q = \frac{e^2}{4k_B T} \int DOS(E) \cdot \text{sech}^2\left(\frac{E-\mu}{2k_B T}\right) dE$
2. 数值溢出：当$|E-\mu|$远大于$2k_B T$时，$\cosh\left(\frac{E-\mu}{2k_B T}\right)$会指数级增长导致溢出，直接取倒数会出现数值不稳定

解决方案 #

核心修正点 #

1. 定义所有物理常量，确保单位一致
2. 使用插值获取任意能量对应的态密度（DOS）
3. 用稳定的数学表达式计算$\text{sech}^2(z)$，避免数值溢出：
   $$\text{sech}^2(z) = \frac{4e^{{-2|z|}}{(1+e}{-2|z|})^2}$$
   当$z$绝对值很大时，该表达式不会产生溢出
4. 覆盖完整的能量数据区间进行积分，或合理截断（因为$\text{sech}^2(z)$在$|z|>5$时已趋近于0）

完整可运行代码 #

import numpy as np
from scipy.integrate import quad

# 定义物理常量
e = 1.602e-19       # 元电荷 (C)
kb = 1.381e-23      # 玻尔兹曼常数 (J/K)
T = 300             # 温度 (K)
mu = 0.5            # 化学势

# 加载E和DOS数据（直接使用提供的数值）
E_data = np.array([
    -18.68526, -17.05825, -15.43124, -13.80423, -12.17723, -10.55022,
    -8.92321, -7.2962, -5.6692, -4.04219, -2.41518, -0.78817, 0.83883,
    2.46584, 4.09285, 5.71986, 7.34686, 8.97387, 10.60088, 12.22789,
    13.85489, 15.4819, 17.10891, 18.73592
])
DOS_data = np.array([
    11.5535, 10.53234, 14.4706, 10.63194, 12.04803, 9.51133, 12.07561,
    25.97326, 26.43385, 21.97081, 14.23236, 5.39077, 3.17652, 13.11387,
    12.40104, 13.84468, 19.24209, 34.36795, 30.42492, 26.10606, 32.17047,
    34.89188, 30.86395, 0.02738
])

# 定义稳定的sech²函数，避免溢出
def sech_squared(z):
    abs_z = np.abs(z)
    exp_term = np.exp(-2 * abs_z)
    return 4 * exp_term / (1 + exp_term) ** 2

# 积分被积函数
def integrand(E):
    # 插值获取当前能量对应的DOS
    dos = np.interp(E, E_data, DOS_data)
    z = (E - mu) / (2 * kb * T)
    return dos * sech_squared(z)

# 计算量子电容
prefactor = e**2 / (4 * kb * T)
# 积分区间取数据的最小和最大能量
E_min = E_data.min()
E_max = E_data.max()
integral_result, error = quad(integrand, E_min, E_max)
Cq = prefactor * integral_result

print(f"量子电容: {Cq:.3e} F/m²")
print(f"积分误差: {error:.3e}")

代码解释 #

1. 数据加载：直接将提供的E和DOS数据转为numpy数组，保证数值计算效率
2. 稳定的sech²实现：通过指数函数的负项避免大数值溢出，保证计算稳定性
3. 插值处理：使用np.interp在离散DOS数据间插值，得到任意能量点的态密度
4. 积分区间：覆盖所有数据点范围，确保积分的完整性；也可根据sech²的衰减特性截断区间（如mu±10*2*kb*T），提升计算效率

从Excel加载数据的替代代码片段 #

如果需要从DOS.xlsx读取数据，替换上述的E_data和DOS_data定义：

import openpyxl

wb = openpyxl.load_workbook(filename='DOS.xlsx', data_only=True)
sheet = wb["Sheet1"]

# 读取所有行（跳过表头）
rows = list(sheet.iter_rows(min_row=2, values_only=True))
E_data = np.array([row[0] for row in rows])
DOS_data = np.array([row[1] for row in rows])

内容的提问来源于stack exchange，提问作者Protima Rani Paul