352 - Topic 07 - Bifurcations

44 days ago by Professor352

John Travis

Mississippi College

This worksheet deals with simple autonomous differential equations of the form $y' = f(y)$.  In sage, one can obtain an entire slope field using built-in functions but this page is focused on creating the "phase line" for autonomous DEs.  The user is allowed to utlize one parameter ($C$ is reserved for that use), determine which value(s) for $C$ may yield possible bifurcations and exhibit these possible bifurcations visually.

Finally, a complete bifurcation diagram is created showing the nature of solutions for varying values of the parameter.

First, one will need to enter the differential equation by specifying $f(y)$ below.  To see a different result throughout this worksheet, one should change $f(y)$ here.
reset() var('x,y,C') def f(y,C): # enter an autonomous function of y only # f = y*(1-y)-C # this one works well. # f = (y^2-C)*(y^2-4) # f = y^2-C*y+1 # this one works ok. f = C*y - C*y^3 # this one works ok. But gives a bifurcation at C = 0 as the slider cell below illustrates. Doesn't give the color map well. # f = C*y # f = y^6-2*y^3+C # f = (y-C)*(y-2*C)*(y-3*C) return f 
       
f(y,C) 
        
-C*y^3 + C*y
-C*y^3 + C*y

Equilibria occur when a solution $y(t)$ for the DE remains constant.  If so, then

$y'=0$ and hence $f(y) = 0$. 

Since once can determine whether an equilibrium is a sink or source using the sign of $f _y(y)$, the nature of the collection of equilibria will only change (under reasonable assumptions) when $f _y(y)$ is also zero.
# Using the Theorem to determine the location of bifurcation values fp = diff(f(y,C),y,1) # Need to set this up properly. Just playing right now. eqns = [ f(y,C)==0, fp==0 ] print 'System of Equations to solve for C is' show(html('$'+latex(eqns[0])+' \\text{ and } '+latex(eqns[1])+'$')) bif=solve(eqns,(y,C),solution_dict=True,ring=RR) # solve the system of (nonlinear) equations for C if bif==[]: print html("\nSolver did not locate any possible bifurcation locations. You might want to try by hand.") else: print html("\nPossible Bifurcations may occur at the following parameter values:") Clist = [] for value in bif: Clist.append(value[C]) print value[C] Clist.sort() ####################### # Below, some things are collected and set up for later use # if (Clist == []): # bottom = 0 # top = 1 # ev_pts.append(0.5) # else: rangeC = max(max(Clist)-min(Clist),1) startC = min(Clist)-rangeC/2 endC = max(Clist)+rangeC/2 previous = startC C_pts = [startC] for val in Clist: mid = (val+previous)/2 # using the midpoints between equilibria for slope evaluation previous = val C_pts.append(mid) C_pts.append(val) C_pts.append( (endC+previous)/2 ) C_pts.append(endC) def phase_line(C): eq_pts = f(y,C).roots(ring=RR, multiplicities=False) eq_pts.sort() eq_graph = Graphics() ev_pts = [] # If there are no real equilibria, you only need to evaluate f(y) once to determine all trajectory directions if len(eq_pts)==0: bottom = 0 top = 1 ev_pts.append(0.5) else: bottom = min(eq_pts)-1 top = max(eq_pts)+1 eq_list = [] previous = bottom for p in eq_pts: eq_graph += point((C,p),size=40,color="red") eq_graph += text(str(p),(C-0.02,p),rgbcolor='black',fontsize=6,horizontal_alignment='right') eq_list.append(p) mid = (p+previous)/2 # using the midpoints between equilibria for slope evaluation previous = p ev_pts.append(mid) ev_pts.append( (top+previous)/2 ) # Finally, create the phase line using the values determined above bottom = min(min(ev_pts)-1,bottom) slopes = Graphics() for pt in ev_pts: p0 = (C,pt) p1 = (C+.10,pt+f(pt,C)/2) # slopes += line([p0,p1],color="green",marker=">") # slopes x 5 for emphasis if (f(pt,C)>0): slopes += text("Increasing",p0,color="green",fontsize=6,horizontal_alignment='center',rotation=30) else: slopes += text("Decreasing",p0,color="red",fontsize=6,horizontal_alignment='center',rotation=-30) Yaxis = line([(C,bottom),(C,top)],color="blue",thickness=2) Yaxis += text(C.n(4),(C,bottom-0.05),fontsize=8,vertical_alignment='top',rotation=-10) phase_line = Yaxis+slopes+eq_graph return phase_line 
        
System of Equations to solve for C is


Possible Bifurcations may occur at the following parameter values:
0
System of Equations to solve for C is


Possible Bifurcations may occur at the following parameter values:
0

Since the DE is autonomous, we can determine the qualitative nature of all solutions by plotting the "phase line".  Experiment with various values of the parameter and look for dramatic changes as you pass through the possible bifurcation values indicated above.
@interact def _(C=slider(startC,endC,0.01,startC)): show(phase_line(C),xmin=startC,xmax=endC,axes=False) 
       
C 

Click to the left again to hide and once more to show the dynamic interactive window

The complete bifurcation diagram is given below.  This could use a bit of work to make it look prettier automatically when you have widely varying equilibria for various $C$ values.  Oh well, the secret is to only assign problems that look good in what has been written!  :)
Biff_diagram = Graphics() for cval in C_pts: Biff_diagram += phase_line(cval) # A contour plot| can quickly give the general regions for increasing/decreasing # Comment out the next statement if you only want the collection of detailed phase lines. show(Biff_diagram,axes=False) 
        
Traceback (click to the left of this block for traceback)
...
cypari2.handle_error.PariError: zero polynomial in roots
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_70.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("QmlmZl9kaWFncmFtID0gR3JhcGhpY3MoKQpmb3IgY3ZhbCBpbiBDX3B0czoKICAgIEJpZmZfZGlhZ3JhbSArPSBwaGFzZV9saW5lKGN2YWwpCiMgIEEgY29udG91ciBwbG90fCBjYW4gcXVpY2tseSBnaXZlIHRoZSBnZW5lcmFsIHJlZ2lvbnMgZm9yIGluY3JlYXNpbmcvZGVjcmVhc2luZwojICBDb21tZW50IG91dCB0aGUgbmV4dCBzdGF0ZW1lbnQgaWYgeW91IG9ubHkgd2FudCB0aGUgY29sbGVjdGlvbiBvZiBkZXRhaWxlZCBwaGFzZSBsaW5lcy4gIApzaG93KEJpZmZfZGlhZ3JhbSxheGVzPUZhbHNlKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmp0FmzuE/___code___.py", line 4, in <module>
    Biff_diagram += phase_line(cval)
  File "/tmp/tmpv6usl3/___code___.py", line 44, in phase_line
    eq_pts = f(y,C).roots(ring=RR, multiplicities=False)
  File "sage/symbolic/expression.pyx", line 11468, in sage.symbolic.expression.Expression.roots (build/cythonized/sage/symbolic/expression.cpp:58525)
  File "sage/rings/polynomial/polynomial_element.pyx", line 7660, in sage.rings.polynomial.polynomial_element.Polynomial.roots (build/cythonized/sage/rings/polynomial/polynomial_element.c:60330)
  File "cypari2/gen.pyx", line 4118, in cypari2.gen.Gen.polroots
  File "cypari2/handle_error.pyx", line 213, in cypari2.handle_error._pari_err_handle
cypari2.handle_error.PariError: zero polynomial in roots

or with the regions smoothly drawn...
Biff_diagram = Graphics() for cval in C_pts: Biff_diagram += phase_line(cval) # A contour plot can quickly give the general regions for increasing/decreasing # Comment out the next statement if you only want the collection of detailed phase lines. Biff_background = contour_plot(f(y,C),(C,startC,endC),(y,-3,3),contours=[-100,0,100],cmap="spring",fill=True,linewidths=3) show(Biff_diagram+Biff_background,axes=False,figsize=(6,20)) # need to get rid of those ugly borders...trouble finding the specific option to turn them off. # I'd also like to significantly lighted the hue of the background so that the lines stand out. 
        
Traceback (click to the left of this block for traceback)
...
cypari2.handle_error.PariError: zero polynomial in roots
Traceback (most recent call last):    show(Biff_diagram+Biff_background,axes=False,figsize=(6,20))
  File "", line 1, in <module>
    
  File "/tmp/tmpa0LKKr/___code___.py", line 5, in <module>
    Biff_diagram += phase_line(cval)
  File "/tmp/tmpv6usl3/___code___.py", line 44, in phase_line
    eq_pts = f(y,C).roots(ring=RR, multiplicities=False)
  File "sage/symbolic/expression.pyx", line 11468, in sage.symbolic.expression.Expression.roots (build/cythonized/sage/symbolic/expression.cpp:58525)
  File "sage/rings/polynomial/polynomial_element.pyx", line 7660, in sage.rings.polynomial.polynomial_element.Polynomial.roots (build/cythonized/sage/rings/polynomial/polynomial_element.c:60330)
  File "cypari2/gen.pyx", line 4118, in cypari2.gen.Gen.polroots
  File "cypari2/handle_error.pyx", line 213, in cypari2.handle_error._pari_err_handle
cypari2.handle_error.PariError: zero polynomial in roots