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 testDS = Generator.Radau_ODEsystem(DSargs) """ 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?