使用C语言rand()函数计算元素半衰期的代码错误排查

问题背景 #

需要用C语言结合rand()函数估算衰变常数λ=0.079(min⁻¹)的元素半衰期，理论值通过公式ln(2)/λ计算约为8.77分钟，但编写的代码输出为12.71分钟，结果偏差较大。

原代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define lambda  0.079       // Decay constant
#define N0      1000        // Number of unstable nuclei in t=0
#define dt      1e-6        // Time interval

main()
{
    double N = N0;   // Number of nuclei
    double t = 0;       // Time counter
    int i;

    while ( N > N0 / 2 ) {
        N = N - ( rand() % 1000 + 1 ) * lambda * dt;
        t = t + dt;
    }

    printf( "The half-life is %.3f", t );
}

错误分析 #

1. 衰变模型完全错误：
   放射性衰变的随机过程是：每个原子核在微小时间间隔dt内衰变的概率为p = lambda * dt。原代码中直接用(rand()%1000+1)*lambda*dt来减少原子核数量，完全不符合随机衰变的概率逻辑，导致N的衰减速度偏离真实规律。
2. 未初始化随机数生成器：
   使用rand()前没有调用srand()初始化随机种子，每次运行程序都会得到完全相同的随机序列，结果不具备统计意义。
3. 时间步长与逻辑不匹配：
   原代码的循环逻辑是每次dt时间就强制减少N，但没有考虑当前N的数量对衰变总数的影响——衰变的原子核数量应该和当前存在的不稳定原子核数量成正比。

修正后的代码 #

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define lambda  0.079       // Decay constant
#define N0      10000       // 增大样本量提升统计准确性
#define dt      1e-3        // 调整时间步长，平衡精度与效率

int main()
{
    // 初始化随机数种子
    srand(time(NULL));
    
    double N = N0;
    double t = 0.0;

    while (N > N0 / 2) {
        // 计算当前时间步内每个原子核的衰变概率
        double decay_prob = lambda * dt;
        int decayed = 0;

        // 遍历每个现存原子核，判断是否衰变
        for (int i = 0; i < (int)N; i++) {
            // 生成0~1之间的随机数，与衰变概率比较
            double rand_val = (double)rand() / RAND_MAX;
            if (rand_val < decay_prob) {
                decayed++;
            }
        }

        N -= decayed;
        t += dt;
    }

    printf("估算的半衰期为 %.3f 分钟，理论值为 %.3f 分钟\n", t, log(2)/lambda);
    return 0;
}

代码说明 #

1. 随机种子初始化：通过srand(time(NULL))用当前时间作为随机种子，确保每次运行的随机序列不同，结果具备统计随机性。
2. 正确的衰变逻辑：
   • 计算每个原子核在dt内的衰变概率decay_prob = lambda * dt
   • 对每个现存的原子核，生成0到1之间的随机数，若随机数小于衰变概率，则判定该原子核衰变，统计衰变总数后从N中减去。
3. 参数优化：
   • 增大N0到10000，提升样本量，减少统计误差；
   • 调整dt为1e-3分钟，既保证模拟精度，又避免过小的时间步导致程序运行过慢。
4. 输出理论值对比：直接在输出中打印理论计算值，方便对比验证结果。

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