Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
| Download
Views: 491
# 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