var('t')
var('R F')
var('R0 F0')
var('a b c d')
## CONSTANTS
a = 0.04 # rabbit birthrate
b = 0.0005 # probability of predation per encounter
c = 0.1*0.0005 # probability of predation per encounter, rabbit -> fox conversion efficiency
d = 0.2 # death rate of foxes
## INITIAL CONDITIONS
R0 = 5000 #initial number of rabbits
F0 = 200 #initial number of foxes
## DIFFERENTIAL EQUATIONS
de_R = (diff(R,t) == a*R - b*R*F)
de_F = (diff(F,t) == c*R*F - d*F)
## CALCULATION PARAMETERS
end_points = 500
stepsize = 1.0
steps = end_points/stepsize
ics = [0,R0,F0]
sol = []
des = [de_R.rhs(), de_F.rhs()]
## NUMERICAL SOLUTION OF DIFFERENTIAL EQUATIONS
sol = desolve_system_rk4(des,[R,F],ics,ivar=t,end_points=end_points,step=stepsize)
|
|
## Clean up to graph
sol_t=[]
sol_R=[]
sol_F=[]
for i in range(steps):
sol_t.append(sol[i][0])
sol_R.append(sol[i][1])
sol_F.append(sol[i][2])
a = plot([])
a += line(zip(sol_R,sol_F))
a += plot_vector_field((des[0], des[1]), (R,0,7000), (F,0,225),xmin=1500,color='orange')
a.axes_labels(['Rabbits','Foxes'])
show(a)
|
|
q = list_plot(zip(sol_t,sol_F), color='orange')
q = q + list_plot(zip(sol_t,sol_R), color='green')
show(q)
|
|
