Mon code C++ fonctionne parfaitement sur le compilateur MS C++ mais me donne NaN sur le compilateur g ++. Pourquoi?

Question

Je simule un modèle biologique impliquant plusieurs (27) équations différentielles ordinaires rigides utilisant C++. Mon programme s'exécute parfaitement sous compilateur d'expression MS C++ 2010 mais échoue sous le compilateur g ++ (NetBeans 6.8, Ubuntu 10.04 LTS). Le problème est que certaines variables deviennent NaN. Ci-dessous sont les valeurs de la Vm variable après chaque étape du programme sous le compilateur g ++:

-59,4 -59,3993 -59,6081 100,081 34,6378 -50392,8 nan nan nan nan nan nan nan nan nan nan nan nan nan nan Nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan

Et voici la sortie du même code, sans aucun changement dans le compilateur de la MS C:

-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 - 

J'ai seulement utilisé les bibliothèques "cmath" et "fstream".

Où est le problème? Le code dans les deux scénarios est exactement le même.

Edit 1:

OK les gars, voici le code complet:

#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; 
} 

Answer 1

EDIT: Bien sûr, ces quatre lignes vont à la fin du tableau Y1 et apparemment en g ++ ils marchent sur YPRIME. Lorsque j'ai changé la déclaration de Y1 en Y1[33] pour faire de la place j'ai obtenu -59.1481 comme résultat avec g ++.

Y1[29] = I_Catotm; 
Y1[30] = I_Natotm; 
Y1[31] = I_Cltotm; 
Y1[32] = I_Ktotm; 

réponse originale:

Si vous avez la possibilité d'exécuter les deux versions dans un débogueur, votre meilleur pari est d'exécuter les deux versions jusqu'au point où ils sont d'accord encore (59.3993), puis les deux débogueurs côte-à-côte, observant la sortie de chaque instruction pour voir où ils divergent. Une fois que vous voyez où les résultats sont différents, il devrait être beaucoup plus clair ce que le code fait différemment.

Une fois, j'ai utilisé cette approche pour trouver un bogue dans une refactorisation/optimisation majeure d'une bibliothèque mathématique interne.

Answer 2

Windows et Linux définir différents paramètres par défaut en virgule flottante. < fenv.h> peut vous aider à déboguer ceci en permettant au code de lancer au lieu de vous silencieusement. Il existe des API spécifiques à Windows qui contrôlent la configuration de la FPU, et je pense que par défaut, ils ne sont pas configurés pour la conformité IEEE-754.

Sous Windows, vous devez appeler le _controlfp si vous souhaitez utiliser le mode standard IEEE-754.

Plus de détails sur this page.

Answer 3

Je suppose que votre code vous donne une réponse avec votre compilateur MS C, mais je suis sceptique qu'il fonctionne parfaitement. NaNs (Not A Number) sont le résultat de fonctions de calcul hors de leurs plages. Puisque vous ne mentionnez pas comment vous résolvez votre système rigide, je n'ai aucune idée si vous prenez le journal d'un nombre négatif, ou quelque autre arithmétique d'écureuil, mais je suis sûr que vous faites la même arithmétique dans votre code de Windows. Ce n'est pas parce que le code tente de s'embrouiller que cela fonctionne correctement. Je suggérerais que vous commenciez à regarder où les deux programmes commencent à diverger, et vous trouverez probablement une opération discutable où le compilateur MS retournant 0 ou Inf dans le cas où g ++ renvoie NaN.

Answer 4

La solution à votre problème n'est pas facile. La variable que vous enregistrez (Y[0]) évolue à la ligne 425: Y[i] = Y[i] + TSTEP * YPRIME[i]. La dérivation YPRIME[i] est calculée en FCN(). Puisque nous regardons seulement Y[0], nous sommes intéressés par la valeur de V_mprime qui est calculée à la ligne 636 en tant que -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). Maintenant, la divergence de MSVC peut commencer dans l'une des variables de l'expression ci-dessus :-) Je vous suggère de mettre une instruction de sortie de débogage juste au-dessus de la ligne 636, et d'imprimer toutes les variables impliquées. Ensuite, exécutez le programme sur les deux plates-formes et rapportez les deux sorties afin que nous puissions nous concentrer sur les différences.