%% 19. Periodic ODEs

%%
% <latex>
% \setcounter{page}{241}
% </latex>

%%
% <latex>
% Suppose we have an ODE defined for
% $t\in (-\infty,\infty)$ whose coefficients are $2\pi$-periodic,
% or in short, a {\bf periodic ODE.}\ \ Will the solutions
% be $2\pi$-periodic too?
% </latex>

%%
% <latex>
% Almost any experiment will show you that
% they don't {\em have} to be.  For example, 
% consider the problem
% $$ y' + (1+\cos(t)) y = 1, \quad y(0) = 0. \eqno (19.1) $$
% By Theorem 11.1, a unique solution exists for all $t$, negative
% and positive.  This solution is
% not periodic, however, as we see by plotting it
% for $t\in [\kern .3pt 0,10\pi]$.
% </latex>
ODEformats
L = chebop(0,10*pi); L.op = @(t,y) diff(y) + (1+cos(t))*y;
L.lbc = 0; y = L\1; plot(y,CO,ivp), ylim([0 3])
title('Fig.~19.1.~~Nonperiodic solution to (19.1)',FS,11)
pilabels = {'0','2\pi','4\pi','6\pi','8\pi','10\pi'};
set(gca,FS,8,XT,pi*(0:2:10),XTL,pilabels), set(gca,YT,0:3)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Although this solution isn't periodic, it is very close to
% periodic after an initial phase.  It would seem
% that there exists a truly periodic solution to this ODE that
% will oscillate between extremes of about $0.5$ and $2.4$.
% (From now on in this chapter, for simplicity, {\bf periodic} usually means
% $2\pi$-periodic.)
% Chebfun will find this solution if you
% specify the boundary condition \verb|'periodic'|.
% </latex>
L = chebop(0,2*pi); L.op = @(t,y) diff(y) + (1+cos(t))*y;
L.bc = 'periodic'; y = L\1; plot(y,CO,bvp), ylim([0 3])
title('Fig.~19.2.~~Periodic solution to (19.1)',FS,11)
pilabels2 = {'0','\pi/2','\pi','3\pi/2','2\pi'};
set(gca,FS,8,XT,pi*(0:.5:2),XTL,pilabels2), set(gca,YT,0:3)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The minimum and maximum are as expected,
% </latex>
ymin = min(y), ymax = max(y)

%%
% <latex>
% \noindent
% and the value at $t=0$ is a little bit above the minimum,
% </latex>
y0 = y(0)

%%
% <latex>
% \noindent
% If we solve the IVP on $[\kern .3pt 0,10\pi]$ with this value
% specified as the initial condition, the periodic behavior is seen.
% </latex>
L = chebop(0,10*pi); L.op = @(t,y) diff(y) + (1+cos(t))*y;
L.lbc = y0; y10pi = L\1; plot(y10pi,CO,ivp), ylim([0 3])
title('Fig.~19.3.~~More periods of the same solution',FS,11)
set(gca,FS,8,XT,pi*(0:2:10),XTL,pilabels), set(gca,YT,0:3)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Encouraged by this experiment, let us try a second-order periodic equation,
% $$ y'' + (1+\cos(t)) y = 1, \quad y(0) = y'(0) = 0. \eqno (19.2) $$
% Again the first solution we find is nonperiodic.
% </latex>
L = chebop(0,10*pi); L.op = @(t,y) diff(y,2) + (1+cos(t))*y;
L.lbc = [4;0]; y = L\1; plot(y,CO,ivp), ylim([-50 50])
title('Fig.~19.4.~~Nonperiodic solution to (19.2)',FS,11)
set(gca,FS,8,XT,pi*(0:2:10),XTL,pilabels)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% This time, the curve is not settling down to a periodic form, but
% oscillating and growing exponentially.
% Nevertheless, again there exists a periodic solution.
% </latex>

L = chebop(0,2*pi); L.op = @(t,y) diff(y,2) + (1+cos(t))*y;
L.bc = 'periodic'; y = L\1; plot(y,CO,bvp), ylim([-4 4])
title('Fig.~19.5.~~Periodic solution to (19.2)',FS,11)
set(gca,FS,8,XT,pi*(0:.5:2),XTL,pilabels2)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% As before, we can recover periodic
% oscillations over $[\kern .3pt 0,10\pi]$ if we set the
% initial conditions appropriately.
% </latex>

%%
% <latex>
% Nonlinear periodic ODE\kern .5pt s often have periodic solutions too.
% For example, the equation
% $$ y'' + (1+\cos(t)) (y + (y/4)^3) = 1 \eqno (19.3) $$
% is a nonlinear variant of (19.2).  Here we compute a periodic
% solution and immediately plot it on $[\kern .3pt 0,10\pi]$.
% </latex>

L = chebop(0,2*pi);
L.op = @(t,y) diff(y,2) + (1+cos(t))*(y+(y/4)^3);
L.bc = 'periodic'; y = L\1;
y10pi = chebfun(@(t) y(t),[0 10*pi],'trig'); plot(y10pi,CO,bvpnl)
title('Fig.~19.6.~~Periodic solution to (19.3)',FS,11)
set(gca,FS,8,XT,pi*(0:2:10),XTL,pilabels), ylim([-4 4])
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Of course, as usual, the behavior of nonlinear problems can vary
% greatly.  If $(y/4)^3$ is changed to $(y/3)^3$ in (19.3), experiments
% indicate that there is no solution.
% For another example, in Chapter 1 we considered the van der Pol
% equation (1.2),
% $$ 0.3y''-(1-y^2)y'+y = 0. $$
% Like all autonomous equations, this ODE is not just $2\pi$-periodic, 
% but $T$-periodic for any period $T$.  For any $T$, it has the $T$-periodic
% solution $y = 0$.  For $T=2\pi$, and for most other values of $T$,
% this is the only $T$-periodic solution.
% For the special value $T\approx 4.0725$, however,
% there is another nontrivial $T$-periodic solution
% corresponding to the limit cycle plotted in Chapter~1.
% In fact, this solution is one of an infinite family of
% distinct $T$-periodic solutions, since it could be shifted in~$t$ by
% any constant $\Delta t$.
% </latex>

%%
% <latex>
% Let us take stock of the situation.  We have seen that
% a periodic ODE may or may not have periodic solutions, and
% that it will almost certainly have nonperiodic ones.
% The remainder of this chapter is organized around three questions for
% periodic ODE\kern .5pt s:
% \par\medskip
% 1. Does there exist a periodic solution?
% \par\smallskip
% 2. If so, is it unique?
% \par\smallskip
% 3. What can we say about nonperiodic solutions?
% \label{question3}
% \par\medskip
% \noindent
% These questions have rich histories going back to
% the nineteenth century, and
% as usual with ODE\kern .5pt s, there are applications involving
% both time and space.  For example, 
% periodicity in time is associated
% with the dynamics of the solar
% system or of rotating or oscillating
% machinery, and periodicity in space is associated 
% with propagation of sound and light waves in crystals.
% We will focus mainly on linear problems, with brief comments
% on nonlinear ones at the end.  With regret, we shall mainly
% look just at questions (1) and (2), though the Application is
% of type (3).
% </latex>

%%
% <latex>
% Before turning to second-order periodic equations, the traditional
% focus in this subject area, let us consider the first-order case.
% Following the notation of Theorem 2.3,
% we start with the scalar linear problem
% $$ y' - a(t) y = g(t) \eqno (19.4) $$
% where $a$ and $g$ are assumed to be
% continuous and $2\pi$-periodic.  We
% can answer questions (1)--(2) above without much difficulty.
% First we note that by Theorem~10.1, there is a unique solution 
% $y(t)$ of (19.4) associated with any initial
% condition $y(0) = y_0$.  If $y(2\pi)\ne y(0)$, then by
% definition, $y(t)$ is not $2\pi$-periodic.  If $y(2\pi) = y(0)$, on
% the other hand, then since the coefficients are periodic, the
% solution on $[2\pi,4\pi]$ will be the same as on $[\kern .3pt 0,2\pi]$,
% the solution on $[4\pi,6\pi]$ will be the same as on $[2\pi,4\pi]$, and
% so on, implying that the solution on all
% of $(-\infty,\infty)$ is $2\pi$-periodic.  Thus the question
% of $2\pi$-periodicity reduces to the question of whether
% $y(2\pi) = y(0)$.
% </latex>

%% 
% One special case of (19.4) is the situation where 
% $a(t) = 0$,
% $$ y' = g(t). \eqno (19.5) $$
% Here the ODE reduces to an integral,
% $$ y(t) = y(0) + \int_0^t g(s)\kern .5pt ds , \eqno (19.6) $$
% implying that the solution for any initial condition
% $y(0) = y_0$ will be periodic if and only if
% $$ \int_0^{2\pi} g(s)\kern .5pt ds = 0 .  \eqno (19.7) $$
% Thus the answer to question (1) is that there is a periodic
% solution if and only if (19.7) holds, and the answer to (2)
% is that it is never unique since
% we could always add a constant.

%% 
% <latex>
% The more substantial special case of (19.4) is the homogeneous
% case $g(t) = 0$,
% $$ y' -a(t) y = 0.  \eqno (19.8) $$
% Here, following Theorem 2.2 (separation of
% variables), the solution is given by
% $$ y(t) = e^{h(t)} y(0), \eqno (19.9) $$
% where $h(t)$ is defined by
% $$ h(t) = \int_0^t a(s)\kern .5pt ds . \eqno (19.10) $$
% In particular, we have
% $$ y(2\pi ) = e^{h(2\pi )} y(0), \eqno (19.11) $$
% which implies that $y(2\pi) = y(0)$ if and only if
% $$ \left(1-e^{h(2\pi)}\right) y(0) = 0. $$
% This equation holds if and only if
% $y(0)$ is zero or $\exp(h(2\pi))=1$, that is,
% $$ h(2\pi) = 2\pi i \times \hbox{integer}. \eqno (19.12) $$
% If (19.12) holds, (19.8) is said to be {\bf critical}; 
% otherwise it is {\bf noncritical}.
% Thus we see that if
% (19.8) is noncritical, then (19.8) has the unique periodic solution
% $y(t) = 0\kern 1pt$: the
% answers to (1) and (2) are yes and yes.  If
% it is critical, then any value $y(0)$ in (19.9)
% gives a $2\pi$-periodic solution, so the answers are yes and no.
% </latex>

%%
% <latex>
% We shall now see that the behavior of the inhomogeneous problem (19.4)
% depends on the same condition (19.12) of criticality.  
% By Theorem 2.3, the solution can be written
% $$ y(t) = e^{h(t)} y(0) + \kern -2pt \int_0^t e^{h(t)-h(s)} g(s)\kern .5pt ds,
% \eqno (19.13) $$
% where $h(t)$ is still defined by (19.10).
% Thus the value at $t = 2\pi$ is 
% $$ y(2\pi) = e^{h(2\pi)} y(0) + \kern -2pt \int_0^{2\pi}
% e^{h(2\pi)-h(s)} g(s)\kern .5pt ds, \eqno (19.14) $$
% so the solution will be $2\pi$-periodic if and only if
% $$ \left(1 - e^{h(2\pi)}\right) y(0) =
% \int_0^{2\pi} e^{h(2\pi)-h(s)} g(s)\kern .5pt ds. \eqno (19.15)  $$
% This leads to the following theorem encompassing the
% two special cases just considered.
% </latex>

%%
% <latex>
% \smallskip
% {\em
% {\bf Theorem 19.1.  First-order scalar linear periodic ODE\kern .5pt s
% (\textsf{\textbf{FLaShi}}).}
% Consider the equation $$ y' - a(t) y = g(t), \eqno (19.16) $$
% where $a(t)$ and $g(t)$ are continuous and $2\pi$-periodic.
% If the equation is 
% noncritical in the sense that\/ $(19.12)$ does not hold,
% then there exists a unique $2\pi$-periodic solution~$y$.
% If it is critical, then there exist infinitely many
% $2\pi$-periodic solutions if the right-hand
% side of $(19.15)$ is zero, and no $2\pi$-periodic solutions if
% it is nonzero.
% }
% \smallskip
% </latex>

%%
% <latex>
% \def\qed{\vrule height4pt width3pt depth4pt}
% {\em Proof.}  If the equation is noncritical, we may divide
% (19.15) by $1-e^{h(2\pi)}$ to get
% $$ y(0) = \left.\int_0^{2\pi} e^{h(2\pi)-h(s)} g(s)\kern .5pt ds\right/
% \kern -2pt  \left( 1 - e^{h(2\pi)}\right) , $$
% producing a $2\pi$-periodic solution (19.13) as claimed.
% Uniqueness follows since this condition on $y(0)$ determines
% $y(t)$ uniquely.
% If the equation is critical, then the left-hand side of (19.15) is
% zero regardless of the value of $y(0)$.  If the right-hand side
% is zero, then each choice of $y(0)$ gives a distinct $2\pi$-periodic
% solution, whereas if the right-hand side is nonzero, no choice of
% $y(0)$ gives a periodic solution. \qed
% \smallskip
% </latex>

%%
% <latex>
% The last few pages have considered scalar problems.  For
% a system of equations, the pattern is similar.  Extending
% Chapter~14, consider the equation
% $$ {\bf y}' - {\bf A}(t) {\bf y} = {\bf g}(t), \eqno (19.17) $$
% where ${\bf y}(t)$ for each $t$ is an $n$-vector,
% ${\bf g}(t)$ is a given $2\pi$-periodic continuous
% $n$-vector function of $t$,
% and ${\bf A}(t)$ is a given $2\pi$-periodic continuous $n\times n$ matrix
% function of $t$.
% The general solution of (19.17) is
% $$ {\bf y}(t)={\bf Y}(t){\bf c}+{\bf Y}(t)\int_0^t{\bf Y}^{-1}(s)\kern .5pt {\bf g}(s)\kern .5pt ds
% \eqno (19.18) $$
% for $n$-vectors ${\bf c}$, where
% ${\bf Y}(t)$ is the {\bf fundamental matrix}, the $n\times n$
% matrix function whose columns are the linearly independent solutions
% of the homogeneous problem
% ${\bf y}' - {\bf A}(t) {\bf y} = {\bf 0}$,
% ${\bf y}^{(1)}(t),\dots,{\bf y}^{(n)}(t)$
% with ${\bf Y}(0) = {\bf I}$, the $n\times n$ identity.
% The values at $t=0$ and $t = 2\pi$ are
% $$ {\bf y}(0)={\bf c} $$
% and
% $$ {\bf y}(2\pi)={\bf Y}(2\pi){\bf c}+{\bf Y}(2\pi)\int_0^{2\pi}
% {\bf Y}^{-1}(s)\kern .5pt {\bf g}(s)\kern .5pt ds, $$
% so the solution will be $2\pi$-periodic if and only if
% $$ \left({\bf I}-{\bf Y}(2\pi)\right) {\bf c} =
% {\bf Y}(2\pi)\int_0^{2\pi} {\bf Y}^{-1}(s)\kern .5pt {\bf g}(s)\kern .5pt ds. 
% \eqno (19.19) $$
% The criticality condition will be the condition that
% ${\bf I}-{\bf Y}(2\pi)$ is singular, 
% or equivalently,
% $$ {\bf Y}(2\pi) \hbox{ has an eigenvalue equal to\/ } 1. \eqno (19.20) $$
% We say that (19.17) is {\bf critical} if this condition
% holds and {\bf noncritical} otherwise, i.e., all eigenvalues
% of ${\bf Y}(2\pi)$ are different from $1$.
% Theorem 19.1 generalizes as follows.
% </latex>

%%
% <latex>
% \smallskip
% {\em
% {\bf Theorem 19.2.  First-order linear periodic systems of ODE\kern .5pt s
% (\textsf{\textbf{FLashi}}).}
% Consider the $n$-dimensional system
% $$ {\bf y}' - {\bf A}(t) {\bf y} = {\bf g}(t), \eqno (19.21) $$
% where ${\bf A}(t)$ and ${\bf g}(t)$ are continuous and $2\pi$-periodic.
% If the equation is noncritical as defined by $(19.20)$,
% then there exists a unique $2\pi$-periodic solution~$\bf y$.
% If it is critical, then there exist infinitely many
% $2\pi$-periodic solutions if
% $\int_0^{2\pi} {\bf Y}^{-1}(s)\kern .5pt {\bf g}(s)\kern .5pt ds$ 
% is in the range of the matrix
% ${\bf I}-({\bf Y}(2\pi))^{-1} $, and no $2\pi$-periodic solutions if
% it is not in this range.
% }
% \smallskip
% </latex>

%%
% <latex>
% \def\qed{\vrule height4pt width3pt depth4pt}
% {\em Proof.}  If (19.21) is noncritical, we may 
% multiply (19.19) on the left by $({\bf I}-{\bf Y}(2\pi))^{-1}$ to get
% $$ {\bf c} = ({\bf I} - {\bf Y}(2\pi))^{-1}\kern 1pt {\bf Y}(2\pi) 
% \int_0^{2\pi} {\bf Y}^{-1}(s)\kern .5pt {\bf g}(s)\kern .5pt ds, $$ 
% that is,
% $$ {\bf c} = (({\bf Y}(2\pi))^{-1}-{\bf I})^{-1}
% \int_0^{2\pi} {\bf Y}^{-1}(s)\kern .5pt {\bf g}(s)\kern .5pt ds, $$ 
% and this vector corresponds to a $2\pi$-periodic solution
% (19.18) as claimed.  Uniqueness follows since
% a given $\bf c$ corresponds to one and only one solution
% ${\bf y}(t) = {\bf Y}(t) {\bf c}$.
% If the equation is critical, then (19.19) is a singular system
% of equations, having infinitely many solutions if the
% right-hand side is in the range of
% ${\bf I}-{\bf Y}(2\pi)$ and no solutions if it is not in that range. For the
% right-hand side of (19.19) to be in the range of ${\bf I}-{\bf Y}(2\pi)$ is the
% same as for $\int_0^{2\pi} {\bf Y}^{-1}(s)\kern .5pt {\bf g}(s)\kern .5pt ds$ 
% to be in the range of $({\bf Y}(2\pi))^{-1}
% - {\bf I}$. \qed
% \smallskip
% </latex>

%%
% <latex>
% Having written down these results for a first-order system
% of ODE\kern .5pt s, let us immediately specialize them to the
% case that has had the most attention.
% {\bf Hill's equation} is the scalar, homogeneous, second-order ODE
% $$ y'' + f(t)y = 0 , \eqno (19.22) $$
% where $f$ is real and $2\pi$-periodic.\footnote{In the literature
% the period is often taken to be $\pi$ instead of $2\pi$.}
% Although one could introduce a nonzero right-hand side,
% the equation is conventionally considered in this homogeneous form.
% Note that there is no first-order term, so there is no damping.
% {\bf Mathieu's equation} is the special case
% $$ y'' + (a + b\cos(t))y = 0 , \eqno (19.23) $$
% which is the simplest instance of Hill's equation in which
% the periodic coefficient is not simply a constant.  It corresponds to
% taking the first nonconstant term in a Fourier series of
% $f$ and is thus a natural first step in the study of (19.22).
% Again the equation is usually considered in this homogeneous form.
% </latex>

%%
% <latex>
% As always, (19.22) is equivalent to a first-order system of
% two variables via the identification
% ${\bf y} = (y_1^{},y_2^{})^T = (y,y')^T$, which transforms it to
% the equations
% $$ y_1' = y_2^{}, \quad y_2' + f(t) y_1^{} = 0. \eqno (19.24) $$
% For this system of ODE\kern .5pt s,
% ${\bf Y}(t)$ takes the form
% $$ {\bf Y}(t) = \pmatrix{u(t) & v(t) \cr\noalign{\vskip 3pt}
% u'(t) & v'(t)} , $$
% where $u$ and $v$ are the unique solutions of the IVPs
% $$ u'' + f(t)u = 0 , \quad u(0) = 1, ~ u'(0) = 0 $$
% and 
% $$ v'' + f(t)v = 0 , \quad v(0) = 0, ~ v'(0) = 1. $$
% If we write the vector
% $\bf c$ as ${\bf c} = (a,b)^T$, the solution (19.18) (here
% with ${\bf g}={\bf 0}$) becomes $y(t)=au(t) + bv(t)$.
% The solution will be $2\pi$-periodic if and only if
% $y(2\pi) = y(0)$ and $y'(2\pi) = y'(0)$, that is,
% $$ {\bf D} \pmatrix{a\cr\noalign{\vskip 3pt}b} = 
% \pmatrix{a\cr\noalign{\vskip 3pt}b} , $$
% with ${\bf D}$ defined by
% $$ {\bf D} = 
% \pmatrix{u(2\pi) & v(2\pi) \cr\noalign{\vskip 3pt} u'(2\pi) & v'(2\pi)} .
% \eqno (19.25) $$
% The criticality condition will be
% $$ {\bf D} \hbox{ has an eigenvalue equal to\/ } 1. \eqno (19.26) $$
% By the usual reasoning we see that if ${\bf D}$ has
% no eigenvalue equal to $1$, then (19.22) has a unique periodic solution,
% namely $y=0$, corresponding to $a=b=0$,
% whereas if ${\bf D}$ has an eigenvalue equal to $1$, then 
% (19.22) has infinitely many periodic solutions, corresponding to a
% nonzero eigenvector $(a,b)^T$.
% </latex>

%%
% We summarize these results in a theorem.

%%
% <latex>
% \smallskip
% {\em
% {\bf Theorem 19.3.  Hill's equation
% (\textsf{\textbf{fLaSHi}}).}
% Consider the Hill equation
% $$ y'' + f(t)y = 0 , \eqno (19.27) $$
% where $f(t)$ is real, continuous, and $2\pi$-periodic. 
% If the equation is noncritical as defined by $(19.25)$--$\kern 1pt (19.26)$,
% then there exists a unique $2\pi$-periodic solution $y$, namely
% $y(t) = 0$.
% If it is critical, then there exist infinitely many
% $2\pi$-periodic solutions.
% }
% \smallskip
% </latex>

%%
% <latex>
% The equations (19.22)--(19.27) are homogeneous, meaning
% that in the generic case, the only solution
% is the zero solution.  One reason such equations are
% interesting nonetheless is that we are effectively speaking here
% of eigenvalue problems.  In the Mathieu equation
% (19.23), for example, $-a$ plays the
% role of an eigenvalue of the operator $y \mapsto y'' + b\cos(t)
% \kern .5pt y$ with periodic boundary conditions,
% For $b=0$, the nonzero eigenvalues are $-1, -4, -9, \dots,$ and each
% is a double eigenvalue because of the translational symmetry.
% For nonzero $b$ the symmetry is broken and the eigenvalues become
% simple.  This has dynamical consequences, a hint of which appears
% in Exercise 19.10.
% </latex>

%%
% <latex>
% We close this chapter with a few more remarks concerning
% periodic solutions of autonomous nonlinear ODE\kern .5pt s, which are
% $T$-periodic for any $T$.
% For a second-order scalar problem of this type, or a system of
% two first-order equations, a phase plane  analysis may
% quickly give an understanding of what periodic solutions may exist.
% For example, the phase plane portrait of the damped linear
% oscillator (8.4) makes it clear that $y=0$ will be the
% unique $T$-periodic solution of this equation for any $T$.  For
% another example, the phase plane of the nonlinear pendulum
% equation (9.7),  $y'' =-\sin(y),$ explains why this problem
% has nonzero $T$-periodic solutions for all $T$ in the range
% $(2\pi, \infty)$ (Exercise 19.5).
% Similarly, the Lotka--Volterra system (10.3)
% has nonzero $T$-periodic solutions for a range of
% values of $T$ (Exercise 19.8).  The van der Pol oscillator (1.2),
% with its unique limit cycle, has nonzero $T$-periodic
% solutions only when $T$ takes exactly the right value, as mentioned
% below (19.3) (or an integer multiple of this value).
% </latex>

%%
% <latex>
% Suppose we have an autonomous problem like the van der Pol 
% oscillator with a limit cycle, corresponding to a nonzero
% $T$-periodic solution for a particular value of $T$.
% As mentioned earlier, this solution will not be unique, since
% $t$ can be shifted by any amount $\Delta t$, but one can make
% it unique, or make its multiplicity finite rather than infinite,
% by imposing one or more {\bf anchor conditions.}
% For example, here is the van der Pol equation with 
% the anchor condition $y(0) = 1$ and an unknown period $T$:
% $$ 0.3y''-(1-y^2)y'+y = 0, ~~ y(0) = 1, ~~
% y \hbox{ $T$-periodic}. \eqno (19.28) $$
% From a phase diagram such as Fig.~9.11, it is
% clear that there will be two solutions to this problem (for
% minimal $T$), depending
% on where in the oscillation the anchor condition is reached, and
% the two solutions will have $y'(0)>0$ and $y'(0) < 0$.
% </latex>

%%
% <latex>
% Although (19.28) has the form of a $T$-periodic equation for
% unknown $T$, it can be converted to a $2\pi$-periodic problem
% if we transform the independent variable to 
% $s = 2\pi t/T\in [\kern .3pt 0,2\pi]$:
% $$ 0.3{T^2\over 4\pi^2}{d^{\kern .7pt 2} y\over d s^2} -
% {T\over 2\pi}(1-y^2){dy\over ds} +y = 0, ~~ 
% y ~ 2\pi\hbox{-periodic with } y(0) = 1. \eqno (19.29) $$
% This is a $2\pi$-periodic problem with an unknown parameter, but
% a trick that is sometimes used is to think of $T$ not as a number
% but as a solution to the trivial differential equation
% $$ {dT\over ds} = 0, ~~
% T \hbox{ is } 2\pi\hbox{-periodic}. \eqno (19.30) $$
% The unknown parameter has become an unknown constant value of
% one component of a solution of a system of ODE\kern .5pt s, namely
% (19.29)--(19.30)!  This formulation
% is not really simpler in any genuine way, but it has the
% advantage that one can now attempt to solve the problem simply
% by applying a software tool
% for solving nonlinear ODE\kern .5pt s.
% </latex>

%%
% <latex>
% On p.~\pageref{question3} the question (3) was raised of nonperiodic
% solutions of periodic ODE\kern .5pt s.  This is a very important subject which is
% dealt with by {\em Floquet theory.}  For a scalar problem, the main
% result asserts that every solution can be written as a complex
% exponential $\exp(iax)$ for some $a$, real or complex, times a
% periodic function, a pattern readily seen in Fig.~19.4.
% For systems of equations $a$ becomes a matrix $A$.
% We shall not discuss this further except in the following application.
% </latex>

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: band gaps and forbidden frequencies}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% The material world is made of atomic and molecular structures
% that are often periodic, at least on
% a microscopic scale, and electromagnetic
% waves propagate through these structures in distinctive ways.
% The resulting field of X-ray crystallography has been the basic tool
% by which structure of all kinds of molecules are determined,
% as recognized by Nobel Prizes to von Laue
% (1914), Bragg and Bragg (1915), Watson and Crick (1962),
% Hodgkin (1964), and others.\footnote{William
% and Lawrence Bragg were father and son, and when the
% prize was announced, 25-year-old Lawrence was serving in the trenches
% in World War I.\ \ Incredibly, he remained in the trenches for
% a further year after winning the award.
% Nobel Prize winners get better treatment nowadays.}  With
% the advent of quantum mechanics it was seen that electrons
% too are associated with waves, whose distinctive modes of propagation
% though periodic structures determine whether a material is
% an insulator, a conductor, or a semiconductor.
% </latex>

%%
% <latex>
% Here we shall give an indication of a fundamental 
% phenomenon of wave propagation in crystals that was investigated
% by Felix Bloch, another Nobel Prize winner, in the
% beginning of the quantum era.
% Let us return to the Schr\"odinger equation as given in
% (6.11), now (just for convenience)
% with periodic boundary conditions,
% $$ -y''/2 +  V(x) y = \lambda y, \quad x\in [\kern .3pt 0,d\kern .5pt ], ~
% y(0) = y(d\kern .5pt ) ,~
% y'(0) = y'(d\kern .5pt ) .~ \eqno (19.31) $$
% The important new feature is that we will take the 
% potential function $V(x)$ to be periodic, consisting
% of a sequence of spikes with spacing $1$,
% $$ V(x) = 60 (\cos(\pi x))^{16}. \eqno (19.32) $$
% (This is a variant of the {\em Kronig--Penney
% model\/} of a 1D crystal.)
% To keep the example simple we set $d=8$, so $V$ extends
% over just 8 periods.
% Here is a plot of the potential, in black, and the
% first 33 eigenfunctions, each one raised up a distance equal
% to its eigenvalue as in Chapter~6.
% </latex>
d = 8; L = chebop(0,d); L.bc = 'periodic';
V = chebfun('60*cos(pi*x)^16',[0 d],'trig');
L.op = @(x,y) -diff(y,2)/2 + V*y;
neigs = 33; [W,D] = eigs(L,neigs);
W = simplify(W,1e-4); e = sort(diag(D));
for k = 1:neigs, W(k) = e(k)+3*W(k); end, plot(W,LW,1), hold on
xlabel('$x$',FS,10,IN,LT), ylabel('eigenfunctions',FS,9,IN,LT)
plot(V,'k',LW,1.4), hold off, ylim([0 100])
title('Fig.~19.7.~~Periodic potential and eigenfunctions',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% The striking thing about this plot is that the 
% lower eigenfunctions are grouped into bands.  At the bottom, hard
% to distinguish, are eight eigenfunctions with eigenvalue $\lambda\approx 5$.
% Then there are eight eigenfunctions less sharply focussed with 
% $\lambda \approx 20$ and another eight broadly distributed with
% $\lambda \approx 45$.  The final 9 eigenfunctions shown correspond to
% values $\lambda$ greater than the maximum value of the
% potential $V(x)$.  Here is a plot.
plot(e,'.k',MS,8), ylim([0 100])
title('Fig.~19.8.~~Eigenvalues',FS,11)
xlabel('$j$',FS,10), ylabel('eigenvalue $j$',FS,9)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The eigenvalues
% fall in groups of 8.  If we had taken thousand or millions
% of periods, these would have fallen
% so close together as to produce effectively a continuum
% of eigenvalues, more or less like this:
% </latex>
for j = 0:3
  plot(j*d+[1 d],e(j*d+[1 d]),'-k',LW,2), hold on
end
title('Fig.~19.9.~~Schematic for a larger crystal',FS,11)
xlabel('$j$',FS,10), ylabel('eigenvalue $j$',FS,9), hold off
set(gca,XTL,{' ',' ',' ',' ',' ',' ',' ',' '})
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% The implications of this kind of structure
% for physics are enormous.  A periodic medium may have
% {\em forbidden frequencies\/} at which no wave can propagate.
% To a condensed matter physicist, these represent an
% {\em energy gap\/} or {\em band gap\/} in the spectrum, and the gaps
% determine how electrons can propagate.
% Roughly speaking, a material is an insulator if it has an energy
% gap, a conductor if there is no gap, and a semiconductor if there
% is a gap but it is very small.
% One place to learn more is in the classic textbook
% {\em Introduction to Solid State Physics} 
% by Kittel, first published in 1966, where it is pointed out
% that the difference in
% conductivity between an insulator and a conductor (a conductor
% of the usual sort, that is, not a superconductor) may
% be a factor as high as $10^{32}$.
% </latex>

%%
% <latex>
% Where do band gaps come from?  To understand this
% we can begin by recalling what the eigenfunctions of the
% Schr\"odinger problem would look like if the potential
% function $V(x)$ were constant rather than periodic.  They would
% simply be sines and cosines, $\sin(kx)$ and $\cos(kx)$,
% for each value of the wave number~$k$.  In a structure with
% period $L$, however, interference effects occur when
% $k$ is close to an integer multiple of $\pi /L$, corresponding
% to a wavelength close to one of the values $2L, 2L/2, 2L/3,\dots .$
% The interfering waves may lead to reflection rather
% than propagation, and in  X-ray crystallography this is known as
% as {\em Bragg reflection.}
% </latex>

%%
% <latex>
% To see the effect concretely, in our computed example,
% let us return to Fig.~19.7 and
% zoom in on the highest eigenfunction of the first
% band, eigenfunction $8$, and the lowest eigenfunction of
% the second band, eigenfunction $9$.
% </latex>
plot(W(:,8:9),LW,2), hold on
for k = 8:9
  w = W{k}; h1 = plot(w,LW,1.6); hold on
  h2 = plot([0 8],mean(w)*[1 1],'--',LW,0.7);
  set(h2,CO,get(h1,CO))
end
xlabel('$x$',FS,10,IN,LT), ylabel('eigenfunctions',FS,9,IN,LT)
plot(V,'k',LW,0.7), hold off, ylim([0 24])
title(['Fig.~19.10.~~Eigenfunctions 8 and 9, ' ...
       'straddling the first band gap'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Both eigenfunctions have the same wavelength, namely 2,
% and if the potential function were constant they would just
% be translates of one another corresponding to exactly the same
% eigenvalue.  However, the periodic potential
% has broken the symmetry between the sines and the cosines.
% Eigenfunction $8$ is concentrated between the spikes and is little
% affected by them, with an eigenvalue not so different from what it
% would be if the potential were zero.  Eigenfunction $9$, on the other
% hand, is concentrated within the spikes, making the potential weigh
% heavily so that the eigenvalue is much higher.
% </latex>

%%
% <latex>
% \smallskip
% {\sc History.}  The study of periodic ODE\kern .5pt s dates to the
% 1880s with the work of Hill, Floquet, and Mathieu.
% George William Hill (1838--1914), unusually for
% mathematicians of that era, was an American, and he was
% concerned with astronomical calculations related to the
% three-body and four-body problems.
% The general theory is due to
% Gaston Floquet (1847--1920) in ``Sur les
% \'equations diff\'erentielles lin\'eaires \`a coefficients
% p\'eriodiques'' (you don't have to speak French to
% translate this one\kern .3pt !), which appeared in {\em Annals of the
% \'Ecole Normale Sup\'erieure,} v.~12, pp.~47--88, 1883.
% Physicists often give credit to Bloch.
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.}  
% An extensive treatment of
% periodic problems can be found in R. Grimshaw, {\em Nonlinear Ordinary
% Differential Equations,} CRC Press, 1991.
% This book includes a chapter on linear periodic problems,
% a chapter on periodic solutions of autonomous nonlinear problems,
% and further chapters on more advanced related topics.
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 19.}  If the coefficients of
% an ODE are\/ $T$-periodic for some\/ $T$, then there
% may or may not exist\/ $T$-periodic solutions.
% If the ODE is linear, there is an associated criticality condition.
% If the equation is noncritical, it has a unique solution, whereas
% if it is critical, it has no solutions or infinitely many solutions
% depending on the inhomogeneous forcing data.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% {\em \underline{Exercise $19.1$}.  Periodic solutions via shooting.}
% Suppose we seek a $2\pi$-periodic solution to
% $y' = -y^3 + \sin(t)$.
% {\em (a)} Produce a plot showing the solutions emanating
% from initial values $y(0) = -1.5,-1.3,\dots, 1.3,1.5$.
% Based on this plot, estimate by eye the value of $y(0)$ that
% gives a periodic solution.
% {\em (\kern .7pt b)} Prove that if $y(0)\in [-1,1]$, then $y(2\pi)\in [-1,1]$.
% {\em (c)} Look up the Brouwer fixed point theorem and state it.
% Use it to prove that this ODE has a periodic solution.
% {\em (d)} Compute this solution with Chebfun using
% the \verb|'periodic'| flag.  You will need to provide an initial
% guess with roughly the form of the true solution.
% What is $y(0)\kern .7pt?$
% \par
% {\em \underline{Exercise $19.2$}.  An unstable variant.}
% Now consider the equation of Exercise 19.1 with the other sign,
% $w' = w^3 + \sin(t)$.
% {\em (a)} Explain why shooting is much more difficult in this
% case with reference to a computer plot showing solutions emanating
% from several initial values $w(0)$.  Instead of the whole interval
% $[\kern .3pt 0,2\pi]$ you may use a shorter interval such as 
% $[\kern .3pt 0,0.5]$ or $[\kern .3pt 0,1]$.
% {\em (\kern .7pt b)}  On the other hand, the plots suggest an easy way
% to reduce this problem to that of Exercise 19.1.  Explain this
% reduction and use it to derive an analytic formula for $w$ in terms
% of the solution $y$ of Exercise 19.1.
% \par
% {\em Exercise $19.3$.  Criticality depends on the interval.}
% {\em (a)} Explain how the theorems of this chapter allow for
% the possibility that a $2\pi$-periodic linear ODE might have
% a unique $2\pi$-periodic solution but not a unique
% $4\pi$-periodic solution.
% {\em (\kern .7pt b)} Find an example of this behavior.
% \par
% {\em \underline{Exercise $19.4$}.  Two electrons and a nucleus.}
% As in Exercise 13.5, consider the idealized problem of two electrons
% of mass 1 and charge $-1$ orbiting a nucleus of
% charge $+2$ fixed at the origin of the $x$-$y$ plane.  Let $z(t)$
% be the position of one electron represented with the
% usual complex variable
% $z = x+iy$, and suppose the configuration is symmetric
% so that the other electron is at $\overline z(t)$.
% {\em (a)} Write down the ODE governing the evolution of $z(t)$
% assuming an inverse-square electrostatic force law with constant $1$.
% {\em (\kern .7pt b)} Plot the trajectories $z(t)$, $t\in [\kern .3pt 0,40\kern .3pt]$
% corresponding to the initial position $z(0) = ia$ and $z'(0) = 1$
% for $a = 1$ and $2$.  
% {\em (c)} Find a value of $a$ (to 3 digits of accuracy or more)
% that gives a periodic solution.  (Use
% periodic boundary conditions in Chebfun or not, as you prefer.)
% What is the period?
% \par
% {\em Exercise $19.5$.  Nonlinear pendulum.}
% Consider $T$-periodic solutions of the nonlinear pendulum
% equation $y'' = -\sin(y)$ satisfying $y(0) = 0$.  Explain why there
% is one such solution for $T\in (0,2\pi]$, three such solutions
% for $T\in (2\pi, 4\pi]$, five such solutions
% for $T\in (4\pi, 6\pi]$, and so on.
% \par
% {\em \underline{Exercise $19.6$}.  Logistic equation with periodic harvesting.}
% As a generalization of Exercise 3.15, consider the equation
% $y' = (1-y/Y)y - \sin(t)^2$, where $Y$ is a positive constant.
% {\em (a)} Solve the equation with Chebfun with $Y=5$ for $t\in [\kern .3pt 0,15]$
% and make a plot of trajectories 
% from initial values $y(0) = 0.5,1,\dots, 8$.
% {\em (\kern .7pt b)}  Use Chebfun to find the oscillatory periodic solution
% that the curves are approaching.  Use the {\tt mean} command to find
% the mean value of this solution.  What would the mean be (figure this
% part out analytically) if $\sin(t)^2$ were replaced by its average
% value $1/2\kern .5pt $?
% {\em (c)}  Find and plot another periodic solution to this equation.
% Again compare the mean to what it would be 
% if $\sin(t)^2$ were replaced by $1/2$.
% \par
% {\em \underline{Exercise $19.7$}.  Periodic caffeine intake.}
% Figure 2.8 showed the caffeine concentration in the bloodstream
% for a drinker enjoying three cups of coffee in a 24-hour period.  Suppose the
% drinker has the same three cups at the same times every day.  Compute the
% periodic solution, taking the time interval
% to be $[-2,22\kern .3pt]$.  How does
% the maximum caffeine concentration for the periodic solution compare
% to that in Fig.~2.8?
% \par
% {\em \underline{Exercise $19.8$}.  Lotka--Volterra system.}
% {\em (a)} In Chapter 10 the period $T$ of the Lotka--Volterra system (10.6)
% was computed for initial populations $(u_0^{},v_0^{})$ = $(1,1)$ of rabbits
% and foxes, respectively.  Now let $u_0^{}$ vary from $0.2$ to $4$ and
% plot the dependence of $T$ on $u_0^{}$.
% {\em (\kern .7pt b)} For $u_0\to 0.2$, determine $T$ analytically.
% \par
% {\em \underline{Exercise $19.9$}.  Tokieda's teacup.} 
% Place a teacup on a table with the handle positioned in a direction
% we shall call north.  Practice tapping gently on the rim with a spoon
% so that you hear a clean tone.  Note that if you tap in the N, E, S,
% or W positions, you hear one tone, whereas if you tap in the NE, SE, SW,
% or NW positions, you hear a higher tone, often about one semitone higher (a factor
% of $2^{1/2}$).  (There is a video of Tadashi Tokieda demonstrating
% this effect at \verb|https://youtu.be/MfzNJE4CK_s|.)
% To explain this, we can imagine that the cup is a ring
% oscillating in certain eigenmodes. {\em (a)} First, to model a simple
% ring with no handle to break the symmetry, compute the
% first five eigenvalues and eigenfunctions 
% of the differential operator $y''$ with
% periodic boundary conditions.  Note that eigenvalues $2n$ and $2n+1$ are
% equal for each $n\ge 1$ because of rotational symmetry.
% {\em (\kern .7pt b)} Now add the handle --- an extra mass at one point
% along the ring --- by considering the differential operator $y''/m(x)$,
% where $m(x) = 1 + 0.6\exp(-20(x-\pi)^2)$.  Again compute the first
% five eigenvalues and eigenfunctions.  Because of the added mass, all
% the frequencies will now be lower.  Because of the broken symmetry,
% there will no longer be any degeneracies.  It is modes 4 and 5 that we
% mainly hear when tapping the cup,
% corresponding to the ring alternately getting taller/thinner and
% shorter/fatter.  Measure the ratio of frequency 5 to frequency 4 and
% discuss the result.
% Plot the eigenfunctions and explain why they are positioned where
% they are. {\em (c)} What connection do you see with the Application of
% this chapter?  (The next exercise is also related.)
% \par
% {\em \underline{Exercise $19.10$}.  Mathieu equation.} 
% Consider the Mathieu eigenvalue problem
% $y'' + b\cos(t)\kern .3pt y = -\lambda y$
% with periodic boundary conditions on $[\kern .3pt 0,2\pi]$.
% For the rest of this exercise fix $b = 1$.
% {\em (a)} Determine the smallest six eigenvalues $\lambda_1^{} < \cdots 
% < \lambda_6^{}$
% and plot the associated eigenfunctions.  Comment on the relationships
% with the case $b=0$.
% {\em (\kern .7pt b)} Now solve
% $y'' + (a+b\cos(t))y = \exp(\sin(t))$ for $t\in[\kern .3pt 0,
% 100\kern .3pt \pi]$ with initial data $y(0)=y'(0) = 0$
% for four choices of $a$:
% $(\lambda_2^{}+\lambda_3^{})/2,\dots,$
% $(\lambda_5^{}+\lambda_6^{})/2;$ use {\tt L.maxnorm = 0}.
% Plot the results and comment on the
% differences between these curves.
% </latex>
