Java中计算边界含x的数值二重积分的问题咨询

问题原因 #

你当前的实现预设了z方向（内积分）的上下限是全局固定的，所以提前生成了固定尺寸的二维网格采样，自然不支持z上下限随x变化的场景。

具体修改方案 #

1. 重构积分计算逻辑，去掉预先生成二维数组的步骤

不需要提前把所有(x,z)点的被积函数值都算出来存着，改为遍历每个x采样点时，单独计算当前x对应的z积分上下限，再对该x做独立的z方向辛普森积分，得到对应x的一维积分结果后，再参与x方向的辛普森积分计算即可。

2. 修复被积函数的语法错误

你当前写的integrand方法括号位置错误，-0.5只乘了二次型的第一项，剩下的交叉项和z的二次项都没有被乘到-0.5，会导致计算结果完全错误，需要把整个二次型的表达式都括起来。

3. 适配动态区间的辛普森采样要求

辛普森1/3规则要求每个积分方向的采样点数为奇数（即区间分割数为偶数），对于每个x对应的不同z积分区间，需要单独计算采样点数，同时保证点数为奇数，避免精度损失。

修改后的代码示例 #

// 补充两个方法定义z的积分上下限，根据你实际的边界表达式修改即可
private double getLowerBoundZ(double x) {
    // 示例：假设z下限是2*x，替换为你实际的边界
    return 2 * x;
}

private double getUpperBoundZ(double x) {
    // 示例：假设z上限是5*x + 1，替换为你实际的边界
    return 5 * x + 1;
}

// 修复后的被积函数
private double integrand(double x, double z) {
    double[][] mat = matrix.getData();
    return Math.exp(-0.5 * (
            mat[0][0] * Math.pow(x, 2) +
            mat[0][1] * z * x +
            mat[1][0] * z * x +
            mat[1][1] * Math.pow(z, 2)
    ));
}

public double calculate() {
    int nx = (int) ((upperBoundX - lowerBoundX) / stepSizeX + 1);
    // 保证x方向采样点数为奇数
    if (nx % 2 == 0) nx += 1;
    double adjustedStepX = (upperBoundX - lowerBoundX) / (nx - 1);
    double answer = 0;

    for (int i = 0; i < nx; ++i) {
        double currentX = lowerBoundX + i * adjustedStepX;
        // 获取当前x对应的z积分上下限
        double lbZ = getLowerBoundZ(currentX);
        double ubZ = getUpperBoundZ(currentX);
        // 计算当前x对应的z方向采样点数，保证为奇数
        int nz = (int) ((ubZ - lbZ) / stepSizeY + 1);
        if (nz % 2 == 0) nz += 1;
        double adjustedStepZ = (ubZ - lbZ) / (nz - 1);
        
        // 对当前x做z方向辛普森积分
        double integralZ = 0;
        for (int j = 0; j < nz; ++j) {
            double currentZ = lbZ + j * adjustedStepZ;
            double val = integrand(currentX, currentZ);
            if (j == 0 || j == nz - 1) {
                integralZ += val;
            } else if (j % 2 == 0) {
                integralZ += 2 * val;
            } else {
                integralZ += 4 * val;
            }
        }
        integralZ *= (adjustedStepZ / 3);
        
        // 累加进入x方向辛普森积分
        if (i == 0 || i == nx - 1) {
            answer += integralZ;
        } else if (i % 2 == 0) {
            answer += 2 * integralZ;
        } else {
            answer += 4 * integralZ;
        }
    }
    answer *= (adjustedStepX / 3);
    return answer;
}

如果需要更高的计算精度，也可以对每个积分区间做自适应辛普森分割，不需要固定全局步长，能在被积函数变化剧烈的区域自动加密采样。

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