我正在使用 C++ 模拟一个涉及几个 (27) 个刚性常微分方程的生物模型。我的程序在 MS C++ 2010 表达式编译器下完美运行，但在 g++ 编译器（NetBeans 6.8、Ubuntu 10.04 LTS）下失败。问题是一些变量变成了NaN。以下是g++编译器下程序每一步运行后变量
Vm下面是相同代码在 MS C++ 编译器下未经任何更改的输出：-59.4 -59.3993 -59.6081 100.081 34.6378 -50392.8 南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南南楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠楠纳米
-59.4 -59.3993 -59.3986 -59.3979 -59.3972 -59.3966 -59.3959 -59.3952 -59.3946 -59.3939 -59.3933 -59.3926 -59.392 -59.3914 -59.3907 -59.3901 -59.3895 -59.3889 -59.3883 -59.3877 -59.3871 -59.3865 -59.3859 -59.3853 -59.3847 -59.3841 -59.3836 -59.383 -59.3824 -59.3819 -59.3813 -59.3808 -59.3802 -59.3797 -59.3791 -59.3786 -59.3781 -59.3775 -59.377 -
我只使用了“cmath”和“fstream”库。
问题出在哪里？两种场景的代码完全相同。
编辑1：
好的大家，这是完整的代码：
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
void FCN(void);
const int TMAX = 10000; //[ms] simulation time
const int TSTEP = 1;
const int MXSTEP = TMAX / TSTEP;
int ISTEPPRINT = MXSTEP / 5000; //time step for storing on disc
int ISTEP = 0;
const double R = 8341.0;
const double temp = 293.0;
const double F = 96487.0;
const double RT_F = R*temp / F;
const double z_K = 1;
const double z_Na = 1;
const double z_Ca = 2;
const double z_Cl = -1;
const double N_Av = 6.022e23;
double Ca_o = 2.0;
double Na_o = 140.0;
double Cl_o = 129.0;
double K_o = 5;
double NE = 0;
double NO = 0;
double cGMP = 0; //[mM] [cGMP]i
double cGMPprime = 0; //var
double IP3 = 0; //[mM] [IP3]i
double IP3d = 0; //var
double IP3prime = 0; //var
double DAG = 0; //[mM]
double DAGprime = 0; //var
double Ca_u = 0.66;
double Ca_r = 0.57;
double Ca_i = 68e-6;
double Na_i = 8.4;
double K_i = 140;
double Cl_i = 59.4;
double V_m = -59.4;
double V_mprime;
double Na_iprime;
double K_iprime;
double Cl_iprime;
double Ca_iprime;
double vol_i = 1;
double I_Natotm;
double I_Ktotm;
double I_Cltotm;
double I_Catotm; //[pA] total membrane Ca current
//Reversal potentials
double E_Ca; //[mV]
double E_Na; //[mV]
double E_K; //[mV]
double E_Cl;
//Membrane capacitance and area
double C_m = 25.0;
double A_m = C_m / 1e6;
//Voltage dependent calcium current I_CaL
double I_CaL;
double P_CaL = 1.88e-5;
double d_L;
double d_Lprime;
double d_Lbar;
double tau_d_L;
double f_L;
double f_Lbar;
double tau_f_L;
double f_Lprime;
//Delayed rectifier current I_K
double I_K;
double g_K = 1.35;
double p_K;
double p_Kbar;
double V_1_2 = -11.0;
double k = 15.0;
double tau_p_K;
double p_Kprime;
double q_1 = 1;
double q_2 = 1;
double q_bar;
double q_1prime;
double q_2prime;
double Pmin_NSC = 0.4344;
double Po_NSC;
double PNa_NSC = (5.11e-7);
double PCa_NSC = (5.11e-7)*4.54;
double PK_NSC = (5.11e-7)*1.06; //
double d_NSCmin = 0.0244;
double K_NSC = 3.0e-3;
double INa_NSC;
double ICa_NSC;
double IK_NSC;
double I_NSC;
//KATP current I_KATP
double I_KATP; //[pA] background K current
double g_KATP = 0.067; //[nS] max. background K current conductance
//Inward rectifier current I_K_i
double I_K_i; //[pA]
double g_maxK_i; //[nS] max. slope conductance
const double G_K_i = 0; // inward rectifier constant
const double n_K_i = 0.5; // inward rectifier constant
//Calcium-activated potassium current I_KCa
double I_KCa;
double i1_KCa;
double P_BKCa = 3.9e-13;
double N_BKCa = 6.6e6;
double P_KCa;
double p_obar;
double V_1_2_KCa;
double p_f;
double p_s;
double p_fprime;
double p_sprime;
double tau_pf = 0.84;
double tau_ps = 35.9;
double dV_1_2_KCa_NO = 46.3;
double R_NO;
double dV_1_2_KCa_cGMP = 76;
double R_cGMP;
double k_leak = 1;
double R_00;
double R_01 = 0.9955;
double R_10 = 0.0033;
double R_11 = 4.0e-6;
double R_01prime;
double R_10prime;
double R_11prime;
const double Kr1 = 2500.0;
const double Kr2 = 1.05;
const double K_r1 = 0.0076;
const double K_r2 = 0.084;
double I_up;
const double I_upbar = 3.34 * (k_leak + 1);
const double K_mup = 0.001;
double I_tr;
const double vol_u = 0.07;
double tau_tr = 1000.0;
double I_rel;
const double vol_r = 0.007;
const double tau_rel = 0.0333; //[ms]
const double R_leak = 1.07e-5 * (k_leak); ////equal to R_10^2 during concentration clamp
// time constant of SR release
double Ca_uprime; // dCa_u/dt
double Ca_rprime; // dCa_r/dt
double S_CM;
const double K_d = 2.6e-4;
const double S_CMbar = 0.1;
double CaCM;
const double K_dB = 5.298e-4;
const double B_Fbar = 0.1;
const double vol_Ca = 0.7;
const double CSQNbar = 15;
const double K_CSQN = 0.8;
double I_PMCA;
double I_PMCAbar = 5.37;
double K_mPMCA = 170e-6;
double I_NaK; ////[pA] Na/K pump
double I_NaKbar = 2.3083;
const double K_mK = 1.6;
const double K_mNa = 22;
const double Q_10_NaK = 1.87;
double I_NCX;
const double gamma2 = 0.45; //
double g_NCX = 0.000487; //[nS]
double d_NCX = 0.0003; //
double Fi_F; //
double Fi_R; //
double I_NaKCl_Cl; //[pA]
double L_cotr = 1.79e-8;
double I_M = 0; //[pA]
double I_MCa = 0;
double I_MNa = 0; //[pA] Na component
double I_MK = 0;
double I_SOC; //[pA]
double I_SOCCa;
double I_SOCNa; //[pA] Na component
const double g_SOCCa = 0.0083; //[nS]
const double g_SOCNa = 0.0575; //[nS]
const double H_SOC = 1;
const double K_SOC = 0.0001;
const double tau_SOC = 100;
double P_SOCbar;
double P_SOC = 0;
double P_SOCprime;
//Chloride currents
double I_Cl;
double g_Cl = 0.23;
double alpha_Cl;
double P_Cl;
//Stimulation current
double I_stim = 0; //[pA]
//IP3 receptor
double I_IP3;
double I_IP3bar = 2880e-6; //[1/ms]
double K_IP3 = 0.12e-3;
double K_actIP3 = 0.17e-3;
double K_inhIP3 = 0.1e-3; //[mM]
double h_IP3;
double k_onIP3 = 1.4;
double h_IP3prime;
double R_TG = 2e4;
double K_1G = 0.01;
double K_2G = 0.2;
double k_rG = 1.75e-7;
double k_pG = 0.1e-3;
double k_eG = 6e-6;
double ksi_G = 0.85;
double G_TG = 1e5;
double k_degG = 1.25e-3;
double k_aG = 0.17e-3;
double k_dG = 1.5e-3;
double PIP2_T = 5e7;
double r_rG = 0.015e-3;
double K_cG = 0.4e-3;
double alpha_G = 2.781e-8;
double vol_IP3 = vol_i;
double gamma_G = N_Av*vol_IP3 * 1e-15;
double RS_G = R_TG*ksi_G;
double RS_PG = 0;
double PIP2; //
double r_hG;
double G;
double delta_G; //
double RS_Gprime;
double RS_PGprime;
double Gprime;
double PIP2prime;
double rho_rG;
//cGMP formation
double k1sGC = 2e3; //[1/mM/ms]
double k_1sGC = 15e-3; //[1/ms]
double k2sGC = 0.64e-5; //[1/ms]
double k_2sGC = 0.1e-6; //[1/ms]
double k3sGC = 4.2; //[1/mM/ms]
double kDsGC = 0.4e-3;
double kDact_deactsGC = 0.1e-3; //[1/ms]
double V_cGMP = 0; //
double V_cGMPprime;
double V_cGMPmax = 0.1 * 1.26e-6; //[mM/ms]
double V_cGMPbar;
double B5sGC = k2sGC / k3sGC;
double A0sGC = ((k_1sGC + k2sGC) * kDsGC + k_1sGC*k_2sGC) / (k1sGC*k3sGC);
double A1sGC = ((k1sGC + k3sGC) * kDsGC + (k2sGC + k_2sGC) * k1sGC) / (k1sGC*k3sGC);
double kpde_cGMP = 0.0695e-3; //[1/ms]
double tausGC;
const int N = 27;
double Y[N], Y1[N], YPRIME[N];
int N1 = 33;
double T = 0;
int main(void) {
ofstream fileY, fileY1, fileT;
// initial conditions SMC
//ICaL
d_Lbar = 1.0 / (1 + exp(-(V_m) / 8.3));
d_L = d_Lbar;
f_Lbar = 1.0 / (1 + exp((V_m + 42.0) / 9.1));
f_L = f_Lbar;
//IKCa
R_NO = NO / (NO + 200e-6);
R_cGMP = pow(cGMP, 2) / (pow(cGMP, 2) + pow(0.55e-3, 2));
V_1_2_KCa = -41.7 * log10(Ca_i) - 128.2 - dV_1_2_KCa_NO * R_NO - dV_1_2_KCa_cGMP*R_cGMP;
p_obar = 1 / (1 + exp(-(V_m - V_1_2_KCa) / 18.25));
p_f = p_obar;
p_s = p_obar;
//I_K
p_Kbar = 1 / (1 + exp(-(V_m - V_1_2) / k));
p_K = p_Kbar;
q_bar = 1.0 / (1 + exp((V_m + 40) / 14));
q_1 = q_bar;
q_2 = q_bar;
//IP3 receptor
h_IP3 = K_inhIP3 / (Ca_i + K_inhIP3);
//norepinephrine receptor
PIP2 = PIP2_T - (1 + k_degG / r_rG) * gamma_G*IP3;
r_hG = k_degG * gamma_G * IP3 / PIP2;
G = (K_cG + Ca_i) / (alpha_G * Ca_i) * r_hG;
delta_G = k_dG * G / (k_aG * (G_TG - G));
Y[0] = V_m;
Y[1] = d_L;
Y[2] = f_L;
Y[3] = p_K;
Y[4] = q_1;
Y[5] = p_f;
Y[6] = p_s;
Y[7] = R_01;
Y[8] = R_10;
Y[9] = R_11;
Y[10] = Ca_u;
Y[11] = Ca_r;
Y[12] = Ca_i;
Y[13] = Na_i;
Y[14] = K_i;
Y[15] = q_2;
Y[16] = P_SOC;
Y[17] = Cl_i;
Y[18] = h_IP3;
Y[19] = RS_G;
Y[20] = RS_PG;
Y[21] = G;
Y[22] = IP3;
Y[23] = PIP2;
Y[24] = DAG;
Y[25] = V_cGMP;
Y[26] = cGMP;
ISTEP = -1;
T = 0.0 - TSTEP;
fileY.open("Y.txt");
fileY1.open("Y1.txt");
fileT.open("T.txt");
for (;;) {
ISTEP = ISTEP + 1;
T = T + TSTEP;
//Norepinephrine
if (T > 10000) NE = 1e-3; //NE [mM] beginning of stimulation
if (T > 70000) NE = 0; //end of stimulation
//Nitric oxide
//IF (T>30000) NO = 1D-3 //NO [mM]
//IF (T>70000) NO = 0
//Extracellular potassium
//IF (T>10000) K_o = 30
//IF (T>70000) K_o = 5
//Current
//IF (T>10000) I_stim = 5 //I_stim [pA] current injection
//IF (T>40000) I_stim = -5
//IF (T>70000) I_stim = 0
// For the time being, I just interested in Y[0] values (which is V_m actually)
fileY << Y[0];
fileY << "\t";
if ((ISTEP % ISTEPPRINT) == 0) {
// for (int i=0; i< N; i++) {
// fileY << Y[i];
// fileY << "\t";
// }
// fileY << endl;
// for (int i=0; i< N1; i++) {
// fileY1 << Y1[i];
// fileY1 << "\t";
// }
// fileY1 << endl;
//
//
//
// fileT << T;
// fileT << "\t";
}
FCN();
for (int i = 0; i < N; i++) {
Y[i] = Y[i] + TSTEP * YPRIME[i];
}
// disp(Yconcat(1))
if (ISTEP == MXSTEP)
break;
}
cout << "It is done!" << endl;
cout << Y[0] << endl;
fileY.close();
fileY1.close();
fileT.close();
return 0;
}
void FCN(void) {
V_m = Y[0];
d_L = Y[1];
f_L = Y[2];
p_K = Y[3];
q_1 = Y[4];
p_f = Y[5];
p_s = Y[6];
R_01 = Y[7];
R_10 = Y[8];
R_11 = Y[9];
Ca_u = Y[10];
Ca_r = Y[11];
Ca_i = Y[12];
Na_i = Y[13];
K_i = Y[14];
q_2 = Y[15];
P_SOC = Y[16];
Cl_i = Y[17];
h_IP3 = Y[18];
RS_G = Y[19];
RS_PG = Y[20];
G = Y[21];
IP3 = Y[22];
PIP2 = Y[23];
DAG = Y[24];
V_cGMP = Y[25];
cGMP = Y[26];
//-------------------------------------- Model equations ---------------------------------------------
//cGMP formation
V_cGMPbar = V_cGMPmax * (B5sGC * NO + pow(NO, 2)) / (A0sGC + A1sGC * NO + pow(NO, 2));
if ((V_cGMPbar - V_cGMP) >= 0) {
tausGC = 1 / (k3sGC * NO + kDact_deactsGC); //kDact_deactsGC different from original kDsGC to uncouple Km from time constants
} else {
tausGC = 1 / (kDact_deactsGC + k_2sGC);
}
V_cGMPprime = (V_cGMPbar - V_cGMP) / tausGC;
cGMPprime = V_cGMP - kpde_cGMP * cGMP * cGMP / (1e-3 + cGMP);
//Norepinephrine receptor
RS_Gprime = (k_rG * ksi_G * R_TG - (k_rG + k_pG * NE / (K_1G + NE)) * RS_G - k_rG * RS_PG);
RS_PGprime = NE * (k_pG * RS_G / (K_1G + NE) - k_eG * RS_PG / (K_2G + NE));
rho_rG = NE * RS_G / (ksi_G * R_TG * (K_1G + NE));
Gprime = k_aG * (delta_G + rho_rG)*(G_TG - G) - k_dG*G;
r_hG = alpha_G * Ca_i / (K_cG + Ca_i) * G;
IP3prime = r_hG / gamma_G * PIP2 - k_degG*IP3;
PIP2prime = -(r_hG + r_rG) * PIP2 - r_rG * gamma_G * IP3 + r_rG*PIP2_T;
DAGprime = r_hG / gamma_G * PIP2 - k_degG*DAG;
//Reversal potentials
E_Ca = RT_F / z_Ca * log(Ca_o / Ca_i);
E_Na = RT_F * log(Na_o / Na_i);
E_K = RT_F * log(K_o / K_i);
E_Cl = RT_F / z_Cl * log(Cl_o / Cl_i);
//Voltage dependent calcium current I_CaL
tau_d_L = 2.5 * exp(-pow((V_m + 40) / 30, 2)) + 1.15;
d_Lbar = 1.0 / (1 + exp(-(V_m) / 8.3));
d_Lprime = (d_Lbar - d_L) / tau_d_L;
f_Lbar = 1.0 / (1 + exp((V_m + 42.0) / 9.1));
tau_f_L = 65 * exp(-pow((V_m + 35) / 25, 2)) + 45;
f_Lprime = (f_Lbar - f_L) / tau_f_L;
if (V_m == 0) {
I_CaL = d_L * f_L * P_CaL * A_m * 1e6 * z_Ca * F * (Ca_i - Ca_o); //[pA]
} else {
I_CaL = d_L * f_L * P_CaL * A_m * 1e6 * V_m * pow(z_Ca * F, 2) / (R * temp)*(Ca_o - Ca_i * exp(V_m * z_Ca / (RT_F))) / (1 - exp(V_m * z_Ca / (RT_F))); //[pA]
}
//Delayed rectifier current I_K
p_Kbar = 1 / (1 + exp(-(V_m - V_1_2) / k));
tau_p_K = 61.49 * exp(-0.0268 * V_m);
p_Kprime = (p_Kbar - p_K) / tau_p_K;
q_bar = 1.0 / (1 + exp((V_m + 40) / 14));
q_1prime = (q_bar - q_1) / 371;
q_2prime = (q_bar - q_2) / 2884;
I_K = 1 * g_K * p_K * (0.45 * q_1 + 0.55 * q_2) * (V_m - E_K);
//Alpha-adrenoceptor-activated nonselective cation channel NSC
Po_NSC = Pmin_NSC + (1 - Pmin_NSC) / (1 + exp(-(V_m - 47.12) / 24.24));
if (V_m == 0) {
INa_NSC = 1 * (DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PNa_NSC * A_m * 1e6 * F * (Na_i - Na_o);
ICa_NSC = 1 * (0 * DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PCa_NSC * A_m * 1e6 * z_Ca * F * (Ca_i - Ca_o);
IK_NSC = 1 * (DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PK_NSC * A_m * 1e6 * F * (K_i - K_o);
} else {
INa_NSC = 1 * (DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PNa_NSC * A_m * 1e6 * V_m * pow(F, 2) / (R * temp)*(Na_o - Na_i * exp(V_m / RT_F)) / (1 - exp(V_m / RT_F));
ICa_NSC = 1 * (0 * DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PCa_NSC * A_m * 1e6 * V_m * pow(z_Ca * F, 2) / (R * temp)*(Ca_o - Ca_i * exp(V_m * z_Ca / RT_F)) / (1 - exp(V_m * z_Ca / RT_F));
IK_NSC = 1 * (DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PK_NSC * A_m * 1e6 * V_m * pow(F, 2) / (R * temp)*(K_o - K_i * exp(V_m / RT_F)) / (1 - exp(V_m / RT_F));
}
I_NSC = ICa_NSC + INa_NSC + IK_NSC;
//KATP current I_KATP
I_KATP = g_KATP * (V_m - E_K);
//Inward rectifier current I_K_i
g_maxK_i = G_K_i * pow(K_o, n_K_i);
I_K_i = g_maxK_i * (V_m - E_K) / (1 + exp((V_m - E_K) / 28.89));
//Calcium-activated potassium current I_KCa
if (V_m == 0) {
i1_KCa = 1e6 * P_BKCa * F * (K_i - K_o); //[pA]
} else {
i1_KCa = 1e6 * P_BKCa * V_m * F / RT_F * (K_o - K_i * exp(V_m / RT_F)) / (1 - exp(V_m / RT_F)); //[pA]
} //Mistry and Garland 1998
R_NO = NO / (NO + 200e-6);
R_cGMP = pow(cGMP, 2) / (pow(cGMP, 2) + pow(1.5e-3, 2));
V_1_2_KCa = -41.7 * log10(Ca_i) - 128.2 - dV_1_2_KCa_NO * R_NO - dV_1_2_KCa_cGMP*R_cGMP;
p_obar = 1 / (1 + exp(-(V_m - V_1_2_KCa) / 18.25));
p_fprime = (p_obar - p_f) / tau_pf;
p_sprime = (p_obar - p_s) / tau_ps;
P_KCa = 0.17 * p_f + 0.83 * p_s;
I_KCa = A_m * N_BKCa * i1_KCa * P_KCa;
//Store operated non-selective cation channel
P_SOCbar = 1 / (1 + pow(Ca_u / K_SOC, H_SOC));
P_SOCprime = (P_SOCbar - P_SOC) / tau_SOC;
I_SOCCa = 1 * P_SOC * g_SOCCa * (V_m - E_Ca);
I_SOCNa = 1 * P_SOC * g_SOCNa * (V_m - E_Na);
I_SOC = I_SOCCa + I_SOCNa;
//Chloride currents
alpha_Cl = pow(cGMP, 3.3) / (pow(cGMP, 3.3) + pow(6.4e-3, 3.3));
P_Cl = pow(Ca_i, 2) / (pow(Ca_i, 2) + pow(365e-6, 2)) * 0.0132 + pow(Ca_i, 2) / (pow(Ca_i, 2) + pow(400e-6 * (1 - alpha_Cl * 0.9), 2)) * alpha_Cl;
I_Cl = P_Cl * g_Cl * C_m * (V_m - E_Cl);
//IP3 receptor current
h_IP3prime = k_onIP3 * (K_inhIP3 - (Ca_i + K_inhIP3) * h_IP3);
I_IP3 = I_IP3bar * pow(IP3 / (IP3 + K_IP3) * Ca_i / (Ca_i + K_actIP3) * h_IP3, 3)*(Ca_u - Ca_i) * z_Ca * F*vol_Ca;
//Calcium-induced calcium release
R_00 = 1 - R_01 - R_10 - R_11;
R_10prime = Kr1 * pow(Ca_i, 2) * R_00 - (K_r1 + Kr2 * Ca_i) * R_10 + K_r2*R_11;
R_11prime = Kr2 * Ca_i * R_10 - (K_r1 + K_r2) * R_11 + Kr1 * pow(Ca_i, 2) * R_01;
R_01prime = Kr2 * Ca_i * R_00 + K_r1 * R_11 - (K_r2 + Kr1 * pow(Ca_i, 2)) * R_01;
I_up = I_upbar * Ca_i / (Ca_i + K_mup);
I_tr = (Ca_u - Ca_r) * (2 * F * vol_u) / tau_tr;
I_rel = (pow(R_10, 2) + R_leak) * (Ca_r - Ca_i) * (2 * F * vol_r) / tau_rel;
Ca_uprime = (I_up - I_tr - I_IP3) / (2 * F * vol_u);
Ca_rprime = (I_tr - I_rel) / (2 * F * vol_r) / (1 + CSQNbar * K_CSQN / pow((K_CSQN + Ca_r), 2));
//Ca buffering and cytosolic material balance
S_CM = S_CMbar * K_d / (K_d + Ca_i);
CaCM = S_CMbar - S_CM;
I_PMCA = I_PMCAbar * Ca_i / (Ca_i + K_mPMCA);
//NaK pump
I_NaK = pow(Q_10_NaK, ((temp - 309.15) / 10)) * C_m * I_NaKbar * ((pow(K_o, 1.1)) / (pow(K_o, 1.1) + (pow(K_mK, 1.1)))
*(pow(Na_i, 1.7)) / ((pow(Na_i, 1.7))+(pow(K_mNa, 1.7)))) *(V_m + 150) / (V_m + 200);
Fi_F = exp(gamma2 * V_m * F / (R * temp));
Fi_R = exp((gamma2 - 1) * V_m * F / (R * temp));
I_NCX = 1 * (1 + 0.55 * cGMP / (cGMP + (45e-3))) * g_NCX * (pow(Na_i, 3) * Ca_o * Fi_F - pow(Na_o, 3) * Ca_i * Fi_R) / (1 + d_NCX * (pow(Na_o, 3) * Ca_i + pow(Na_i, 3) * Ca_o));
I_NaKCl_Cl = (1 + 7 / 2 * cGMP / (cGMP + 6.4e-3))*(-A_m * L_cotr * R * temp * z_Cl * F * log(Na_o / Na_i * K_o / K_i * pow(Cl_o / Cl_i, 2)));
I_Catotm = I_SOCCa + I_CaL - 2 * I_NCX + I_PMCA + ICa_NSC + I_MCa;
Ca_iprime = -(I_Catotm + I_up - I_rel - I_IP3) / (2 * F * vol_Ca) / (1 + S_CMbar * K_d / (pow(K_d + Ca_i, 2)) + B_Fbar * K_dB / (pow(K_dB + Ca_i, 2)));
I_Natotm = -0.5 * I_NaKCl_Cl + I_SOCNa + 3 * I_NaK + 3 * I_NCX + INa_NSC + I_MNa;
Na_iprime = -(I_Natotm) / (F * vol_i);
I_Ktotm = -0.5 * I_NaKCl_Cl + I_K + I_KCa + I_K_i + IK_NSC + I_KATP - 2 * I_NaK + I_MK;
K_iprime = -(I_Ktotm) / (F * vol_i);
I_Cltotm = I_NaKCl_Cl + I_Cl;
Cl_iprime = -(I_Cltotm) / (z_Cl * F * vol_i);
//Transmembrane potential
V_mprime = -1 / C_m * (-I_stim + I_Cl + I_SOC + I_CaL + I_K + I_KCa + I_K_i + I_M + I_NCX + I_NaK + I_PMCA + I_NSC + I_KATP);
//YPRIME = zeros(1, N);
YPRIME[0] = V_mprime;
YPRIME[1] = d_Lprime;
YPRIME[2] = f_Lprime;
YPRIME[3] = p_Kprime;
YPRIME[4] = q_1prime;
YPRIME[5] = p_fprime;
YPRIME[6] = p_sprime;
YPRIME[7] = R_01prime;
YPRIME[8] = R_10prime;
YPRIME[9] = R_11prime;
YPRIME[10] = Ca_uprime;
YPRIME[11] = Ca_rprime;
YPRIME[12] = Ca_iprime;
YPRIME[13] = Na_iprime;
YPRIME[14] = K_iprime;
YPRIME[15] = q_2prime;
YPRIME[16] = P_SOCprime;
YPRIME[17] = Cl_iprime;
YPRIME[18] = h_IP3prime;
YPRIME[19] = RS_Gprime;
YPRIME[20] = RS_PGprime;
YPRIME[21] = Gprime;
YPRIME[22] = IP3prime;
YPRIME[23] = PIP2prime;
YPRIME[24] = DAGprime;
YPRIME[25] = V_cGMPprime;
YPRIME[26] = cGMPprime;
//Non state variables
Y1[0] = I_CaL;
Y1[1] = I_K;
Y1[2] = I_K_i;
Y1[3] = I_NSC;
Y1[4] = I_KCa;
Y1[5] = I_up;
Y1[6] = I_rel;
Y1[7] = I_PMCA;
Y1[8] = I_NCX;
Y1[9] = I_NaK;
Y1[10] = INa_NSC;
Y1[11] = ICa_NSC;
Y1[12] = IK_NSC;
Y1[13] = I_SOCCa;
Y1[14] = I_SOCNa;
Y1[15] = I_Cl;
Y1[16] = I_NaKCl_Cl;
Y1[17] = I_IP3;
Y1[18] = I_tr;
Y1[19] = NE;
Y1[20] = I_KATP;
Y1[21] = I_stim;
Y1[22] = V_cGMPbar;
Y1[23] = NO;
Y1[29] = I_Catotm;
Y1[30] = I_Natotm;
Y1[31] = I_Cltotm;
Y1[32] = I_Ktotm;
}
我建议您开始查看两个程序开始分歧的地方，您可能会发现一个有问题的操作，即在 g++ 返回 NaN 的情况下 MS 编译器返回 0 或 Inf。
可以通过允许代码抛出而不是默默地通知您来帮助您调试此问题。有一些特定于 Windows 的 API 可以控制 FPU 的配置，我认为默认情况下它们并未设置为符合 IEEE-754 标准。
_controlfp
。更多详情请参见
此页面。
Y1
数组的末尾，显然在 g++ 中它们走在 YPRIME 的顶部。当我将
Y1的声明更改为
Y1[33]以腾出空间时，我得到了 -59.1481 作为 g++ 的结果。
Y1[29] = I_Catotm;
Y1[30] = I_Natotm;
Y1[31] = I_Cltotm;
Y1[32] = I_Ktotm;
原答案：
如果您能够在调试器中运行两个版本，那么最好的选择是运行两个版本，直到它们仍然一致 (59.3993)，然后并排单步执行两个调试器，观察每个语句的输出看看他们的分歧在哪里。一旦您看到结果不同的地方，就应该更清楚代码的不同之处。
我曾经使用这种方法在内部数学库的主要重构/优化中发现错误。
Y[0]
) 在第 425 行演变：
Y[i] = Y[i] + TSTEP * YPRIME[i]。导数
YPRIME[i]在
FCN()中计算。由于我们只关注
Y[0]，因此我们对
V_mprime的值感兴趣 在第 636 行计算为
-1 / C_m * (-I_stim + I_Cl + I_SOC + I_CaL + I_K + I_KCa + I_K_i + I_M + I_NCX + I_NaK + I_PMCA + I_NSC + I_KATP)。现在，与 MSVC 的分歧可能始于上面表达式中的任何变量:-) 我建议你在第 636 行上方放置一条调试输出语句，并打印所有涉及的变量。然后在两个平台上运行程序并报告两个输出，以便我们可以关注差异。