%% 6. Eigenvalues of linear BVPs

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

%%
% The first BVP of the last chapter, equation (5.1), can be written
% $$ y'' + b\kern .3pt y = 0, \quad x \in [\kern .3pt 0,60\kern.3pt], ~ y(0) = 1, ~y(60)=0
% \eqno (6.1) $$
% with $b=1$, and the solution had amplitude
% about $3.3$.  Look what happens
% if $b$ is reduced by just $1\%$, to $0.99$.
ODEformats, b = 0.99; L = chebop(0,60); L.op = @(x,y) diff(y,2) + b*y;
L.lbc = 1; L.rbc = 0; y = L\0; plot(y,CO,bvp)
title(['Fig.~6.1.~~Solution of (6.1) for $b=0.99$, ' ...
      'near an eigenvalue'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The amplitude has increased to $111.3\kern .5pt!$\ \ With
% $b = 0.9897$, it is about $16606$, and
% with $b = (19\pi/60)^2\approx 0.98970199,$
% it reaches ``infinity'': the BVP has no solution.
% </latex>

%%
% What is happening is that the solution depends
% on a division by a number that is close to zero.
% In this chapter we will explore the mathematics of situations
% like this.  For a simple problem like (6.1), the algebra is elementary,
% but the implications, including the phenomenon
% of resonance for time-dependent partial differential equations
% (see Chapter 22), are far-reaching.

%%
% Let's begin with the elementary algebra.  From equation
% (4.5), we know that
% the general solution of (6.1) is
% $$ y(x) = A \sin(kx) + B\cos(kx) $$
% with $k = \sqrt b$.
% The boundary conditions in (6.1) correspond to the equations
% $$ B = 1, \quad A\sin(60\sqrt b\kern 1pt)
% + B\cos(60\sqrt b\kern 1pt) = 0. $$
% This can be regarded as a $2\times 2$ matrix problem for
% the unknowns $A$ and $B$, an especially easy one 
% since the first equation is trivial.  The solution is
% $$ B = 1, \quad A = - \cot(60\sqrt b\kern 1pt) , $$
% where $\cot x = \cos x / \sin x$ as always.
% Here is the value of $A$ for $b=0.99$.
b = 0.99; A = -cot(60*sqrt(b))

%%
% <latex>
% \noindent
% And here it is for $b=0.9897$.
% </latex>
b = 0.9897; A = -cot(60*sqrt(b))

%%
% <latex>
% \noindent
% With $b = (19\pi/60)^2$, we get $60\sqrt b = 19\pi$,
% and the cotangent is infinite.  Likewise it would be
% infinite with $b = (\kern 1pt j\pi/60)^2$ for any positive integer $j$.
% </latex>

%%
% <latex>
% We are seeing here the phenomenon of eigenvalues of a
% differential operator with boundary conditions. 
% To frame the matter more generally,
% let $L$ be the linear operator $L: y\mapsto y''$ on
% $[\kern .3pt 0,60\kern .3pt]$ and let us consider the problem
% $$ Ly = y'' = \lambda\kern .3pt y + f, \quad x \in [\kern .3pt 0,60\kern.3pt],
% ~ y(0) = \alpha , ~y(60)=\beta , \eqno (6.2) $$
% where $\lambda$ is a given number and $f$ is a function of $x$.
% It is clear that regardless of the choices of $f$, $\alpha$,
% and $\beta$, there can be no unique solution of (6.2)
% if there exists a nonzero function $v(t)$ that satisfies
% the homogeneous equation
% $$ Lv = v''= \lambda v , \quad v(0) = v(60) = 0. \eqno (6.3) $$
% The reason is that any multiple of $v$ could be added to $y$,
% and the result would still satisfy (6.2).
% Such a function is called an {\bf eigenfunction}
% of $L$, and $\lambda$ is the corresponding
% {\bf eigenvalue.}  Another way to say it is that an eigenvalue
% of a linear operator $L$ is a number $\lambda$ such that
% the operator $L-\lambda I$, acting on functions
% satisfying the homogeneous boundary conditions,
% has a nontrivial nullspace.  (Here $I$ denotes the
% identity operator, mapping a function $y$ to itself.)
% An eigenfunction $v$ of $L$
% is a nonzero function in this nullspace.
% Note that any nonzero multiple of an eigenfunction
% is also an eigenfunction for the same $\lambda$.
% </latex>

%%
% Let's do an experiment to scan systematically for eigenvalues of the
% operator $L$ of (6.2).  Since $\lambda$ is the
% negative of $b$ in (6.1), the eigenvalues will be negative numbers, and
% the smallest will be quite close to zero.  Here we arbitrarily pick
% boundary conditions $y(0)=e$ and $y(60)=\pi$
% and the forcing function $f(x) = \exp(-(x-10)^2)$, and we
% solve the ODE successively with $\lambda = -.0300, -.02998, 
% -.02996,\dots, 0$.  For each value of $\lambda$ we plot a
% dot representing $\max_{x\in [\kern .3pt 0,60\kern.3pt]} |y(x)|$.
f = chebfun('exp(-(x-10)^2)',[0 60]); L.lbc = exp(1); L.rbc = pi;
for lambda = -0.03:.0002:0
  L.op = @(x,y) diff(y,2) - lambda*y; y = L\f;
  plot(lambda,norm(y,inf),'.k',MS,7), hold on, drawnow
end
axis([-.03 0 0 400]), hold off, xlabel('$\lambda$',FS,10,IN,LT)
title(['Fig.~6.2.~~$\max_x |y(x)|$ for solutions to (6.2)' ...
       ' with various values of $\lambda$'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The evidence is clear: there are three eigenvalues of $L$
% in this range, lying near $-0.003$, $-0.01$, and $-0.025$.
% Here are plots of the solutions $y(x)$ for these values of $\lambda$.
% </latex>
f = chebfun('exp(-(x-10)^2)',[0 60]); llam = [-.003 -.01 -.025];
for j = 1:3
  lambda = llam(j); L.op = @(x,y) diff(y,2) - lambda*y;
  subplot(1,3,j), plot(L\f,CO,bvp), axis([0 60 -115 115])
  text(30,95,['$\lambda = ' num2str(lambda) '$'],FS,9,HA,CT)
end
subplot(1,3,2)
title('Fig.~6.3.~~Three large solutions of (6.2)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The amplitudes of all three functions are on the order of $100$, and we could
% have made them larger by estimating the eigenvalues more
% accurately.  Just as the three values of $\lambda$ approximate the
% first three eigenvalues of $L$ (i.e., the three
% smallest eigenvalues in absolute
% value), the three functions approximate multiples of the corresponding
% eigenfunctions.  Note that these curves have 1, 2, and 3
% humps, respectively, separated by 0, 1, and 2 interior zeros, called
% {\bf nodes}.  The
% property that the $j$th eigenfunction of $L$ has exactly $j-1$ nodes
% is typical of eigenfunctions for a wide range of differential
% operators.\footnote{In particular this holds for the
% self-adjoint equation $(\kern .6pt p(x)y')' + q(x) y = \lambda\kern .3pt y$
% on an interval $[a,b]$ with homogeneous Dirichlet, Neumann,
% or Robin boundary conditions, where $p'$ and $q$ are continuous
% real functions and $p(x) > 0$.  Such problems are the
% subject of Sturm--Liouville theory.}
% </latex>

%%
% Eigenvalues and eigenfunctions are much too important
% for us to have to find them by 
% ``scanning values of $\lambda$''!
% Chebfun computes these quantities directly
% with the |eigs| command.  We can find the first five eigenvalues of $L$
% like this.
L.op = @(x,y) diff(y,2); sort(eigs(L,5),'descend')'

%%
% <latex>
% \noindent
% To get eigenfunctions as well as eigenvalues we can type
% </latex>
[V,D] = eigs(L,5); llam = diag(D);
[llam,ii] = sort(llam,'descend'); V = V(:,ii);

%%
% <latex>
% \noindent
% A plot of the first three eigenfunctions shows convincing approximations
% of multiples of the curves above.\footnote{As mentioned
% earlier, eigenfunctions are only defined
% up to multiplicative constants.  Those returned by Chebfun are
% scaled to have square-integral equal to 1, a
% normalization that determines their
% shape apart from an arbitrary choice of sign.}
% \label{chap5foot}
% </latex>
for j = 1:3
  subplot(1,3,j), plot(V(:,j),CO,bvp), axis([0 60 -.3 .3])
  text(30,.24,['$\lambda = ' num2str(llam(j)) '$'],HA,CT,FS,9)
end
subplot(1,3,2)
title('Fig.~6.4.~~First three eigenfunctions of (6.3)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% From these images it appears that the $j$th
% eigenfunction consists of the sine function scaled so
% that exactly $j$ lobes, i.e.\ $j/2$ periods, fit in the
% interval $[\kern .3pt 0,60\kern .3pt]$,
% $$ v_j^{}(t) = \sin(\kern 1pt j\pi t/60) , \quad j = 1,2,\dots .
% \eqno (6.4) $$
% Substituting this trial solution in (6.3) confirms that it is
% indeed an eigenfunction.  The associated eigenvalue comes out as
% $$ \lambda_j^{} = -(\kern 1pt j\pi/60)^2,
% \quad j = 1,2,\dots,  \eqno (6.5) $$
% and this formula enables us to confirm the five numbers computed above,
% </latex>
-((1:5)*pi/60).^2

%%
% <latex>
% The eigenfunctions of $L$ are {\bf orthogonal}, which
% means that the number
% $$ a_{ij}^{} = \int_0^{60} v_i^{} (x) v_j^{}(x) dx \eqno (6.6) $$
% is zero whenever $i\ne j$.\footnote{Eigenvectors of real symmetric matrices
% are also orthogonal when they correspond to distinct
% eigenvalues.  The proofs are closely related.}
% To see why, consider two arbitrary functions 
% $v_i^{}$ and $v_j^{}$ satisfying the homogeneous boundary conditions.
% Using integration by parts, we compute
% $$ \int_0^{60}\kern -2pt  (L v_i^{}(x)) v_j^{}(x) dx
% = \int_0^{60}\kern -2pt  v_i'' (x) v_j^{}(x) dx
% = -\kern -2pt \int_0^{60}\kern -2pt  v_i' (x) v_j'(x) dx, $$
% where we have discarded boundary terms that are
% zero because of the boundary conditions on $v_j^{}$.
% A second integration by parts exploiting the zero boundary
% conditions on $v_i^{}$ shows that this is equal to
% $$ \int_0^{60}\kern -2pt  v_i^{} (x) v_j''(x) dx
% = \int_0^{60}\kern -2pt  v_i^{}(x) (L v_j^{}(x))dx. $$
% Thus we have
% $$ \int_0^{60} (Lv_i^{})v_j^{}dx = \int_0^{60}
% \kern -2pt v_i^{} (L v_j^{}) dx. \eqno (6.7) $$
% An operator satisfying this property for arbitrary functions 
% $v_i^{}$ and $v_j^{}$ is said to be {\bf self-adjoint}.  Now suppose
% $v_i^{}$ and $v_j^{}$ are eigenfunctions of a self-adjoint
% operator $L$ corresponding to
% eigenvalues $\lambda_i^{}$ and $\lambda_j^{}$ with
% $\lambda_i^{}\ne \lambda_j^{}$.  Then (6.7) becomes
% $$ \lambda_i^{} \kern -2pt \int_0^{60} \kern -2pt v_i^{} v_j^{} dx =
% \lambda_j^{}\kern -2pt  \int_0^{60} \kern -2pt v_i^{} v_j^{} dx, $$
% and if $\lambda_i^{}\ne \lambda_j^{}$, this implies that
% $v_i^{}$ and $v_j^{}$ must be orthogonal.
% We shall repeat this calculation in greater generality
% for Theorem~7.1 in the next chapter.
% </latex>

%%
% <latex>
% Let us verify orthogonality numerically for the operator $L$ 
% of (6.2).  The object $\bf V$ computed above by
% {\tt eigs} is a {\em quasimatrix\/} of dimensions $\infty\times 5$,
% whose five ``columns'' are the functions $v_1^{}, \dots,v_5^{}$.
% The object ${\bf V}^T {\bf V}$ is the product of the $5\times \infty$ transpose
% of ${\bf V}$ with ${\bf V}$ itself, that is, the $5\times 5$ matrix of entries
% $a_{ij}^{}$ defined
% by (6.6), $1\le i,j \le 5$, known as a {\em Gram matrix}.
% Orthogonality of the
% eigenfunctions is confirmed by the fact
% that this matrix has zero entries off the diagonal.
% </latex>
format short, A = V'*V

%%
% <latex>
% \noindent
% The fact that the entries on the diagonal are equal
% to $1$, making $\bf A$ the identity matrix,
% indicates that the eigenfunctions are not just orthogonal but
% {\bf orthonormal}, a consequence of the
% normalization mentioned in the footnote on p.~\pageref{chap5foot}.
% </latex>

%%
% <latex>
% Eigenvalue problems are associated with
% homogeneous boundary conditions, because they are all about
% identifying nontrivial solutions of the homogeneous problem.
% It does not make sense, for example, to specify $v(0) = 1$
% for an eigenvalue problem.  A variation that does make
% sense, however, is to change a homogeneous boundary 
% condition from Dirichlet to Neumann form.  For example, here is
% what happens if we change the right-hand boundary
% condition of (6.3) from $v(60)=0$ to $v'(60) = 0$.
% </latex>
L.rbc = 'neumann'; sort(eigs(L,5),'descend')'

%%
% <latex>
% \noindent
% The eigenfunctions have zero slope at the right boundary.
% </latex>
[V,D] = eigs(L,5); llam = diag(D);
[llam,ii] = sort(llam,'descend'); V = V(:,ii);
for j = 1:3
  subplot(1,3,j), plot(V(:,j),CO,bvp), axis([0 60 -.3 .3])
  text(30,.24,['$\lambda = ' num2str(llam(j)) '$'],HA,CT,FS,9)
end
subplot(1,3,2)
title('Fig.~6.5.~~Eigenfunctions with a Neumann BC',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The images reveal $1/2$, $3/2$, $5/2,\dots$ 
% lobes of the sine function scaled to $[\kern .3pt 0,60\kern .3pt]$, suggesting
% the general formula
% $$ v_j^{}(x) = \sin((\kern 1pt j-\textstyle{1\over 2})\pi x/60) ,
% \quad j = 1,2,3,\dots . \eqno (6.8) $$
% As before, substitution confirms that these are eigenfunctions, with
% eigenvalues
% $$ \lambda_j^{} = -((\kern 1pt j-\textstyle{1\over 2})\pi/60)^2,
% \quad j = 1,2,\dots.  \eqno (6.9) $$
% The numbers match properly,
% </latex>
-((0.5:4.5)*pi/60).^2

%%
% <latex>
% \noindent
% Again the eigenfunctions are orthonormal,
% </latex>
V = chebfun(V); A = V'*V

%%
% Eigenfunctions of differential operators are not always orthogonal.
% For example, here is an advection-diffusion problem adapted
% from equation (4.9):
% $$ Ly = y'' +0.1 y' , \quad x\in [\kern .3pt 0,60\kern .3pt], ~y(0)=1, ~ y'(0) = 0.
% \eqno (6.10) $$
% Unlike the other operators of this chapter, this one is not self-adjoint.
% Here are the first four eigenfunctions.
% Note that the nonzero boundary value specified in (6.10)
% is ignored by Chebfun in solving the eigenvalue problem, which
% is defined in terms of homogeneous boundary conditions as always.
L = chebop(0,60);
L.op = @(x,y) diff(y,2) + 0.1*diff(y); L.lbc = 1; L.rbc = 0;
[V,D] = eigs(L,5); llam = diag(D);
[llam,ii] = sort(llam,'descend'); V = V(:,ii); clf
for j = 1:4
  subplot(2,2,j), plot(V(:,j),CO,bvp), axis([0 60 -.5 .5])
  text(57,.30,['$\lambda = ' num2str(llam(j)) '$'],FS,9,HA,RT)
end
subplot(2,2,1)
title('\kern 2.5in Fig.~6.6.~~Eigenfunctions of (6.10)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The matrix of inner products confirms that
% the eigenfunctions are nonorthogonal.
% </latex>
V = chebfun(V); A = V'*V

%%
% <latex>
% \noindent
% Here are some higher eigenfunctions, corresponding to
% eigenvalues $\lambda_8$, $\lambda_{16}$, $\lambda_{24}$,
% and $\lambda_{32}$.  The higher the
% eigenvalue, the more oscillations.
% </latex>
[V,D] = eigs(L,32); llam = diag(D);
[llam,ii] = sort(llam,'descend'); V = V(:,ii);
for j = 1:4
  subplot(2,2,j), plot(V(:,8*j),CO,bvp), axis([0 60 -.5 .5])
  text(57,.30,['$\lambda = ' num2str(llam(8*j)) '$'],FS,9,HA,RT)
end
subplot(2,2,1)
title('\kern 2.45in Fig.~6.7.~~Some higher eigenfunctions',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: eigenstates of the Schr\"odinger equation}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% One of the biggest scientific discoveries of all time was
% the Schr\"odinger equation of quantum mechanics,
% in 1926, which in principle
% reduces much of the subject of chemistry to self-adjoint eigenvalue calculations.
% It is from this date that eigenvalues and eigenfunctions became
% a universally known tool in the physical sciences.
% </latex>

%%
% <latex>
% The steady-state (time-reduced) 1D Schr\"odinger
% equation on a finite or infinite
% interval $[-d,d\kern .5pt]$ takes the form
% $$ -h^2 y'' +  V(x) y = \lambda\kern .3pt y, \quad x\in [-d,d\kern .5pt ], ~
% y(-d\kern .5pt ) = y(d\kern .5pt ) = 0, \eqno (6.11) $$
% where $V(x)$ is a fixed potential function and
% $h$ is Planck's constant.
% Equation (6.11) describes what states a particle may occupy,
% in the quantum theory, if it is confined to the given interval 
% and subject to an energy potential $V(x)$ there.
% (More precisely, $h^2$ in (6.11) should be
% $h^2/8\pi^2 \mu$, where $h$ is the true
% Planck constant and $\mu$ is the mass of the particle.)
% The eigenvalue $\lambda$ corresponds to the energy of a particle
% in that state, and its probability density of being
% at each point $x\in [-d,d\kern .5pt]$ is equal to 
% $|y(x)|^2$, where $y$ is the corresponding normalized
% eigenfunction.\footnote{This is an ODE problem
% because there is just one particle in one dimension.
% For $n$ particles in three dimensions,
% it would be a PDE in $3n$ independent variables, far too difficult
% in most cases to be solved numerically.  That is why we say
% that Schr\"odinger's equation reduces chemistry to mathematical
% calculations ``in principle.''
% In practice, chemists have developed powerful methods for
% approximating the equations for multiple-particle systems,
% including those that led to the Nobel Prize for Kohn and
% Pople in 1998.}
% </latex>

%%
% <latex>
% The simplest choice of $V(x)$ would be a constant; the
% corresponding eigenfunctions are sines and cosines.  Next
% simplest is a parabola such as $V(x) = x^2$, corresponding
% to a linear restoring force as in a Hookean spring; this is
% the {\em quantum harmonic oscillator.}
% Here we use Chebfun to compute the first ten eigenstates
% of this problem with $h=0.1$, taking $d=3$ as a good
% approximation to $d=\infty$.  Following a format standard in physics,
% we show the potential function as a heavy black curve and
% plot each eigenfunction raised up by an amount corresponding
% to its eigenvalue.  Higher curves correspond to higher
% energies.
% </latex>
clf, x = chebfun('x',[-3 3]); V = x^2; h = 0.1;
L = chebop(-3,3); L.bc = 0;
L.op = @(x,y) -h^2*diff(y,2) + V*y;
[W,D] = eigs(L,10); diag(D)'
for k = 1:10, plot(D(k,k)+0.06*W{k},LW,1), hold on, end
xlabel('$x$',FS,10,IN,LT), ylabel('eigenfunctions',FS,9,IN,LT)
title(['Fig.~6.8.~~First ten eigenstates of ' ...
      '(6.11) for harmonic oscillator'],FS,11)
plot(V,'k',LW,1), axis([-1.5 1.5 0 2.3]), hold off
%% 
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% This problem can
% be solved exactly with Hermite polynomials (for $d=\infty$);
% the eigenvalues are $1,3,5,\dots$ times
% $h$, hence $0.1, 0.3, 0.5,\dots$ for our choice $h=0.1$.
% </latex>

%%
% The energy levels of systems like this
% determine the states of atomic and molecular systems and
% of the photons they emit or absorb.  A photon emitted by
% a radiating system will have a wavelength corresponding to
% the difference between two energies, and much of what we
% know about objects outside our solar system, including
% the expanding universe and the Big Bang, comes from
% analysing such wavelengths.
% Closer to home, the glowing color of a ``neon'' fabric comes
% from the phenomenon of fluorescence, in which the fabric
% absorbs light at ultraviolet wavelengths and radiates it away
% at a wavelength corresponding to the difference between
% two quantum energy eigenstates.

%%
% <latex>
% Chebfun has a command {\tt quantumstates} for
% automating such calculations (which sets $h=0.1$ by default).
% We illustrate this here with a modified function~$V$
% by adding a small peak in the middle, a
% {\em potential barrier,}
% making $V$ what is known as a {\em double-well potential.} 
% </latex>
V = x^2 + 1.5*exp(-(x/.25)^4);
[U,D] = quantumstates(V); axis([-1.5 1.5 0 2.3])
eigenvalues = diag(D)'
title('Fig.~6.9.~~Eigenstates of a double well potential',FS,11)
xlabel('$x$',FS,10,IN,LT), ylabel('eigenfunctions',FS,9,IN,LT)
%% 
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% The effect on the lower eigenvalues is dramatic if you
% look at the numbers printed above: the lower ones now
% fall in {\em nearly degenerate pairs.}  This is because the potential
% barrier nearly decouples
% the left and right halves, so that it is approximately as if
% we had two identical, uncoupled copies of a single-well
% eigenvalue problem. 
% The lowest eigenvalue
% $0.4436$ corresponds to an even eigenfunction, whereas the next
% lowest, only slightly higher at $0.4459$, corresponds to an
% odd eigenfunction that is almost equal to the same function multiplied
% by $\hbox{sign\kern .5pt}(x)$.  The third and fourth
% eigenvalues and eigenfunctions have similar behavior, but
% not as pronounced; at this energy level the two
% wells are less perfectly decoupled.
% This kind of close but inexact agreement of eigenvalues
% as a result of approximate
% but imperfect symmetry, known as {\em line splitting,}
% is of tremendous importance in physics and led to Nobel Prizes for
% the Zeeman effect in 1902 (line splitting in a 
% magnetic field) and the Stark effect in 1919 (line splitting in
% an electric field).
% </latex>

%%
% <latex>
% Further remarkable things happen if we change the
% potential of Figure~6.9 slightly, moving
% the barrier to the right a distance $0.01$.
% The eigenvalues, still in nearly degenerate pairs, do
% not change very much.  The lower eigenfunctions, however,
% change completely.  Now they are {\em localized,}
% with the lowest eigenfunction being concentrated
% in the bigger well, the
% next-lowest in the narrower well, and so on.
% In Figure 6.10 we plot just these two lowest
% eigenfunctions in the two cases of the symmetric-well
% potential and the slightly asymmetric one.
% </latex>

V2 = x^2 + 1.5*exp(-((x-.01)/.25)^4);
[U2,D] = quantumstates(V2,'noplot');
eigenvalues = diag(D)'
subplot(1,2,1), plot(U(1:2)), axis([-1.4 1.4 -1.9 1.9])
text(.2,-.7,'symmetric:',FS,9)
text(.2,-1.05,'mode 1 even',FS,9)
text(.2,-1.4,'mode 2 odd',FS,9)
subplot(1,2,2), plot(U2(1:2)), axis([-1.4 1.4 -1.9 1.9])
text(.2,-.7,'asymmetric:',FS,9)
text(.2,-1.05,'mode 1 left',FS,9)
text(.2,-1.4,'mode 2 right',FS,9)
subplot(1,2,1)
title('\kern 2.4in Fig.~6.10.~~Eigenfunction localization',FS,11)
%% 
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% More advanced phenomena of localization are related
% to transparency of glass, insulating properties of
% ceramics, and further Nobel Prizes.  In fact, it can
% be argued that eigenvalue problems
% are related to the Physics Nobel Prizes of 
% 1902 (Zeeman),
% 1903 (Becquerel, Curie \& Curie),
% 1907 (Michelson),
% 1911 (Wien),
% 1917 (Barkla),
% 1918 (Planck),
% 1919 (Stark),
% 1921 (Einstein),
% 1922 (Bohr),
% 1923 (Millikan),
% 1924 (Siegbahn),
% 1927 (Compton),
% 1929 (de Broglie),
% 1930 (Raman),
% 1932 (Heisenberg),
% 1933 (Schr\"odinger \& Dirac),
% 1943 (Stern),
% 1944 (Rabi),
% 1945 (Pauli),
% 1952 (Bloch \& Purcell),
% 1954 (Born),
% 1955 (Lamb),
% 1956 (Schockley, Bardeen \& Brattain),
% 1961 (M\"ossbauer),
% 1962 (Landau), 
% 1963 (Wigner, Mayer \& Jensen),
% 1964 (Townes, Basov \& Prokhorov),
% 1965 (Tomonago, Schwinger \& Feynman),
% 1966 (Kastler)$,\dots.$
% </latex>

%%
% <latex>
% \smallskip
% {\sc History.}  The starting ideas of eigenvalue analysis
% go back about 200 years to Fourier, Poisson, Sturm, and Liouville.
% Sylvester and Cayley were diagonalizing
% matrices in the 1850s, and by the 1880s the analysis of
% eigenfunctions for vibrating membranes
% and the theory of sound was well established.  The terms
% ``eigenvalue'' and ``spectral theory,'' however, seem to have come
% later, coined by Hilbert around the turn of the 20th century.
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.}  Although eigenvalues are everywhere
% in the mathematical sciences, it is hard to find
% discussions of {\em why} they matter so much.  What is it about
% the seemingly arbitrary algebraic property $L y = \lambda\kern .3pt y$ that
% demands our attention?  Such a discussion can be found in Chapter 1 of
% Trefethen and Embree,
% {\em Spectra and Pseudospectra: The Behavior of Nonnormal Matrices
% and Operators,} Princeton U. Press, 2005, where four answers to
% the question are proposed.
% \smallskip
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 6.}  
% A second-order scalar linear ODE\/ $Ly = y'' + a(x) y' + b(x) y = f(x)$
% has a two-dimensional affine space
% of solutions, assuming $a$, $b$, and $f$ are continuous or
% piecewise continuous.  For an IVP\/ $Ly=f$
% with two initial conditions, it follows that
% there exists a unique solution.  For a BVP with two boundary
% conditions, however, there exists a unique solution if and only if
% there is no nonzero solution of the homogeneous problem $Lv = 0$ with
% homogeneous boundary conditions.  Such a solution is
% called an eigenfunction of $L$ with eigenvalue $0$. 
% More generally, if\/ $\lambda$ is a number and $v$ is a nonzero
% solution of\/ $Lv=\lambda v$ with homogeneous boundary
% conditions, then $v$ is an eigenfunction of\/ $L$ with
% eigenvalue~$\lambda$. \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>
% \small\parindent=0pt\parskip=1pt
% {\em \underline{Exercise $6.1$}.  Generalized eigenfunctions and an approximate
% delta function.}  Consider the generalized eigenvalue problem
% $y''(x) + \lambda F_\varepsilon^{}(x) y(x) = 0$ with $y(-1) = y(1) = 0$,
% where $F_\varepsilon^{}(x)$ is $(2\kern .3pt\varepsilon)^{-1}$ for
% $|x|<\varepsilon$ and $0$ otherwise.  The first $k$ eigenvalues 
% $D(1,1),\dots,D(k,k)$ and
% eigenfunctions $V\{1\},\dots, V\{k\}$ can be obtained by executing
% {\tt [V,D] = eigs(L,M,k)} with
% {\tt x = chebfun('x'),}
% \verb|F = (abs(x)<ep)/(2*ep),|
% \verb|L = chebop(@(x,y) diff(y,2),[-1 1],0,0),|
% \verb|M = chebop(@(x,y) -F(x)*y).|  
% {\em (a)} Plot the first eight eigenfunctions with $\varepsilon = 0.5$
% and also with $\varepsilon = 0.1$.  For plotting
% consistency, normalize each one by multiplying
% it by $-1$ if it is negative at $x = 0.75$.
% {\em (\kern .7pt b)} Which eigenfunctions are even and which are odd?  How many
% zeros do they have in $(-1,1)$?
% \par
% {\em Exercise $6.2$.  Even and odd eigenfunctions.} 
% {\em (a)} Prove that if $u$ is a function on $[-1,1]$ that is 
% neither even nor odd, then $u(x)$ and $u(-x)$ are linearly independent.
% {\em (\kern .7pt b)} 
% Let $L$ be an ``even'' linear operator 
% on $[-1,1]$ in the sense that if $L$ maps $u(x)$ to $v(x)$,
% then it maps $u(-x)$ to $v(-x)$.
% Prove that if $v(x)$ is an eigenfunction of $L$ with
% eigenvalue~$\lambda$,
% then $v(-x)$ is an eigenfunction of $L$ with the same $\lambda$.
% {\em (c)} Assuming the eigenvalues of~$L$ are all of multiplicity one,
% show that every eigenfunction of $L$ is either even or odd.
% \par
% {\em \underline{Exercise $6.3$}.  Flute and clarinet.} 
% Let $p(x,t)$ denote the deviation from atmospheric pressure at time $t$ at a distance
% $x$ along a flute or clarinet of length $D$ that is playing its lowest note with all its
% keys closed.  A not unreasonable
% approximation is that $p$ satisfies the partial differential equation
% $\partial^2 p/\partial t^2 = c^2 \partial^2p/\partial x^2$, where $c$ is the
% speed of sound.  A further not unreasonable approximation is that the note will
% consist of a superposition of waves $p(x,t) = \sin(\omega_j^{} t) v_j^{}(x)$,  where
% $v_j^{}(x)$ is the $j$th eigenfunction of
% the differential operator $Ly = c^2y''$ and the corresponding eigenvalue is
% $-\omega_j^2$.
% {\em (a)} For a flute, the instrument is open to the atmosphere at both ends and the
% appropriate boundary conditions are $y(0) = y(D) = 0$.  Assuming
% $c = 340$ m/s and $D = 0.66$ m, what would you expect the first three values
% $\omega_1^{}, \omega_2^{},
% \omega_3^{}$ to be?  What frequency
% in Hertz does $\omega_1^{}$ correspond to? 
% What note of the scale is this closest to?
% (Middle C is 261.6 Hz and the A above it is 440 Hz.)
% {\em (\kern .7pt b)} For a clarinet, it is the velocity rather than the pressure deviation
% that is zero at the mouthpiece, and this leads to 
% boundary conditions $y'(0) = y(D) = 0$.
% Answer the same questions as before, now taking $D = 0.60$ m.
% {\em (c)} The flute has an {\em octave key} that opens a hole
% about half-way along, which
% raises certain notes by an octave.  Explain
% in a general way how this works, based on attenuation of
% certain eigenmodes.  (The full details are
% nonlinear, but do not worry about this.)
% {\em (d)} The clarinet has a {\em register key} that opens a hole
% about one-third of the way along, which raises
% certain notes by an interval of a twelfth,
% i.e., a frequency ratio of approximately $3$.
% Adapt your explanation of part {\em (a)} to explain
% why the hole is located where it is and why the frequencies are tripled.
% \par
% {\em \underline{Exercise $6.4$}. Exponential weight.}
% Consider the exponentially weighted eigenvalue problem
% $y'' = -\lambda\kern .3pt  e^x y$ on $[\kern .3pt 0,10\kern .3pt]$
% with $y(0) = y(10) = 0$.  Compute the lowest eight
% eigenvalues and plot the corresponding eigenfunctions, multiplying each one
% by $-1$ if necessary to ensure $y'(0)>0$.  (This is actually
% a generalized eigenvalue problem, as in Exercise 6.1,
% so you will need to define
% two chebops and call e.g.\ \verb|eigs(L,M,8)|.)
% Which is the first eigenfunction
% of this problem that is not positive on $(0,5\kern .3pt]\kern .5pt?$
% \par
% {\em \underline{Exercise $6.5$}. Davies's complex harmonic oscillator.}
% {\em (a)} Use Chebfun to compute the first four
% eigenvalues of the problem
% $-y''+i\kern.3pt x^2y = \lambda\kern .3pt y$ on $(-\infty,\infty)$ with
% $u(\pm \infty) = 0$.  For the computation it will suffice to replace
% $\infty$ by $8$.
% {\em (\kern .7pt b)} The eigenvalues of the usual harmonic oscillator, 
% without the factor $i$, are
% well known to be $1,3,5,7,\dots .$  Use this fact and a suitable change
% of variables to explain the result of part {\em (a)} analytically.
% \par
% {\em \underline{Exercise $6.6$}.  Robin boundary condition.}  
% {\em (a)} Use Chebfun to compute the first six eigenvalues of
% the problem $y'' = - \lambda\kern .3pt y$ on $[\kern .3pt 0,\pi]$ with
% $y(0) = 0$ and $y(\pi) = y'(\pi)$.  Plot the eigenfunctions.
% {\em (\kern .7pt b)} One of the eigenfunctions is quite different from the
% others.  Explain this by working out the form of the eigenfunctions
% analytically.  
% </latex>
