Ok, I'm stumped. I have a relatively system that has a Hopf point. I can compute the Hopf point but cannot find the limit cycles past this point.
Here's the code. I get the error "cannot compute the tangent vector" and have fiddled with the tolerances. Can someone point me to some help here on what is happening? (it gets stuck at PC['LC1'].forward())
import PyDSTool as dst
import numpy as np
from matplotlib import pyplot as plt
DSargs.tdomain = [0,5] # set the range of integration.
ode = dst.Generator.Vode_ODEsystem(DSargs) # an instance of the 'Generator' class.
traj = ode.compute ('braking') # integrate ODE
pts = traj.sample(dt=0.01) # Data for plotting
PyPlot commands
plt.plot(pts['t'], pts['x'])
plt.xlabel('time') # Axes labels
plt.ylabel('Lambda') # ...
plt.ylim([0,1]) # Range of the y axis
plt.title(ode.name) # Figure title from model name
Prepare the system to start close to a steady state
ode.set(pars = {'xbar': .075} ) # Lower bound of the control parameter 'i'
ode.set(ics = {'x': .075,'y': 2.67 } ) # Close to one of the steady states present for i=-220
PC = dst.ContClass(ode) # Set up continuation class
PCargs = dst.args(name='EQ1', type='EP-C') # 'EP-C' stands for Equilibrium Point Curve. The branch will be labeled 'EQ1'.
PCargs.freepars = ['xbar'] # control parameter(s) (it should be among those specified in DSargs.pars)
PCargs.MaxNumPoints = 30 # The following 3 parameters are set after trial-and-error
PCargs.MaxStepSize = 1
PCargs.MinStepSize = 1e-5
PCargs.StepSize = .001
PCargs.LocBifPoints = 'all' # detect limit points / saddle-node bifurcations
PCargs.SaveEigen = True # to tell unstable from stable branches
PC.newCurve(PCargs)
PC['EQ1'].forward()
PC['EQ1'].display(['xbar','x'], stability=True, figure=2) # stable and unstable branches as solid and dashed curves, resp.
ok. I moved past the solution error with success but now have even more of a basic question on plotting. I have a simple system with one Hopf point and I can continue now successfully past this point. I can compute the cycles.
What I want to do is to find several of the limit cycles for different values of the parameter, save these and then plot them all on one figure. How simple? I'm confused by the examples and could use just a few wee words of "do this..." or "look at this example."
In particular how can I (a) set where I want to compute the limit cycle (for a fixed and specified parameter) and (b) how can I save the computed cycles and (c) how can I put them all on one plot?
I have a (d) while I'm at it. I would like to compute the L_infty norm (ie max) of the limit cycle - what is the actual representation that I can grab hold of to compute the function?
Thanks - here's the new code that actually works...
Clas
import PyDSTool as dst
import numpy as np
from matplotlib import pyplot as plt
DSargs.tdomain = [0,5] # set the range of integration.
ode = dst.Generator.Vode_ODEsystem(DSargs) # an instance of the 'Generator' class.
traj = ode.compute ('braking') # integrate ODE
pts = traj.sample(dt=0.01) # Data for plotting
PyPlot commands
plt.plot(pts['t'], pts['x'])
plt.xlabel('time') # Axes labels
plt.ylabel('Lambda') # ...
plt.ylim([0,1]) # Range of the y axis
plt.title(ode.name) # Figure title from model name
Prepare the system to start close to a steady state
ode.set(pars = {'xbar': .075} ) # Lower bound of the control parameter 'i'
ode.set(ics = {'x': .075,'y': 2.67 } ) # Close to one of the steady states present for i=-220
PC = dst.ContClass(ode) # Set up continuation class
PCargs = dst.args(name='EQ1', type='EP-C') # 'EP-C' stands for Equilibrium Point Curve. The branch will be labeled 'EQ1'.
PCargs.freepars = ['xbar'] # control parameter(s) (it should be among those specified in DSargs.pars)
PCargs.MaxNumPoints = 30 # The following 3 parameters are set after trial-and-error
PCargs.MaxStepSize = 1
PCargs.MinStepSize = 1e-5
PCargs.StepSize = .001
PCargs.LocBifPoints = 'all' # detect limit points / saddle-node bifurcations
PCargs.SaveEigen = True # to tell unstable from stable branches
PC.newCurve(PCargs)
PC['EQ1'].forward()
PC['EQ1'].display(['xbar','x'], stability=True, figure=2) # stable and unstable branches as solid and dashed curves, resp.
Ok, I'm stumped. I have a relatively system that has a Hopf point. I can compute the Hopf point but cannot find the limit cycles past this point.
Here's the code. I get the error "cannot compute the tangent vector" and have fiddled with the tolerances. Can someone point me to some help here on what is happening? (it gets stuck at PC['LC1'].forward())
import PyDSTool as dst
import numpy as np
from matplotlib import pyplot as plt
we must give a name
DSargs = dst.args(name='Controller')
parameters
DSargs.pars = {'om':80,'k': 1,'xbar':.075,'Tmax':1000,'Tmin':0,'vr2':23.991.4,'vr1':0.650821,'vr3': 0.328354,'Fz':1e5,'m':61201.5,'J':10.025,'r':0.463463,'press':3000,'Tmax':107.940967366519183000,'Scale':10000}
DSargs.varspecs = {'x': '-((1-x)/(Jom))( (r + (J/(rm))(1-x))(Fz(vr1(1-exp(-vr2x))-xvr3)) - yScale)',
'y': '(1/Scale)k(1/(Jom))(x - xbar)(yScale - Tmax)(yScale - Tmin)' }
initial conditions
DSargs.ics = {'x': 0.2, 'y': 2.5 }
DSargs.tdomain = [0,5] # set the range of integration.
ode = dst.Generator.Vode_ODEsystem(DSargs) # an instance of the 'Generator' class.
traj = ode.compute ('braking') # integrate ODE
pts = traj.sample(dt=0.01) # Data for plotting
PyPlot commands
plt.plot(pts['t'], pts['x'])
plt.xlabel('time') # Axes labels
plt.ylabel('Lambda') # ...
plt.ylim([0,1]) # Range of the y axis
plt.title(ode.name) # Figure title from model name
Prepare the system to start close to a steady state
ode.set(pars = {'xbar': .075} ) # Lower bound of the control parameter 'i'
ode.set(ics = {'x': .075,'y': 2.67 } ) # Close to one of the steady states present for i=-220
PC = dst.ContClass(ode) # Set up continuation class
PCargs = dst.args(name='EQ1', type='EP-C') # 'EP-C' stands for Equilibrium Point Curve. The branch will be labeled 'EQ1'.
PCargs.freepars = ['xbar'] # control parameter(s) (it should be among those specified in DSargs.pars)
PCargs.MaxNumPoints = 30 # The following 3 parameters are set after trial-and-error
PCargs.MaxStepSize = 1
PCargs.MinStepSize = 1e-5
PCargs.StepSize = .001
PCargs.LocBifPoints = 'all' # detect limit points / saddle-node bifurcations
PCargs.SaveEigen = True # to tell unstable from stable branches
PC.newCurve(PCargs)
PC['EQ1'].forward()
PC['EQ1'].display(['xbar','x'], stability=True, figure=2) # stable and unstable branches as solid and dashed curves, resp.
plt.show()
PCargs = dst.args(name='LC1', type='LP-C')
PCargs.initpoint = 'EQ1:H1'
PCargs.freepars = ['xbar']
PCargs.StepSize = 1.e-5
PCargs.MinStepSize = 1.e-10
PCargs.NumCollocation = 4
PCargs.NumIntervals = 5
PCargs.LocBifPoints = 'all'
PCargs.FuncTol = 1e-10
PCargs.VarTol = 1e-7
PCargs.TestTol = 1e-6
PCargs.MaxNumPoints = 10
PCargs.SolutionMeasures = 'all'
PCargs.verbosity=4
PCargs.NumSPOut=1000
PC.newCurve(PCargs)
PC['LC1'].forward()
PC['LC1'].display(['xbar','x'])
ok. I moved past the solution error with success but now have even more of a basic question on plotting. I have a simple system with one Hopf point and I can continue now successfully past this point. I can compute the cycles.
What I want to do is to find several of the limit cycles for different values of the parameter, save these and then plot them all on one figure. How simple? I'm confused by the examples and could use just a few wee words of "do this..." or "look at this example."
In particular how can I (a) set where I want to compute the limit cycle (for a fixed and specified parameter) and (b) how can I save the computed cycles and (c) how can I put them all on one plot?
I have a (d) while I'm at it. I would like to compute the L_infty norm (ie max) of the limit cycle - what is the actual representation that I can grab hold of to compute the function?
Thanks - here's the new code that actually works...
Clas
import PyDSTool as dst
import numpy as np
from matplotlib import pyplot as plt
we must give a name
DSargs = dst.args(name='Controller')
parameters
DSargs.pars = {'om':80,'k': 1,'xbar':.075,'Tmax':1000,'Tmin':0,'vr2':23.991.4,'vr1':0.650821,'vr3': 0.328354,'Fz':1e5,'m':61201.5,'J':10.025,'r':0.463463,'press':3000,'Tmax':107.940967366519183000,'Scale':10000}
DSargs.varspecs = {'x': '-((1-x)/(Jom))( (r + (J/(rm))(1-x))(Fz(vr1(1-exp(-vr2x))-xvr3)) - yScale)',
'y': '(1/Scale)k(1/(Jom))(x - xbar)(yScale - Tmax)(yScale - Tmin)' }
initial conditions
DSargs.ics = {'x': 0.2, 'y': 2.5 }
DSargs.tdomain = [0,5] # set the range of integration.
ode = dst.Generator.Vode_ODEsystem(DSargs) # an instance of the 'Generator' class.
traj = ode.compute ('braking') # integrate ODE
pts = traj.sample(dt=0.01) # Data for plotting
PyPlot commands
plt.plot(pts['t'], pts['x'])
plt.xlabel('time') # Axes labels
plt.ylabel('Lambda') # ...
plt.ylim([0,1]) # Range of the y axis
plt.title(ode.name) # Figure title from model name
Prepare the system to start close to a steady state
ode.set(pars = {'xbar': .075} ) # Lower bound of the control parameter 'i'
ode.set(ics = {'x': .075,'y': 2.67 } ) # Close to one of the steady states present for i=-220
PC = dst.ContClass(ode) # Set up continuation class
PCargs = dst.args(name='EQ1', type='EP-C') # 'EP-C' stands for Equilibrium Point Curve. The branch will be labeled 'EQ1'.
PCargs.freepars = ['xbar'] # control parameter(s) (it should be among those specified in DSargs.pars)
PCargs.MaxNumPoints = 30 # The following 3 parameters are set after trial-and-error
PCargs.MaxStepSize = 1
PCargs.MinStepSize = 1e-5
PCargs.StepSize = .001
PCargs.LocBifPoints = 'all' # detect limit points / saddle-node bifurcations
PCargs.SaveEigen = True # to tell unstable from stable branches
PC.newCurve(PCargs)
PC['EQ1'].forward()
PC['EQ1'].display(['xbar','x'], stability=True, figure=2) # stable and unstable branches as solid and dashed curves, resp.
plt.show()
PCargs = dst.args(name='LC1', type='LC-C')
PCargs.initpoint = 'EQ1:H1'
PCargs.freepars = ['xbar']
PCargs.Corrector='MoorePenrose'
PCargs.StepSize = 1.e-5
PCargs.MaxStepSize=.01
PCargs.MinStepSize = 1.e-5
PCargs.NumCollocation = 5
PCargs.NumIntervals = 10
PCargs.LocBifPoints = 'all'
PCargs.FuncTol = 1e-7
PCargs.VarTol = 1e-7
PCargs.TestTol = 1e-7
PCargs.MaxNumPoints = 50
PCargs.SolutionMeasures = 'all'
PCargs.verbosity=4
PCargs.NumSPOut=1000
PC.newCurve(PCargs)
PC['LC1'].forward()
PC['LC1'].update({'MaxNumPoints': 50, 'NumSPOut': 1000})
PC['LC1'].forward()
PC['LC1'].display(['xbar','x'],figure=3)
PC.plot.toggleAll('off',bytype=['P','RG','MX'])
PC.plot.toggleCycles('on')
PC['LC1'].plot_cycles(normalized=True,coords=('x','y'), method='highlight',figure=4)
plt.show()
so a bit of google'ing gives me this page...http://www2.gsu.edu/~matrhc/CodeTopics.html
with this information...
**PyCont 0.3.1:
Two plotting methods were added to plot_cycles: 'stack' and 'highlight'.
Added continuation argument 'StopAtPoints', allowing computation to stop
at specified special points.
Added domain checking through introduction of a new special point labeled
'B' for 'Boundary'. Note that if 'B' is active (through specification in
LocBifPoints), domain information is saved along the curve with the labels 'inside' and 'outside'.
Added continuation argument 'description' allowing the user to give details
on the specific curve computation. When the info() method is called from the curve class, the description will be displayed.
Added argument SPOut to the LimitCycle class, allowing the user to stop at
specified values of variables or parameters.**
How about humming a few more words (or examples...) - this seems to have the functionality that I want...
Clas