蒙特卡洛模拟C代码输出-nan值：原因排查及含义解析

什么是-nan？ #

-nan是**无效数值（Not a Number）**的变体，代表计算过程中出现了无意义的操作，常见场景包括：

• 0除以0
• 对负数进行开平方运算
• 无穷大与无穷大相减
• 未初始化的垃圾值参与数值计算

代码中的核心问题 #

1. 数组未初始化 #

你定义的num_1、num_2、num_3、num_a四个局部数组默认存储的是垃圾值，直接执行num_1[time] += cover1相当于给随机垃圾值累加，最终导致均值计算出现无效数值。

2. 模拟状态未重置 #

外层MCSTEPS循环中，每次蒙特卡洛步都应该重新初始化网格和粒子计数器（particle1、particle2等），否则所有步骤都会在同一个网格上持续模拟，不仅不符合独立迭代的要求，还可能导致粒子数变为负数，引发异常计算。

3. 标准差计算逻辑错误 #

CalcStdDev函数中，你错误地将方差进行平方运算（Square(sum/(double)steps)），但标准差的正确计算是方差的平方根，而非方差的平方。

4. 近邻检测不完整 #

模拟要求检查目标位点的四个最近邻，但你的代码仅检测了右侧邻居，模拟逻辑不符合需求，会导致粒子生成/捕获的规则执行错误，进而影响数值合理性。

修复方案 #

步骤1：初始化数组 #

在定义数组时直接将所有元素初始化为0：

double num_1[MCSTEPS] = {0};
double num_2[MCSTEPS] = {0};
double num_3[MCSTEPS] = {0};
double num_a[MCSTEPS] = {0};

步骤2：重置模拟状态 #

将网格初始化和计数器重置放到外层MCSTEPS循环内部，确保每次蒙特卡洛步都是独立的模拟：

for(int time = 0; time < MCSTEPS; time++ ) {
    // 每次MC步重置网格与计数器
    grid_init(grid);
    particle1 = 0;
    particle2 = 0;
    particle3 = 0;
    availcells = MAX_X * MAX_Y;
    fullcells = 0;

    for(int iter = 0; iter < ITERATIONS; iter++){
        // 原有迭代逻辑...
    }
}

步骤3：修正标准差计算 #

修改CalcStdDev函数，使用平方根计算标准差：

double CalcStdDev( double number[], int steps, double * pMean )
{
    double sum = 0;
    double stddev = 0;
    
    for (int k = 0; k < steps; k++) {
        sum += Square( number[k] - *pMean );
    }
    
    // 计算方差的平方根得到标准差
    stddev = sqrt( sum / (double) steps );

    return stddev;
}

步骤4：完善四近邻检测 #

实现覆盖四个方向的近邻检测函数，并替换原有单方向逻辑：

// 检查目标位点的四个近邻是否存在指定粒子
int has_target_neighbor(int x, int y, gstate target, gstate grid[MAX_Y][MAX_X]){
    int dx[] = {-1, 1, 0, 0}; // 左、右、上、下偏移
    int dy[] = {0, 0, -1, 1};
    
    for(int d=0; d<4; d++){
        int nx = get_limited_coord(x + dx[d], MAX_X);
        int ny = get_limited_coord(y + dy[d], MAX_Y);
        if(grid[ny][nx] == target){
            return 1;
        }
    }
    return 0;
}

在粒子选择逻辑中使用该函数：

if(r <= Y){//选择粒子1
    if(has_target_neighbor(i, j, P2_OCCUPIED, grid)){
        // 遍历近邻，找到并移除一个粒子2（逻辑可根据需求调整）
        int dx[] = {-1, 1, 0, 0};
        int dy[] = {0, 0, -1, 1};
        for(int d=0; d<4; d++){
            int nx = get_limited_coord(i + dx[d], MAX_X);
            int ny = get_limited_coord(j + dy[d], MAX_Y);
            if(grid[ny][nx] == P2_OCCUPIED){
                grid[ny][nx] = S_EMPTY;
                particle3++;
                particle2--;
                fullcells--;
                availcells++;
                break;
            }
        }
    } else { 
        grid[j][i] = P1_OCCUPIED;
        particle1++;
        availcells--;
        fullcells++;
    }   
}

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