Limit Cycle Fitzhugh Nagumo

Help
Greg
2012-07-04
2013-11-04
  • Greg

    Greg - 2012-07-04

    Hi,

    I'm trying to do a simple bifurcation analysis to
    the standard FN model. Here is my code:

    from PyDSTool import *
    from PyDSTool.Toolbox import phaseplane as pp

    #==========SET UP SYSTEM==========================                                            

    pars = {'b' : 1,'bet':1,'eps' : 0.01,'gam' : 2}
    icdict={'x1':-0.3,'x2':-0.7}

    # Set up model                                                                                
    x1eq = '(x1-x1**3-x2)/eps'
    x2eq = 'gam*x1-bet*x2 + b'

    DSargs = args(name='fn')
    DSargs.pars = pars
    DSargs.varspecs = {'x1': x1eq, 'x2': x2eq}
    DSargs.ics = icdict         #initial conditions                                               
    DSargs.tdomain =      # set the range of integration.                                   
    DSargs.xdomain = {'x1': , 'x2': }
    DSargs.fnspecs = {'Jacobian': (,
                                    """[,                               
                                        ]""")}
    ode = Generator.Vode_ODEsystem(DSargs)

    #============PHASEPLANE=================================                                      

    figure(4)
    pp.plot_PP_vf(ode, 'x1', 'x2', scale_exp=-1.2,N=20)

    fp_coord = pp.find_fixedpoints(ode, n=4, eps=1e-8)
    #fp = pp.fixedpoint_2D(ode, Point(fp_coord), eps=1e-8) # not working :(                       

    nulls_x, nulls_y = pp.find_nullclines(ode, 'x1', 'x2', n=3, eps=1e-8, max_step=0.1,           \
                      fps=)

    #======BIFURCATIONS===========================================                                

    PC = ContClass(ode)

    PCargs = args(name='EQ1', type='EP-C')
    bpar='gam'
    PCargs.freepars =
    PCargs.StepSize = 1e-2
    PCargs.MaxNumPoints = 50
    PCargs.MaxStepSize = 1e-1
    PCargs.LocBifPoints = 'all'
    PCargs.verbosity = 2
    PCargs.SaveEigen    = True                       # to tell unstable from stable branches      
    PC.newCurve(PCargs)

    print 'Computing curve…'
    start = clock()
    PC.forward()
    PC.backward()
    print 'done in %.3f seconds!' % (clock()-start)

    PC.display((bpar,'x2'),stability=True, figure=2)

    PCargs = args(name='LC1', type='LC-C')
    PCargs.freepars =
    PCargs.initpoint = 'EQ1:H1'
    PCargs.StepSize = 1e-2
    PCargs.MaxNumPoints = 250
    PCargs.MaxStepSize = 1e-1
    PCargs.LocBifPoints = 'LPC'
    PCargs.verbosity = 2
    PCargs.SolutionMeasures = 'avg'
    #PCargs.SaveEigen    = True                       # to tell unstable from stable branc        
    PC.newCurve(PCargs)

    It detects the Hopfbifurcation, but the limit cycle continuation fails with the
    following error:

    $ python bif1.py
    Computing curve…
    Maximum test function iterations reached.

    Checking…
      |q| = 1.000000
      |<p,q> - 1| =  1.11086097034e-16
      |Aq - iwq| = 0.000000
      |A*p + iwp| = 0.000000

    H Point found

    0 :
    x2  =  -0.384885697323
    gam  =  2.41077655768
    x1  =  -0.574456265102

    Eigenvalues =

         (-0.000000,15.494440)
         (-0.000000,-15.494440)

    w = 15.494439516254179
    l1 = -25.671139504469721

    B Point found

    0 :
    x2  =  0.899999982327
    gam  =  0.0768795307873
    x1  =  -1.30073657641

    done in 0.410 seconds!
    /home/whir/Dropbox/p53/py/FN/auto_temp/fn_vf.c: In function ‘vfieldfunc’:
    /home/whir/Dropbox/p53/py/FN/auto_temp/fn_vf.c:51:16: error: invalid type argument of unary ‘*’ (have ‘int’)
    Traceback (most recent call last):
      File "bif1.py", line 94, in <module>
        PC.newCurve(PCargs)
      File "/home/whir/progs/PyDSTool/PyCont/ContClass.py", line 221, in newCurve
        self.loadAutoMod()
      File "/home/whir/progs/PyDSTool/PyCont/ContClass.py", line 432, in loadAutoMod
        self.compileAutoLib()
      File "/home/whir/progs/PyDSTool/PyCont/ContClass.py", line 717, in compileAutoLib
        raise RuntimeError
    RuntimeError

    Has anyone an idea what goes wrong?

    Thx in advance,
    Gregor

     
  • Greg

    Greg - 2012-07-04

    ..the PyCont tests like  "PyCont_MorrisLecar_TypeII.py" all go fine!

     
  • Rob Clewley

    Rob Clewley - 2012-07-04

    Yes, I suspect it's your use of ** for power in your system's RHS. You should use the pow(x,p) function instead. This is explained in the web pages about functional specifications for your systems. There's even a string conversion function somewhere that will take care of that automatically in the future, e.g. if you have a large number of legacy strings to convert

     
  • Greg

    Greg - 2012-07-05

    Hey,

    yes thank you very much! I thought it was valid because the fixed point
    search went so well.. I want to say that I really appreciate this project,
    I find it very promising and hopefully it  can replace xppaut in a convenient way,
    thanks a lot! I try to spread it here in Berlin in my scientific community!

     

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks