%% 17. Bifurcation

%%
% <latex>
% \setcounter{page}{212}
% \parindent=.6cm
% \parskip=-1em
% </latex>

%%
% As parameters vary, the solutions of an
% ODE may change in structural properties.  This is the
% subject of bifurcation theory.

%%
% <latex>
% For a prototypical example, consider a marble resting
% on a surface whose height as a function of horizontal position $y$ is
% $$ h(y) = -c \kern .5pt y^2,  \eqno (17.1) $$
% where $c$ is a constant.
% For any value of $c$, the marble is in equilibrium at
% the fixed point $y=0$: there
% are no net forces on it.  But it
% is clear that the equilibrium will be stable
% if $c<0$, neutrally stable if $c=0$, and unstable if $c>0$.\footnote{As
% mentioned in the first footnote of Chapter~15, there is some
% inconsistency between formal and more casual uses of terms related to
% stability.  Our use of ``stable'' and ``unstable'' here matches
% the formal definitions given in that chapter.  We did not
% define ``neutrally stable,'' which here corresponds to a
% case that is unstable, but with unstable trajectories growing just
% algebraically, not exponentially.}  The critical value $c=0$ is the
% {\bf bifurcation point} at which the behavior switches from
% one class to the other.
% </latex>
ODEformats, y = chebfun('y',[-1.3 1.3]);
marble = .14i + .12*exp(pi*1i*(0:60)/30);
s = []; s{1} = ' '; s{2} = ' neutrally '; s{3} = ' un';
for k = 1:3
  c = -.4+.2*k; surface = -c*y^2;
  subplot(1,3,k), plot(surface,'k',LW,1.2), hold on
  fill(real(marble),imag(marble),[0 0 0]), hold off
  axis([-1.3 1.3 -1.2 1.4]), axis square, axis off
  text(0,1.2,[sprintf('$c = %3.1f$,',c) s{k} 'stable'],FS,10,HA,CT)
end
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Let us examine the ODE implicit in these pictures.  
% We make the physics as elementary as possible, regarding
% the marble simply as a point mass moving along a 
% surface in a gravitational field with constant $1$.
% (Thus this marble has no internal angular momentum.)
% Since the slope of the surface at position $y$ is
% $h'(y) = -2\kern .3pt c\kern .3pt y$, the
% ODE for its position $y(t)$ as a function of time
% $t$ is
% $$ y'' = 2\kern .3pt c\kern .3pt y. \eqno (17.2) $$
% Following the pattern of Chapter~15,
% we plot some typical trajectories in a neighborhood
% of the fixed point $(0,0)$ in the $(y,y')$ phase plane.
% </latex>
clf, th = (pi/6)*(1:12)+.000001; u0 = cos(th); v0 = sin(th);
cc = chebfun('exp(1i*pi*x)'); L = chebop(0,2.5);
for j = 1:3
  subplot(1,3,j), plot(0,0,'.k',MS,6), hold on, plot(cc,'k',LW,.5)
  c = -.4+.2*j; L.op = @(t,y) diff(y,2) - 2*c*y;
  for k = 1:12
    L.lbc = [u0(k); v0(k)] + 1e-4*[1;1];
    y = L\0; arrowplot(y,diff(y),CO,ivp,MS,3)
  end
  axis([-2 2 -2 2]), axis square, hold off
  set(gca,XT,-2:2:2,YT,-2:2:2)
  title([sprintf('$c = %3.1f$,',c) s{j} 'stable'],FS,11)
end
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The three images show a center, a degenerate
% node, and a saddle point, corresponding to $c=-0.2,$
% $0,$ and $0.2.$
% The Jacobian matrices are
% $$ {\bf J}_*^{} = \pmatrix{0 & 1\cr -0.4 & 0}, ~
% \pmatrix{0 & 1\cr 0 & 0}, ~
% \pmatrix{0 & 1\cr 0.4 & 0}, $$
% with eigenvalues
% $$ \{\lambda_1,\lambda_2\} =
% \{di, -di\}, ~ \{0,0\}, ~ \{d,-d\}, \quad  d = \sqrt{0.4}\kern 1pt . $$
% </latex>

%%
% <latex>
% As an example of the interpretation of these results, what is
% the physics of the third problem of Fig.~17.1,
% the unstable configuration with the marble at the top of a hill?  
% From the phase plot we can see the answer.  Almost any nonzero
% initial position in the phase plane will lead to divergence to
% infinity in a northeasterly or southwesterly direction: to be
% precise, in the direction of a positive or negative multiple
% of the eigenvector $(1, \sqrt{0.4}\kern 1pt)$ of
% ${\bf J}_*^{}$.  The exception is the special situation in which the marble
% begins with position $(y,y')$ exactly equal to a multiple of
% the other eigenvector, $(1,-\sqrt{0.4}\kern 1pt )$.\ \ In
% this case the marble is moving up
% the hill with just the right amount
% of energy to reach the top with velocity zero at $t=\infty$. 
% (In general, the set of initial points that converge
% to a fixed point of a dynamical system is called the
% {\bf stable manifold} of that point.)
% Of course, this is a physically unstable situation, which would be
% undone by the slightest perturbation.
% </latex>

%%
% <latex>
% For problems depending on a parameter, it is common
% to draw a {\bf bifurcation diagram} indicating the
% dependence of some measure of a fixed point --- an
% equilibrium --- on the parameter.  For the marble on the surface,
% the bifurcation diagram is very simple.
% </latex>
clf, plot([-2 0],[0 0],'-k',[2 0],[0 0],'--k')
title('Fig.~17.3.~~Bifurcation diagram for (17.2)',FS,11)
xlabel('$c$',FS,10,IN,LT), ylabel('$y_*^{}$',FS,10,IN,LT)
text(-0.5,0.12,'stable',FS,10,HA,CT)
text(0.6,0.12,'unstable',FS,10,HA,CT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% For each value of $c$ there is a single fixed point
% $y_*^{}$, and the diagram plots the dependence of $y_*^{}$
% on $c$.\ \ For this very simple problem, the dependence is
% trivial since the fixed point is $y_*^{}=0$
% for all $c$.  For $c<0$ the equilibrium is stable, and 
% the curve is shown solid.  For $c>0$ it is unstable, and
% the curve is shown dashed.  In general, dashed curves
% in bifurcation diagrams correspond to unstable steady states
% of a system, which one would not ordinarily expect to observe in an
% experiment.
% </latex>

%%
% <latex>
% There is an interesting way to trace parts of a bifurcation
% diagram dynamically: set up a time-dependent problem in which
% the parameter of interest varies slowly in time, slowly enough that at each
% $t$, the behavior is approximately that of a 
% constant-coefficient system.   (In quantum physics this would
% be called an {\em adiabatic transition\/} from one parameter value
% to another.)
% For example, suppose we modify (17.2) to the equation
% $$ y'' = 2\kern .3pt c(t)y, \quad c(t) = -2+t/150. \eqno (17.3) $$
% Here is the solution with $y(0)=0.02$, $y'(0)=0$.
% </latex>
L = chebop(0,600); L.lbc = [0.02;0]; L.maxnorm = 0.9;
L.op = @(t,y) diff(y,2) - 2*(-2+t/150)*y;
y = L\0; plot(y,CO,ivp)
title(['Fig.~17.4.~~Solution of (17.3) ' ...
    'with slowly varying parameter $c$'],FS,11)
xlabel('$t$',FS,10,IN,LT), ylabel('$y$',FS,10,IN,LT)
axis([0 600 -1 1])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% For $t<300$ we see quasi-steady oscillation around
% $y=0$, with frequency decreasing to $0$ as $t$ approaches
% $300$.  After this, the character of the problem changes
% and an exponential explosion begins. There is no steady
% state; the orbit is diverging to $\infty$.
% </latex>

%%
% By varying features of the marble problem, we can begin to
% explore the rich world of bifurcation theory.  
% A good way to start is to introduce a quartic term in (17.1),
% $$ h(y) = -c\kern .3pt y^2 + y^4.  \eqno (17.4) $$
% Now the surfaces look like this.
y = chebfun('y',[-1.4 1.4]);
for k = 1:3
  c = -2+k; surface = -c*y^2 + y^4;
  subplot(1,3,k), plot(surface,'k',LW,1.2), hold on
  fill(real(marble),imag(marble),[0 0 0]), hold off
  axis([-1.3 1.3 -1.2 1.4]), axis square, axis off
  text(0,1.7,[sprintf('$c = %2.0f$,',c) s{k} 'stable'],FS,9,HA,CT)
end
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The slope is $h'(y) = -2\kern .3pt c\kern .3pt y + 4y^3$, and
% (17.2) becomes
% $$ y'' = 2\kern .3pt c\kern .3pt y - 4y^3. \eqno (17.5) $$
% Here are some trajectories in the $(y,y')$ phase plane.
% </latex>
N = chebop(0,2.5);
for j = 1:3
  subplot(1,3,j), plot(0,0,'.k',MS,6), hold on
  c = -2+j; N.op = @(t,y) diff(y,2) - 2*c*y + 4*y^3;
  if c>0, plot(sqrt(c/2)*[-1 1],[0 0],'.k',MS,6), end
  for k = 1:12
    N.lbc = [u0(k); v0(k)];
    y = N\0; arrowplot(y,diff(y),LW,.5,CO,ivpnl,MS,3)
  end
  axis([-2 2 -2 2]), axis square, hold off, set(gca,YT,-2:2:2)
  title([sprintf('$c = %3.1f$,',c) s{j} 'stable'],FS,11)
end
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The bifurcation diagram for this problem is drawn below.
% This is called a {\bf pitchfork bifurcation,} or more fully a
% {\bf supercritical pitchfork bifurcation}.
% For $c>0$, we se that the marble has stable rest positions
% at $x = \pm \sqrt{c/2}$.
% </latex>
clf, plot([-2 0],[0 0],'-k',[2 0],[0 0],'--k')
ystar = chebfun('y'); c = 2*ystar^2; 
hold on, plot(c,ystar,'k'), hold off
title(['Fig.~17.7.~~Bifurcation diagram for (17.5): ' ...
       'supercritical pitchfork'],FS,11)
xlabel('$c$',FS,10,IN,LT), ylabel('$y_*^{}$',FS,10,IN,LT)
text(1.08,0.62,'stable',FS,10)
text(1.08,0.12,'unstable',FS,10)
%%
% <latex> \vskip 1.02em </latex>

%%
% Let us again consider a variable coefficient problem with a
% slowly changing coefficient $c(t)$, 
% $$ y'' = 2\kern .3pt c(t)y - 4y^3, \quad c(t) = -2+t/150. \eqno (17.6) $$
% This time we see a continuous transition from oscillation about
% $0$ to oscillation about a nonzero value, beautifully
% matching one branch of the bifurcation diagram.
N = chebop(0,600); N.lbc = [0.02;0];
N.op = @(t,y) diff(y,2) - 2*(-2+t/150)*y + 4*y^3;
y = N\0; plot(y,CO,ivpnl,LW,.5)
title(['Fig.~17.8.~~Solution of (17.6) ' ...
    'with slowly varying parameter $c$'],FS,11)
xlabel('$t$',FS,10,IN,LT), ylabel('$y$',FS,10,IN,LT)
axis([0 600 -1 1])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The solution settles on the upper rather than the lower branch
% for no particular reason, a phenomenon
% known as {\em symmetry breaking} (discussed
% further in Chapters~18 and~22).  If we alter the initial conditions, it might as
% easily find the lower branch.  Here is an example of
% that behavior, where the only difference is that the initial amplitude
% $y(0)$ has been increased from $0.02$ to $0.05$.
% </latex>
N.lbc = [0.05;0]; y = N\0; plot(y,CO,ivpnl,LW,.5)
title(['Fig.~17.9.~~Same, with adjustment ' ...
       'of initial condition'],FS,11)
xlabel('$t$',FS,10,IN,LT), ylabel('$y$',FS,10,IN,LT)
axis([0 600 -1 1])
%%
% <latex> \vskip 1.02em </latex>

%%
% In the curves above, the departure from the steady state
% has been seeded by a non-zero initial condition.
% Another approach is to start from a zero condition but
% introduce perturbations along the way.  Here we carry
% out such an experiment, with the perturbation consisting of
% $0.003$ times a smooth random function $f(t)$ of the kind described in
% Chapter 12,
% $$ y'' = 2\kern .3pt c(t)y - 4y^3 + 0.003f(t),
% \quad c(t) = -2+t/150. \eqno (17.7) $$
% The details differ, but the overall
% behavior is as before, with the trajectory again happening to find
% the lower branch of the bifurcation diagram.
N.lbc = [0;0]; rng(2); rhs = .003*randnfun(1,[0 600],'norm');
y = N\rhs; plot(y,CO,ivpnl,LW,.5)
title(['Fig.~17.10.~~Solution of (17.7), ' ...
       'now driven by random perturbations'],FS,11)
xlabel('$t$',FS,10,IN,LT), ylabel('$y$',FS,10,IN,LT)
axis([0 600 -1 1])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% A conspicuous feature of the last three plots is persistent
% oscillations, an effect we saw in Chapter 8 in
% the discussion of resonance.  To construct problems that are not
% so influenced by past history, we can introduce a
% damping term in the ODE.\ \ Specifically, 
% from now on we shall add the term $-0.2\kern .3pt y'$, so that
% (17.5), for example, becomes
% $$ y'' = 2\kern .3pt c\kern .3pt y - 4y^3 - 0.2\kern .3pt y'. \eqno (17.8) $$
% The trajectories in the phase plane change their shapes.
% Now they always spiral into fixed
% points, just as the marble in Fig.~17.5 will eventually
% come to rest if there is friction.
% </latex>
N = chebop(0,5);
th = (pi/2)*(1:4)+.000001; u0 = cos(th); v0 = sin(th);
for j = 1:3
  subplot(1,3,j), plot(0,0,'.k',MS,6), hold on
  c = -2+j; N.op = @(t,y) diff(y,2) - 2*c*y + 4*y^3 + .2*diff(y);
  if c>0, plot(sqrt(c/2)*[-1 1],[0 0],'.k',MS,6), end
  for k = 1:4
    N.lbc = [u0(k); v0(k)];
    y = N\0; arrowplot(y,diff(y),LW,.5,CO,ivpnl,MS,3,YS,.4)
  end
  axis([-2 2 -2 2]), axis square, hold off
  set(gca,YT,-2:2:2)
  title([sprintf('$c = %3.1f$,',c) s{j} 'stable'],FS,11)
end
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Here is a trajectory for the damped time-dependent problem
% $$ y'' = 2\kern .3pt c(t)y - 4y^3 - 0.2\kern .3pt y' + 0.003f(t),
% \quad c(t) = -2+t/150, \eqno (17.9) $$
% with the same random forcing function $0.003f(t)$ as before.
% </latex>
N = chebop(0,600); N.lbc = [0;0];
N.op = @(t,y) diff(y,2) - 2*(-2+t/150)*y + 4*y^3 + .2*diff(y);
y = N\rhs; clf, plot(y,CO,ivpnl,LW,.5)
title('Fig.~17.12.~~Solution of (17.9), with damping',FS,11)
xlabel('$t$',FS,10,IN,LT), ylabel('$y$',FS,10,IN,LT)
axis([0 600 -1 1])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \vskip 1em
% \begin{center}
% *~~~*~~~*
% \end{center}
% \vskip -1em
% </latex>

%%
% <latex>
% Let us review what we have done so far in this
% chapter.  We have looked at two dynamical
% systems: first the quadratic problem (17.1)--(17.2),
% then the quartic problem (17.4)--(17.5).
% For each one we displayed four figures showing the physical
% interpretation by a marble
% on a surface, the phase plane, a bifurcation diagram,
% and a trajectory for a problem with slowly varying $c(t)$.\ \ In
% the second case we actually showed four such trajectories:
% two for different initial conditions, then one driven by
% a small random forcing term $0.003f(t)$, then a fourth
% of the same kind but with damping included in the equation.
% </latex>

%%
% <latex>
% We are now going to follow the same pattern for
% one final kind of bifurcation: a {\bf subcritical
% pitchfork bifurcation.}  The simplest starting point of
% this discussion could be equations (17.4) and (17.8) again, but
% with a sign change on the quartic term:
% $$ h(y) = -c\kern .3pt y^2 - y^4  \eqno (17.10) $$
% and
% $$ y'' = 2\kern .3pt c\kern .3pt y + 4y^3 - 0.2\kern .3pt y'. \eqno (17.11) $$
% This gives a new marble diagram,
% </latex>
y = chebfun('y',[-1.4 1.4]);
for k = 1:3
  c = -2+k; surface = -c*y^2 - y^4;
  subplot(1,3,k), plot(surface,'k',LW,1.2), hold on
  fill(real(marble),imag(marble),[0 0 0]), hold off
  axis([-1.3 1.3 -1.4 1.2]), axis square, axis off
  text(0,0.9,[sprintf('$c = %2.0f$,',c) s{k} 'stable'],FS,9,HA,CT)
end
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% and a new bifurcation diagram,
% </latex>
clf, plot([-2 0],[0 0],'-k',[2 0],[0 0],'--k')
ystar = chebfun('y'); c = -2*ystar^2; 
hold on, plot(c,ystar,'--k'), hold off
title(['Fig.~17.14.~~Bifurcation diagram for (17.11): ' ...
       'subcritical pitchfork'],FS,11)
xlabel('$c$',FS,10,IN,LT), ylabel('$y_*^{}$',FS,10,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% However, let us immediately modify the problem
% further to make it more realistic, and more interesting, by
% adding a stabilizing $y^6$ term:
% $$ h(y) = -c\kern .3pt y^2 - y^4 + y^6  \eqno (17.12) $$
% and
% $$ y'' = 2\kern .3pt c\kern .3pt y + 4y^3 - 6y^5 - 0.2\kern .3pt y'.
% \eqno (17.13) $$
% Now the marble diagram looks like this.
% </latex>
y = chebfun('y',[-1.4 1.4]);
for k = 1:3
  c = -.4+.1*k; surface = -c*y^2 - y^4 + y^6;
  subplot(1,3,k), plot(surface,'k',LW,1.2), hold on
  fill(real(marble),.5*imag(marble),[0 0 0]), hold off
  axis([-1.3 1.3 -.7 .6]), axis square, axis off
  text(0,0.7,[sprintf('$c = %3.1f$,',c) s{k} 'stable'],FS,9,HA,CT)
end
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The bifurcation diagram becomes more complicated.
% </latex>
clf, plot([-2 0],[0 0],'-k',[2 0],[0 0],'--k')
ystar = chebfun('y',[-1 1]/sqrt(3)); c = -2*ystar^2 + 3*ystar^4; 
hold on, plot(c,ystar,'--k'), axis([-2 2 -1.2 1.2])
ystar = chebfun('y',[1/sqrt(3) 1.1]); c = -2*ystar^2 + 3*ystar^4; 
plot(c,ystar,'k'), plot(c,-ystar,'k'), hold off
title(['Fig.~17.16.~~Bifurcation diagram for (17.13): ' ...
       'subcritical pitchfork'],FS,11)
xlabel('$c$',FS,10,IN,LT), ylabel('$y_*^{}$',FS,10,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% New dynamical possibilities are implied by this figure.
% One implication is that if $c$ increases slowly through $c=0$, one may
% expect a sudden {\em jump transition} to an amplitude 
% $y \approx \pm \sqrt{2/3}$.  A simulation confirms this prediction.
% </latex>
N.lbc = [0;0];
N.op = @(t,y) diff(y,2)-2*(-2+t/150)*y-4*y^3+6*y^5+.2*diff(y);
y = N\rhs; clf, plot(y,CO,ivpnl,LW,.5)
title('Fig.~17.17.~~Jump transition',FS,11)
xlabel('$t$',FS,10,IN,LT), ylabel('$y$',FS,10,IN,LT)
axis([0 600 -1.2 1.2])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent 
% Another implication is that if
% we now let $c$ {\em decrease} through $0$, we may observe
% the effect known as {\bf hysteresis}, in which the
% jump back to the initial state
% occurs at a parameter value different from the initial one.
% </latex>
y0 = y(600); yp0 = deriv(y,600); N.lbc = [y0;yp0];
N.op = @(t,y) diff(y,2)-2*(2-t/150)*y-4*y^3+6*y^5+.2*diff(y);
y = N\rhs; hold on, treverse = chebfun('600-t',[0 600]);
plot(treverse,y,'r',LW,.5), hold off
title('Fig.~17.18.~~Hysteresis',FS,11)
xlabel('$t$',FS,10,IN,LT), ylabel('$y$',FS,10,IN,LT)
axis([0 600 -1.2 1.2])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: FitzHugh--Nagumo equations of
% neural signals}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% One of the great achievements of
% twentieth century science was the model of propagation of
% neural signals by Alan Hodgkin and Andrew Huxley in 1952,
% which won them the 1963 Nobel Prize for Physiology or Medicine.
% In the words of J. D. Murray in his book
% {\em Mathematical Biology I\/} (Springer, 2002), from which this
% discussion is adapted,
% ``The theory of neuron firing and
% propagation of nerve action potentials is one of the major successes of
% real mathematical biology.''
% </latex>

%%
% <latex>
% Hodgkin and Huxley worked 
% on the neuron of the giant squid, expressing the current out of the
% axon of the neuron in terms of
% the oscillations of an activator--inhibitor system.
% The full description requires a PDE, but the model is interesting even
% when reduced to an ODE by ignoring variation along the
% length of the axon.  FitzHugh and Nagumo
% subsequently simplified the equations even further to a two-variable ODE
% system now known as the {\em FitzHugh--Nagumo equations,}
% $$ 
% v' = v - {1\over 3}v^3 - w + I, \quad w' = 0.08(v-0.8w+0.7).
% \eqno (17.14)
% $$
% Here $v(t)$ represents the electric potential across the membrane of the
% neuron, which is experimentally accessible. The $w(t)$ component is an
% approximation of the current through the membrane due to movement of
% ions. Finally, $I$ is a current applied to stimulate the cell
% experimentally.   This will be our bifurcation parameter.
% </latex>

%%
% <latex>
% In view of the plus and minus signs in (17.14), we can say that
% $v$ is self-activating but inhibited by $w$. Conversely, $w$
% is activated by $v$ but otherwise decays (down to the level
% $w=.7/.8$). This a bit like the
% interaction of rabbits and foxes in the Lotka--Volterra model, but with
% different details. While~$v$ and~$w$ tend to grow and decrease in
% opposition to one another, they are also given some inhomogeneous growth.
% If~$I$, the external forcing for $v$, is too small, then~$v$ does not
% recover fast enough to excite the system. However, when $I$ passes a
% threshold~$I_1$, the neuron begins to fire, with $v$ reaching peaks and
% troughs repeatedly. If~$I$ continues to increase, it passes another
% threshold $I_2$ past which $w$ is pushed into a higher steady state that
% again shuts down the activation of $v$.
% </latex>

%%
% We can see all three of these behaviors
% by allowing $I$ to grow slowly in time, with a small random
% perturbation included to make the system noisy.
% We plot just the component $w$; the behavior of $v$ is similar.
N = chebop([0 1000]); t = chebfun('t',[0 1000]);
I = t/500 + .01*randnfun(2,[0 1000],'norm');
N.op = @(t,v,w) [ diff(v)-(v-v^3/3-w)
                  diff(w)-0.08*(v-0.8*w+0.7) ];
N.lbc = @(v,w) [v;w];
[v,w] = N\[I; 0*I]; plot(w,CO,ivpnl), ylim([-4 4])
title(['Fig.~17.19.~~Firing of a neuron in ' ...
       '(17.14) as $I$ is increased'],FS,11,IN,LT)
xlabel('$t$',FS,10,IN,LT), ylabel('$w$',FS,10,IN,LT)
text(103,-.96,'neuron switches on',FS,9)
text(592,-.3,'neuron switches off',FS,9)
%%
% <latex> \vskip 1.02em </latex>
%%
% <latex>
% \noindent
% This image reveals that there must be a pair of bifurcations as
% the parameter~$I$ is increased.  At $t\approx 180$, with
% $I\approx 0.36$, the oscillation switches on, and at $t\approx 700$,
% with $I\approx 1.4$, it switches off again.  Bifurcations like
% this where an oscillation turns on or off are known as
% {\bf Hopf bifurcations,} and they are characterized by a
% pair of eigenvalues of
% the Jacobian matrix crossing the imaginary axis to move into or
% out of the right half of the complex plane.\footnote{There are
% many examples of Hopf bifurcation in our lives, such as that
% annoying vibration that sets in
% when your car passes a particular speed.}
% </latex>

%%
% We can confirm these guesses as follows.  A fixed point
% $v_*^{}, w_*^{}$ of (17.14) is characterized by the conditions
% $$ v_*^{} - {1\over 3}v_*^3 - w_*^{} + I = 0,
% \quad v_*^{}-0.8w_*^{}+0.7 = 0, $$
% and eliminating $w_*^{}$ with the aid of the second equation gives
% $$ I = -v_*^{} + {v_*^3\over 3} + {v_*^{}+0.7\over 0.8}. \eqno (17.15) $$
% By inverting this equation we get a plot $v_*^{}$ as a function of $I$.
v = chebfun('v',[-2,2]); I = -v + v^3/3 + (v+.7)/.8;
vstar = inv(I); plot(vstar,'k')
title(['Fig.~17.20.~~Fixed point $v_*^{}$ as a ' ...
       'function of $I$'],FS,11,IN,LT), xlim([0 3])
xlabel('$I$',FS,10,IN,LT), ylabel('$v_*^{}$',FS,10,IN,LT)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% By differentiating the two equations of (17.14) with respect
% to $v_*^{}$ and $w_*^{}$ we find the corresponding Jacobian matrix,
% which depends only on $v_*^{}$ since $w$ enters the equations
% linearly:
% $$ {\bf J}_*^{} =
% \pmatrix{1 - v_*^2 & -1 \cr\noalign{\vskip 4pt} .08& -.064}. \eqno (17.16) $$
% We now plot the maximal
% real part of the eigenvalues of ${\bf J}_*^{}$,
% the {\em spectral abscissa} (introduced
% in Chapter~14), as a function of $I$.
% </latex>
J = @(v) [1-v^2 -1; .08 -.064];
abscissa = chebfun(@(I) max(real(eig(J(vstar(I))))),[0 3], ...
   'splitting','on'); 
plot(abscissa,'m'), hold on
title(['Fig.~17.21.~~Spectral abscissa of ${\bf J}_*^{}$ ' ...
       'with bifurcation points'],FS,11,IN,LT)
xlabel('$I$',FS,10,IN,LT), ylabel('spectral abscissa',FS,9,IN,LT)
Ibifurc = roots(abscissa);
plot(Ibifurc,0*Ibifurc,'.k',MS,14), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The two bifurcation points are shown as
% black dots.  They match our earlier estimates nicely.
% </latex>
Ibifurc

%%
% <latex>
% You may wonder why the plot of the spectral abscissa is
% not smooth, but breaks into four regions separated by
% discontinuities of the derivative.
% Actually, this effect appeared already in Fig.~14.7.
% The discontinuities correspond to transitions between
% regimes where the dominant eigenvalues of ${\bf J}_*^{}$ 
% are a complex
% conjugate pair (the second and the fourth parts of the
% curve) and regimes where there is a single dominant real eigenvalue
% (the first and third).  As predicted, when the spectral
% abscissa passes through zero for this problem,
% it is with a pair of complex
% conjugate eigenvalues, signalling a Hopf bifurcation:
% </latex>
eig(J(vstar(Ibifurc(1))))

%%
% <latex>
% \noindent
% The magnitude of the imaginary part of these eigenvalues, $0.2755$,
% gives the frequency of the linearized oscillation.  This corresponds to a 
% period of $2\pi/0.2755 \approx 22.8$, slightly less than
% the period observed in Fig.~17.19.
% </latex>

%%
% <latex>
% Here are plots of trajectories in the $v$-$w$ plane
% for $I=0.31$, below the bifurcation, and $I=0.34$, above.
% </latex>
dom = [0 75]; N = chebop(dom);
N.op = @(t,v,w) [ diff(v)-(v-v^3/3-w)
                  diff(w)-0.08*(v-0.8*w+0.7) ];
I = chebfun('0.31+0*t',dom); vs = vstar(0.31); ws = (vs+.7)/.8;
subplot(1,2,1), plot(vs,ws,'.k',MS,6), hold on
N.lbc = @(v,w) [v-vs;w-ws-.2];
[v,w] = N\[I; 0*I]; arrowplot(v,w,LW,.5,CO,ivpnl,MS,4,YS,1.5), hold off
axis([-1.5 -0.7 -.45 -.05]), axis square
text(-1.1,-.11,'$I=0.31$: sink',FS,9,IN,LT,HA,CT)
set(gca,XT,-1.4:.2:-.8,YT,-.4:.2:-.2)
xlabel('$v$',FS,10), ylabel('$w$',FS,10)
I = chebfun('0.34+0*t',dom); vs = vstar(0.34); ws = (vs+.7)/.8;
pp = get(gca,'pos'); pp(1) = pp(1)+.025; set(gca,'pos',pp)
subplot(1,2,2), plot(vs,ws,'.k',MS,6), hold on
N.lbc = @(v,w) [v-vs;w-ws-.2];
[v,w] = N\[I; 0*I]; arrowplot(v,w,LW,.5,CO,ivpnl,MS,4,YS,3), hold off
axis([-2.8 2.8 -.8 2]), axis square
text(0,1.6,'$I=0.34$: limit cycle',FS,9,IN,LT,HA,CT)
set(gca,XT,-2:2,YT,0:2), xlabel('$v$',FS,10), ylabel('$w$',FS,10)
pp = get(gca,'pos'); pp(1) = pp(1)-.025; set(gca,'pos',pp)
%%
% <latex> \vskip 1.02em </latex>

%%
% The equations above are simple enough that although
% our calculation has been numerical, it is not hard to
% do it exactly.  The trace of ${\bf J}_*^{}$, the sum of
% its diagonal entries, is $0.936-v_*^2$, and this will be
% equal to the sum of the eigenvalues.  When the eigenvalues
% form a conjugate pair, both have the same real part, and
% thus the Hopf bifurcation will occur when the trace
% passes through zero, that is, $v_*^{} = \pm \sqrt{0.936}\kern 1pt$.
% Thus by (17.15), the two bifurcation points correspond to
% the parameter values
% $$ I_c = \pm a + {\pm a^3\over 3} + {\pm a+0.7\over 0.8},
% \quad a = \sqrt{0.936}\kern 1pt. $$

%%
% <latex>
% {\sc History.} 
% The observation that a system may change stability
% when a parameter passes through a critical value
% was well known to 19th-century applied mathematicians
% such as Helmholtz, Kelvin, and Rayleigh.
% The technical use of the term ``bifurcation'' may be
% due to Poincar\'e in 1885.  Further
% foundations of the mathematical theory of bifurcations
% were laid by Andronov and Pontryagin in the Soviet Union in
% the 1930s.  
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.}  For a marvelously rich discussion
% of all kinds of bifurcation effects see Seydel, 
% {\em Practical Bifurcation and Stability Analysis,} 3rd ed.,
% Springer, 2009.  The original edition appeared in 1994.
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 17.}  Bifurcation refers to the
% change of the qualitative nature of a solution to a nonlinear problem
% as a parameter varies through a critical value.  Typically what
% changes at the critical value is the stability structure of
% a fixed point of the system. For example, if an eigenvalue crosses
% into the right half-plane, then a stable solution may become
% unstable, so that trajectories jump instead to other solutions.
% If it is a complex conjugate pair of eigenvalues that crosses into
% the right half-plane, we have a Hopf bifurcation, leading to the
% onset of oscillations.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% {\em \underline{Exercise $17.1$}.  Van der Pol equation.}
% Consider the fixed point $y = y' = 0$ of the
% van der Pol equation $y''+y-\mu(1-y^2)y'=0$, as
% in (9.6).  {\em (a)} Determine the Jacobian matrix and its eigenvalues
% analytically, and show that there is a Hopf bifurcation as
% $\mu$ passes from negative to positive.
% {\em (\kern .7pt b)} A different bifurcation takes place
% as $\mu$ passes through $-2$ or $2$.  Explain. 
% \par
% {\em \underline{Exercise $17.2$}.  FitzHugh--Nagumo experiment and noise.}
% Rerun the experiment of Fig.~17.19 without
% the random perturbation term.  What is the change
% in the output?  Can you explain why?
% \par
% {\em \underline{Exercise $17.3$}.  Slowly-varying logistic map.}
% This book deals with continuous-time processes, also known as {\em flows,}
% defined by equations like $y' = f(y)$.  Every topic we consider has
% an analogue for discrete-time processes, also known as {\em maps,}
% defined by equations like $y_{n+1}^{} = f(y_n^{})$.
% Execute the code
% \begin{verbatim}
%     r = linspace(2,4,npts); y = 0*r; y(1) = .5;
%     for j = 1:npts-1, y(j+1) = r(j)*y(j)*(1-y(j)); end
%     plot(r,y,'.','markersize',3)
% \end{verbatim}
% for \verb|npts =| $10^2, 10^3, \dots , 10^7$ and show the plots that result.
% Write down five careful
% sentences describing what you think these plots
% reveal in the light of the discussion of this chapter.
% Make specific reference to relevant figure numbers.
% \par
% {\em Exercise $17.4$.  Bifurcation diagram for the nonlinear pendulum.}
% Consider the nonlinear pendulum equation (9.6), $y'' = -\sin (y)$, on the
% interval $t\in [\kern .3pt 0,d\kern .5pt]$
% with boundary conditions $y(0) = y(d\kern .5pt ) = 0$.  
% {\em (a)} Making reference to the
% phase plane diagram, explain why for $d\le \pi$ this problem has a unique
% solution and describe that solution.
% {\em (\kern .7pt b)}~Explain why there is a bifurcation at $d=\pi$ and describe
% the two new solutions that appear for $d>\pi$.  
% {\em (c)} Explain why there is another bifurcation at $d=2\pi$ and describe
% the two new solutions that appear there.  Meanwhile what has
% happened to the three solutions already present?
% {\em (d)} Exactly how many
% solutions will there be for $d=100$, and what will they look like?
% \par
% {\em Exercise $17.5$.  Subcritical pitchfork marble.}
% Consider the ODE (17.13) for a marble on a surface with $c=-0.2$ sketched in
% Fig.~17.15, whose bifurcation diagram was shown in Fig.~17.16.
% Consider the BVP for this problem $y(0) = y(T\kern .5pt )$ = 0, where
% $T>0$ is a fixed constant.  Assume in what follows that
% we are only interested in solutions with $y(t)\ge 0$ for all
% $t\in [\kern .3pt 0,T\kern .5pt ]$.
% {\em (a)} Describe physically, in terms of the marble rolling
% along the surface, three structurally different types of solutions
% that may exist for this problem (for appropriate values of $T\kern .5pt $).
% {\em (\kern .7pt b)} Draw a sketch of the phase plane and interpret your
% three solutions in this context.
% \par
% {\em \underline{Exercise $17.6$}.  Emergence of a limit cycle.}
% The system $u' = -v + u(\mu-u^2-v^2)$,
% $v' = u + v(\mu-u^2-v^2)$, where $\mu$ is a real parameter, has
% a fixed point $(u,v) = (0,0)$.
% {\em (a)} Compute solutions for $t\in [\kern .3pt 0,20\kern .3pt]$
% starting from $(u,v) = (2,0)$ for $\mu = -1, -0.5, \dots, 1$.  Plot
% each curve in the $(u,v)$ plane with axes equal.  What dependence on
% $\mu$ do you observe?
% {\em (\kern .7pt b)} Analyze the stability of the fixed point and discuss the
% bifurcation situation.
% </latex>
