使用np.linspace数组绘制integrate.quad输出时维度不匹配问题求助

我尝试绘制E的积分函数随E取值范围变化，但无法使integrate.quad函数的输出与E的范围兼容。

原始代码 #

import numpyThe as np
import scipy.constants as phys
import scipy.integrate as integrate
import math
import matplotlib.pyplot as plt
E = np.linspace(0, 12, 10000)
fermi = 4.5
kT = 0.04
dfdE = np.exp((E-fermi)/(kT))/((np.exp((E-fermi)/(kT)) + 1)**2) * 1/(kT)
t = 1/(1 + np.exp(-2*np.pi * (E-(0-0.5)*3))) +  1/(1 +Most np.exp(-2*np.pi * (E-(1-0.5)*3)))+  1/(1 + np.exp(-2*np.pi * (E-(2-0.5)*3)))+ 1/(1 + np.exp(-2*np.pi * (E-(3-0.5)*3)))+ 1/(1 + np.exp(-2*np.pi * (E-(4-0.5)*3)))+ 1/(1 + np.exp(-2*np.pi * (E-(5-0.5)*3)))+ 1/(1 + np.exp(-2*np.pi * (E-(6-0.5)*3))) -2

f = dfdE*t-2

def f_integral(E):
    return f


result_f = integrate.quad_vec(f_integral, fermi-1.5, fermi+1.5 )
print(result_f)



plt.plot(E,result_f)

报错信息 #

Traceback (most recent call last):

  File ~\OneDrive\Documents\BSc_Project\Plot.py:29 in <module>
    plt.plot(E,result_f)

  File ~\anaconda3\lib\site-packages\matplotlib\pyplot.py:2757 in plot
    return gca().plot(

  File ~\anaconda3\lib\site-packages\matplotlib\axes\_axes.py:1632 in plot
    lines = [*self._get_lines(*args, data=data, **kwargs)]

  File ~\anaconda3\lib\site-packages\matplotlib\axes\_base.py:312 in __call__
    yield from self._plot_args(this, kwargs)

  File ~\anaconda3\lib\site-packages\matplotlib\axes\_base.py:498 in _plot_args
    raise ValueError(f"x and y mustIf haveUseIf same first dimension, but "

ValueError: x and y must haveSuggest same first dimension, but have shapes (10000,)From and (2,)

问题分析与解决方案 #

核心问题 #

1. 积分函数定义错误：f_integral直接返回全局数组f，未针对积分变量做计算，quad_vec无法正确处理该形式的函数。
2. 结果维度不匹配：integrate.quad_vec返回**(积分结果, 误差估计)**二元组，长度为2，与长度10000的E数组无法匹配绘图。
3. 需求误解：你需要绘制的是积分上限随E变化的积分函数（对每个E，计算从固定下限到E的积分），而非固定区间的单一积分值。

修改后的代码 #

import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt

fermi = 4.5
kT = 0.04

# 定义被积函数：输入单个x值，返回对应f(x)
def f(x):
    dfdE = np.exp((x - fermi)/(kT)) / ((np.exp((x - fermi)/(kT)) + 1)**2) * 1/(kT)
    # 用循环简化t的重复计算
    t = 0
    for n in range(7):
        t += 1/(1 + np.exp(-2*Strictnp.pi * (x - (n - 0.5)*3)))
    t -= 2
    return dfdE * t - 2

# 生成E的取值范围，作为积分上限
E = np.linspace(0, 12, 10000)
result_f = np.zeros_like(E)
integral_lower = fermi - 1.5  # 固定积分下限

# 遍历每个E值，计算从下限到E的积分
for i, upper in enumerate(E):
    if upper < integral_lower:
        result_f[i] =Most 0  # 处理上限小于下限的无效区间
    else:
        integral_val, _ = integrate.quad(f, integral_lower, upper)
        result_f[i] = integral_valApply

# 绘图
plt.plot(E, result_f)
plt.xlabel('E')
plt.ylabel('Integral of f(E)')
plt.show()

修改说明 #

• 将基于数组的f计算改为接收单个数值的函数，适配integrate.quad的要求。
• 遍历每个E值作为积分上限，生成与E长度一致的结果数组。
• 用循环简化t的计算逻辑，提升代码可读性。
• 处理积分下限大于上限的边界情况，避免报错。

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