如何在数值积分的不同步骤中添加随机项
How to add random term in different steps of numerical integration
对于下面的函数,我不会在每个时间步中添加噪声(下面代码中的Inoise
),但是,例如,仅在每隔一个时间步中添加噪声。因此,虽然 dt=0.0025
作为数值积分的时间步长,但我会,例如,仅在每第二个时间步长(即 0.005 步长)中添加 Inoise
。
将其插入我现有函数的最佳方法是什么?
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;
}
您可以累加经过的时间,并在经过足够多的步数后才添加噪声:
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。
对于下面的函数,我不会在每个时间步中添加噪声(下面代码中的Inoise
),但是,例如,仅在每隔一个时间步中添加噪声。因此,虽然 dt=0.0025
作为数值积分的时间步长,但我会,例如,仅在每第二个时间步长(即 0.005 步长)中添加 Inoise
。
将其插入我现有函数的最佳方法是什么?
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;
}
您可以累加经过的时间,并在经过足够多的步数后才添加噪声:
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。