%% 20. Boundary and interior layers
% <latex>
% \setcounter{page}{253}
% </latex>

%%
% <latex>
% We have seen examples of boundary and interior layers
% already in Chapters~4, 7, and~18.  This is a topic where
% the use of a tool like Chebfun is particularly compelling.
% With a computation or two, one quickly sees where
% layers lie and how their thickness depends on parameters.
% </latex>

%%
% <latex>
% To start the discussion, here is
% a reprise of the example of equation (5.3),
% $$ y'' = y , \quad x \in [\kern .3pt 0,40\kern .3pt],
% ~y(0) = y(40) = 1.  \eqno (20.1) $$
% </latex>
ODEformats
L = chebop(0,40); L.op = @(x,y) diff(y,2) - y; L.lbc = 1;
L.rbc = 1; y = L\0; plot(y,CO,bvp), axis([0 40 -1.5 1.5])
title('Fig.~20.1.~~Simplest example of boundary layers (20.1)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Just by looking at the image, one sees immediately at some
% level what is going on: the solution ``wants to be zero'' and
% more or less achieves this aim, apart from regions
% near the boundaries,
% where it goes through rapid transitions to satisfy the
% boundary conditions.  What does it mean that it ``wants to be
% zero''?  This ODE is easy to analyze: the general solution
% is $A\exp(-x) + B \exp(x)$.  This implies
% that all components grow exponentially
% either as $x$ increases or as $x$ decreases.  Therefore, the
% only way $y$ can deviate much from zero without deviating
% {\em hugely\/} from zero is for this to happen near a boundary.
% </latex>

%%
% <latex>
% Problems with boundary or interior layers usually contain
% a large or a small parameter.  In (20.1), the parameter is the length
% of the interval, 40.  To move to a more standardized
% formulation, let us
% transplant the equation to the interval $[-1,1]$
% with a small constant $\varepsilon$ multiplying
% the highest-order derivative:
% $$ \varepsilon \kern .3pt y'' = y , \quad x \in [-1,1],
% ~y(-1) = y(1) = 1.  \eqno (20.2) $$
% The rescaling multiplies
% $y''$ by $(20)^2$, so the right choice if we wish
% to match (20.1) is $\kern .5pt \varepsilon = (20)^{-2}$.
% As mentioned in Chapter~7, a problem like this with a small
% parameter multiplying the highest derivative is called a
% {\bf singular perturbation problem.}
% Here is the solution to (20.2), looking the
% same as for (20.1).
% </latex>
L = chebop(-1,1); ep = (1/20)^2;
L.op = @(x,y) ep*diff(y,2) - y; L.lbc = 1; L.rbc = 1;
y = L\0; plot(y,CO,bvp), axis([-1 1 -1.5 1.5])
title('Fig.~20.2.~~Same problem rescaled to $[-1,1]$, (20.2)',FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%

%%
% This example is so simple that we can quickly work out 
% exact formulas.
% The general solution to the rescaled problem (20.2) is
% $A\exp(-\varepsilon^{-1/2}x) +
% B \exp(\varepsilon^{-1/2}x)$, or
% as we may rewrite it with redefinitions of $A$ and $B$,
% $$ y(x) = A\exp(-\varepsilon^{-1/2}(1+x))
% + B \exp(-\varepsilon^{-1/2}(1-x)). \eqno (20.3) $$
% The boundary conditions correspond to the equations
% $$ A + B\exp(-2\kern .5pt\varepsilon^{-1/2}) =
% B + A\exp(-2\kern .5pt\varepsilon^{-1/2}) = 1, $$
% with solution
% $$ A =  B = {1-\exp(-2\kern .5pt\varepsilon^{-1/2})\over
% 1-\exp(-4\kern .5pt\varepsilon^{-1/2})} = 
% 1 + O(\exp(-2\kern .5pt\varepsilon^{-1/2})). $$
% Therefore (20.3) can be written
% $$ y(x) = \exp(-\varepsilon^{-1/2}(1+x)) +
% \exp(-\varepsilon^{-1/2}(1-x)) + O(\exp(-3\kern.3pt\varepsilon^{-1/2})). 
% \eqno (20.4) $$
% At $x=0$, for example, we have
% $$ y(0) = 2 \exp(-\varepsilon^{-1/2})
% + O(\exp(-3\kern .5pt \varepsilon^{-1/2})), $$
% which for $\varepsilon = (20)^{-2}$ becomes
% $$ y(0) = 2\exp(-20) + O(\exp(-60)). $$
% We confirm this result numerically:
y(0), 2*exp(-20)

%%
% <latex>
% Although this example is very simple, it illustrates a general idea.
% Consider the approximation (20.4), a sum of two terms.
% The first term consists of an approximation
% at the left boundary, and the second
% at the right.  Together, they provide an 
% accurate approximation to the solution throughout $[-1,1]$.
% Such combinations are the business
% of {\bf boundary layer analysis,} which constructs
% approximate solutions to all kinds of problems
% by combining different approximations in different
% regions.  This is an advanced art,
% applicable for nonlinear as well as linear problems.
% In the remainder of this chapter, we will take some
% introductory steps into this subject for linear equations.
% </latex>

%%
% <latex>
% The starting point of boundary layer analysis is this principle:
% \begin{enumerate}
% \item {\em Outside a boundary or interior
% layer, terms involving\/ $\varepsilon$ are negligible.}
% \end{enumerate}
% For example, for the equation (20.2), outside any boundary or interior layers,
% we can expect
% $$ y = 0 $$
% to hold approximately.  This equation applicable outside
% layers is called the {\bf outer equation}.
% </latex>

%%
% <latex>
% If (20.2) is modified to
% include a further term, the outer equation will adjust accordingly.  For
% example, suppose we consider
% $$ \varepsilon \kern .3pt y'' = y + 2x , \quad x \in [-1,1],
% ~y(-1) = y(1) = 1.  \eqno (20.5) $$
% Now the outer equation is
% $$ y = -2x, $$
% and the solution changes shape to look like $y = -2x$ away
% from $x = \pm 1$.
% </latex>
L.op = @(x,y) ep*diff(y,2) - y - 2*x;
y = L\0; plot(y,CO,bvp), axis([-1 1 -2 2])
title('Fig.~20.3.~~Modified problem (20.5)',FS,11)
%%
% <latex> \vskip 1.02em </latex>
%%

%%
% <latex>
% To proceed further with boundary layer analysis we add
% a pair of principles applicable inside a transition layer
% (i.e., a boundary or interior layer):
% \begin{enumerate}
% \setcounter{enumi}{1}
% \item
% {\em Inside a layer, a low order derivative
% is negligible compared with a higher order one.}
% \item
% {\em Inside a layer, a variable coefficient can be
% approximated by a constant, at least if it is locally nonzero.}
% \end{enumerate}
% Let us see how these ideas work out in a
% slightly more complicated example, equation (7.8),
% $$ \varepsilon \kern .3pt y''  + x\kern .1pt y' + x\kern .1pt y = 0, 
% \quad x\in [-2,2\kern .2pt ], ~ y(-2) = -4, ~
% y(2) = 2.  \eqno (20.6) $$
% With $\varepsilon = 0.001$, the
% solution shows an interior layer at $x=0$, as we saw already
% in Fig.~7.6.
% </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), ylim([-6 18])
title('Fig.~20.4.~~Interior layer (20.6)',FS,11), hold on
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% For a boundary layer analysis of this
% result, applying principle (1), we replace
% $\varepsilon$ by zero in (20.6) and get the outer equation
% $$ x\kern .1pt y' + x\kern .1pt y = 0 . \eqno (20.7) $$
% Thus we immediately have a prediction as to the behavior
% of a solution to (20.6).  In every region,
% (i) $y''$ will
% be very large (of order $\varepsilon^{-1}$),
% or (ii) $x$ will be close to 
% zero, or (iii) $y(x)$ will approximate $Ce^{-x}$ for some $C$.
% </latex>

%%
% From the picture it is clear that case (iii) applies for 
% $|x|\gg 0$, with outer solutions $-4\exp(-(x+2))$ on the left
% and $2\exp(-(x-2))$ on the right.  Here are thick dashed lines
% marking these approximations:
fleft = chebfun('-4*exp(-(x+2))',[-2,.2]);
h1 = plot(fleft,'--r',LW,2);
fright = chebfun('2*exp(-(x-2))',[-.1,2]);
h2 = plot(fright,'--r',LW,2); ylim([-6 18])
title('Fig.~20.5.~~Outer solutions on left and right',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% To explain the rest of the solution, we
% apply principle (2).  In a transition layer,
% we can expect $y''$ to be much larger than $y'$ and $y'$ to be
% much larger than $y$.   An approximation of (20.6) in a layer
% is accordingly
% $$ \varepsilon \kern .3pt y''  + x\kern .1pt y' = 0. \eqno (20.8) $$
% First, let us use this equation to see why there is no
% boundary layer at $x=2$.  Applying
% principle (3), we approximate (20.8) further by
% $$ \varepsilon \kern .3pt y''  + 2\kern .3pt y' = 0. \eqno (20.9) $$
% The general solution to this equation is
% $A + B\exp(-2\kern .3pt \varepsilon^{-1}x)$, and the
% minus sign in the second of these terms implies that 
% it decreases rather than increases with increasing $x$.
% Similarly, at $x=-2,$
% $$ \varepsilon \kern .3pt y''  - 2\kern .3pt y' = 0, \eqno (20.10) $$
% with general solution $A + B\exp(2\kern .3pt \varepsilon^{-1}x)$,
% and again the sign is such that no boundary layer is possible.
% </latex>

%%
% <latex>
% Only $x=0$ remains a candidate for a transition layer, matching
% the shape of the solution in the plot.
% Let us accordingly consider (20.8) for $x\approx 0$.  Here we
% rescale variables to see what is happening.   We define
% $$ s = x\kern .5pt \varepsilon^{-1/2}, \quad u(s) = y(x), $$
% which leads to
% $$ y''(x) = \varepsilon^{-1} u''(s), \quad
% y'(x) = \varepsilon^{-1/2} u'(s), $$
% whereupon (20.8) becomes
% $$ u'' + s u' = 0, \quad u(-\infty) =
% -4\kern .3pt e^{-2}, ~ u(\infty) = 2\kern .3pt e^2 . \eqno (20.11) $$
% One solution of the ODE part of this equation
% is $u(s) = A$ for any constant~$A$, and the other can
% be obtained by setting $w=u'$, yielding the equation
% $w' + sw = 0$.  Separation of variables gives
% $w(s) = B\exp(-s^2/2)$ for a constant $B$, or after
% integration and redefinition of $B$,
% $u(s) = B \hbox{\kern  1pt erf} (s/\sqrt 2)$, where
% $\kern .5pt\hbox{erf}\kern .5pt$ is the error function.  Thus we have
% $$ u(s) = A + B\hbox{\kern  1pt erf}\left({s\over \sqrt 2 \kern 1pt}\right ) , \quad 
% A = e^2 - 2e^{-2}, ~ B = e^2+2e^{-2}. $$
% Here we superimpose on the plot a new dashed line corresponding to
% this approximation.
% </latex>
set(h1,CO,[1 .7 .7]), set(h2,CO,[1 .7 .7])
A = exp(2)-2*exp(-2); B = exp(2)+2*exp(-2);
ep = .001; s = @(x) x/sqrt(ep);
fmiddle = chebfun(@(x) A + B*erf(s(x)/sqrt(2)),[-.8 .8]);
plot(fmiddle,'--r',LW,2), hold off
title('Fig.~20.6.~~Approximation to the transition layer',FS,11)
%%
% <latex> \vskip .2em </latex>

%%
% <latex>
% This completes our study of (20.6).  Now let us look
% at the same equation, except with a sign change on the
% second derivative term:
% $$ -\varepsilon \kern .3pt y''  + x\kern .1pt y' + x\kern .1pt y = 0,
% \quad x\in [-2,2\kern .2pt ], ~ y(-2) = -4, ~
% y(2) = 2.  \eqno (20.12) $$
% The solution looks completely different.
% Instead of an interior layer, we
% have boundary layers at each end, but they are so thin
% as to be almost invisible in the plot.
% </latex>
L.op = @(x,y) -.001*diff(y,2) + x*diff(y) + x*y;
y = L\0; plot(y,CO,bvp), ylim([-5,5])
title('Fig.~20.7.~~Sign change (20.12)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% An experienced eye immediately
% conjectures that the boundary layers are probably of width
% $O(\varepsilon)$ rather than $O(\varepsilon^{1/2})$, as a
% closeup at the left boundary confirms.
% </latex>
plot(y{-2,-1.98},CO,bvp), ylim([-5,5])
title('Fig.~20.8.~~Closeup near left boundary',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% The scaling of the widths of boundary and interior layers
% is of fundamental importance in applications.
% Let us illustrate this matter further with a constant-coefficient
% equation, an advection-diffusion problem,
% $$ \varepsilon \kern .3pt y'' + y' + y = 0 , \quad
% x \in [\kern .3pt 0,1],~ y(0)=y(1) = -1. \eqno (20.13) $$
% Here are the solutions for $\varepsilon = 0.1,$ $0.02,$
% and $0.004$.
L = @(ep) chebop(@(x,y) ep*diff(y,2)+diff(y)+y,[0,1],'dirichlet');
y = @(ep) L(ep)\(-1);
for ep = [.1 .02 .004]
  plot(y(ep),CO,bvp), ylim([0 2]), hold on
end
title(['Fig.~20.9.~~Solutions to (20.13) with ' ...
  '$\varepsilon = 0.1, 0.02, 0.004$'],FS,11), hold off
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Let us informally define the width of a boundary layer to be 
% the distance from~$0$ to the value of~$x$ at which
% the solution first reaches $y(x)=0.5$ (compare Exercise~18.2).
% The widths for four values of $\varepsilon$
% show convincing $O(\varepsilon)$ dependence.
% </latex>
width = @(ep) min(roots(y(ep)-0.5));
disp('   epsilon   boundary layer width')
for ep = 10.^(-1:-1:-4)
  fprintf('%10.5f      %10.7f\n',ep,width(ep))
end

%%
% <latex>
% Here is another variation on the theme of (20.6) and (20.12),
% with the left-hand boundary condition made positive to enhance the
% effect:
% $$ 0.001y''  + x\kern .1pt y' - y = 0, 
% \quad x\in [-2,2\kern .2pt], ~ y(-2) = 4, ~
% y(2) = 2.  \eqno (20.14) $$
% The transition here is called a {\em corner layer}.
% </latex>
L = chebop(-2,2); L.lbc = 4; L.rbc = 2;
L.op = @(x,y) .001*diff(y,2) + x*diff(y) - y;
y = L\0; plot(y,CO,bvp), ylim([-1 5])
title('Fig.~20.10.~~Corner layer (20.14)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Almost anything is possible with interior and boundary
% layers, and we give this example as just another illustration.
% </latex>

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: why is New York hotter than San Francisco?} \\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% New York is hot in July, around 9 degrees hotter than
% San Francisco ($25^\circ$C vs.\ $16^\circ$ are typical
% 24-hour averages).
% Yet in January, New York is 9 degrees colder than
% San Francisco ($1^\circ$ vs.\ $10^\circ$).
% What's going on?  Why is the weather so much more moderate
% on the west coast of the USA than the east coast?
% We shall explain the effect in words first,
% then explore the mathematics and the connection with boundary layers.
% </latex>

%%
% The first key fact is that in the middle latitudes of Earth, prevailing
% winds blow from west to east.
% This means that San Francisco's weather comes from
% over the ocean, whereas New York's comes from over the land.  

%%
% <latex>
% The second key fact is that ocean is liquid, and land is solid.  
% When the summer sun beats down on the land, it heats up, but only
% down to a depth of a few meters, because the heat transfer is limited
% by conduction, which is slow.  This lets the surface
% get very hot.  Water, on the other hand, transfers heat
% by the much faster mechanism of convection --- fluid motion.
% The stirring isn't enough to bring
% the summer heat down to the great depths, but it certainly brings it deeper
% than a few meters.  The heat is spread over a much greater
% volume, so that the surface of the ocean never gets
% anywhere near as hot as the surface of the land, and San Franciscans 
% can get by without air conditioning.
% </latex>

%%
% <latex>
% To model this effect mathematically, we shall be very rough,
% since the point is just to understand the main mechanism.
% Let us imagine that the earth is 100 meters deep, with a temperature
% $T$ fixed at $13^\circ $C at the bottom (the average of the figures 
% reported above for both New York and San Francisco).
% Our depth variable will be $x$, in units of meters,
% running from $-100$ to $0$.
% At the earth's surface, let us imagine that heat flows in at
% the rate of 100 watts per square meter in July and flows out
% at the same rate in January, with a sinusoidal dependence on $t$ in-between,
% $$ \hbox{incoming heat in } W/m^2 = 100 \exp(2\pi it) , $$
% where $t$ is measured in years.
% (We use a complex exponential for the usual reason of convenience
% explained in Chapter~4; the physical temperature is the real part.)
% To get to an ODE, 
% we start with a PDE, the heat equation,
% $$ \rho \kern .5pt c \kern .5pt {\partial T\over \partial t} =
% \left({\hbox{year}\over\hbox{sec}}\right) k
% \kern .7pt {\partial^2 T\over \partial x^2}, \eqno (20.15) $$
% where $\rho$ denotes density, $c$ is heat capacity, and $k$ is
% conductivity.
% (PDE\kern .5pt s are considered systematically in Chapter 22.)
% We measure $\rho$ in kilograms per cubic meter,
% $c$ in Joules per kilogram per degree, and $k$ in watts per meter per degree.
% The quotient $(\hbox{year}/\hbox{sec})$, which is equal to
% about $3.1\times 10^7$, is included to reconcile our
% choice of units of years for $t$ and watts (Joules per second)
% in the definition of $k$.
% With the not too unreasonable choices
% $$ \rho = 3000, \quad c = 2000 , $$
% (20.15) reduces to 
% $$ {\partial T\over \partial t} =
% 5.3 \kern .5pt k \kern 1pt {\partial^2 T\over \partial x^2}.
% \eqno (20.16) $$
% Meanwhile the boundary condition at $x = -100$ is
% $$ T(-100,t) = 13 , $$
% and our assumption concerning
% the heat influx gives the boundary condition at $x=0$,
% $$ {\partial T\over \partial x}(0,t) = 100 \kern .5pt
% k^{-1}\exp(2\pi it). $$
% We now reduce the problem to an ODE by separating variables, assuming
% the solution takes the form
% $$ T(x,t) = 13 + e^{2\pi it} (y(x)-13) $$
% for some function $y$.  The PDE (20.16) becomes
% $$ 2\pi i (y-13) = 5.3 \kern .3pt k \kern .5pt y'', \eqno (20.17)  $$
% and the boundary conditions become
% $$ y(-100) = 13, \quad y'(0) = 100\kern .5pt k ^{-1}. \eqno (20.18) $$
% </latex>

%%
% <latex>
% We are ready to compute a solution.  Taking $k=3$ as a
% reasonable value of the conductivity, here is the
% temperature distribution we get from (20.17)--(20.18).
% </latex>
L = chebop(-100,0);
k = 3; L.op = @(y) 5.3*k*diff(y,2) -2i*pi*(y-13);
L.lbc = 13; L.rbc = @(y) diff(y)-100/k;
y = L\0; plot(real(y),CO,bvp), grid on, ylim([0 60])
title('Fig.~20.11.~~July temperature in the earth',FS,11)
xlabel('depth $x$',FS,10), ylabel('temperature (${}^\circ$C)',FS,9)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Note how hot it is in the boundary layer at the surface!  New York's summers
% begin to make sense.  Note also the dip below $13^\circ$C in the
% temperature at a depth of about 3 meters.  This effect is genuine,
% reflecting a bit of wavelike behavior in the conduction of
% heat.\footnote{The thick stone walls of Italian
% villas exploit this effect to keep cool by day and warm by night.
% A thickness of $3\hbox{m}/20 = 15$cm ought to be of the right
% order of magnitude, since
% days are about $20^2$ times shorter than years and the
% effect in question involves a square root.}
% </latex>

%%
% As for mild San
% Francisco, we are not about to go into details of the physics of
% convection in the upper ocean.  To get the general idea, however, let's
% imagine that ocean is like land, but with 100 times greater
% conductivity.  The heat now reaches much deeper,
% giving a much more moderate surface temperature.
k = 300; L.op = @(y) 5.3*k*diff(y,2) -2i*pi*(y-13);
L.rbc = @(y) diff(y)-100/k;
y = L\0; plot(real(y),CO,bvp), grid on, ylim([0 60])
title('Fig.~20.12.~~July temperature in the ocean',FS,11)
xlabel('depth $x$',FS,10), ylabel('temperature (${}^\circ$C)',FS,9)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% The details of our models have hardly been precise,
% but the main message of
% Figs.~20.11 and~20.12 is genuine.  New York is hotter than
% San Francisco in the summer because its weather is controlled by
% a thinner boundary layer.
% </latex>

%%
% <latex>
% \smallskip
% {\sc History.}
% Boundary-layer analysis began with an epochal
% contribution by the German fluid mechanician
% Ludwig Prandtl.  An important feature of fluid mechanics is viscosity,
% that is, friction, but viscous effects are often insignificant
% away from solid boundaries.
% Prandtl realized that one could exploit this effect mathematically
% by dropping certain terms from the Navier--Stokes partial differential
% equations of fluid mechanics near a wall and dropping
% other terms away from the wall.  These simplifications led to
% major advances in engineering, making
% it possible to solve problems that had previously been
% inaccessible.  Prandtl presented these methods
% in 1904, and ever since then,
% boundary-layer analysis has been recognized as one of the
% fundamental tools of applied mathematics.
% </latex>

%%
% <latex>
% \smallskip
% {\sc Our favorite reference.}  Generations of
% applied mathematics graduate
% students at MIT have learned asymptotics from
% the exceptionally rich book by Bender and 
% Orszag, {\em Advanced Mathematical Methods for Scientists
% and Engineers,} McGraw-Hill, 1978.
% \smallskip
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 20.} Singular perturbation problems, where
% a small parameter\/ $\varepsilon$
% multiplies the highest derivative, typically have
% solutions featuring rapid transitions known as boundary
% or interior layers.  The development of approximate solutions
% for such problems starts from three principles.
% (1) Outside a layer, terms involving\/ $\varepsilon$ are negligible.
% (2) Inside a layer, a low order derivative
% is negligible compared with a higher order one.
% (3) Inside a layer, a variable coefficient can be
% approximated by a constant, at least if it is locally nonzero.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>\smallskip\small\parskip=1pt\parindent=0pt
% \newpage
% {\em \underline{Exercise $20.1$}.  Corner layer.}
% {\em (a)}  Plot the solution of $\varepsilon \kern .1pt y''
% + x\kern .1pt y' - y = 0$, $y(-1) = 1$, $y(1) = 1$, $\varepsilon = 0.001$
% (compare Fig.~20.10).
% Which two terms are balanced to determine the
% outer solution, outside the transition region $x\approx 0\kern .5pt$?
% Write down the outer equation and solve it to explain
% the behavior of $y$ in this region.
% {\em (\kern .7pt b)} Which two terms are balanced for $x\approx 0\kern .5pt$?
% Write down the inner equation and show how it can capture the shape
% of $y(x)$ in this region.
% {\em (c)}
% Plot $y$ again, but now on a semilogy scale and zooming
% in with {\tt axis([-.1 .1 1e-4 1])}.
% On the same plot, superimpose curves corresponding to the other
% two terms of the ODE, $\varepsilon\kern.1pt y''$ and $x\kern .1pt y'$.
% For what range of values of $x$ are the two terms that are balanced for
% the outer solution at least 10 times bigger than the omitted term?
% What is the answer to the same question for the inner solution?
% \par
% {\em \underline{Exercise $20.2$}.  Cusp.}
% {\em (a)} If we change the equation of the previous exercise by just
% one coefficient, $\varepsilon \kern .1pt y''
% + x\kern .1pt y' - y/2 = 0$, the solution changes in an interesting way.
% Plot the solution for $\varepsilon = 10^{-3}$.  
% Find the outer equation and explain the shape of the solution.
% {\em (\kern .7pt b)} Use {\tt tic} and {\tt toc} to measure the
% computing time for this problem with $\varepsilon = 10^{-3}$, 
% $10^{-4}$, and $10^{-5}$.  The reason for the slow-down is that
% Chebfun is utilizing a Chebyshev grid of thousands of points, requiring
% the solution of a matrix problem of dimension in the thousands. 
% Change the chebop domain to {\tt L.domain = [-1 0 1]}
% and repeat the same timings; also plot the solution in the case
% $\varepsilon = 10^{-6}$.  Now Chebfun has introduced a breakpoint
% at $x=0$, so that separate Chebyshev grids are utilized on either side, 
% and this makes the matrices smaller.  (Although the user can
% introduce breakpoints like this at fixed locations, Chebfun
% does not offer adaptive gridding capabilities for ODE BVPs.)
% \par
% {\em \underline{Exercise $20.3$}.  Variable oscillations.}
% Solve $\varepsilon \kern .3pt y'' + (x^2-1) y = 0$, $x\in[-2,2\kern .2pt]$
% with $y(-2) = 1$ and $y(2) = 2$ for
% $\varepsilon =10^{-5}$ and plot $y$.
% Explain the nature of this solution and derive an approximate
% formula for the wavelength of the oscillations in the oscillatory
% regions.  (Compare Fig.~7.5.)
% Note that for this problem, principle (1) of our general method of
% boundary layer analysis does not apply because if the
% term involving $\varepsilon$ is dropped, there is just one
% other term remaining, with nothing to balance it.
% \par
% {\em \underline{Exercise $20.4$}.  A double boundary layer.}
% {\em (a)} Solve $\varepsilon^3 \kern .1pt y'' + x^3y' +
% (x^3-\varepsilon) y = 0$, $x\in[\kern .3pt 0,1]$
% with $y(0) = 2$ and $y(1) = 1$ for
% $\varepsilon =0.005$ and plot the resulting solution.
% This problem has two boundary layers at $x=0$, one of width
% $O(\varepsilon)$ and the other of width $O(\varepsilon^{1/2})$.
% Find a way to show this graphically based on solutions computed
% for various values of~$\varepsilon$.
% {\em (\kern .7pt b)}  A good global approximation for $y$ is
% $y(x) \approx
% 2\kern .2pt \exp(-x/\varepsilon) + e\kern .3pt 
% [\kern .3pt \exp(-\varepsilon/2\kern .2pt x^2)
% + \exp(-x) - 1\kern .3pt]$.  What is the maximum error
% of this approximation for $\varepsilon = 10^{-1},$
% $10^{-2}$, and $10^{-3}\kern.5pt$?
% \par
% {\em \underline{Exercise $20.5$}.  Exponential
% ill-conditioning and pseudo-nonuniqueness.}
% {\em (a)}
% Solve $\varepsilon\kern .1pt  y'' - x\kern .1pt y' +y= 1$, $x\in[-1,1]$
% with $y(-1) = 0$ and $y(1) = 0$ for
% $\varepsilon = 1/8, 1/16, 1/32, \dots$ until you get a solution
% that looks completely unlike the others.  Mathematically,
% this solution is not correct.  However, it is a 
% {\em pseudo-solution} in the sense that the residual 
% $\varepsilon\kern .1pt  y'' - x\kern .1pt y' + y-1$
% is very small.  Confirm this by computing
% $\max_{x\in[-1,1]}^{} |\kern .3pt \varepsilon \kern .1pt
% y''-x\kern .1pt y' + y - 1|$ for each of
% your values of\/ $\varepsilon$.
% {\em (\kern .7pt b)}
% Explain this effect by boundary layer analysis
% as follows.  Find a formula for outer solutions to this problem and
% note that it matches the observed incorrect solution.
% Now find formulas for inner
% solutions corresponding to boundary layers at 
% $x\approx \pm 1$ and note that
% these provide enough boundary conditions to match any choice of
% the outer solution.
% Thus although the linear operator $\cal L$ that
% maps right-hand side functions to solutions of this BVP is
% mathematically nonsingular, hence has no null function, it
% has {\em pseudo null functions} with exponentially small residual.
% This ODE BVP is exponentially ill-conditioned, and we may think of it as
% exponentially close to being underdetermined.
% (The reason this has happened is because of the sign change in the
% coefficient of $y'$ in this equation at $x=0$.  This is related
% to the theories of {\em exponential dichotomy} and {\em pseudospectra}.
% See Chapters 10 and 11 of Trefethen and Embree, {\em Spectra and 
% Pseudospectra,} Princeton, 2005.)
% \par
% {\em \underline{Exercise $20.6$}.  Exponential ill-conditioning
% of the adjoint problem.}
% Show that the adjoint of the operator $\cal L$
% of Exercise 20.5 is ${\cal L}^* \kern -3pt: v\mapsto
% \varepsilon\kern .1pt  v'' + x\kern .1pt v' + 2v$ with boundary
% conditions $v(\pm 1) = 0$.
% Solve the problem ${\cal L}^*\kern -1pt v = 1$ numerically
% for $\varepsilon = 1/8, 1/16, 1/32, \dots.$ What goes wrong now when
% when $\varepsilon$ gets sufficiently small?
% </latex>
