时间戳

首页 > 解决方案 > 如何在数值积分的不同步骤中添加随机项

c++ - 如何在数值积分的不同步骤中添加随机项

问题描述

对于下面的函数，我不会Inoise在每个时间步中添加噪声（在下面的代码中），但例如，仅在每个第二个时间步中。因此，虽然dt=0.0025用作数值积分的时间步长，但例如，我Inoise只会在每第二个时间步长（即 0.005 步）中添加。

将其插入现有功能的最佳方法是什么？

    runs = 1000;
    t_end = 5;
    dt = 0.0025;
    t_steps = t_end/dt;

    for(int j=0; j<runs; j++){

        double vT = v0;
        double mT = m0;
        double hT = h0;
        double nT = n0;

        for(int i=0; i<t_steps; i++){

            double IStim = 0.0;
            if ((delay / dt <= (double)i) && ((double)i <= (delay + duration) / dt))
               IStim = I;

            mT = (mT + dt * alphaM(vT)) / (1.0 + dt * (alphaM(vT) + betaM(vT)));
            hT = (hT + dt * alphaH(vT)) / (1.0 + dt * (alphaH(vT) + betaH(vT)));
            nT = (nT + dt * alphaN(vT)) / (1.0 + dt * (alphaN(vT) + betaN(vT)));

            const double iNa = gNa * pow(mT, 3.0) * hT * (vT - vNa);
            const double iK = gK * pow(nT, 4.0) * (vT - vK);
            const double iL = gL * (vT-vL);
            const double Inoise =  (doubleRand() * knoise * sqrt(gNa * A));
            const double IIon = ((iNa + iK + iL) * A) + Inoise;

            vT += ((-IIon + IStim) / C) * dt;
            voltage[i] = vT;

            if(vT > 60.0) {
                count++;
                break;
            }
        }
    }
    return count;
}

标签： c++numerical-integration

解决方案

您可以累积经过的时间，并且只有在经过足够的步骤后才添加噪音：

double elapsedTime = 0;
double INoiseThreshold = 0.005;

for(int j=0; j<runs; j++){

    //...

    for(int i=0; i<t_steps; i++){

        //...

        double Inoise = 0;
        elapsedTime += dt;

        if(elapsedTime >= INoiseThreshold){
            Inoise = (doubleRand() * knoise * sqrt(gNa * A));
            elapsedTime = 0;
        }
        const double IIon = ((iNa + iK + iL) * A) + Inoise;

        //...
    }
}
return count;

  • 您可以检查它们的差异是否在一个小的 epsilon 内，以允许舍入误差，而不是直接比较浮点数。

  • 与其使 Inoise 的值取决于条件，不如使 IIon 公式中的存在取决于例如：

    const double IIon = ((iNa + iK + iL) * A) + (elapsedTime >= INoiseThreshold) ? Inoise : 0;

    只要记住一旦超过阈值就重置elapsedTime。

推荐阅读