Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Catching fixedpoints

Help
Greg
2012-11-07
2013-11-04
  • Greg
    Greg
    2012-11-07

    Hello,

    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')
    DSargs.pars = pars
    DSargs.varspecs = {'P': Peq,  'X': Xeq, 'M' : Meq}
    DSargs.ics = icdict         #initial conditions
    DSargs.tdomain = [0,50]     # set the range of integration.
    #===================Domain==========================================
    DSargs.xdomain = { 'X': [-100, 15.],'M' : [-10,15],'P' : [-10,15]}
    #===================================================================
    #==========================================0
    ode = Generator.Vode_ODEsystem(DSargs)
    #==========================================0
    #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,
    Greg

     
  • Rob Clewley
    Rob Clewley
    2012-11-07

    Hi,

    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!

    -Rob