SharedLS 30A F17 NOW / Lac Operon Interactive.sagewsOpen in CoCalc
# 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