#predator-prey equations

#continuous system

#x = population of flounder
#y = population of striped-bass

#dxdt = growth rate of flounder (x) over time (t)
#dydt = growth rate of striped-bass (y) over time (t)

#alpha = the coefficient of the growth rate of the flounders when they survive from natural resources.

#beta  = the coefficient of growth rate of the flounders in the presence of striped-bass.

#theta = the coefficient of the growth rate when the striped-bass have no flounders to eat and ends up dying off.

#gamma = the coefficient of the growth rate of the striped-bass when they eat the flounders for survival.

#Constraints
# alpha*x   > 0  ; assuming there is unlimited natural resources to feed the flounders, flounders can have no limit to their population (carrying capacity)
# beta*x*y  < 0  ; due to the hunger of the striped-bass, the striped-bass will be consuming the flounders until the flounders are near extincted
# theta*y   < 0  ; assuming there is only two species in the pond, the strip-bass will die off in an exponential rate as soon as they cannot find any more flounder to eat
# gamma*x*y > 0  ; the population of the strip-bass will be increased whenever the flonders are still available

#inital values:
alpha = 0.5
beta  = 0.02
theta = 0.4
gamma = 0.004
#(The data above is given by a Postdoctoral Research Associate Courtney Davis from Utah)
#link: http://www.math.utah.edu/~davis/1180FILES/lab2.pdf

#Model:
dxdt(x,y) = alpha*x - beta*x*y
dydt(x,y) = -theta*y + gamma*x*y

solve([dxdt==0,dydt==0],x,y) #eq points (0,0) and (100,25)

def f1(x,y):
    return alpha*x - beta*x*y

def f2(x,y):
    return -theta*y + gamma*x*y

def simulation(delT,N,x10,y10):
    x1= [x10]
    y1= [y10]
    for i in range(N):
        # most recent values
        x1old= x1[-1]
        y1old= y1[-1]
        # next set of values
        x1new= x1old + delT*f1(x1old,y1old)
        y1new= y1old + delT*f2(x1old,y1old)
        # append values to lists
        x1.append(x1new)
        y1.append(y1new)
    return [x1,y1]
s = simulation(0.01,10000,100,50)
u = simulation(0.0115,10000,100,50)
p = simulation(0.0085,10000,100,50)

FL = list_plot((s[0]),axes_labels = ['$time(years)$', '$population$'],legend_label='flounder',color = "green")
SB = list_plot((s[1]),axes_labels = ['$time(years)$', '$population$'],legend_label='striped-bass',color = "pink")
show(FL+SB)
#If the population of folunders decreased, as you would guess the population of the striped-bass would also go down, becasue they wouldn't have nothing to eat and would end up dying.

#when there's less striped-bass in the pond, the flounder's growth rate will increase.

g1a = list_plot(s[0], axes_labels = ['$time(years)$', '$population$'],legend_label='delT = 0.01',color = "green")  #normal delT  0.01
g1b = list_plot(s[1], axes_labels = ['$time(years)$', '$population$'],legend_label='delT = 0.01',color = "pink")  #normal delT  0.01
g2a = list_plot(u[0], axes_labels = ['$time(years)$', '$population$'],legend_label='delT = 0.0115',color = "darkgreen") #bigger delT  0.0115
g2b = list_plot(u[1], axes_labels = ['$time(years)$', '$population$'],legend_label='delT = 0.0115',color = "orange") #bigger delT  0.0115
g3a = list_plot(p[0], axes_labels = ['$time(years)$', '$population$'],legend_label='delT = 0.0085',color = "lightgreen") #smaller delT 0.0085
g3b = list_plot(p[1], axes_labels = ['$time(years)$', '$population$'],legend_label='delT = 0.0085',color = "lightpink") #smaller delT 0.0085

show(g1a+g2a+g3a,axes_labels = ['$time (years)$', '$Flounder$ $Population$'])
show(g1b+g2b+g3b,axes_labels = ['$time (years)$', '$Striped-bass$ $Population$'])

#As the time step gets smaller, it takes a longer time for both prey(flounder) and predator(striped-bass) to reach the equilibrium point.
#As the time step gets bigger, it takes shorter time for both prey(flounder) and predator(striped-bass) to reach the equilibrium point
#It's because the action is still taking effect.
#The situation is either the flounders are still reprodcucing back to the stable level after being overeaten by the striped-bass or the striped-bass is still dying becuase of no more preys to feed on.

#finding Eigenvalues 
f = (dxdt,dydt)
A = jacobian(f,(x,y))
show(A)

A(0,0).eigenvalues() #The real parts are not all negative ; unstable eq point
A(100,25).eigenvalues() #The real parts are not all negative ; unstable eq point
abs(n(1/5*I*sqrt(5)))
abs(n(-1/5*I*sqrt(5)))
n(-2/5)

p1 = implicit_plot(dxdt==100,(x,0,650),(y,0,100))
p2 = implicit_plot(dydt==25,(x,0,650),(y,0,100), color = 'red' )
p3 = plot_vector_field((dxdt,dydt),(x,0,650),(y,0,100))
show(p1+p2+p3,aspect_ratio = 3,axes_labels = ['$flounder$', '$striped-bass$'])