352 - Topic 05 - Equilibria and Phase Lines -- Sage

Solving Separable Differential Equations

John Travis

Mississippi College

Blanchard, Devaney and Hall, 4th edition.

This worksheet deals with simple autonomous differential equations of the form $y' = f(y)$. That is, there are no explicit references to the variable $t$ in f(t,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. However, since changes in $t$ will not affect $y'$, slope fields will have constant slopes along horizontal lines. Hence, a phase line represents the slopes on a vertical slice of this autonomous DE slope field.

First, one will need to let this worksheet know what particular the differential equation $y' = f(y)$ is to be investigated by specifying $f(y)$ below. To see a different result throughout this worksheet, one should change $f(y)$ here and then be certain to reevaluate the cell.

# Worksheet still under development. Use at your own risk! 8-o

%auto def f(y): # f = y*(1-y) f = (y^2)*(y^2-4) # f = y^2-y+1 # f = y - y^3 return f

Equilibria occur when a solution $y(t)$ for the DE remains constant. If so, then $y'=0$ and hence $f(y) = 0$.

Also, one can determine whether an equilibrium is a sink or source using the sign of $f'(y)$.

# Determining the location and nature of equilibria # We can use the sign of the derivative of f to indicate sink, source, or possible node fp = diff(f(y),y,1) # Need to set this up properly. Just playing right now. eqn = f(y)==0 html('\nSolving \n'+'$'+latex(eqn)+'$') equil = solve(eqn,y,solution_dict=True) if equil==[]: print html("\nSolver did not locate any possible equilibria. You might want to try by hand.") else: print html("\nEquilibria occur at the following parameter values:") Eq_list = [] for value in equil: equ = value[y] if sgn(fp(y=equ))==1: print equ,'Source' elif sgn(fp(y=equ))==-1: print equ,'Sink' else: print equ,'Node?' Eq_list.append(value[y]) Eq_list.sort()

Solving {\left(y^{2} - 4\right)} y^{2} = 0

Equilibria occur at the following parameter values:


-2 Sink
2 Source
0 Node?

Solving {\left(y^{2} - 4\right)} y^{2} = 0

Equilibria occur at the following parameter values:


-2 Sink
2 Source
0 Node?

rangeC = max(max(Eq_list)-min(Eq_list),1) startC = min(Eq_list)-rangeC/2 endC = max(Eq_list)+rangeC/2 previous = startC C_pts = [startC] for val in Eq_list: 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(): 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_list)==0: bottom = 0 top = 1 ev_pts.append(0.5) else: bottom = min(Eq_list)-1 top = max(Eq_list)+1 eq_list = [] previous = bottom for p in Eq_list: eq_graph += point((0,p),size=40,color="red") eq_graph += text(str(p),(-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 = (0,pt) p1 = (0+.10,pt+f(pt)/2) # slopes += line([p0,p1],color="green",marker=">") # slopes x 5 for emphasis if (f(pt)>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([(0,bottom),(0,top)],color="blue",thickness=2) Yaxis += text(0,(0,bottom-0.05),fontsize=8,vertical_alignment='top',rotation=-10) phase_line = Yaxis+slopes+eq_graph return phase_line

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.

show(phase_line(),xmin=startC,xmax=endC,axes=False)

or with the regions smoothly drawn...

For the fun of it, let's let sage attempt to get the exact solution for the input differential equation.

y=function('y',t) de = diff(y,1)-f(y) de_soln = desolve(de,y,ivar=t) show(de_soln) solve(de_soln,y,t)

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{\log\left(y\left(t\right) - 2\right) y\left(t\right) - \log\left(y\left(t\right) + 2\right) y\left(t\right) + 4}{16 \, y\left(t\right)} = c + t

([y(t) == 4/(16*c + 16*t - log(y(t) - 2) + log(y(t) + 2))], [1])

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{\log\left(y\left(t\right) - 2\right) y\left(t\right) - \log\left(y\left(t\right) + 2\right) y\left(t\right) + 4}{16 \, y\left(t\right)} = c + t

([y(t) == 4/(16*c + 16*t - log(y(t) - 2) + log(y(t) + 2))], [1])

352 - Topic 05 - Equilibria and Phase Lines

3404 days ago by Professor352