Analysis of the May model with SAGE

Dimensionless the May model reads:

\[\frac{dx}{d\tau} = (1-x)\, x - \alpha \frac{x y}{d+x}\]

\[\frac{dy}{d\tau} = \beta \left(1- \frac{y}{x}\right) y\]

where we have defined the following dimensionless parameters:

\[\tau = r t, \quad \alpha = \frac{b_0}{h r}, \quad d = \frac{D}{K} , \quad \beta = \frac{s}{r}\]

This scaling procedure leads to a set of two differential equations with only three parameters. The dimensionless time is scaled to the reproducinng rate \(r\) of prey. The parameter \(\beta\) is the relation of the reproducing rate of predators \(s\) to the reproducing rate of prey \(r\). If \(\beta \lt 1\) then the reproducting tempo of predators is smaller than prey and therefore the prey population can survive. On the other hand, f \(\beta \gt 1\) then the reproducting tempo of predators is greater than for prey and therefore the prey population can become extinct (can fail to survive). But because the set of two equations is nonlinear, such simplified considerations cannot be fully true. We have to use precise mathematical methods to check this.

STATIONARY STATES

Stationary states are determined by a set of two algebraic equations:

\[(1-x)\, x - \alpha \frac{x y}{d+x} = 0\]

\[\beta \left(1- \frac{y}{x}\right) y = 0\]

One stationary state can easily be found:

\[x_1=1, \quad y_1 =0\]

In theis state there are no predators and the state of prey is the same as in the Verhulst system. We should cjeck whether this stse is stable.

From the second equation we find that \(y=x\) solves this equation. We put it in the first equation and then the other stationary states are determined by the equations:

\[y=x , \quad \quad (1-x) - \alpha \frac{ y}{d+x} = 0\]

From the second equation we obtain the condition:

\[x^2 + (\alpha + d -1) x -d =0\]

It is a trinomial square and we take into account only non-negative solutions of this equation, \(x \ge 0\). Discriminant is always positive, namely,

\[\Delta = (\alpha + d -1)^2 + 4d \gt 0\]

So, the second stationary state is

\[x_2= y_2 = \frac{1}{2} \left[- (\alpha + d -1) + \sqrt{\Delta}\right]\]

Let us note that this state does not depend on the parameter \(\beta\), i.e. on the birth rate of prey and predators.

STABILITY OF STATIONARY STATES

1. In the first step, we calculate the Jacobi matrix:

\[\begin{split}\begin{aligned} \quad \quad \quad J = \begin{bmatrix}\frac{ \partial f}{\partial x}& \frac{\partial f}{\partial y}\\ \frac{\partial g}{\partial x}& \frac{\partial g}{\partial y} \end{bmatrix} = \begin{bmatrix}1-2x-\alpha y \frac{d}{(x+d)^2}& -\alpha \frac{x}{x+d} \\ \frac{\beta y^2}{x^2}& \beta - \frac{2\beta y}{x} \end{bmatrix} \end{aligned}\end{split}\]

2. In the second step, we determine the eigenvalues of the Jacobi matrix \(|J -\lambda I|=0\):

1. For the stationary state \((1, 0)\) one gets:

\[\begin{split}\begin{aligned} \quad \quad \quad J(1, 0) = \begin{bmatrix}-1& - \frac{\alpha}{1+d} \\ 0& \beta \end{bmatrix} \end{aligned}\end{split}\]

and the eigenvalues are: \(\lambda_{1} = -1, \quad \lambda_{2} = \beta > 0\). So, this state is not stable because one eigenvalue is positive, \(\lambda_{2} \gt 0\). A small perturbation will displace the system from this state.

2. For the second stationary state the problem of stability is more complicated because the Jacobi matrix takes the form

\[\begin{split}\begin{aligned} \quad \quad \quad J(x_2, y_2) = \begin{bmatrix}x_2\left[ \frac{\alpha x_2}{(x_2+d)^2} -1\right]& - \frac{\alpha x_2}{x_2+d} \\ \beta & -\beta \end{bmatrix} \end{aligned}\end{split}\]

We used the equation for the stationary state to transform the term \(J_{11}\) to the form as in this matrix. One can try to determine the eigenvalues of this matrix. But it is not good way. What we want to know is the sign of a real part (positive or negative) of both eigenvalues \((\lambda_1, \lambda_2)\). Because this matrix is \(2 \times 2\), we get a trinomial square for \(\lambda\). The conditions for its negative part reads:

\[\lambda_1 + \lambda_2 \lt 0 \quad \mbox{and} \quad \lambda_1 \; \lambda_2 \gt 0, \quad \quad \mbox{it means that} \quad \mbox{Tr} \, J \lt 0, \quad \quad \mbox{det} \,J \gt 0\]

WORK:

Show that for arbitrary (positive) values of parameters \(\alpha, \beta, d\), the second condition \(\mbox{det} \,J \gt 0\) is always fulfiled.

The first condition for the stability of the state \((x_2, y_2)\) takes the form:

\[b \gt x_2\left[ \frac{\alpha x_2}{(x_1+d)^2} -1\right] = \phi (\alpha, d)\]

Because \(x_2\) depends on 2 parameters \(\alpha\) and \(d\), the right-hand side defines the surface in the 3-dimensional space.

---

var('a b d x y')
ode_lotka=[x*(1-x)-(a*x*y)/(x+d),b*y*(1-y/x)]
show(ode_lotka)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-{\left(x - 1\right)} x - \frac{a x y}{d + x}, -b y {\left(\frac{y}{x} - 1\right)}\right]$$

\[\left[-{\left(x - 1\right)} x - \frac{a x y}{d + x}, -{\left(\frac{y}{x} - 1\right)} b y\right]\]

J = jacobian(ode_lotka,[x,y])

detJ = det(J)

detJ.variables()

(a, b, d, x, y)

detJ.subs({a:1})

(b*(y/x - 1) + b*y/x)*(2*x + y/(d + x) - x*y/(d + x)^2 - 1) + b*y^2/((d + x)*x)

for sol_ in solve(ode_lotka,[x,y]):
    show(sol_)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[x = 0, y = 0\right]$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[x = 1, y = 0\right]$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[x = -d, y = 0\right]$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[x = -\frac{1}{2} \, a - \frac{1}{2} \, d - \frac{1}{2} \, \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} + \frac{1}{2}, y = -\frac{1}{2} \, a - \frac{1}{2} \, d - \frac{1}{2} \, \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} + \frac{1}{2}\right]$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[x = -\frac{1}{2} \, a - \frac{1}{2} \, d + \frac{1}{2} \, \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} + \frac{1}{2}, y = -\frac{1}{2} \, a - \frac{1}{2} \, d + \frac{1}{2} \, \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} + \frac{1}{2}\right]$$

sols = solve(ode_lotka,[x,y])

We check if solutions are in first quadrant of coordinate system. We need it for positive \(a\) and \(d\).

sol_ = sols[4]
region_plot(sol_[0].rhs()>=0 and sol_[1].rhs()>=0,(a,-2,2),(d,-2,2),incol='red',plot_points=400,figsize=4)

/usr/local/lib/SageMath82/local/lib/python2.7/site-packages/sage/plot/contour_plot.py:1515: RuntimeWarning: invalid value encountered in less
  neg_indices = (xy_data_arrays < 0).all(axis=0)

We check if \(\det J\) is positive al all parametres. If \(b>0\) then \(\frac{\det J}{b}\) depends only on \((a,d)\):

(detJ.subs(sol_)/b).full_simplify().variables()

(a, d)

thus we can check numerically using region plot if it takes positive values for all \((a,d)\) in the first quadrant of parameter space.

region_plot((detJ.subs(sol_)/b).full_simplify()>0,(a,-2,2),(d,-2,2),incol='red',figsize=4,plot_points=400)

trJ = J.trace()

One can solve equation \(tr J==0\) for \(b\) and the inequality becomes: \(b< \mbox{trJ_b}\).

trJ_b = trJ.subs(sol_).solve(b)[0].rhs()

trJ_b.variables()

(a, d)

implicit_plot3d(trJ.subs(sol_),(a,1e-3,1.41),(d,1e-3,.61),(b,0,1.5)).show(viewer='tachyon')

plt = plot3d(trJ_b,(a,1e-3,1),(d,1e-3,1))
plt += plot3d(0,(a,0,1),(d,0,1),color='red')
plt.show(viewer='tachyon')

implicit_plot(trJ_b,(a,0,2),(d,0,1))

plt_cond_trJ = region_plot( trJ_b>0,(a,0,2),(d,+1e-4,2),plot_points=400,
                          axes_labels=[r'$a$',r'$d$'])
plt_cond_trJ

var('X')
expr = trJ_b.subs((a^2 + 2*a*d + d^2 - 2*a + 2*d + 1) == X)
expr.full_simplify()

1/2*(4*a^3 + 2*(a - 1)*d^2 - 10*a^2 + 2*(3*a^2 + a - 2)*d - (3*a^2 - d^2 + X - 4*a + 1)*sqrt(X) + 8*a - 2)/(a^2 + d^2 - sqrt(X)*(a - d - 1) - 2*a + 2*d + 1)

expr.subs([sqrt(X)])

1/2*(4*a^3 + (2*a + sqrt(X) - 2)*d^2 - a^2*(3*sqrt(X) + 10) + 2*(3*a^2 + a - 2)*d + 4*a*(sqrt(X) + 2) - X^(3/2) - sqrt(X) - 2)/(a^2 + d^2 - a*(sqrt(X) + 2) + d*(sqrt(X) + 2) + sqrt(X) + 1)

var('a,d,b,x,y,t')
ode_lotka=vector([x*(1-x)-(a*x*y)/(x+d),b*y*(1-y/x)]);

b_in = 0.12
d_in = 0.15
a_in = 1.2
p = {a:a_in, d:d_in, b:b_in}

b_in<trJ_b.subs(p).n()

True

plt_cond_trJ + point( (a_in,d_in),color='red',size=30,zorder=10,plot_points=400)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-49-00ec68891fc8> in <module>()
----> 1 plt_cond_trJ + point( (a_in,d_in),color='red',size=Integer(30),zorder=Integer(10),plot_points=Integer(400))

NameError: name 'plt_cond_trJ' is not defined

ode_lotka.subs(p).list()

[-(x - 1)*x - 1.20000000000000*x*y/(x + 0.150000000000000),
 -0.120000000000000*y*(y/x - 1)]

ode_lotka

(-(x - 1)*x - a*x*y/(d + x), -b*y*(y/x - 1))

def phase_may(p, xylims  = [(x,1e-4,1),(y,1e-4,1)]):
    pkts = solve(ode_lotka.subs(p).list(),[x,y])

    colors = ['red','green','blue']



    plt_fixpoints = point([vector([x,y]).subs(pkt) for pkt in pkts if pkt[0].rhs()>=0 and pkt[1].rhs()>=0],color='red',size=30)


    plt_null = sum([implicit_plot(eq,*xylims,color=c) for eq,c in zip(ode_lotka.subs(p),colors)])
    plt_phase = plot_vector_field(ode_lotka.subs(p).normalized(),*xylims,plot_points=20)
    plt_phase += plt_fixpoints+plt_null

    ith_1st = [i for i,pkt in enumerate(pkts) if pkt[0].rhs()>0 and pkt[1].rhs()>0 ][0]

    Jp = jacobian(ode_lotka.subs(p),[x,y]).subs(pkts[ith_1st])

    has_LC = False

    if Jp.trace()>0:
        print("fix point is unstable")
        has_LC = True

    print('detJ>trJ^2/4?',Jp.det().n()>Jp.trace().n()^2/4)

    plt_trajs=Graphics()

    ics_lst = []
    for dy in srange(0.01,0.71,0.05):
        ics_lst.append( vector((x,y+dy)).subs(pkts[ith_1st]) )

    ics_lst = [ (random(),random()) for i in range(50)]

    for ics in ics_lst:
        T = srange(0,50,0.01)
        traj = desolve_odeint(ode_lotka.subs(p),ics,T,[x,y])

        plt_trajs += line(traj,thickness=0.3)

    if has_LC:
        T = srange(0,350,0.2)
        traj_lc = desolve_odeint(ode_lotka.subs(p),(1,1.23),T,[x,y])[-205:,:]
        plt_traj_lc = line(traj_lc,color='red',zorder=11,thickness=3)


    plt_phase += plt_trajs
    if has_LC:
        plt_phase += plt_traj_lc
    return plt_phase,pkts[ith_1st],Jp.trace(),Jp.det()

phase_may(p)[0]

fix point is unstable
('detJ>trJ^2/4?', True)

@interact
def _(a_in = slider(0,2,0.01,default=1.0),\
        b_in = slider(0,2,0.01,default=0.1),\
        d_in = slider(0,2,0.01,default=0.1) ):

    p = {a:a_in,d:d_in,b:b_in}
    plt,pkt,tr,det_ = phase_may(p)
    plt_dettr = plot(x^2/4, (x,-3,3),ymin=-1,ymax=1,legend_label=r'$\detJ=\frac{\mathrm{Tr}J^2}{4}$')
    plt_dettr += point((tr,det_),color='red',size=30,axes_labels=['trJ','detJ'],legend_label="")
    #plt.show(figsize=6)
    #plt_dettr.show(figsize=5)
    graphics_array([[plt,plt_dettr]]).show(figsize=[10,5])

var('a,d,b,x,y,t')
ode_lotka=vector([x*(1-x)-(a*x*y)/(x+d),b*y*(1-y/x)]);

p = {a:1.0, d:0.1, b:1}

sol = solve(list(ode_lotka.subs(p)),[x,y])[4]

J = jacobian(ode_lotka.subs(p),[x,y]).subs(sol).n()
show(J)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 0.262515743197908 & -0.729843788128358 \\ 1.00000000000000 & -1.00000000000000 \end{array}\right)$$

J.trace(),

(-0.737484256802092,)

J.det()

0.467328044930449

J.trace()^2/4

0.135970757257733

xylims  = [(x,1e-4,1),(y,1e-4,1)]

phase_may(p,xylims=[(x,+1e-4,1),(y,+1e-4,1)])[0]

('detJ>trJ^2/4?', True)

eps = 1e-3
solve(ode_lotka.subs(p).subs(y==eps)*vector([0,1])>=0,x)

[[x < 0], [x >= (1/1000)]]

solve(ode_lotka.subs(p).subs(x==0.00001)*vector([1,0])>0,x)

[[y < (10001/100001)]]

J

[ 0.262515743197908 -0.729843788128358]
[ 0.100000000000000 -0.100000000000000]

@interact
def May_phase(a_in = slider(0,2,0.01,default=1.0),\
        b_in = slider(0,2,0.01,default=0.1),\
        d_in = slider(0,2,0.01,default=0.1) ):
    p = {a:a_in,d:d_in,b:b_in}
    ode_lotka_num =[i.subs(p) for i in ode_lotka]
    pkt_osob=solve(ode_lotka_num,x,y, solution_dict=True)
    x_osobliwy,y_osobliwy=0,0
    plt_pkt=[]
    for n_pkt,pkt in enumerate(pkt_osob):
       x_osobliwy,y_osobliwy=pkt[x].n(),pkt[y].n()
       plt_pkt.append(point([x_osobliwy,y_osobliwy],size=30,color='red') )
       JJ=jacobian(ode_lotka_num,[x,y])
       JJ0=JJ.subs({x:x_osobliwy+1e-8,y:y_osobliwy+1e-8})
       print n_pkt+1,":",x_osobliwy.n(digits=3),y_osobliwy.n(digits=3),vector(JJ0.eigenvalues()).n(digits=3)
       if pkt[x]>0 and pkt[y]>0 :
           print "Czy pkt. jest stabilny?",bool(b_in>f(a_in,d_in))
           print "detJ>(TrJ/2)^2",(detJ.subs(sol_)-(trJ.subs(sol_)/2)^2).subs(p).n()
    plt1 = plot_vector_field(vector(ode_lotka_num)/vector(ode_lotka_num).norm(),(x,-0.1,2),(y,-0.1,2))
    #plt2a = implicit_plot(ode_lotka_num[0],(x,-0.10,2),(y,-0.10,2),color='green')
    plt2a = plot(solve(ode_lotka_num[0],y)[0].rhs(),(x,-0.10,2),ymin=-0.10,ymax=2, color='green')
    show(ode_lotka_num)
    plt2b = implicit_plot(ode_lotka_num[1],(x,-0.10,2),(y,-.010,2),color='blue',xmin=-0.03)
    T = srange(0,123,0.1)
    sol1 = desolve_odeint(vector(ode_lotka_num), [0.82,0.85], T, [x,y])
    plt_solution = list_plot(sol1.tolist(), plotjoined=True, color='brown')
    show( sum(plt_pkt)+plt1+plt2a+plt2b+plt_solution )

var('Delta',latex_name='\sqrt{\Delta}')
subDelta = sqrt(a^2 + 2*(a + 1)*d + d^2 - 2*a + 1)==Delta

sol_[1].rhs().subs(subDelta)

1/2*Delta - 1/2*a - 1/2*d + 1/2

detJ.subs(sol_).subs(subDelta).expand().full_simplify()

-((Delta + 2*a + 2)*b*d^2 + b*d^3 - (Delta^2 - a^2 + 2*a - 1)*b*d - (Delta^3 + Delta*a^2 + 2*Delta^2 - 2*(Delta^2 + Delta)*a + Delta)*b)/(Delta^2 - 2*(Delta + 1)*a + a^2 + 2*(Delta - a + 1)*d + d^2 + 2*Delta + 1)

point1 = sols[0]
point2 = sols[1]
point3 = sols[4]
#point3 = map(lambda x:x.subs(subDelta),point3)

point3

[x == -1/2*a - 1/2*d + 1/2*sqrt(a^2 + 2*(a + 1)*d + d^2 - 2*a + 1) + 1/2,
 y == -1/2*a - 1/2*d + 1/2*sqrt(a^2 + 2*(a + 1)*d + d^2 - 2*a + 1) + 1/2]

trJ.subs(point3).simplify().numerator()

a^3 - a^2*b + a^2*d + 2*a*b*d + a*d^2 - b*d^2 + d^3 - 3*sqrt(a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*a^2 + 2*sqrt(a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*a*b - 2*sqrt(a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*b*d + sqrt(a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*d^2 + 3*(a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*a - 4*a^2 - (a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*b + 2*a*b - (a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*d - 6*a*d - 2*b*d - (a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)^(3/2) + 4*sqrt(a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*a - 2*sqrt(a^2 + 2*a*d + d^2 - 2*a + 2*d + 1)*b + 5*a - b - 3*d - sqrt(a^2 + 2*a*d + d^2 - 2*a + 2*d + 1) - 2

detJ.subs(point3).simplify().numerator().show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}-a^{2} b d - 2 \, a b d^{2} - b d^{3} + \sqrt{a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1} a^{2} b - \sqrt{a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1} b d^{2} - 2 \, {\left(a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1\right)} a b + {\left(a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1\right)} b d + 2 \, a b d - 2 \, b d^{2} + {\left(a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1\right)}^{\frac{3}{2}} b - 2 \, \sqrt{a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1} a b + 2 \, {\left(a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1\right)} b - b d + \sqrt{a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1} b$$

2d

\[\det J = \left(\frac{\mathrm{Tr} J}{2}\right)^2\]

(detJ.subs(sol_)-(trJ.subs(sol_)/2)^2).subs({a:1,d:1,b:.1}).n()

0.00357530887414820

point3

[x == -1/2*a - 1/2*d + 1/2*sqrt(a^2 + 2*(a + 1)*d + d^2 - 2*a + 1) + 1/2,
 y == -1/2*a - 1/2*d + 1/2*sqrt(a^2 + 2*(a + 1)*d + d^2 - 2*a + 1) + 1/2]

var('x0')

x0

trJ.subs({x:x0,y:x0}).factor().numerator().collect(b).show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}-a d x_{0} - 2 \, d^{2} x_{0} - 4 \, d x_{0}^{2} - 2 \, x_{0}^{3} - {\left(d^{2} + 2 \, d x_{0} + x_{0}^{2}\right)} b + d^{2} + 2 \, d x_{0} + x_{0}^{2}$$

y_z_pierwszego=solve(ode_lotka[0],y,solution_dict=True)[0]
drugie=ode_lotka[1].subs(y_z_pierwszego)
show(drugie)
show(solve(drugie,x,solution_dict=True)[0])
x_0=x.subs(solve(drugie,x,solution_dict=True)[1])
y_0=y_z_pierwszego[y].subs({x:x_0}).expand()
show(x_0)
show( y_0 )

$$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left({\left(d - 1\right)} x + x^{2} - d\right)} b {\left(\frac{{\left(d - 1\right)} x + x^{2} - d}{a x} + 1\right)}}{a}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left\{x : -\frac{1}{2} \, a - \frac{1}{2} \, d - \frac{1}{2} \, \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} + \frac{1}{2}\right\}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{2} \, a - \frac{1}{2} \, d + \frac{1}{2} \, \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} + \frac{1}{2}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{2} \, a - \frac{1}{2} \, d + \frac{1}{2} \, \sqrt{a^{2} + 2 \, a d + d^{2} - 2 \, a + 2 \, d + 1} + \frac{1}{2}$$

ode_lotka[0].diff(x).show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}-2 \, x - \frac{a y}{d + x} + \frac{a x y}{{\left(d + x\right)}^{2}} + 1$$

\[-\frac{a y}{d + x} + \frac{a x y}{{\left(d + x\right)}^{2}} - 2 \, x + 1\]

JJ = jacobian(ode_lotka,[x,y])
show(JJ)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} -2 \, x - \frac{a y}{d + x} + \frac{a x y}{{\left(d + x\right)}^{2}} + 1 & -\frac{a x}{d + x} \\ \frac{b y^{2}}{x^{2}} & -b {\left(\frac{y}{x} - 1\right)} - \frac{b y}{x} \end{array}\right)$$

\[\begin{split}\begin{aligned} \left(\begin{array}{rr} -\frac{a y}{d + x} + \frac{a x y}{{\left(d + x\right)}^{2}} - 2 \, x + 1 & -\frac{a x}{d + x} \\ \frac{b y^{2}}{x^{2}} & -{\left(\frac{y}{x} - 1\right)} b - \frac{b y}{x} \end{array}\right) \end{aligned}\end{split}\]

#mamy x0=y0 ;-)
var('x0')
JJ0 = JJ.subs({x:x0,y:x0})

show(JJ0)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} -\frac{a x_{0}}{d + x_{0}} - 2 \, x_{0} + \frac{a x_{0}^{2}}{{\left(d + x_{0}\right)}^{2}} + 1 & -\frac{a x_{0}}{d + x_{0}} \\ b & -b \end{array}\right)$$

\[\begin{split}\begin{aligned} \left(\begin{array}{rr} -\frac{a x_{0}}{d + x_{0}} + \frac{a x_{0}^{2}}{{\left(d + x_{0}\right)}^{2}} - 2 \, x_{0} + 1 & -\frac{a x_{0}}{d + x_{0}} \\ b & -b \end{array}\right) \end{aligned}\end{split}\]

show(JJ0.trace())

$$\newcommand{\Bold}[1]{\mathbf{#1}}-b - \frac{a x_{0}}{d + x_{0}} - 2 \, x_{0} + \frac{a x_{0}^{2}}{{\left(d + x_{0}\right)}^{2}} + 1$$

\[-\frac{a x_{0}}{d + x_{0}} + \frac{a x_{0}^{2}}{{\left(d + x_{0}\right)}^{2}} - b - 2 \, x_{0} + 1\]

\[b = x_1\left[ \frac{\alpha x_1}{(x_1+d)^2} -1\right]\]

trJ.subs({x:x0,y:x0}).operands()

[-b, -a*x0/(d + x0), -2*x0, a*x0^2/(d + x0)^2, 1]

expr_murray = x0*(a*x0/(x0-d)^2-1)
expr_murray.show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}x_{0} {\left(\frac{a x_{0}}{{\left(d - x_{0}\right)}^{2}} - 1\right)}$$

\[{\left(\frac{a x_{0}}{{\left(d - x_{0}\right)}^{2}} - 1\right)} x_{0}\]

show( JJ0.trace().subs({x0:x_0})+b )

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(a + d - \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} - 1\right)}^{2} a}{{\left(a - d - \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} - 1\right)}^{2}} - \frac{{\left(a + d - \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} - 1\right)} a}{a - d - \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} - 1} + a + d - \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1}$$

\[\begin{split}\begin{aligned} -\frac{{\left(a + d - \sqrt{2 \, {\left(a + 1\right)} d + a^{2} + d^{2} - 2 \, a + 1} - 1\right)} a}{a - d - \sqrt{2 \, {\left(a + 1\right)} d + a^{2} + d^{2} - 2 \, a + 1} - 1} \\ + \frac{{\left(a + d - \sqrt{2 \, {\left(a + 1\right)} d + a^{2} + d^{2} - 2 \, a + 1} - 1\right)}^{2} a}{{\left(a - d - \sqrt{2 \, {\left(a + 1\right)} d + a^{2} + d^{2} - 2 \, a + 1} - 1\right)}^{2}} \\ + a + d - \sqrt{2 \, {\left(a + 1\right)} d + a^{2} + d^{2} - 2 \, a + 1} \end{aligned}\end{split}\]

p={a:1.23,d:1.01}
show( JJ0.trace().subs({x0:x_0}).subs(p) )
expr_murray.subs({x0:x_0}).subs(p)

$$\newcommand{\Bold}[1]{\mathbf{#1}}-b - 0.404054289657954$$

1.35696399470668

\[\begin{split}\begin{aligned} -b - 0.404054289657954 \\ 1.35696399470668 \end{aligned}\end{split}\]

var('a,d,b,x,y,t')
ode_lotka=[x*(1-x)-(a*x*y)/(x+d),b*y*(1-y/x)];
#Murray eq. 3.28
f(a,d)=(a-sqrt(  (1-a-d)^2+4*d) )*(1+a+d-sqrt((1-a-d)^2+4*d))/(2*a)
@interact
def myf(a_in = slider(0,2,0.01,default=1.0),b_in = slider(0,2,0.01,default=0.1),d_in = slider(0,2,0.01,default=0.1) ):
    p={a:a_in,d:d_in,b:b_in}
    ode_lotka_num=[i.subs(p) for i in ode_lotka]
    pkt_osob=solve(ode_lotka_num,x,y, solution_dict=True)
    x_osobliwy,y_osobliwy=0,0
    plt_pkt=[]
    for n_pkt,pkt in enumerate(pkt_osob):
       x_osobliwy,y_osobliwy=pkt[x].n(),pkt[y].n()
       plt_pkt.append(point([x_osobliwy,y_osobliwy],size=30,color='red') )
       JJ=jacobian(ode_lotka_num,[x,y])
       JJ0=JJ.subs({x:x_osobliwy+1e-8,y:y_osobliwy+1e-8})
       print n_pkt+1,":",x_osobliwy.n(digits=3),y_osobliwy.n(digits=3),vector(JJ0.eigenvalues()).n(digits=3)
       if pkt[x]>0 and pkt[y]>0 :
           print "Czy pkt. jest stabilny?",bool(b_in>f(a_in,d_in))
           print "detJ>(TrJ/2)^2",(detJ.subs(sol_)-(trJ.subs(sol_)/2)^2).subs(p).n()
    plt1 = plot_vector_field(vector(ode_lotka_num)/vector(ode_lotka_num).norm(),(x,-0.1,2),(y,-0.1,2))
    #plt2a = implicit_plot(ode_lotka_num[0],(x,-0.10,2),(y,-0.10,2),color='green')
    plt2a = plot(solve(ode_lotka_num[0],y)[0].rhs(),(x,-0.10,2),ymin=-0.10,ymax=2, color='green')
    show(ode_lotka_num)
    plt2b = implicit_plot(ode_lotka_num[1],(x,-0.10,2),(y,-.010,2),color='blue',xmin=-0.03)
    T = srange(0,123,0.1)
    sol1 = desolve_odeint(vector(ode_lotka_num), [0.82,0.85], T, [x,y])
    plt_solution = list_plot(sol1.tolist(), plotjoined=True, color='brown')
    show( sum(plt_pkt)+plt1+plt2a+plt2b+plt_solution )

var('a,d',domain='real')
f(a,d)=(a-sqrt(  (1-a-d)^2+4*d) )*(1+a+d-sqrt((1-a-d)^2+4*d))/(2*a)
show(f)
plt2 = implicit_plot( f(a,d),(d,0,.61),(a,0,2),aspect_ratio=0.3,  figsize=(6, 3), axes_labels=[r'$d$','$a$'] )
plt2

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left( a, d \right) \ {\mapsto} \ \frac{{\left(a + d - \sqrt{{\left(a + d - 1\right)}^{2} + 4 \, d} + 1\right)} {\left(a - \sqrt{{\left(a + d - 1\right)}^{2} + 4 \, d}\right)}}{2 \, a}$$

b_ad = (trJ.subs(point3)).solve(b)[0].rhs().full_simplify()

b_ad.full_simplify()

(2*a^3 + (a - 1)*d^2 - 5*a^2 + (3*a^2 + a - 2)*d - (2*a^2 + (a + 1)*d - 3*a + 1)*sqrt(a^2 + 2*(a + 1)*d + d^2 - 2*a + 1) + 4*a - 1)/(a^2 + d^2 - sqrt(a^2 + 2*(a + 1)*d + d^2 - 2*a + 1)*(a - d - 1) - 2*a + 2*d + 1)

contour_plot(b_ad,(d,1e-4,.61),(a,0,2),aspect_ratio=.31,contours=srange(0,1.1,0.031),cmap='jet',plot_points=432)

b_ad.subs(d==1e-4).show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, a^{3} - 4.99970000000000 \, a^{2} - {\left(2 \, a^{2} - 2.99990000000000 \, a + 1.00010000000000\right)} \sqrt{a^{2} - 1.99980000000000 \, a + 1.00020001000000} + 4.00010001000000 \, a - 1.00020001000000}{a^{2} - \sqrt{a^{2} - 1.99980000000000 \, a + 1.00020001000000} {\left(a - 1.00010000000000\right)} - 2 \, a + 1.00020001000000}$$

plot(b_ad.subs(d==1e-5),(a,0.0,4.98))

\[\left( a, d \right) \ {\mapsto} \ \frac{{\left(a - \sqrt{{\left(a + d - 1\right)}^{2} + 4 \, d}\right)} {\left(a + d - \sqrt{{\left(a + d - 1\right)}^{2} + 4 \, d} + 1\right)}}{2 \, a}\]

|image|

Pod tą płaszczyzną

\[b = x_1\left[ \frac{\alpha x_1}{(x_1+d)^2} -1\right] = \phi (\alpha, d)\]

mamy cykl graniczny:

var('a,d')
f(a,d)=(a-sqrt(  (1-a-d)^2+4*d) )*(1+a+d-sqrt((1-a-d)^2+4*d))/(2*a)
show(f)
sol1=solve(f(a,d)==0, d)
show(sol1)
dm_expr=sol1[1].rhs()
assume(a>0)
f(a,0).full_simplify()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left( a, d \right) \ {\mapsto} \ \frac{{\left(a + d - \sqrt{{\left(a + d - 1\right)}^{2} + 4 \, d} + 1\right)} {\left(a - \sqrt{{\left(a + d - 1\right)}^{2} + 4 \, d}\right)}}{2 \, a}$$

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[d = -a - \sqrt{a^{2} + 4 \, a} - 1, d = -a + \sqrt{a^{2} + 4 \, a} - 1, d = -a + \sqrt{a^{2} + 2 \, {\left(a + 1\right)} d + d^{2} - 2 \, a + 1} - 1\right]$$

1/2*(2*a^2 - a*(2*abs(a - 1) + 1) - abs(a - 1) + 1)/a

forget()
assume(a<1)
f(a,0).simplify()

1/2*(a - sqrt((a - 1)^2) + 1)*(a - sqrt((a - 1)^2))/a

def b2(a,d):
    return (a-np.sqrt(  (1-a-d)**2+4*d) )*(1+a+d-np.sqrt((1-a-d)**2+4*d))/(2*a)
def dm(a):
    return -a + np.sqrt(a**2 + 4*a) - 1
def bm(a):

    if type(a) == np.ndarray:
        out = np.empty_like(a)
        a_lt_1 = a[a<1]
        a_gt_1 = a[a>=1]
        out[a<1] = 1/a_lt_1

        out[a>=1] =  0.5*(a_gt_1 - np.sqrt((a_gt_1 - 1)**2) + 1)*(a_gt_1 - np.sqrt((a_gt_1 - 1)**2))/a_gt_1
        return out
    else:
        if a>1:
            return 1/a
        else:
            return 0.5*(a - np.sqrt((a - 1)**2) + 1)*(a - np.sqrt((a - 1)**2))/a



from mpl_toolkits.mplot3d import Axes3D
import matplotlib
from matplotlib import cm
from matplotlib import pyplot as plt
step = 0.04
maxval = 1.0
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d',azim=134)
x = np.linspace(0.5,6,115)
y = np.linspace(0.00,1.,35)
X,Y = np.meshgrid(x,y)
# transform them to cartesian system
X,Y = X,Y*(np.sqrt(X**2+4*X)-(1.0+X))
#Y[:,i((X[7,0]**2+4*X[7,0])**0.5 - (1+X[7,0]) )
#0.99*d(:,i)*((a1.^2+4.*a1)^0.5 - (1+a1) )
Z = b2(X,Y)
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap=cm.jet)
#ax.plot_surface(X, Y, Z,  cmap=cm.jet)
#ax.plot_wireframe(X, Y, Z)
ax.set_zlim3d(0, 1)
ax.set_xlim3d(0, 3)
ax.set_ylim3d(0, .7)
ax.set_xlabel(r'$a$')
ax.set_ylabel(r'$d$')
ax.set_zlabel(r'$b(a,d)$')
ax.plot(x, dm(x), np.zeros_like(x), color=(.6,.1,.92),linewidth=3)
ax.plot(x,np.zeros_like(x), bm(x),  color='red',linewidth=3)
ax.plot([0],[0],[0])
ax.view_init(elev=35, azim=134)
plt.show()#plt.savefig("1.png")

\[\left[d = -\sqrt{a + 4} \sqrt{a} - a - 1, d = \sqrt{a + 4} \sqrt{a} - a - 1, d = -a + \sqrt{2 \, {\left(a + 1\right)} d + a^{2} + d^{2} - 2 \, a + 1} - 1\right]\]

|image|

Rozwiązania dążące do stabilnego stanu stacjonarnego: \(a \in (0, 0.5), \beta \gt 0, d \gt 0\)

var('x,y')
a, b, d = 0.3, 0.35, 0.1
T = srange(0,30,0.01)
sol2 = desolve_odeint(\
 vector([x*(1-x) - (a*x*y/(x+d)), b*y*(1-y/x)]),\
 [0.2, 0.5],T,[x,y])
line( zip ( T,sol2[:,0]) ,color='green', figsize=(6, 3), legend_label="x")+\
 line( zip ( T,sol2[:,1]) ,color='black',legend_label="y")

|image|

a, b, d = 0.3, 0.35, 0.1
F(x,y)=x*(1-x) - a*x*y/(x+d)
G(x,y)= b*y*(1-y/x)
T = srange(0,30,0.01)
sol1=desolve_odeint(vector([F,G]), [0.2,0.5], T, [x,y])
list_plot(sol1.tolist(), plotjoined=True,  figsize=(6, 3))

|image|

Rozwiązania dążące do stabilnego cyklu granicznego: \(a \gt 0.5 , \beta \gt 0, d \gt 0\)

var('x,y')
a, b, d = 1.3, 0.33, 0.1
T = srange(0,200,0.01)
sol2=desolve_odeint(\
 vector([x*(1-x) - (a*x*y/(x+d)), b*y*(1-y/x)]),\
 [0.2, 0.5],T,[x,y])
line( zip ( T,sol2[:,0]) ,color='green', figsize=(6, 3), legend_label="x")+\
 line( zip ( T,sol2[:,1]) ,color='black',legend_label="y")

|image|

a, b, d = 1.3, 0.33, 0.1
F(x,y)=x*(1-x) - a*x*y/(x+d)
G(x,y)= b*y*(1-y/x)
T = srange(0,250,0.01)
sol1=desolve_odeint(vector([F,G]), [0.2,0.5], T, [x,y])
list_plot(sol1.tolist(), plotjoined=True,  figsize=(6, 3))

|image|