%% 17. Bifurcation %% % % \setcounter{page}{212} % \parindent=.6cm % \parskip=-1em % %% % As parameters vary, the solutions of an % ODE may change in structural properties. This is the % subject of bifurcation theory. %% % % 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. % 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 %% % \vskip 1.02em %% % % 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. % 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 %% % \vskip 1.02em %% % % \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 .$$ % %% % % 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. % %% % % 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. % 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) %% % \vskip 1.02em %% % % \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. % %% % % 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$. % 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]) %% % \vskip 1.02em %% % % \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$. % %% % 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 %% % \vskip 1.02em %% % % \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. % 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 %% % \vskip 1.02em %% % % \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}$. % 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) %% % \vskip 1.02em %% % 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]) %% % \vskip 1.02em %% % % \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$. % 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]) %% % \vskip 1.02em %% % 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]) %% % \vskip 1.02em %% % % 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. % 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 %% % \vskip 1.02em %% % % \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. % 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]) %% % \vskip 1.02em %% % % \vskip 1em % \begin{center} % *~~~*~~~* % \end{center} % \vskip -1em % %% % % 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. % %% % % 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, % 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 %% % \vskip 1.02em %% % % \noindent % and a new bifurcation diagram, % 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) %% % \vskip 1.02em %% % % \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. % 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 %% % \vskip 1.02em %% % % \noindent % The bifurcation diagram becomes more complicated. % 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) %% % \vskip 1.02em %% % % \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. % 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]) %% % \vskip 1.02em %% % % \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. % 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]) %% % \vskip 1.02em %% % % \begin{center} % \hrulefill\\[1pt] % {\sc Application: FitzHugh--Nagumo equations of % neural signals}\\[-3pt] % \hrulefill % \end{center} % %% % % 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.'' % %% % % 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. % %% % % 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$. % %% % 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) %% % \vskip 1.02em %% % % \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.} % %% % 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) %% % \vskip 1.02em %% % % \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$. % 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 %% % \vskip 1.02em %% % % \noindent % The two bifurcation points are shown as % black dots. They match our earlier estimates nicely. % Ibifurc %% % % 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: % eig(J(vstar(Ibifurc(1)))) %% % % \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. % %% % % Here are plots of trajectories in the $v$-$w$ plane % for $I=0.31$, below the bifurcation, and $I=0.34$, above. % 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) %% % \vskip 1.02em %% % 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.$$ %% % % {\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. % %% % % \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. % %% % % \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} % %% % \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. %