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