CoCalc Shared FilesLS 30A F17 NOW / Lac Operon Interactive.sagews
Views : 23
# Lac Operon Interactive

var("x")
@interact
def bifurcation(r=(0.2, 1, 0.01)):
a = 0.001
inflow = plot(((a+x^2)/(1+x^2)), (x,0,6), ymin=0, ymax=1.2, axes_labels = ["X",""], legend_label = "Inflow")
outflow = plot(r*x, (x,0,6), ymin=0, ymax=1.2, color="red", legend_label = "Outflow")
label = text(" r= " + str(r), [1,0.9], color = "black") # Plot input, output
show(inflow + outflow + label)

x
Interact: please open in CoCalc
# Labels equilibria (but doesn't use a) - fancier version below!
var("x")
@interact
def bifurcation(r=(0.2, 1, 0.01)):
assume(x >= 0) # Restricts solve to find only solutions that are real numbers >= 0
assume(x, "real")
eqs = point([0,0],size=1)
soln = solve((x^2)/(1+x^2) - r*x == 0, x, solution_dict = True) # Solve to get equilibria values of r
for s in soln:
eqs += point([s[x].n(30), r*(s[x].n(30))], size=50, color = "black") # Adds dots for equilibria # Note: .n(30) means numerical approximation with 30 bits
inflow = plot(((x^2)/(1+x^2)), (x,0,6), ymin=0, ymax=1.2, axes_labels = ["X",""], legend_label = "Input")
outflow = plot(r*x, (x,0,6), ymin=0, ymax=1.2, color="red", legend_label = "Output")
label = text(" r= " + str(r), [1,0.9], color = "black") # Plot input, output
show(inflow + outflow + label + eqs)

x
Interact: please open in CoCalc
# Alternative demo
def list_classify(funct, var):
s = solve(funct,var)
s = [i.rhs() for i in s]
s.sort()
output = []
for i in srange(0,len(s)):
state = ""
if s[i]== max(s):
upper = s[i]+0.1
else:
upper = (s[i]+s[i+1])/2
if s[i]== min(s):
lower = s[i]-0.1
else:
lower =(s[i]+s[i-1])/2
if funct(upper)<0 and funct(lower)>0:
state = "stable"
elif funct(upper)>0 and funct(lower)<0:
state = "unstable"
else:
state= "semi-stable"
output.append((s[i],state))
return output

def bifurcation_interact(inputs,outputs, var, lowx,highx, lowerbound,upperbound):
rlist=[]
zerolist=[]
lenlist=[]
directionlist =[]
@interact
def display(rval = ("Parameter $r$",slider(lowerbound, upperbound)),showtable = checkbox(False,"Display info"),auto_update = False):
rlist.append(rval)
eqs= list_classify(inputs(var,rval)-outputs(var,rval),var)
s = solve(inputs(var,rval)-outputs(var,rval), var)
s = [i.rhs().n(5) for i in s]
s.sort()
zerolist.append(str(s).strip("[]"))
lenlist.append(len(eqs))
pointsum= Graphics()
for i in eqs:
if i[1]== "stable":
pointsum+=point([i[0],outputs(i[0],rval)], color = "green", size = 40)
elif i[1] =="unstable":
pointsum+=point([i[0],outputs(i[0],rval)], color = "yellow", size = 40)
else:
pointsum+=point([i[0],outputs(i[0],rval)], color = "orange", size = 40)
if len(s)<2:
directionlist.append("Too large - Choose smaller value")
elif len(s)==2:
directionlist.append("At bifurcation point")
elif len(s) >2:
directionlist.append("Too small - Choose larger value")
p = list_plot([0], color = "yellow", legend_label = "Unstable Eq Point")+ list_plot([0], color = "orange", legend_label = "Semi-stable Eq Point")+list_plot([0], color = "green", legend_label = "Stable Eq Point")
p+= plot(inputs(var,rval),(var,lowx,highx), color = "blue", axes_labels = ["X","Rate of Change for X"], legend_color = "black", legend_label= "Input for $r$: "+ str(rval.n(5)))+plot(outputs(var,rval),(var,lowx,highx), color = "red", axes_labels = ["X","Rate of Change for X"], legend_color = "black", legend_label= "Output for $r$: "+ str(rval.n(5)))+pointsum+ list_plot([], color = "green", legend_label = "Stable Eq Point")
show(p)
if showtable:
print table(columns = [["Rvalues"]+rlist, ["Eq Vals"]+zerolist, ["Number of Eqs"]+ lenlist,["Changes"]+directionlist], frame = true, header_row = true)
x,r =var("x r")
assume(x>=0)
assume(x,"real")
a = 0.01
inp(x,r)= (x^2)/(1+x^2) # Doesn't work for (a + x^2)/(1+x^2)
out(x,r)= r*x
bifurcation_interact(inp,out,x,0,4,0.3,0.7)

Interact: please open in CoCalc