What is the mistake in this C implementation of Fitzhugh-Nagumo model?

Question

I have implemented a MATLAB code for Fitzhugh-Nagumo model and got the plots, but when I translated it to C code like below, its not giving me proper output. Specifically neurons are not spiking. Can any one point out the mistake?

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
double rand2(void);
int main()
{
    int n,i;
    double dt,r[120000],g[120000],rr ;// 20min = 20*60/dt = 120000
    double a=0.08,b=0.7,c=0.8,gg;
    //FILE * temp = fopen("fhn.dat", "w");
    srand(time(NULL));
    //T=200*60;
    dt=0.01;
    r[0]=-0.212002413260425;//initial values
    g[0]=-0.869601930608358;

    for(n=1;n<=120000;n++)
    {
    r[n]=(dt)*((a*g[n-1])+b-(c*r[n-1]))+r[n-1];
    rr=rand2(); //uniform random number between 0 and 1
    gg=pow(g[n-1],3); // g[n-1]^3
    g[n]=(dt)*(g[n-1]-gg-r[n-1]+rr)+g[n-1];
    //printf("rr %f\n",rr);
    }

    FILE *gnuplot = popen("gnuplot -persistent", "w");
    fprintf(gnuplot, "set title \"Fitzhugh-Nagumo\" \n");
    fprintf(gnuplot, "plot '-' with lines \n");
    for (i=0;i<=120000-1;i++)
    {
        fprintf(gnuplot, "%g %g\n", r[i],g[i]);
        //fprintf(temp, "%lf %lf \n",r[i], g[i]); //Write the data to a temporary file

    }
    fprintf(gnuplot, "e\n");
    fflush(gnuplot);
    pclose(gnuplot);
    //fclose(temp);
    return 0;
}

double rand2()
{
    return (double)rand() / (double)RAND_MAX ;
}

NOTE : I expect an output like below plot :

UPDATE : Link to MATLAB Code.

But What I get is like :

What debugging have you done? Where in the process do the results diverge? What outputs did you get, and what were you expecting? — jonrsharpe, Aug 30 '15 at 08:41
I think `r[0]`'s initial value should be `-0.212002443260425` instead ;-) Come on man, why do you think we can spot a problem here, and why even post these graphs!? — meaning-matters, Aug 30 '15 at 09:12
To start with, `for(n=1;n<=120000;n++)` results in an out-of-bound access. Newer versions of GCC even produces a helpful warning for this. Bump up your warning levels. — T.C., Aug 30 '15 at 09:16
And in general, start by examining why your "new" version doesn't have anything in the vicinity of -0.2/-0.8. — T.C., Aug 30 '15 at 09:18
@T.C., good point, also answers the question of @meaning-matters;) As the data explicitly starts at those coordinates, the plotting is definitely off, so the computation could even be OK. And @dexterdev, is there a reason why you have a few `floats` thrown in there, like the third power of your `double` `gg`? Doesn't sound right, nor practical. — Andras Deak -- Слава Україні, Aug 30 '15 at 09:45
Sorry there is an error in the initial points , I will edit it. — dexterdev, Aug 30 '15 at 09:51
@dexterdev your initial `g` is still with the wrong sign... Makes me wonder: what is the code you posted if that's not what you use? — Andras Deak -- Слава Україні, Aug 30 '15 at 09:58
@AndrasDeak : I have corrected everything to double. About `g` I was also wondering from yesterday? My code is the same. But the plot gives different. When I checked by seeing the .dat file of `r` and `g` , I noted this. But don't know why? — dexterdev, Aug 30 '15 at 10:07
@4566976 : Uploaded the matlab code link in pastebin. http://pastebin.com/fPdJSZrc — dexterdev, Aug 30 '15 at 11:22
My recommendation: Replace rand in both cases with a identical deterministic random number stream, then run both codes side-to-side in the debugger. — Daniel, Aug 30 '15 at 12:00

score 1 · Accepted Answer · answered Aug 30 '15 at 12:34

1

r[n] is wrongly parenthesized. Replace

r[n]=(dt)*((a*g[n-1])+b-(c*r[n-1]))+r[n-1];

with

r[n]=(dt)*(a*(g[n-1]+b-c*r[n-1]))+r[n-1];

and the result looks like expected.

answered Aug 30 '15 at 12:34

4566976

2,419
1
10
14

Thank you very much. Still the initial stick initial standing as in wrong output. LINK: http://i.imgur.com/2fgBQ6G.png?1?8436 – dexterdev Aug 30 '15 at 12:41
Thanks for the hard work. I will figure it out. Just gave you the output. – dexterdev Aug 30 '15 at 12:45
1

@dexterdev I can't reproduce this on my machine and can't see why this could be, from your posted code. – 4566976 Aug 30 '15 at 13:37

What is the mistake in this C implementation of Fitzhugh-Nagumo model?

1 Answers1