Limit Cycle Continuation

Help
whir
2013-01-09
2013-11-04
  • whir
    whir
    2013-01-09

    Hello,

    I'm trying to plot the limit cycles in the standard Fitzhugh Nagumo model.
    The FP bifurcation analysis goes well, but it fails to actually follow the emerging
    limit cycle oscillations. I only got 'MX' Points. That also happens with other more complicated
    systems, although I could swear that it all was working very well some time ago..

    Here is the code:

    """ FitzHugh Nagumo """

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

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

    pars = {'b' : 1,'bet':1,'eps' : 0.1,'I' : 1}
    icdict={'x1':-1.3,'x2':-1.2}

    # Set up model                                                                                                       
    x1eq = '(x1-pow(x1,3)/3.-x2) + I'
    x2eq = 'eps*(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(1)
    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=[\
    fp_coord])

    # plot the nullclines                                                                                                
    plt.plot(nulls_x, nulls_x, 'b')
    plt.plot(nulls_y, nulls_y, 'g')

    #Traj                                                                                                                

    traj = ode.compute('pol')              # integrate ODE                                                               
    pts  = traj.sample(dt=0.01)                      # Data for plotting                                                 
    plt.plot(pts, pts,'r-')
    plt.plot(fp_coord, fp_coord,'c*',markersize=10)

    figure(2)
    plt.plot(pts, pts,'c-')

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

    # ode.set(pars = {'b' : 0.7,'bet':0.8,'eps' : 0.01,'eps' : 2.0})                                                     
    PC = ContClass(ode)

    PCargs = args(name='EQ1', type='EP-C')
    bpar='I'
    PCargs.freepars =
    PCargs.StepSize = 1e-3
    PCargs.MaxNumPoints = 150
    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.forward()
    PC.backward()
    print 'done in %.3f seconds!' % (clock()-start)

    PCargs = args(name='LC1', type='LC-C')
    PCargs.freepars =
    PCargs.initpoint = 'EQ1:H1'
    PCargs.MinStepSize = 0.0001
    PCargs.MaxStepSize = .001
    PCargs.StepSize = 0.001
    PCargs.MaxNumPoints = 620
    PCargs.LocBifPoints = 'LPC'
    PCargs.verbosity = 2
    PCargs.NumSPOut = 10000
    PCargs.SolutionMeasures = 'all'
    PCargs.SaveEigen    = True                       # to tell unstable from stable branc                                
    PC.newCurve(PCargs)

    print 'Computing LC curve…'
    start = clock()
    PC.backward()
    PC.forward()

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

    show()

    I can not find proper parameters to circumvent the numerical issues!
    Any Ideas?

    Greets,
    Greg

     
  • Drew LaMar
    Drew LaMar
    2013-01-09

    Hi, Greg.  Hmmmm, I was able to get your script to work using my version of PyDSTool.  Here's a link to the resulting bifurcation diagram so you can see if it looks proper based on your memory:

     

    I've made a few changes to PyCont in my local version, mostly in regard to how information is passed back to PyCont from auto, but if you would like access to my version, just let me know and I'll make them available.

    - Drew

     
  • whir
    whir
    2013-01-10

    Hi Drew,

    thanks for your quick reply!

    Assuming that you got my script working without (parameter) alterations,
    I will try to make a fresh install of pydstool.I got some updates from rob a few
    month ago, i.e. phaseplane.py is altered for sure.

    I could also first run an appropriate
    test script.. /tests/PyCont_vanDerPol.py seems to work, although it is hard
    to see for me what they are really doing there.
    This is the output on my machine:
    (uh I have no URL for my image)

    Greets,
    Greg

     
  • whir
    whir
    2013-01-10

    quick update… /tests/PyCont_MorrisLecar_TypeII.py runs flawlessly and plots
    nice limit cycles in the bifurcation diagram!

     
  • whir
    whir
    2013-01-11

    Ok..as far as I figured it out, I had to remove all *auto* related folders and files in the active directory.
    There seemed to be sth. wrong. Thus with that reinitialization of auto, everything goes fine (as remembered) !!

    Maybe there should be some clean up routine implemented?

    Greets,
    Greg

     
  • Drew LaMar
    Drew LaMar
    2013-01-11

    Ah, yes.  I should have thought of that.  On my machine I use a Makefile with the command "make clean" that will remove all compiled dynamic libraries and the auto_temp directory.  It might be nice as you've mentioned to just make this a method in PyCont so that this can be done within Python.

    Glad it works now!

    - Drew