On Jul 4, 2013, at 20:14 , Berend Hasselman wrote: > > On 04-07-2013, at 19:56, peter dalgaard <pda...@gmail.com> wrote: > >> >> On Jul 4, 2013, at 19:11 , Berend Hasselman wrote: >> >>> >>> On 04-07-2013, at 18:42, Jannetta Steyn <janne...@henning.org> wrote: >>> >>>> Hi Ben and others >>>> >>>> I don't quite know how to explain the "doesn't work" in more detail without >>>> any visual aid. >>> >>> You said that R got into an indefinite loop, whatever that maybe. >>> >>>> When you run the two scripts it is easy to see the >>>> difference. MatLab produces a line on x= -55. This is what I expect - a >>>> more or less straight line. R on the other hand the result drops from -55 >>>> (the initial value) to -80 and then goes up to -71.37092. >>>> >>>> The two scripts have exactly the same equations. >>> >>> >>> I don't think so. >>> In the R script you have >>> >>> init = c(v_axon_AB=-55,mNa_axon_AB=1,hNa_axon_AB=0,mK_axon_AB=1) >>> >>> That is not the same as in your Matlab script. To make them the same you >>> should replace the line with >>> >>> init = c(v_axon_AB=-55,mNa_axon_AB=0,hNa_axon_AB=1,mK_axon_AB=0) >>> >>> Using this line gives quite different results. But not the same as Matlab. >>> >>> It seems that you ought to carefully check all your equations. >>> I don't know enough about Matlab/Octave syntax to determine if the results >>> of the two systems are identical. >> >> Also, the line >> >> iLeak_axon_AB <- gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB) >> >> differs from the Matlab code. Changing it gave me all-NA results, though. >> And the with() construct assumes that init and parms are named vectors; >> parms is not, so the values are taken from the global environment. It is not >> clear whether that makes a difference. >> > > Same NA results here. Making parms a named vector makes no difference but is > is certainly necessary to do that. > > In that case there is probably an inconsistency in the file simulate.m near > the end: > The last non empty line before the line out = … reads > > iLeak_axon = gLeak_axon_AB.*(v_axon_AB-ELeak_axon_AB); > > and this doesn't agree with the line > > iLeak_axon = ELeak_axon_AB*(v_axon_AB-ELeak_axon_AB); > > in the function xdot. > > Final question for the OP: if the model is supposed to be dynamic, isn't it > suspicious that Matlab gives a constant result?
Also, in this block gNa_axon_AB=300e-3 gK_axon_AB=52.5-3 gLeak_axon_AB=0.0018e-3 the middle line looks like a typo for 52.5e-3 and the 3rd line looks very low -- googling suggests values like (120, 36, 0.3) mS/cm^3, which would be more consistent with a value around 1e-3. > > Berend > >> >>> >>> Berend >>> >>>> I have even named the >>>> variables the same in the two scripts and copied the equations across to >>>> make sure I haven't made any typos. >>>> >>>> Can one attach images to posts? I'll try. The flatline image is the plot >>>> from MatLab and the other is the plot from R. >>>> >>>> Thanks >>>> Jannetta >>>> >>>> >>>> On 4 July 2013 16:52, Ben Bolker <bbol...@gmail.com> wrote: >>>> >>>>> Berend Hasselman <bhh <at> xs4all.nl> writes: >>>>> >>>>>> >>>>>> >>>>>> On 04-07-2013, at 17:15, Jannetta Steyn <jannetta <at> henning.org> >>>>> wrote: >>>>>> >>>>>>> Hi folks >>>>>>> >>>>>>> I have implemented a model of a neuron using Hodgkin Huxley equations >>>>> in >>>>>>> both R and MatLab. My preference is to work with R but R is not giving >>>>> me >>>>>>> the correct results. I also can't use ode45 as it just seems to go >>>>> into an >>>>>>> indefinite loop. However, the MatLab implementation work fine with >>>>>>> the same >>>>>>> equations, parameters and initial values and ode45. Below is my R and >>>>>>> MatLab implementations. >>>>>>> >>>>>> >>>>>> No problem in running your R file. Have plot. >>>>>> (Mac mini Core2Duo 2.66Ghz running R 3.0.1 Patched >>>>>> (2013-06-19 r62992) on Mac OS X 10.8.4: 16.5 seconds.) >>>>>> >>>>>> Trying to run your Matlab file in Octave. >>>>>> No success yet because of unavailable ode45. >>>>> >>>>> I'm impressed that you (BH) went to the trouble of checking >>>>> on this vague "doesn't work" question. If you want to go farther >>>>> I think you can get ode45 for octave by installing the odepkg >>>>> package: >>>>> >>>>> http://octave.sourceforge.net/odepkg/index.html >>>>> >>>>> ______________________________________________ >>>>> R-help@r-project.org mailing list >>>>> https://stat.ethz.ch/mailman/listinfo/r-help >>>>> PLEASE do read the posting guide >>>>> http://www.R-project.org/posting-guide.html >>>>> and provide commented, minimal, self-contained, reproducible code. >>>>> >>>> >>>> >>>> >>>> -- >>>> >>>> =================================== >>>> Web site: http://www.jannetta.com >>>> Email: janne...@henning.org >>>> =================================== >>>> <Rplot01.png><MatPlot01.png>______________________________________________ >>>> R-help@r-project.org mailing list >>>> https://stat.ethz.ch/mailman/listinfo/r-help >>>> PLEASE do read the posting guide >>>> http://www.R-project.org/posting-guide.html >>>> and provide commented, minimal, self-contained, reproducible code. >>> >>> ______________________________________________ >>> R-help@r-project.org mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help >>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >>> and provide commented, minimal, self-contained, reproducible code. >> >> -- >> Peter Dalgaard, Professor, >> Center for Statistics, Copenhagen Business School >> Solbjerg Plads 3, 2000 Frederiksberg, Denmark >> Phone: (+45)38153501 >> Email: pd....@cbs.dk Priv: pda...@gmail.com >> >> >> >> >> >> >> >> > -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
