TypeError for fixedpoint_2D

Help
whir
2012-07-12
2013-11-04
  • whir
    whir
    2012-07-12

    Hallo again,

    I want to calculate the fixedpoints and their respective stability.
    Here is my code:
    """ EXAMPLE: Two patch Rosenzweig-MacArthur predator-prey model (see )

        When continuue fold curve from ZH point, it detects the ZH point again.  More elegant way to handle
        this?

        Detects DH point

        Drew LaMar, March 2006
    """

    from PyDSTool import *
    from PyDSTool.Toolbox import phaseplane as pp
    from PyDSTool.Toolbox.phaseplane import *

    #==========SET UP SYSTEM==========================
    #k2=4 !
    pars = {'b' : 0.111,'V1' : 0.889,'gam' : 2.1,'k1' : 0.7,'k2' : 1.0,'kgam' : 0.1,'kflux' : 8.1,'mu' : 0.564, 'tau' : 5}

    icdict={'c' : 0.2 , 'n' : 1.05}

    # Set up model
    ceq = 'kflux*mu*n*(b+V1*c/(k1+c)) - gam*c/(kgam + c)'
    neq = '(1-c*c/(k2*k2 + c*c) -n)/tau'

    DSargs = args(name='atri')
    DSargs.pars = pars
    DSargs.varspecs = {'c': ceq, 'n': neq}
    DSargs.ics = icdict         #initial conditions
    DSargs.tdomain =      # set the range of integration.
    DSargs.xdomain = {'c': , 'n': }
    DSargs.fnspecs = {'Jacobian': (,"""[,]""")}

    #==========================================0
    ode = Generator.Vode_ODEsystem(DSargs)
    #==========================================0
    def plotfp(fps):
        for fp in fps:
            plt.plot(fp, fp,'c*',markersize=10)

    #============PHASEPLANE=================================

    figure(4)
    pp.plot_PP_vf(ode, 'c', 'n', scale_exp=-1.2,N=20)

    fp_coord = pp.find_fixedpoints(ode, n=4, eps=1e-8)
    #fp = pp.fixedpoint_2D(ode, Point(fp_coord), eps=1e-8) # not working :(
    print fp_coord
    nulls_x, nulls_y = pp.find_nullclines(ode, 'c', 'n', n=3, eps=1e-8, max_step=0.04,                             fps=[fp_coord])

    # plot the nullclines
    plt.plot(nulls_x, nulls_x, 'b')
    plt.plot(nulls_y, nulls_y, 'g')

    #Traj

    traj = ode.compute('pol')              # integrate ODE
    pts  = traj.sample(dt=0.01)                      # Data for plotting
    plt.plot(pts, pts,'r+')
    plotfp(fp_coord)

    # fps=
    # fps.append(fixedpoint_2D(ode, Point(fp_coord), coords=,
    #                           description='bottom', eps=1e-8))
    # fps.append(fixedpoint_2D(ode, Point(fp_coord), coords=,
    #                           description='middle', eps=1e-8))
    # fps.append(fixedpoint_2D(ode, Point(fp_coord), coords=,
    #                          description='top', eps=1e-8))

    # for fp in fps:
    #     print "F.p. at (%.5f, %.5f) is a %s and has stability %s" % (fp.point,
    #                             fp.point, fp.classification, fp.stability)

    figure(3)

    plt.plot(pts,pts)

    and here is the output:

    Traceback (most recent call last):
      File "bif.py", line 69, in <module>
        description='bottom', eps=1e-8))
      File "/home/gregor/Progs/PyDSTool/Toolbox/phaseplane.py", line 1893, in __init__
        normord=normord, eps=eps)
      File "/home/gregor/Progs/PyDSTool/Toolbox/phaseplane.py", line 1816, in __init__
        self.D = asarray(jac(**jac_test_arg))
    TypeError: jac() got an unexpected keyword argument 'c'

    The strange thing is, when I alter the Jacobian to the (wrong) value c*c, the error disappears,
    but of course I get wrong fixed points. If I leave out the jacobian, there are only 2 fixed points found,
    but there are three.
    I appreciate any help!

    Greets,
    Gregor

     
  • whir
    whir
    2012-07-12

    Oh sorry, I pasted the code where the critical section is commented out, it is the one
    beginning with fps= .

     
  • whir
    whir
    2012-07-17

    Hey,

    I moved on to another system, cause I do some testing an playing with
    pydstools…but thank you very much!

    greets,
    Greg

     
  • whir
    whir
    2012-07-17

    So now everything goes fine with the fixedpoints!
    However, there are convergence problems for the
    bifurcation analysis. I wonder if that is the case because of the
    specific nullcline intersection? One can run the sample code,
    it will draw the phaseplane, I'm interested in bifurcations of
    the middle fixed point! I simply want to shift the 2nd (green) nullcline
    with the paramter 'cc' to the top direction to make the saddle disappear.
    But with the required min stepsize I hardly advance.

    ''' Atri Model '''

    from PyDSTool import *
    from PyDSTool.Toolbox import phaseplane as pp
    from PyDSTool.Toolbox.phaseplane import *

    #==========SET UP SYSTEM==========================
    #k2=4 !
    pars = {'b' : 0.111,'V1' : 0.889,'gam' : 2.1,'k1' : 0.7,'k2' : 1.0,'kgam' : 0.1,'kflux' : 8.1,'mu' : 0.564, 'tau' : 5,'cc' : 1}

    icdict={'c' : 0.2 , 'n' : 1.05}

    # Set up model
    ceq = 'kflux*mu*n*(b+V1*c/(k1+c)) - gam*c/(kgam + c)'
    neq = '(cc-c*c/(k2*k2 + c*c) -n)/tau'

    DSargs = args(name='atri')
    DSargs.pars = pars
    DSargs.varspecs = {'c': ceq, 'n': neq}
    DSargs.ics = icdict         #initial conditions
    DSargs.tdomain =      # set the range of integration.
    DSargs.xdomain = {'c': , 'n': }
    DSargs.fnspecs = {'Jacobian': (,"""[,]""")}

    #==========================================0
    ode = Generator.Vode_ODEsystem(DSargs)
    #==========================================0
    def plotfp(fps):
        for fp in fps:
            plt.plot(fp, fp,'c*',markersize=10)

    #============PHASEPLANE=================================

    figure(4)
    pp.plot_PP_vf(ode, 'c', 'n', scale_exp=-1.2,N=20)

    fp_coord = pp.find_fixedpoints(ode, n=4, eps=1e-12)
    #fp = pp.fixedpoint_2D(ode, Point(fp_coord), eps=1e-8) # not working :(
    print fp_coord
    nulls_x, nulls_y = pp.find_nullclines(ode, 'c', 'n', n=3, eps=1e-8, max_step=0.04,                             fps=[fp_coord])

    # plot the nullclines
    plt.plot(nulls_x, nulls_x, 'b')
    plt.plot(nulls_y, nulls_y, 'g')

    #Traj

    traj = ode.compute('pol')              # integrate ODE
    pts  = traj.sample(dt=0.01)                      # Data for plotting
    plt.plot(pts, pts,'r-')
    plotfp(fp_coord)

    fps=
    fps.append(fixedpoint_2D(ode, Point(fp_coord), coords=,
                              description='bottom', eps=1e-8))
    fps.append(fixedpoint_2D(ode, Point(fp_coord), coords=,
                              description='middle', eps=1e-8))
    fps.append(fixedpoint_2D(ode, Point(fp_coord), coords=,
                             description='top', eps=1e-8))

    for fp in fps:
        print "F.p. at (%.5f, %.5f) is a %s and has stability %s" % (fp.point,
                                fp.point, fp.classification, fp.stability)

    figure(3)

    plt.plot(pts,pts)
    #======BIFURCATIONS===========================================

    ode.set( ics = { 'c': fp_coord,'n': fp_coord } )
    PC = ContClass(ode)

    PCargs = args(name='EQ1', type='EP-C')
    bpar='cc'
    PCargs.freepars =
    PCargs.StepSize = 1e-1
    PCargs.MaxNumPoints = 350
    PCargs.MaxStepSize = 1e-2
    PCargs.MinStepSize = 1e-6
    PCargs.LocBifPoints = 'all'
    PCargs.verbosity = 2
    PCargs.SaveEigen    = True                       # to tell unstable from stable branches
    PC.newCurve(PCargs)

    print 'Computing curve…'
    start = clock()
    PC.forward()
    #PC.forward()
    #PC.backward()

    print 'done in %.3f seconds!' % (clock()-start)

    PC.display((bpar,'n'),stability=True, figure=1)

    show()