### STOCHASTIC
### Initialize some stuff
var('t dt tf Tprint')
var('at m i')
X = [0,2,10] ## holds [species], ICs (B,M,F)
k = [3.5,0.15,0.15,.1,.5,.1] ## holds rate constants
a = [0,0,0,0,0,0] ## holds calc'd rxn rates (once they ARE calculated)
r = [0,0,0,0,0,0] ## holds the random numbers (once they ARE created)
sol = [] ## holds the solution
t = 0 ## time, start
tf = 150 ## time, end
dt = 0.01 ## time, step (observational frequency)
m = len(a) ## the number of reactions under consideration
Tprint = t
### Loop 0: continues until target time ###
while (t<tf):
at = 0 ## 'overall reaction rate'; takes into account all reactions
## calculate the current reaction rates
a[0] = k[0]*X[2]
a[1] = k[1]*X[0]
a[2] = k[2]*X[0]
a[3] = k[3]*X[2]
a[4] = k[4]*X[1]
a[5] = k[5]*X[0]^2
for j in range(0,m):
at = at + a[j] ## calculate the 'overall reaction rate'
r[j] = random() ## generate the needed random numbers
t = t+(1/at)*ln(1/r[0]) ## next rxn occurs in time increment 'tau'
while(t>=Tprint): ## until the next reaction occurs..
sol.append([Tprint,X[0],X[1],X[2]]) ## the concentrations stay the same
Tprint=Tprint+dt ## as time increments
i = 1 ## (which reaction is being investigated)
while (i <= m): ## and when a reaction DOES occur...
sum = 0
for l in range(0,i):
sum = sum + a[l] ## sum together the rates (for normalization)
if (sum > r[1]*at): ## randomly select a value (0, total rate)
if (i == 1): ## if value falls within rate #1, rxn #1 occurs
X[0] = X[0] + 1
X[1] = X[1]
X[2] = X[2]
break
elif (i == 2): ## if value falls within rate #2, rxn #2 occurs
X[0] = X[0] - 1
X[1] = X[1] + 1
X[2] = X[2]
break
elif (i == 3): ## if value falls within rate #3, rxn #3 occurs
X[0] = X[0] - 1
X[1] = X[1]
X[2] = X[2] + 1
break
elif (i == 4): ## if value falls within rate #4, rxn #4 occurs
X[0] = X[0]
X[1] = X[1]
X[2] = X[2] - 1
break
elif (i == 5): ## if value falls within rate #5, rxn #5 occurs
X[0] = X[0]
X[1] = X[1] - 1
X[2] = X[2]
break
elif (i == 6): ## if value falls within rate #6, rxn #6 occurs
X[0] = X[0] - 1
X[1] = X[1]
X[2] = X[2]
break
else:
show("something went wrong...")
i = i + 1
#show(sol)
### Simple Plotting
sol0 = []
sol1 = []
sol2 = []
sol3 = []
for q in range(0,len(sol)):
sol0.append(sol[q][0])
sol1.append(sol[q][1])
sol2.append(sol[q][2])
sol3.append(sol[q][3])
a = list_plot(zip(sol0,sol1), color='mediumpurple')
a += list_plot(zip(sol0,sol2), color='dodgerblue')
a += list_plot(zip(sol0,sol3), color='indianred')
#show(a)
|
|
### DETERMINISTIC
## This sheet simulates the population dynamics in a 'rabbit farm'
## VARIABLES
var('t B M F B0 M0 F0')
var('k1 k2 k3 k4 k5 k6')
var('r1 r2 r3')
## INITIAL CONDITIONS
B0 = 0 ## Baby Rabbits
M0 = 2 ## Male Rabbits
F0 = 10 ## Female Rabbits
## CONSTANTS
k1 = 3.5
k2 = 0.15
k3 = 0.15
k4 = 0.1
k5 = 0.5
k6 = 0.1
## CALCULATION PARAMETERS
ivar = t
dvars = [B,M,F]
ics = [0,B0,M0,F0]
end_points = 150
stepsize = 0.1
steps = end_points/stepsize
## EQUATIONS
r1 = (diff(B,t) == k1*F-k2*B-k3*B-k6*B^2)
r2 = (diff(M,t) == k2*B-k5*M)
r3 = (diff(F,t) == k3*B-k4*F)
## NUMERICAL SOLUTION OF A SERIES OF DIFFERENTIAL EQUATIONS
des = [r1.rhs(), r2.rhs(), r3.rhs()]
sol = desolve_system_rk4(des,dvars,ics,ivar,end_points=end_points,step=stepsize)
## Process the output into a form that can be graphed
sols_1=[]
sols_2=[]
sols_3=[]
for i in range(steps):
sols_1.append([sol[i][0],sol[i][1]])
sols_2.append([sol[i][0],sol[i][2]])
sols_3.append([sol[i][0],sol[i][3]])
################################################
#### Unnecessarily Fancy Plotting Stuff ####
################################################
## Create a plot object
#a = plot([])
## Set the plot parameters
title = "Population Dynamics" ## Graph Title
xmin = 0 ## X minimum
xmax = end_points ## X maximum
ymin = 0 ## Y minimum
ymax = 100 ## Y maximum
## Add a title to the plot
a += text(title,(xmax/1.8,ymax),color='black',fontsize=15)
## Add the desired lines to the plot
a += list_plot(sols_1,color='purple',legend_label='B')
a += list_plot(sols_2,color='blue',legend_label='M')
a += list_plot(sols_3,color='red',legend_label='F')
## For more information on plots in general, evaluate 'plot?'
## For a list of legend options, evaluate 'a.set_legend_options?'
## For a list of Sage predefined colors, evaluate 'sorted(colors)'
## View predefined colors: http://en.wikipedia.org/wiki/X11_color_names
## Additional plot options
#a.set_axes_range(xmin,xmax,ymin,ymax)
a.axes_labels(['time','number'])
a.axes_label_color('grey')
a.set_legend_options(ncol=3,borderaxespad=5,back_color='whitesmoke',fancybox=true)
show(a)
|
|
|
|