Example 3.1 - Applied Nonlinear Dynamics [A. Nayfeh, 1995]

[Bruno Cayres]

Hello, guys.

I'm trying to reproduce the Example 3.1 from the book Applied Nonlinear Dynamics: Finding periodic solutions and bifurcation analysis. However, I stocked. Here is the code:

import PyDSTool as dst
import matplotlib.pyplot as plt
from PyDSTool.Toolbox import phaseplane as pp

pars = {'mu':1,
        'alpha':-1,
        'omega':1,
        'beta':1}

icdict = {'x': 0.5, 'y': 0.1}

xstr = 'mu*x - omega*y + (alpha*x - beta*y)*(x*x + y*y)'
ystr = 'omega*x + mu*y + (beta*x + alpha*y)*(x*x + y*y)'

DSargs = dst.args(name='2D-System')
DSargs.pars = pars
DSargs.varspecs = {'x': xstr, 'y': ystr}
DSargs.ics = icdict
DSargs.tdomain = [0, 15] # set the range of integration.
DSargs.xdomain = {'x': [-3, 3], 'y': [-3, 3]}
DSargs.pdomain = {'mu': [-1, 1]}
DSargs.algparams = {'max_pts': 3000, 'init_step': 0.02, 'stiff':True}
ode = dst.Generator.Vode_ODEsystem(DSargs)

traj = ode.compute('polarization')
pts  = traj.sample(dt=0.01)

plt.plot(pts['t'], pts['x'])
plt.xlabel('time') # Axes labels
plt.ylabel('x')
plt.show()

plt.plot(pts['x'], pts['y'])
plt.xlabel('x') # Axes labels
plt.ylabel('y')
plt.show()

# plot vector field, using a scale exponent to ensure arrows are well spaced
# and sized
pp.plot_PP_vf(ode, 'x', 'y', scale_exp=-1, N=50)

plt.show()


PyCont = dst.ContClass(ode)

PCargs = dst.args(name='EQ1', type='EP-C')
PCargs.freepars = ['mu']
PCargs.StepSize = 0.1
PCargs.MaxNumPoints = 50
PCargs.MaxStepSize = 2
PCargs.MinStepSize = 1e-5
PCargs.LocBifPoints = ['all']
PCargs.verbosity = 2
PCargs.SaveEigen = True
PyCont.newCurve(PCargs)

print('Computing curve...')
start = dst.clock()
PyCont['EQ1'].forward()
PyCont['EQ1'].backward()
print('done in %.3f seconds!' % (dst.clock()-start))

PyCont['EQ1'].display(['mu','y'], stability=True)
plt.show()

For mu>=0, I should've observe a Hopf Point, but it doesn't happen. What did I miss?
Also, the other options as 'H-C1' or H-C2' create this error:

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

Any help would be pretty much appreciated.

Thanks.

[Maurizio De Pitta']

Hi Bruno,
have you tried playing with StepSize, MaxStepSize and MinStepSize? Often it is related to the initialization of the tangent vector for continuation. Indeed this is why often, you can run the same continuation with the same parameters and get different results.

Also see if you can play with the different tolerances values. See this post.

Hope this helps.
If not, Rob and the other folks will come out shortly :)

M

[Bruno Cayres]

Hi, Maurizio.

Yes, I have tried. Also, I checked the suggested post, but still the same result =[

Thanks for your answer.

Bruno.

[Drew LaMar]

Hi, Bruno. I think it's the PCargs.LocBifPoints = ['all'] line. I changed it to PCargs.LocBifPoints = 'all' and it worked.

[Bruno Cayres]

Hi, Drew.

I did what you said and it did not work.
I downloaded the package again. First, it presents a SciPy version error and here I found a solution.

Even overcome these, another error arises:

PyDSTool/PyCont/TestFunc.py", line 381, in __init__
    self.data.P = zeros((n*(n-1)/2, n*(n-1)/2), float)
TypeError: 'float' object cannot be interpreted as an integer

So, I did this:

self.data.P = zeros((n*(n-1) // 2, n*(n-1) // 2), float)

And then, the previous error arises:

PyDSTool/PyCont/TestFunc.py", line 405, in bialtprodeye
    self.data.P[p*(p-1)/2 + q][r*(r-1)/2 + s] = A[p][p] + A[q][q]
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

But this error only occurs when I demand Hopf Point (Or Limit Cycle Curve).

Thanks for your help.

Bruno.

[Bruno Cayres]

Hello, guys.

I already suppressed the current errors and I did found the Hopf bifurcation point at mu=0 as we can see in Figure_1 (Cheers =D)!!!
I know that, from mu>= 0, the stability jumps to a stable limit cycle (or periodic attractor).
So, I'm did this,

PCargs = dst.args(name='SN1', type='LC-C')
PCargs.initpoint    = 'EQ1:H1'
PCargs.freepars     = ['alpha']
PCargs.MaxStepSize  = 0.01
PCargs.StepSize = 0.001
PCargs.LocBifPoints = 'all  '
PCargs.MaxNumPoints = 400
PCargs.verbosity = 5
PyCont.newCurve(PCargs)
# PyCont['SN1'].forward()
PyCont['SN1'].backward()
PyCont['SN1'].display(['mu','y'], stability=True)

and the Figure_2 showed up.

Did I do something wrong? What does "MX1" mean?

I'm already very happy!!!
Thank you, guys =D

[Drew LaMar]

Hi, Bruno. Question: why is your free parameter for the limit cycles alpha? I changed it to mu, the same as the equilibrium point curve, and got the attached diagram. Does this look close to what you are expecting? Also, if you add PCargs.StopAtPoints = 'B' to your equilibrium point curve arguments, the curve will stop at the predetermined boundary points (plus and minus 1).

[Bruno Cayres]

Hi, Drew.

That is it. I used alpha parameter by mistake =/,
Now, I did reproduce the results from literature =D

But, What python version do you, guys, use?
I'm already with SciPy v1.0.0 and Python 3.6.4.

This is probably why some errors popped up to me.

Thank you very much =D