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
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