## Limit Cycle Fitzhugh Nagumo document.SUBSCRIPTION_OPTIONS = { "thing": "thread", "subscribed": false, "url": "subscribe", "icon": { "css": "fa fa-envelope-o" } };

Help
Greg
2012-07-04
2013-11-04
• 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
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?

Gregor

• Greg - 2012-07-04

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

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