## Limit Cycle Continuation

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

• 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
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
2013-01-10

Hi Drew,

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
2013-01-10

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

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