## Limit point continuation on a circle document.SUBSCRIPTION_OPTIONS = { "thing": "thread", "subscribed": false, "url": "subscribe", "icon": { "css": "fa fa-envelope-o" } };

Help
2012-06-04
2013-11-04
• Recently, I tried to continue the roots of a 6th order polynomial with PyDSTool. Although I was able to continue the same system using matcont, I did not succeed to do this using pydstool.

To investigate what went wrong I downscaled the problem to limit point continuation on a circle, while this works without any problem with matcont, I am not able to do this using pydstool.

I wonder what's wrong in my current implementation and what I can do to improve the convergence. The code I use is:

```from PyDSTool import *
from pylab import plt
from numpy import log,exp,arange
y=0.8
R=1.0
pars = { 'y': y, 'R': R}
xval=0.0
icdict = {'x': xval}
# Set up model
xstr = '-(x*x+y*y-R*R)'
DSargs = args(name='BenchmarkPolynomialContinuation')
DSargs.pars = pars
DSargs.varspecs = {'x': xstr}
DSargs.ics = icdict
"""
simulate time trace...
"""
testDS.set(tdata=[0,30])
start = clock()
traj=testDS.compute('test')
print '  ... finished in %.3f seconds.\n' % (clock()-start)
plotData = traj.sample()  # works but not VODE accurate for non-terminal events
yaxislabelstr = 'x'
pylab.ylabel(yaxislabelstr)
pylab.xlabel('t')
vlinex=plot(plotData['t'], plotData['x'])
"""
do continuation:
"""
PyCont = ContClass(testDS)
PCargs = args(name='EQ1', type='EP-C')
icdict={'x': plotData['x'][-1]}
DSargs.ics = icdict
PCargs.freepars = ['y']
PCargs.LocBifPoints = 'all'
PCargs.verbosity = 2
PCargs.SaveEigen = True
PyCont.newCurve(PCargs)
"""
simulate EP continuation
"""
print 'Computing curve...'
start = clock()
PyCont['EQ1'].forward()
print 'done in %.3f seconds!' % (clock()-start)
try:
start = clock()
PyCont['EQ1'].backward()
print 'EQ1 backwards: done in %.3f seconds!' % (clock()-start)
except:
print "\nError occurred in EQ1"
print sys.exc_info()[0], sys.exc_info()[1]

PyCont['EQ1'].display(('y','x'), figure='fig1',stability=True)
PCargs=args()
PCargs.verbosity = 12
PCargs.SaveEigen = True
PCargs.name = 'FoldCurve1'
PCargs.type = 'LP-C'
PCargs.freepars = ['y','R']
PCargs.initpoint = 'EQ1:LP1'
PyCont.newCurve(PCargs)
print 'Computing curve...'
try:
start = clock()
PyCont['FoldCurve1'].forward()
PyCont['FoldCurve1'].backward()
print 'FoldCurve1: done in %.3f seconds!' % (clock()-start)
PyCont['FoldCurve1'].display(('y','R'),figure='FoldCurve1',stability=True)
except:
print "\nError occurred in FoldCurve1"
print sys.exc_info()[0], sys.exc_info()[1]
show()
```

This implementation gives this output on my system (64-bit ubuntu):

… finished in 0.050 seconds.

Warning: Variable dimension must be larger than 1 to detect Hopf points.
Computing curve…
/usr/lib/python2.6/dist-packages/scipy/linalg/decomp.py:755: RuntimeWarning: Diagonal number 1 is exactly zero. Singular matrix.
RuntimeWarning)
Warning: Problem updating border vectors.  Using svd…
/usr/lib/python2.6/dist-packages/scipy/linalg/decomp.py:755: RuntimeWarning: Diagonal number 2 is exactly zero. Singular matrix.
RuntimeWarning)
Maximum test function iterations reached.

False positive!

Warning: Problem updating border vectors.  Using svd…

Checking…
|q| = 1.000000
<p,q> = 1.000000
|Aq| = 0.000000
|transpose(A)p| = 0.000000

# LP Point found

0 :
y  =  -1.0
x  =  -1.08018980052e-11

Eigenvalues =

(0.000000,0.000000)

a = 1.0000000000380076

done in 0.790 seconds!
Warning: Problem updating border vectors.  Using svd…
Maximum test function iterations reached.

False positive!

Warning: Problem updating border vectors.  Using svd…

Checking…
|q| = 1.000000
<p,q> = 1.000000
|Aq| = 0.000000
|transpose(A)p| = 0.000000

# LP Point found

0 :
y  =  -1.0
x  =  1.08018980052e-11

Eigenvalues =

(0.000000,0.000000)

a = 1.0000000000380076

Error occurred in EQ1
<type 'exceptions.KeyError'> 'data'
Computing curve…
_compute:  Log(Condition #) = 11.508194
_compute:  Log(Condition #) = 0.629977
_compute:  Log(Condition #) = 0.629977

Error occurred in FoldCurve1
<type 'exceptions.ValueError'> array must not contain infs or NaNs

Apparently, something goes wrong due to singularity problems in this system, but how do I circumvent those? It would be nice if I could solve the problems in this system, as I think they are related with the continuation of the 6th order polynomial I investigate for my research (there I do not have a ValueError, but I have complaints about not finding the initial point on the curve with a 'Could not find starting point on curve.  Stopping continuation.' message).

Is there a fundamental difference between the algorithms used in matcont and pydstool?

• Drew LaMar - 2012-06-04

Hi, Thomas.  I'm actually in the process of correcting some of the issues in PyCont that you have just witnessed.  In particular, the method of continuing limit points uses bordered matrix methods and by default tries to update things in the numerical algorithm so that root solving is as well-conditioned as possible.  This is currently not working properly.  Your issue might be solved by using the following __init__ code for the AddTestFunction class in TestFunc.py in the corresponding place:

```class AddTestFunction(Function):
"""Only works with testfuncs that don't have PreTestFunc and rely only on sysfunc.
BorderMethod update is set to False for now since updating is not working."""
def __init__(self, C, TF_type):
self.sysfunc = C.sysfunc
self.testfunc = TF_type(self.sysfunc, C, update=False)
self.coords = self.sysfunc.coords
self.params = self.sysfunc.params
Function.__init__(self, (self.sysfunc.n, self.sysfunc.m + self.testfunc.m))
```

These are the figures I get:

Please let me know if the fix works!

- Drew

• Nice work! This solved my problem (both the one I had with this simple circle example as the one with my 6th order polynomial). And at first sight I have similar results as in matcont.