如何用代码求解含嵌套高斯积分的反常二重积分？

先明确：原积分是发散的 #

当$y \to \infty$时，内层积分$\int_0^y e^{-x2}dx$收敛到$\frac{\sqrt{\pi}}{2}$，因此被积函数$y^2 e^{\int_0y e^{-x2}dx}$趋近于$y^2 e^{{\frac{\sqrt{\pi}}{2}}$，而$\int_0}\infty y^2 dy$是发散的，所以原无穷积分不存在有限值。这应该是你用quadgk得不到理想结果的核心原因——积分本身发散，必须截断到有限上限计算近似值。

优化计算的可行方案 #

1. 用MATLAB内置误差函数替代内层积分，提升计算精度 #

内层积分$\int_0^y e^{-x2}dx$可以用误差函数精确表示：
$$\int_0^y e^{-x2}dx = \frac{\sqrt{\pi}}{2} \text{erf}(y)$$
MATLAB的erf是高度优化的内置函数，计算精度远高于多项式拟合，不会出现病态问题。原积分改写为：
$$\int_0^Y y^2 e^{\frac{\sqrt{\pi}}{2} \text{erf}(y)} dy$$
其中$Y$是你选定的截断上限。直接用quadgk计算的示例代码：

Y = 5; % 可根据需求调整截断上限
f = @(y) y.^2 .* exp( (sqrt(pi)/2)*erf(y) );
result = quadgk(f, 0, Y, 'RelTol', 1e-8, 'AbsTol', 1e-10);
disp(result);

2. 自适应选择截断上限，平衡精度与计算量 #

如果不确定选多大的$Y$，可以通过判断被积函数的增长速率来确定：当$y$足够大时，$\text{erf}(y) \approx 1 - \frac{e^{-y2}}{\sqrt{\pi}y}$，此时被积函数近似为$y^2 e^{\frac{\sqrt{\pi}}{2}} \left(1 - \frac{e^{-y2}}{\sqrt{\pi}y}\right)$，后续积分的增量可以近似估算。比如设置一个阈值，当新增的积分增量小于你设定的精度要求时，就停止增大$Y$。示例逻辑：

target_tol = 1e-6;
Y = 1;
prev_result = 0;
while true
    f = @(y) y.^2 .* exp( (sqrt(pi)/2)*erf(y) );
    current_result = quadgk(f, 0, Y);
    delta = current_result - prev_result;
    if abs(delta) < target_tol
        break;
    end
    prev_result = current_result;
    Y = Y + 1; % 每次增大1，也可调整步长
end
disp(['截断上限Y=', num2str(Y), '时，积分结果为', num2str(current_result)]);

3. 变量替换降低被积函数增长速度 #

考虑做变量替换$t = e^{-y}$，把无穷积分转化为有限区间积分，同时减缓被积函数的增长。令$y = -\ln t$，则$dy = -\frac{1}{t}dt$，当$y=0$时$t=1$，$y\to\infty$时$t\to0$，原积分变为：
$$\int_0^1 (-\ln t)^2 e^{\frac{\sqrt{\pi}}{2} \text{erf}(-\ln t)} \cdot \frac{1}{t} dt$$
这个积分的被积函数在$t\to0$时的增长速度会比原函数慢，可能更适合数值积分计算。示例代码：

f = @(t) (-log(t)).^2 .* exp( (sqrt(pi)/2)*erf(-log(t)) ) ./ t;
result = quadgk(f, 0, 1, 'RelTol', 1e-8);
disp(result);

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