%% 7. Variable coefficients and adjoints

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

%%
% <latex>
% In this chapter we focus on second-order problems
% with variable coefficients.  Locally, near a particular
% value of $x$, solutions look
% like what we've seen in Chapters 4--6, featuring
% oscillation or exponential growth and decay depending on the
% coefficients.  Globally, over a wider interval of\/ $x$ values,
% these structures can combine in interesting ways.
% </latex>

%%
% <latex>
% To begin, let us consider the simplest variable coefficient
% of all, one that jumps abruptly from one constant
% to another:
% $$ y'' - 10 \kern 1.3pt \hbox{sign}(x)\kern .5pt  y = 0.  \eqno (7.1) $$
% For $x<0$ the coefficient is $10$, and for $x>0$ it is
% $-10$.  Thus this is a second-order ODE with piecewise
% continuous coefficients, generalizing the first-order
% such problems of Chapter 2, and a solution is defined as a function
% $y$ with a continuous derivative that satisfies the equation except
% at the points of discontinuity.\footnote{One way to justify this condition is to
% think of the discontinuous coefficient ODE as the limit of
% smooth ODE\kern .5pt s with sharper and sharper transitions.  The second derivative
% $y''$ cannot be continuous at $x=0$, since it satisfies (7.1).}
% Here is the solution on $[-10,10\kern .3pt]$ with boundary conditions
% $y(-10) = y(10)=1$.
% </latex>
ODEformats, L = chebop(-10,10); L.op = @(x,y) diff(y,2) - 10*sign(x)*y;
L.lbc = 1; L.rbc = 1; y = L\0; plot(y,CO,bvp), ylim([-2 2])
title('Fig.~7.1.~~Equation with a sign change (7.1)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% There is no mystery about this picture.
% On the left, the ODE is $y'' + 10y = 0$, with oscillatory solutions
% $A\sin(\sqrt{10}\kern .7pt x) + B\cos(\sqrt{10}\kern .7pt x)$.
% On the right, it is $y'' - 10y = 0$, with exponentially
% growing and decaying solutions
% $C\exp(\sqrt{10}\kern .7pt x) + D\exp(-\sqrt{10}\kern .7 pt x)$.
% The solution makes a continuously differentiable
% transition between these behaviors at $x=0$ and matches
% the boundary conditions.
% Note in particular the exponential boundary layers to the
% right of $x=0$ and to the left of $x=10$.
% It is not hard to work out this solution exactly by solving four
% equations for the four unknowns $A, B, C, D$ (Exercise 7.1).

%%
% <latex>
% Abrupt changes of coefficients arise frequently in
% applications.  Equally important in practice, and mathematically
% more interesting, are cases where coefficients change
% continuously, and the prototypical ODE of
% this kind is the {\bf Airy equation},
% $$ y'' - x\kern .1pt y = 0.  \eqno (7.2) $$
% Note that $10 \kern 1pt \hbox{sign}(x)$ has given way to simply $x$.
% Again we have one sign for $x<0$ and another for $x>0$,
% but now with a smooth transition.  Here is the solution
% with the same domain and boundary conditions as before.
% </latex>
L.op = @(x,y) diff(y,2) - x*y; y = L\0; plot(y,CO,bvp)
ylim([-20 20]), title('Fig.~7.2.~~Airy equation (7.2)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Qualitatively speaking, this image is much like the last one.
% For $x>0$, the two solutions of
% the ODE are exponentially growing and decaying, slowly near
% $x=0$ and faster near $x=10$.
% For $x<0$, they are oscillatory, with slow oscillations near $x=0$
% and faster ones near $x=-10$.
% The point $x=0$ where the equation changes
% from one type to another is called a {\bf turning point}.
% A notable difference from the last image
% is that the amplitude is ten times larger. 
% This is because, by chance, we are close to an eigenvalue.
% If the interval had been $[-9.7,9.7]$ instead
% of $[-10,10\kern .3pt]$, the roles would have
% been reversed, with the solution to (7.1) having the
% high amplitude (Exercise 7.3).
% </latex>

%%
% <latex>
% The Airy equation arises in applications in
% optics, quantum mechanics, and other fields in the analysis of
% wave effects near interesting edges, including George
% Airy's original application in 1838 related to the
% colors of the rainbow.\footnote{The Airy equation is also
% a prototype in the theory of asymptotics, where it illustrates
% fundamental ideas related to steepest descent contours, 
% asymptotics-beyond-all-orders, and Stokes lines.  One can also
% think of it as a step toward the subject of {\em multiphysics,}
% which is concerned with problems and computations coupling
% physical media of different kinds.}  
% </latex>

%%
% <latex>
% With almost any ODE, it is generally interesting to reduce
% the magnitude of the coefficient of the highest derivative
% term to see what happens, and specifically, to sharpen
% the separation between local and global behavior.  (An ODE
% problem with a small leading coefficient is called a
% {\bf singular perturbation problem}; see Chapters~18 and~20.)
% Here we replace $y''$ by $0.01y''$.
% </latex>
L.op = @(x,y) 0.01*diff(y,2) - x*y;
y = L\0; plot(y,CO,bvp)
title('Fig.~7.3.~~Airy equation with leading coefficient 0.01',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% On the left, instead of 3 or 4 wavelengths, there are now
% about 35, as can be verified by Chebfun:
% </latex>
length(max(y{-10,0},'local'))

%%
% <latex>
% \noindent
% This increase is by a factor of about $10$, that is, $\sqrt{100}$.
% </latex>

%%
% <latex>
% Like any second-order linear homogeneous ODE,
% the Airy equation has a two-dimensional
% vector space of solutions.  One unique solution
% in this vector space (up to
% a constant factor) has the special property of decaying to
% 0 as $t\to \infty$, and a particular normalization
% of this solution is called the {\em Airy function}
% $\hbox{Ai}(x)$.  The curves we have seen above are almost, but
% not quite exactly, multiples of $\hbox{Ai}(x)$, or
% $\hbox{Ai}(10 \kern 1pt x)$ for the last example.
% As it happens, there is
% an integral representation for $\hbox{Ai}(x)$,
% $$ \hbox{Ai}(x) = {1\over \pi} \int_0^\infty \cos\Big(xs+{s^3\over 3}
% \kern 1pt \Big) ds.  \eqno (7.3) $$
% It can be shown that at $x=0$, $\hbox{Ai}(x)$ and its derivative
% take the values
% $$ \hbox{Ai}(0) = {3^{-2/3}\over \Gamma(2/3)}, \quad
% \hbox{Ai}'(0) = {-3^{-1/3} \over\Gamma(1/3)}, \eqno (7.4) $$
% where $\Gamma$ is the gamma function, and with
% these two conditions the solution to (7.2) is
% specified fully.  The standard choice of
% a second linearly independent solution, which goes by the label
% $\hbox{Bi}(x)$, is the solution of (7.2) satisfying
% $$ \hbox{Bi}(0) = {3^{-1/6}\over \Gamma(2/3)}, \quad
% \hbox{Bi}'(0) = {3^{1/6}\over \Gamma(1/3)}, \eqno (7.5) $$
% which has the integral representation
% $$ \hbox{Bi}(x) = {1\over \pi} \int_0^\infty \big(
% e^{xs - s^3/3} + \sin(xs + s^3/3)\big)\kern .5pt ds. \eqno (7.6) $$
% Whereas $\hbox{Ai}(x)$ approaches $0$ as $x\to \infty$,
% $\hbox{Bi}(x)$ approaches $\infty$.
% As $x\to -\infty$, $\hbox{Ai}$ and $\hbox{Bi}$ have asymptotically
% the same amplitude, with a phase difference of $\pi/2$, so that
% one is approximately zero where the other is
% approximately extremal.
% We can plot these two solutions by solving the ODE
% as a ``middle-value problem'' (this is not a standard term), with
% $y$ and $y'$ prescribed by (7.4) and (7.5) at $x=0$.
% </latex>
L = chebop(-10,10); L.op = @(x,y) diff(y,2) - x*y;
g1 = gamma(1/3); g2 = gamma(2/3);
Ai0 = 3^(-2/3)/g2; Aip0 = -3^(-1/3)/g1;
L.bc = @(x,y) [y(0)-Ai0; feval(diff(y)-Aip0,0)];
Ai = L\0; plot(Ai,CO,green)
Bi0 = 3^(-1/6)/g2; Bip0 = 3^(1/6)/g1;
L.bc = @(x,y) [y(0)-Bi0; feval(diff(y)-Bip0,0)];
Bi = L\0; hold on, plot(Bi,CO,purple), hold off
title('Fig.~7.4.~~Airy functions',FS,11), ylim([-2 2])
text(4.4,0.3,'$\hbox{Ai}(x)$',FS,11)
text(1.5,1.25,'$\hbox{Bi}(x)$',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% In applications, turning points of ODE\kern .5pt s may
% arise wherever the nature of the ``physics'' changes.
% Though the variable coefficient at such a point may not be exactly
% linear, it will usually
% be possible to approximate it locally by a linear function.
% Consequently the Airy function is relevant to many problems
% with turning points.  We can illustrate
% these effects by changing the variable coefficient of (7.2) 
% to $\sin(x)$.  Specifically, consider
% $$ 0.003 y'' - \sin(x) y = 1, \quad x\in [-4\pi, 4\pi], ~
% y(-4\pi) = y(4\pi) = 0. \eqno (7.7) $$
% The solution is charming.
L = chebop(-4*pi,4*pi); L.op = @(x,y) 0.003*diff(y,2)-sin(x)*y;
L.lbc = 0; L.rbc = 0; y = L\1; plot(y,CO,bvp), ylim([-40 40])
title('Fig.~7.5.~~Sinusoidal variable coefficient (7.7)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Though one could hardly expect to find an
% exact formula for such a function, the essence of the matter
% is simple enough: oscillation
% in regions with $\sin(x)<0$, exponential growth and decay
% in regions with $\sin(x)>0$.  Locally near each turning point,
% $y(x)$ is approximately equal to a suitably scaled
% Airy function, and the approximation would get more precise
% if the small coefficient $0.003$ were reduced further.
% </latex>

%%
% <latex>
% Airy-like equations feature
% transitions between exponential and oscillatory
% behavior.  When a first-order term is introduced in a 
% variable-coefficient equation, new possibilities arise.
% For example, consider the BVP
% $$ 0.001y''  + x\kern .1pt y' + x\kern .1pt y = 0,
% \quad x\in [-2,2\kern .2pt ], ~ y(-2) = -4, ~
% y(2) = 2.  \eqno (7.8) $$
% The solution features an {\bf interior layer} at the turning
% point $x=0$, which
% would grow sharper if the coefficient $0.001$ were further reduced.
% </latex>
L = chebop(-2,2); L.op = @(x,y) .001*diff(y,2)+x*diff(y)+x*y;
L.lbc = -4; L.rbc = 2; y = L\0; plot(y,CO,bvp)
title('Fig.~7.6.~~Interior layer (7.8)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% We shall look at this and other problems with boundary
% and interior layers in Chapter~20.
% </latex>

%%
% <latex> 
% Having explored some ODE\kern .5pt s with variable coefficients,
% we now turn to an area of ODE theory that is of particular 
% interest when the coefficients are variable:
% the notion of an {\bf adjoint} and the associated property
% of {\bf self-adjointness}.\footnote{Not much in the rest of
% the book depends on this material.} Many theoretical properties
% are related to adjoints, including whether eigenvalues are
% real and whether eigenfunctions are orthogonal.  Many practical
% matters are tied to adjoints too, 
% notably the behavior of solutions with respect to perturbations
% and the examination of how the outputs of processes governed by
% ODE\kern .5pt s depend on the inputs.
% </latex>

%%
% <latex>
% Let $L$ be a second-order linear differential operator acting on
% functions on $[-1,1]$.  The coefficients may be variable,
% so $L$ takes the form
% $$ L y = a(x) y'' + b(x) y' + c(x) y, $$
% which as usual we will write more simply as
% $$ L y = a\kern .5pt y'' + b\kern .5pt y' + c\kern .5pt y. \eqno (7.9) $$
% Now we consider the inner product of $v$ and $Ly$, where $v$ and $y$ are
% two functions defined on $[-1,1]$.
% For full generality, it is necessary in this
% business to include the possibility of complex functions,
% because even if an operator has real coefficients, its eigenfunctions
% and eigenvalues may be complex.  Thus the
% inner product formula includes a bar over $v$ denoting the complex
% conjugate:
% $$ \int_{-1}^1 \bar v \kern .7pt Ly  = \int_{-1}^1 a\kern .5pt  \bar v\kern .5pt y'' + b 
% \kern .5pt \bar v \kern .5pt y' + c \kern .5pt \bar v\kern .5pt y. $$
% (We omit the $dx$ factors for simplicity.)
% If we integrate the first and second terms by parts, this
% becomes
% $$ \int_{-1}^1 \bar v \kern .7pt Ly =
% \left. \vrule height 10pt width 0pt depth 0pt (a\kern .5pt \bar v
% \kern .5pt y' + b\kern .5pt \bar v \kern .5pt y)\right|_{-1}^1
% + \int_{-1}^1 -(a\kern .5pt
% \bar v)' y' - (b\kern .5pt \bar v)' y 
% + c\kern .5pt \bar v\kern .5pt  y, $$
% and if we integrate the first term in the integral
% by parts a second time, we get
% $$ \int_{-1}^1 \bar v \kern .7pt Ly = 
% \left. \vrule height 10pt width 0pt depth 0pt
% (a\kern .5pt \bar v \kern .5pt y'  -(a\kern .5pt \bar v)'y
% + b\kern .5pt \bar v \kern .5pt y)\right|_{-1}^1
% + \int_{-1}^1 (a\kern .5pt\bar v)''y - (b\kern .5pt \bar v)' y
% + c\kern .5pt \bar v\kern .5pt y. $$
% In other words,
% $$ \int_{-1}^1 \bar v \kern .7pt Ly-
% \overline{L^*\kern -1pt v}\kern .7pt y = J(\kern .5pt \bar v
% \kern .5pt,y), \eqno (7.10) $$
% where $L^*$ is the {\bf formal adjoint} of $L$, defined by
% $$ L^* \kern -1pt v = (\bar a\kern .5pt v)'' - (\bar b\kern .5pt v)' + 
% \bar c\kern .5pt v, \eqno (7.11) $$
% and $J(\kern .5pt \bar v \kern .5pt ,y)$ is the bilinear
% {\bf conjunct} or {\bf concomitant}
% of $\kern .5pt\bar v \kern .5pt$ and $y$,
% $$ J(\kern .3pt\bar v\kern .5pt ,y) =
% \left.\vrule height 10pt width 0pt depth 0pt
% (a\kern .5pt \bar v\kern .5pt y'
% - (a\kern .5pt \bar v)' y + b\kern .5pt \bar v\kern .5pt y)\
% \right|_{-1}^1 . \eqno (7.12) $$
% In words, {\em $\int_{-1}^1 \bar v\kern .7pt Ly$ is the same as
% $\int_{-1}^1 \overline{L^*\kern -1pt v}\kern .7pt y$, apart from
% boundary effects.}
% Equation (7.10) is known as {\bf Green's formula,}
% and it is a univariate prototype of the multivariate
% Green's formulas that one finds
% in the study of partial differential equations.
% </latex>

%%
% <latex>
% Comparing (7.9) and (7.11), we see that $L^*$ differs from
% $L$ in three ways: the functions are conjugated,
% the sign of the odd-order derivative
% term is negated, and the parentheses have
% moved so that $a$ and $b$ are differentiated.
% We can write out (7.11) fully as
% $$ L^*\kern -1pt v = \bar a \kern .5pt v'' + 
% (2\kern .5pt \overline{a'} \kern .3pt -\kern .3pt \bar b\kern .5pt ) v' + 
% (\kern .5pt \overline{a''}\kern .5pt -\kern .5pt \overline{b'}
% \kern .5pt +\kern .5pt \bar c\kern .5pt  )v , \eqno (7.13) $$
% and similarly (7.12) can be expanded as
% $$ J(\kern .5pt \bar v\kern .5pt ,y) = 
% \left.\vrule height 10pt width 0pt depth 0pt
% a(\kern .5pt \bar v\kern .5pt y' - y\kern .5pt
% \overline{v'}\kern .8pt ) + (b-a'\kern .3pt
% )\kern .5pt \bar v \kern .5pt y\right|_{-1}^1 . \eqno (7.14) $$
% By taking the adjoint of (7.13), it is readily verified that
% the formal adjoint of the formal adjoint is the original operator,
% $$ (L^*)^* = L . $$
% An operator is {\bf formally self-adjoint} if $L^* = L$, and
% from (7.13) it is readily verified that if the coefficient functions are
% real, then this is the case if
% and only if $a' = b$.
% </latex>

%%
% <latex>
% Thus, for example, if $Ly = y'' + y' + y$, we have
% $L^*\kern -.5pt v = v'' - v' + v$ and $J(\kern .5pt \bar v \kern .5pt,y) = 
% (\kern .5pt \bar v \kern .5pt y'-y\kern .5pt \overline{v'}
% +\kern .5pt \bar v \kern .5pt y)\left.\vrule height 8pt width 0pt depth 0pt
% \right|_{-1}^1$.
% Clearly $L$ is not formally self-adjoint since
% $L$ and $L^*$ are different.
% Taking the arbitrary choices $v(x) = \exp(x)$ and $y(x) = \hbox{Ai}(x)$, 
% here is the left side of (7.10):
% </latex>
L = chebop(-1,1); L.op = @(x,y) diff(y,2) + diff(y) + y;
x = chebfun('x'); v = exp(x); y = airy(x);
v'*(L*y) - y'*(L'*v)

%%
% <latex>
% \noindent
% A calculation of $J(\kern .5pt \bar v \kern .5pt ,y)$
% confirms that the right side of (7.10) is the same:
% </latex>
Jfun = @(vb,y) vb.*diff(y) - y.*diff(vb) + vb.*y;
J = @(vb,y) feval(Jfun(vb,y),1) - feval(Jfun(vb,y),-1);
J(conj(v),y)

%%
% As a similar example with variable coefficients, for the entirely arbitrary
% operator $Ly = x\kern .1pt y'' + J_0^{}(x) y' + \sec(x) y$, we have
% $L^*\kern -1pt v = x v'' +(2-J_0^{}(x)) v' + (\sec(x) -
% J_0'(x) )v$ and $J(\kern .5pt \bar v \kern .5pt ,y) = 
% x(\kern .5pt \bar v \kern .5pt y'-y
% \kern .5pt \overline{v'}\kern .5pt)+(J_0^{}(x)-1) \kern .5pt \bar v \kern .5pt y\left.\vrule height 8pt width 0pt depth 0pt
% \right|_{-1}^1$.
% Here again is the left side of (7.10),
L.op = @(x,y) x*diff(y,2) + besselj(0,x)*diff(y) + sec(x)*y;
v'*(L*y) - (L'*v)'*y

%%
% <latex>
% \noindent
% and here is the right side,
% </latex>
Jfun = @(vb,y) x.*(vb.*diff(y)-y.*diff(vb)) + ...
               (besselj(0,x)-1).*vb.*y;
J = @(vb,y) feval(Jfun(vb,y),1) - feval(Jfun(vb,y),-1);
J(conj(v),y)

%%
% <latex>
% The word ``formally'' in the expressions ``formal adjoint'' and
% ``formally self-adjoint'' means: ``apart from boundary conditions.''
% To speak of a true adjoint without this qualification, we
% need to consider not just an operator but a BVP.\ \ Let $L$
% together with certain homogeneous boundary conditions define a BVP, which
% we will write with a script letter as $\cal L$.\ \ The
% {\bf adjoint BVP} ${\cal L}^*$ is the BVP in which $L$ is changed to
% $L^*$ and the boundary conditions are changed
% to {\bf adjoint boundary conditions,} which
% are homogeneous boundary conditions on $v$ such that the conjunct
% (7.14) is always zero.  There exists a unique set of adjoint boundary
% conditions with this property, though we shall not prove this.
% </latex>

%%
% For example, suppose a BVP for the operator $L$ of
% (7.9) has boundary conditions $y(-1)=y(1) = 0$.
% Then (7.14) reduces to $J(\kern .5pt \bar v \kern .5pt ,y) = 
% a(1)\kern .5pt \overline{ v(1)} \kern .5pty'(1) - a(-1)\kern .5pt 
% \overline{v(-1)} \kern .5pt y'(-1)$.
% The values of $y'(1)$ and $y'(-1)$ could be anything, so to ensure
% $J(\kern .5pt \bar v \kern .5pt ,y) = 0$, we must assume two conditions:
% $$ \overline{v(-1)}\kern .5pt =0, \quad \overline{v(1)}
% \kern .5pt = 0, $$
% or equivalently,
% $$ v(-1) =0, \quad v(1) = 0. $$
% These are the adjoint boundary conditions.

%% 
% <latex>
% For another example,
% if the boundary conditions are $y(-1) = y'(1) = 0$, (7.14) 
% reduces to $J(\kern .5pt \bar v\kern .5pt,y) =
% -a(-1)\kern .5pt \overline{v(-1)}\kern .5pt y'(-1) -a(1)y(1)
% \kern .5pt \overline{v'(1)}\kern .5pt
% + [\kern .5pt b(1)-a'(1)]\kern .7pt \overline{v(1)}\kern .5pt y(1) = 0,$ and this will be zero for all
% $y$ if the following two conditions hold:
% $$ \overline{v(-1)}=0, \quad a(1) \overline{v'(1)}\kern .5pt +
% [b(1)-a'(1)]\kern .5pt \overline{v(1)} \kern .5pt = 0,  $$
% or equivalently
% $$ v(-1)=0, \quad \overline{a(1)}\kern .5pt  v'(1) +
% [\kern .5pt \overline{b(1)}\kern .5pt - \kern .5pt 
% \overline{a'(1)}\kern .5pt ]\kern .3pt v(1) \kern .5pt = 0.  $$
% So these are the adjoint boundary conditions.
% In Chebfun, where each chebop may have boundary conditions attached, 
% the computed adjoint of a chebop is not just the formal adjoint
% but the adjoint BVP, with adjoint boundary conditions.
% </latex>

%%
% <latex>
% Suppose $\cal L$ is a BVP and ${\cal L}^*$ is its
% adjoint.  Then for any functions
% $y$ and $v$ satisfying the BCs of $\cal L$ and 
% ${\cal L}^*$, respectively, the right-hand side of (7.10) will vanish, giving
% $$ \int_{-1}^1 \kern .5pt \bar v \kern .7pt Ly = 
% \int_{-1}^1 \overline{L^*v} \kern .7pt y . \eqno (7.15) $$
% A BVP is {\bf self-adjoint} if it is the same as its
% adjoint BVP, i.e.,
% $$ \int_{-1}^1 \bar v\kern .7pt  Ly= \int_{-1}^1
% \overline{L v}\kern .7pt y, \eqno (7.16)  $$
% again for any $v$ and $y$ satisfying the appropriate boundary conditions (which
% are now the same for both $v$ and $y$).   From this identity we
% can derive the following conclusions.
% </latex>

%%
% <latex>
% \smallskip
% {\em
% {\bf Theorem 7.1.  Self-adjoint eigenproblems
% (\textsf{\textbf{fLaSHi}}).}
% Let $\cal L$ be a self-adjoint
% linear homogeneous BVP.\ \ Then
% the eigenvalues of\/ ${\cal L}$ are real, and
% the eigenfunctions corresponding to distinct
% eigenvalues are orthogonal, i.e., if $y$ and $v$
% are eigenfunctions with eigenvalues $\lambda\ne \mu$, 
% then $\int_{-1}^1 \bar v\kern .5pt y^{} = 0$.
% }
% \smallskip
% </latex>

%%
% <latex>
% \def\qed{\vrule height4pt width3pt depth4pt}
% {\em Proof.}
% Let $y$ be an eigenfunction of $\cal L$ with
% eigenvalue $\lambda$.  Then (7.16) gives
% $$ \lambda \int_{-1}^1 \bar y \kern .5pt y
% = \int_{-1}^1 \kern .5pt \bar y \kern .7pt Ly = 
% \int_{-1}^1 \overline{Ly} \kern .7pt y =
% \bar\lambda \int_{-1}^1  \bar y \kern .5pt y . $$
% Since $y$ is nonzero (part of the definition of
% an eigenfunction), the integral of $\bar y \kern .5pt y$ is nonzero,
% and thus this identity implies $\lambda = \bar\lambda$,
% which establishes the first claim.  Now let $v$ be another
% eigenfunction of $\cal L$ with eigenvalue $\mu\ne \lambda$.
% Then another calculation of the same pattern gives
% $$ \lambda \int_{-1}^1 \bar v \kern .5pt y
% = \int_{-1}^1 \kern .5pt \bar v \kern .7pt Ly = 
% \int_{-1}^1 \overline{Lv} \kern .7pt y =
% \mu \int_{-1}^1  \bar v \kern .5pt y . $$
% Since $\mu \ne \lambda$, this establishes the claim of orthogonality.
% \qed
% \smallskip
% </latex>

%%
% If a self-adjoint operator $\cal L$ has real coefficients, then
% the eigenfunctions can be taken to be real,
% so the complex conjugates in the calculations above can be dropped.

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: adjoints and optimization} \\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% Whenever something has to be perturbed, optimized, or designed
% and there is a differential equation involved, there's a good chance
% there may be a role for the adjoint.  Here we'll give just a hint
% of how this may work.  An interesting paper from which to learn
% more is M. B. Giles and N. A. Pierce, An introduction to the adjoint
% approach to design, {\em Flow, Turbulence and Combustion} 65
% (2000), pp.~393--415.
% </latex>

%%
% We show the idea with an
% artificial example.  Suppose $y$ is the solution of
% $$ y'' + \sin(x/10) y' + y = f(x) = \exp(-(x-k)^2), \quad 
% y(0 ) = y(100) = 0 \eqno (7.17) $$
% for some integer $k$ between $1$ and $100$, and we want to know,
% which choice of $k$ will 
% maximize the integral $I(f) = \int_0^{100} y(x)\kern .7pt dx \kern
% .7pt ?$ 
% You can think of this as a prototype of a problem of finding the maximal
% lift on an airfoil in various flow conditions.

%%
% The obvious approach is to simply solve all 100 of these BVPs.
% However, this will take a while, as we can see by solving
% just the first 10 of them.
L = chebop(@(x,y) diff(y,2) + sin(x/10)*diff(y) + y, [0 100]);
L.lbc = 0; L.rbc = 0;
x = chebfun('x',[0 100]); format compact, format short, tic
for k = 1:10
  f = exp(-(x-k)^2); y = L\f; disp([k sum(y)])
end
toc

%%
% <latex>
% \noindent
% Use of the adjoint enables us to solve one BVP rather
% than 100 of them.  The key observation is that since $I(f)$
% is just a single number that depends linearly on $f$ (i.e.,
% $I$ is a linear functional), it 
% must be possible to write $I$ in the form 
% $$ I(f) = \int_0^{100} v \kern .2pt f $$
% for some function $v$ defined on $[\kern .5pt 0,100\kern .5pt]$.
% (This property of linear functionals
% is called the {\em Riesz representation theorem.}
% In this application all the functions are real, so we dispense
% with complex conjugates.)
% In fact, if $g$ is the constant function $1$ on
% $[\kern .5pt 0,100\kern .5pt]$, then
% $v$ is the solution of the adjoint BVP ${\cal L}^*v = g$.
% To derive this characterization, we note that
% $$ I(f) = \int_0^{100} g\kern .5pt  y . $$
% Therefore the two equations
% $$ {\cal L} y = f, \quad {\cal L}^* v = g  \eqno (7.18) $$
% imply 
% $$ I(f) = \int_0^{100} g\kern .5pt
% y = \int_0^{100} [{\cal L}^* v] y = \int_0^{100} v \kern .5pt
% {\cal L} \kern .7pt y = \int_0^{100} v \kern .5pt f . 
% \eqno (7.19) $$
% By taking advantage of this formula, instead of solving
% 100 ODE BVPs, we can now evaluate 100 inner products.
% First we confirm that the first 10 of these give
% the same numbers as before.
% </latex>
tic, v = L'\1;
for k = 1:10
  f = exp(-(x-k)^2); disp([k v'*f])
end

%%
% <latex>
% \noindent
% Next we evaluate all 100 inner products and
% plot the maximal solution, which turns out to correspond to $k=32$.
% </latex>
tic, v = L'\1; d = [];
for k = 1:100
  f = exp(-(x-k)^2); d(k) = v'*f;
end
[maxint, kmax] = max(d), toc
f = exp(-(x-kmax)^2); y = L\f;
plot(y,LW,1,CO,bvp), xlabel('$x$',FS,10), ylabel('$y(x)$',FS,10)
title(['Fig.~7.7.~~The solution of (7.17) with $k=32$, ',...
       'maximizing $I(f)$'],FS,11)
%%
% <latex> \vskip 1.02 em </latex>

%%
% <latex>
% The function $v$ has an interpretation: it quantifies
% the effect on $I(f)$ of the value of $f$
% at each $x\in [\kern .5pt 0,100\kern.5pt]$.
% Consequently $v$ is also 
% the {\em most efficient shape} that $f$ may take,
% measured in the 2-norm,
% if our goal is to maximize $I(f)$.  That is, suppose we
% want to find a function $f$ with $\|f\|_2^{} = 1$ such that $I(f)$ 
% is maximized.
% Applying the Cauchy--Schwarz inequality to (7.19) gives
% $$ I(f) \le  \|v\|_2^{} \kern .5pt \|f\|_2^{} , $$ 
% with equality if $f$ is a multiple of $v$.  So the maximal
% value of $I(f)$ for $\|f\|_2^{} = 1$ is achieved
% by the choice $f = v/\|v\|_2^{}$, and the value
% is equal to $\|v\|_2^{}$,
% </latex>
maxval = norm(v)

%%
% <latex>
% \noindent
% A plot of $v$ shows that $I(f)$ is
% most sensitive to the values of $f(x)$ for $x\approx 32$, as is
% consistent with Figure 7.7.
% </latex>
plot(v,LW,1,'k')
xlabel('$x$',FS,10), ylabel('$v(x)$',FS,10)
title(['Fig.~7.8.~~The adjoint solution $v$ ',...
       'governs sensitivity of $I(f)$'],FS,11)

%%
% <latex>
% \medskip
% \smallskip
% {\sc History.}  
% George Biddell Airy, who used Airy functions in analysis of the
% colors of the rainbow in 1838, was one of the leading astronomers
% and mathematical physicists of the 19th century.  As Astronomer
% Royal in Britain during 1835--1881, nearly half a
% century, he was largely responsible for the adoption of the longitude
% line through Greenwich as the earth's prime meridian.
% \smallskip
% </latex>

%%
% <latex>
% {\sc Our favorite reference.}  Though there are a hundred
% functional analysis books that discuss adjoints in the abstract,
% it is hard to find more down-to-earth discussions of adjoint
% differential equations in the style of this book.
% One source we like is Max D. Gunzburger, {\em Perspectives in
% Flow Control and Optimization,} SIAM, 2003.  As promised in
% the title, this book gives fascinating perspectives on the
% meaning of adjoints --- even in the opening 8 pages.
% \smallskip
% </latex>


%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 7.}  
% Algebraically, second-order
% linear ODE\kern .5pt s with variable coefficients have the same structure
% as those with constant coefficients: there is a two-dimensional
% space of solutions, and an IVP can always be solved but
% a BVP will have a unique solution if and only if it does not
% correspond to an eigenvalue of the associated homogeneous problem.
% Locally, solutions approximate solutions of ODE\kern .5pt s
% with constant coefficients, especially when the coefficient of
% the second derivative term is small.
% Interesting transitions occur at turning points, where a coefficient
% changes sign; the prototype is $x=0$ for the Airy equation
% $y'' - x\kern .1pt y = 0$.  A linear ODE BVP $\cal L$ has an adjoint
% ${\cal L}^*$ consisting of the formal adjoint operator
% together with adjoint boundary conditions.  If $\cal L$ is
% self-adjoint, its eigenvalues are real and its eigenfunctions
% corresponding to distinct eigenvalues are orthogonal.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>
% \small\parindent=0pt\parskip=1pt
% {\em Exercise $7.1$.  Exact solution of the jump problem $(7.1)$.}
% {\em (a)} Determine the exact solution of (7.1) with $y(-10)=y(10) = 1$.
% {\em (\kern .7pt b)} Determine the eigenvalues of the coresponding
% linear operator with boundary conditions $y(-10) = y(10) = 0$.
% \par
% {\em Exercise $7.2$.  An ODE of Euler.}
% Consider the ODE $y'' + (a/t)y' + (b/t^2) y = 0$.  Show
% that it can be reduced to a simpler problem by the substitution
% $\tau = \log(t)$ and use this method to find the general solution.
% \par
% {\em \underline{Exercise $7.3$}.  Dependence on interval length.} 
% {\em (a)} What is the value of the ratio of $\max_x^{} |y(x)|$ for (7.2) divided
% by the same quantity for (7.1)? {\em (\kern .7pt b)} What happens
% to this ratio if the interval is changed
% from $[-10,10\kern .5pt]$ to $[-9.7,9.7\kern .3pt]\kern 1pt$?
% \par
% {\em \underline{Exercise $7.4$}.  Intersection of $\hbox{\rm Ai}(x)$ and
% $\hbox{\rm Bi}(x)$.}   Use Chebfun {\tt roots} to calculate the
% value $x_0^{}$ closest to $0$ for which 
% $\hbox{Ai}(x_0^{}) = \hbox{Bi}(x_0^{})$.  In Matlab you can calculate
% $\hbox{Ai}(x)$ and $\hbox{Bi}(x)$ with \verb|airy(x)|
% and \verb|airy(2,x)|, respectively.
% \par
% {\em \underline{Exercise $7.5$}.  Airy functions with
% many oscillations.}  Solve the Airy equation
% on the interval $[-D,0\kern .3pt]$ with boundary
% conditions $y(-D) = 1$, $y(0)=0$ for
% $D = 10, 20, 40, \dots, 320$.  Plot the solution in
% each case and measure how long it takes Chebfun to
% solve the differential equation.  For comparison, measure how long
% it takes to construct a chebfun for
% $\hbox{Ai}(x)$ on $[-32,0\kern .3pt ]$ directly from Matlab's {\tt airy}
% command.
% \par
% {\em \underline{Exercise $7.6$}.  Sine oscillations.} 
% The solution of (7.7) shows seven maxima
% in $[-\pi,0\kern .5pt]$ in Fig.~7.5.  What does the
% count become if the coefficient $0.003$ 
% is reduced by a factor of $16\kern .7pt$?
% \par
% {\em Exercise $7.7$.  Taylor series solution of Airy equation.} 
% Exercise 2.10 mentioned a classic analytical solution method for
% linear ODE\kern .5pt s, the use of Taylor series.  Suppose
% $y(t) = \sum_{k=0}^\infty a_k^{} t^k$ is a convergent Taylor series
% solution for (7.2) near $t=0$.
% {\em (a)}  Show that $a_2^{} = a_5^{} = a_8^{} = \cdots = 0$. 
% {\em (\kern .7pt b)}  Give a formula for $a_3^{}, a_6^{}, a_9^{},\dots$
% in terms of $a_0^{}$.
% {\em (c)}  Give a formula for $a_4^{}, a_7^{}, a_{10}^{},\dots$
% in terms of $a_1^{}$.
% {\em (d)}  Show that the radius of convergence of the series is
% $\infty$, implying that Airy functions are analytic throughout
% the complex $t\kern .5pt$-plane.
% \par
% {\em \underline{Exercise $7.8$}.  Eigenvalues of the Airy operator.} 
% {\em (a)} Calculate the roots of ${\rm Ai}(x)$ in $[-10,10\kern .3pt]$.
% {\em (\kern .7pt b)} Calculate the first six eigenvalues
% of the Airy operator $L: y \mapsto y'' - x\kern .1pt y$ with
% homogeneous boundary conditions on $[\kern .3pt 0,\infty)$.
% (It suffices to replace $\infty$ by 15.)
% {\em (c)} Explain this coincidence with an
% analytical verification.  To get the idea it
% may help to plot the eigenfunctions.
% \par
% {\em \underline{Exercise $7.9$}.  Legendre equation and polynomials.} 
% If $n$ is a nonnegative integer, the solution to the BVP
% $[(1-x^2)y']' + n(n+1) y = 0,$ $y(-1) = (-1)^n,$ 
% $y(1) = 1$ has a special property: it is a polynomial, known
% as the {\em Legendre polynomial\/} of degree~$n,$ for which
% the standard notation is $P_n^{}$.
% {\em (a)} Compute $P_{10}^{}$ and $P_{40}^{}$ by solving the
% BVP and plot them.  You can confirm that your calculation
% is successful by comparing these functions with
% {\tt legpoly(10)} and {\tt legpoly(40)}.
% {\em (\kern .7pt b)} The Legendre polynomials are orthogonal:
% $\int_{-1}^1 P_m^{}(x)P_n^{}(x) \kern .5pt dx$ is equal
% to $2/(2n+1)$ if $m=n$ and $0$ otherwise.
% Confirm this numerically by computing the
% inner products of $P_{10}^{}$ with itself,  
% $P_{10}^{}$ with $P_{40}^{}$, 
% and $P_{40}^{}$ with itself.  (In Chebfun, the
% inner product of {\tt u} and {\tt v} is \verb|u'*v|.)
% \par
% {\em \underline{Exercise $7.10$}.  An oscillatory pulse.} 
% {\em (a)} Plot the solution of $0.01y''-0.01y'-x\kern .1pt y=1$ for $x\in [-10,10\kern .5pt]$
% with $y(-10) = y(10)=0$. Determine (i) the maximum value, (ii) the
% width of the highest peak measured at half its height, and (iii) the
% distance over which the oscillations decay by a factor of 10.
% {\em (\kern .7pt b)} Produce the same plots and numbers in the three alternative cases where
% the $y''$ coefficient is reduced to $0.005$, or the $y'$ coefficient,
% or both.
% \par
% {\em Exercise $7.11$.  Condition for self-adjointness.} 
% {\em (a)}
% Use (7.13) to verify that if the coefficient functions are real,
% then the operator $L$ of (7.9) is formally
% self-adjoint if and only if $a' = b$.
% {\em (\kern .7pt b)} What is the analogous condition if
% the coefficients are complex?
% \par
% {\em \underline{Exercise $7.12$}.~~Self-adjoint and nonself-adjoint.}
% Let $\cal L$ be the BVP $e^xy''+e^x y' + x\kern .2pt y=0$
% with boundary conditions
% $y(-1)=y'(1)=0$.
% {\em (a)} Show analytically that $\cal L$ is self-adjoint.
% {\em (\kern .7pt b)} Compute
% the first six eigenvalues and eigenfunctions of $\cal L\kern .3pt$,
% showing that the
% eigenvalues are real and the eigenfunctions are numerically orthogonal.
% {\em (c)}
% Show that if the ODE is changed to $e^x y''+ y' + x\kern .2pt y=0$, the eigenvalues
% remain real but the eigenfunctions are no longer orthogonal, and verify
% that the eigenvalues of $\cal L$ are the same as those of ${\cal L}^*$.
% {\em (d)}
% On the other hand, it can be shown by an adaptation of the argument of
% Theorem~7.1 that the $j\kern .3pt$th eigenfunction of ${\cal L}^*$ is orthogonal
% to the $k\kern .3pt$th eigenfunction of $\cal L$ for $j \ne k$.  Verify
% this numerically.
% \par
% {\em \underline{Exercise $7.13$}.~~Optimizations with the adjoint.}
% {\em (a)} Repeat the optimization illustrated in the Application, but maximizing
% $I(f)=\int_0^{100} x\kern .4pt y(x)\kern .7pt dx$. 
% {\em (\kern .7pt b)} Repeat again, but maximizing
% $I(f)=\int_{30}^{50} y(x)\kern .7pt  dx$.
% </latex>
