Catching fixedpoints

  • Greg

    Greg - 2012-11-07


    I'm a bit confused with the behavior of the following program:

    from PyDSTool import *
    from PyDSTool.Toolbox import phaseplane as pp
    #==========SET UP SYSTEM==========================
    pars = {'a_P' : 0.9, 'gam' : .4 , 'b_P' : 0.2 , 'a_XP' : 0,
            'a_PM' : 1, 'b_M' : 0.2, 'b_XM' : 0.3,'k_MP' : 0.01,
             'b_X' : 0.3, 'a_PX' : 1.}
    icdict={'P': 0.4, 'X': 0, 'M' : 0.3}
    # Set up model
    Peq = 'a_P - gam*M*P/(k_MP + P) - b_P*P + a_XP*X' 
    Meq = 'a_PM*P - b_M*M - b_XM*X*M'
    Xeq = 'a_PX*P - b_X*X' #X protein
    DSargs = args(name='mod1') = pars
    DSargs.varspecs = {'P': Peq,  'X': Xeq, 'M' : Meq}
    DSargs.ics = icdict         #initial conditions
    DSargs.tdomain = [0,50]     # set the range of integration.
    DSargs.xdomain = { 'X': [-100, 15.],'M' : [-10,15],'P' : [-10,15]}
    ode = Generator.Vode_ODEsystem(DSargs)
    #ode.set(algparams   =   {'max_pts': 500000})
    fp_coord = pp.find_fixedpoints(ode, n=10, eps=1e-8)
    print 'Nr. of FP:', len(fp_coord)
    print fp_coord

    If I run it with that configuration I get the output:

    Nr. of FP: 3
    ({'P': -0.34769564640811784, 'M': 2.3541360552183699, 'X': -1.1589854880270594}, {'P': 2.6474727845698838, 'M': 0.92976227864800787, 'X': 8.8249092818996129}, {'P': -0.0097771381617756583, 'M': -0.051398333866360696, 'X': -0.032590460539252193})

    But when I make the domain actually smaller, by setting  the xdomain for 'X':  I only
    get 2 FP:

    Nr. of FP: 2
    ({'P': -0.34769564640809114, 'M': 2.3541360552183415, 'X': -1.1589854880269703}, {'P': 2.6474727845698833, 'M': 0.92976227864800809, 'X': 8.8249092818996111})

    Note that the missed FP is also contained in the smaller domain. I came to this when I wondered what
    the parameter 'n' from the 'find_fixedpoints' method is doing. In fact, when I set n<10 the 3rd FP is also
    missed with the larger domain setting. I assured with mathematica that there are indeed 3FP and they
    are also numerically correct with the 1st setting. So why gives a larger domain better results for FB finding,
    and what controls the parameter 'n' ?

    Thanks in advance and greetings,

  • Rob Clewley

    Rob Clewley - 2012-11-07


    I'm afraid this is fairly typical of using fsolve to find roots of equations numerically. It's not necessarily the larger domain that makes the FP found, but more likely the effect it had on the mesh of initial starting points in the space for fsolve. It probably helped put a point close enough to the 3rd FP for fsolve to converge there. The trouble is that fsolve (and other root finders of this kind) have a tendency to be attracted to some roots more than others, for technical numerical reasons I won't go into here. In lieu of a whole better way to implement finding fixed points, I provided the n parameter to increase the mesh resolution of starting guesses. It's hard in general to predict how big n should be, and making it large by default can be very inefficient for analysis of many systems. I am sure there are smarter ways to find fixed points and if anyone would like to contribute code, or if I find the time later, I will try to improve the algorithm!



Log in to post a comment.