dupont - 2013-05-28

Hi,

I am trying to move from Matcont to pydstool and I am not sure about the results from pydstool. I have been using the same equations and the pydstool seem to be missing some Generalized Hopf bifircation points. There should be two additional GH points at the x-axis values roughly -1 and 3 respectively. Has anyone gotten any similar issue?

Thank you for your help.

Here is my code

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
#! /opt/local/bin/python2.6
# -*- coding: utf-8 -*-
"""
Created on Wed May 22 15:03:55 2013
"""
class bcolors:
    HEADER = '\033[95m'
    OKBLUE = '\033[94m'
    OKGREEN = '\033[92m'
    WARNING = '\033[93m'
    FAIL = '\033[91m'
    ENDC = '\033[0m'
    def disable(self):
        self.HEADER = ''
        self.OKBLUE = ''
        self.OKGREEN = ''
        self.WARNING = ''
        self.FAIL = ''
        self.ENDC = ''
        
from PyDSTool import *
from scipy.special import *
import sys
pars = {'Jee': 10.0,
        'Jei': -12.0,
        'Jie': 10.0,
        'Jii': -10.0,
        'tauE': 0.375,
        'tauI': 1.0,
        'Ii': 1.0,
        'Ee': 2.5,
        'alpha': 4.0}
auxfndict = {'S': (['v'], '1.0/(1.0 + exp(-v))'), \
           'L': (['v'], 'log(1.0+exp(v))'), \
             'Q': (['v'],'log(1.0+exp(v*10.0-.0))/10.20')   
            }
   
icdict = {'E': 0.0,
          'I': 0.0}
Estr = '(-E+S(Jee*E+Jei*I+Ee))/tauE'
Istr = '(-I+S(Jie*E+Jii*I+Ii))/tauI'   
###############################################################################
# parameter
###############################################################################
DSargs = args(name='ISN97')
DSargs.pars = pars
DSargs.varspecs = {'E': Estr, 'I': Istr}
DSargs.fnspecs = auxfndict
DSargs.ics = icdict
###############################################################################
# find IC
###############################################################################
DSargs.tdomain = [0,200]
DSargs.pdomain = {'Ii': [-4,7], 'Ee': [-2, 6]}
DSargs.algparams = {'init_step' :0.01, 'strictopt':False}
testDS = Generator.Vode_ODEsystem(DSargs)
print 'Integrating...'
start = clock()
testtraj = testDS.compute('testDS')
print '  ... finished in %.3f seconds.\n' % (clock()-start)
if 1:
    plotData=testtraj.sample(dt=0.01)
    mline=plot(plotData['t'],plotData['E'])
    show()
###############################################################################
# perform coninuation
###############################################################################
# Set up continuation class
PyCont = ContClass(testDS)
PCargs = args(name='EQ1', type='EP-C')
PCargs.freepars = ['Ii']
PCargs.StepSize = 1e-4
PCargs.MinStepSize = 1e-4
PCargs.MaxNumPoints = 2000
PCargs.MaxStepSize = 5e-2
PCargs.LocBifPoints = 'all'
PCargs.verbosity = 2
PCargs.SaveEigen = True
PCargs.pdomain = {'Ii': [-4,7], 'Ee': [-4, 6]}
PyCont.newCurve(PCargs)
print bcolors.WARNING + '--> Computing EQ1...'+ bcolors.ENDC 
start = clock()
PyCont['EQ1'].forward()
PyCont['EQ1'].backward()
print 'done in %.3f seconds!' % (clock()-start)
if 1:
    PyCont.display(('Ee','I'),stability=True);
    show()
#sys.exit()
###############################################################################
# courbe des Hopfs
###############################################################################
print '--> Hopf Curve ...'
PCargs = args(name='HO1', type='H-C2')
PCargs.initpoint = 'EQ1:H1'
PCargs.freepars = ['Ee', 'Ii']
PCargs.MaxNumPoints = 100
PCargs.LocBifPoints = 'all'
PCargs.verbosity = 2
PyCont.newCurve(PCargs)
print 'Computing Hopf curve...'
start = clock()
PyCont['HO1'].forward()
print '----> backward'
PyCont['HO1'].backward()
print 'done in %.3f seconds!' % (clock()-start)
PyCont['HO1'].display(('Ee','Ii'),stability=True);
###############################################################################
# courbe des Hopfs
###############################################################################
print '--> Hopf Curve #2 ...'
PCargs = args(name='HO2', type='H-C2')
PCargs.initpoint = 'EQ1:H2'
PCargs.freepars = ['Ee', 'Ii']
PCargs.MaxNumPoints = 200
PCargs.LocBifPoints = 'all'
PCargs.verbosity = 2
PyCont.newCurve(PCargs)
print 'Computing Hopf curve...'
start = clock()
PyCont['HO2'].forward()
print '----> backward'
PyCont['HO2'].backward()
print 'done in %.3f seconds!' % (clock()-start)
PyCont['HO2'].display(('Ee','Ii'),stability=True);
PyCont['HO1'].display(('Ee','Ii'),stability=True);grid();axis([-2,6,-6,6])
###############################################################################
# courbe des LP
###############################################################################
print '--> LP Curve ...'
PCargs = args(name='LP1', type='LP-C')
PCargs.initpoint = 'HO1:BT1'
PCargs.freepars = ['Ee', 'Ii']
PCargs.MaxNumPoints = 100
PCargs.LocBifPoints = 'all'
PCargs.verbosity = 2
PyCont.newCurve(PCargs)
print 'Computing Hopf curve...'
start = clock()
PyCont['LP1'].forward()
PyCont['LP1'].backward()
print 'done in %.3f seconds!' % (clock()-start)
PyCont['HO2'].display(('Ee','Ii'),stability=True);
PyCont['HO1'].display(('Ee','Ii'),stability=True);
PyCont['LP1'].display(('Ee','Ii'),stability=True);
show()
print '--> LP Curve ...'
PCargs = args(name='LP2', type='LP-C')
PCargs.initpoint = 'HO2:BT2'
PCargs.freepars = ['Ee', 'Ii']
PCargs.MaxNumPoints = 100
PCargs.LocBifPoints = 'all'
PCargs.verbosity = 2
PyCont.newCurve(PCargs)
print 'Computing Hopf curve...'
start = clock()
PyCont['LP1'].forward()
PyCont['LP1'].backward()
print 'done in %.3f seconds!' % (clock()-start)
###############################################################################
# Display
###############################################################################
PyCont['HO2'].display(('Ee','Ii'),stability=True);
PyCont['HO1'].display(('Ee','Ii'),stability=True);
PyCont['LP1'].display(('Ee','Ii'),stability=True);
PyCont['LP2'].display(('Ee','Ii'),stability=True);
show()