Limit point continuation on a circle

Help
2012-06-04
2013-11-04
  • Thomas Van Vaerenbergh

    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?

     
  • Drew LaMar

    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

     
  • Thomas Van Vaerenbergh

    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.

    Thanks a lot for your fast and accurate reply!

     

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks