嘿，我懂你复现杨氏双缝实验时卡在数值积分这一步的苦恼——虽然问题根源于物理实验，但核心是Matlab里的无穷区间积分实现难题，咱们一步步来拆解。

你需要计算的是 q(xp, yp) = ∫∫(-∞到+∞) f(x,y,xp,yp) dx dy，这类带无穷区间的二重数值积分，其实不用硬啃“无穷”这个坎，结合物理特性和Matlab工具就能搞定：

• 先把无穷区间“砍”成有限范围
  杨氏双缝对应的干涉/衍射函数，在远离缝的位置衰减极快，完全可以根据函数特性截断到一个合理的有限区间。比如先分析f(x,y,xp,yp)的衰减规律，找到一个阈值（比如函数值小于1e-8时就可忽略），确定积分上下限[-L, L]，L的大小可以结合你的实验参数（缝宽、波长、接收屏距离等）来定，这样既保证精度又提升计算效率。

• 选对Matlab的积分工具

  • 如果你的f是光滑函数，优先用integral2，它专门处理二重数值积分，用法示例：

    % 先定义四变量函数f，把xp、yp作为额外参数传入
    f = @(x,y,xp,yp) ... % 这里填入你的f(x,y,xp,yp)具体表达式
    % 对给定的xp、yp计算积分
    q = integral2(@(x,y) f(x,y,xp,yp), -L, L, -L, L);
    

  • 如果函数存在奇点或者振荡剧烈，integral2精度不够的话，可以试试嵌套integral，先对x积分再对y积分：

    q = integral(@(y) integral(@(x) f(x,y,xp,yp), -L, L), -L, L);
    

• 批量计算xp、yp网格点的q值
  要是你需要计算一系列(xp, yp)点的结果，别用循环逐个跑，效率太低。可以把xp和yp做成网格矩阵，用arrayfun批量处理：

  xp_range = linspace(-5, 5, 100); % 示例xp范围和点数
  yp_range = linspace(-5, 5, 100);
  [Xp, Yp] = meshgrid(xp_range, yp_range);
  Q = arrayfun(@(xp,yp) integral2(@(x,y) f(x,y,xp,yp), -L, L, -L, L), Xp, Yp);
  

• 验证积分结果的可靠性
  调整L的大小，看看当L增大到某个值后，q的结果不再明显变化，这就说明截断区间是合理的。另外，如果杨氏双缝的场景有解析解，也可以拿数值结果和解析解对比，验证积分的正确性。

要是你的f函数包含复指数、三角函数这类复杂形式，记得用Matlab对应的内置函数（比如exp、sin），同时注意物理量的单位统一，避免因单位混乱导致计算错误。

内容的提问来源于stack exchange，提问作者A.Maine