%% 22. Time-dependent PDEs
%%
%
% \setcounter{page}{278}
%
%%
%
% Many ODE BVPs that arise in applications
% are related to time-dependent PDE\kern .5pt s.\ \ In the simplest case
% the solution $y(x)$ of the BVP is a {\bf time-independent}
% or {\bf steady solution}
% of the PDE, that is, a solution $y(x,t)$ of the PDE that is
% independent of $t$. Often a solution of
% a PDE may converge to time-independent form as
% $t\to\infty$.
%
%%
%
% \smallskip
% {\em Laplace and heat equations.}
% The archetypical pair of equations related in this way
% are the Laplace equation, an ODE BVP when restricted to one
% space dimension, and the heat equation, a PDE IBVP
% ({\em initial boundary-value problem}). Specifically,
% the {\bf 1D Laplace equation} is the equation
% $y'' = 0$, that is, $d^{\kern .7pt 2}y/dx^2 = 0$. In this chapter,
% since time derivatives will also come into play,
% we will write this in the notation $y_{xx}^{} = 0$, where $y_{xx}^{}$
% is an abbreviation for $\partial^2y/\partial x^2.$\ \ For
% example, here is the 1D Laplace equation with
% homogeneous Dirichlet boundary conditions on $[-\pi,\pi]$:
% $$ y_{xx}^{} = 0, \quad x\in [-\pi,\pi], ~ y(\pm\pi) = 0.
% \eqno (22.1) $$
% We hardly need a computer to see that the solution is
% $y(x) = 0.$
%
%%
%
% Now let us consider how (22.1) may arise in
% describing steady solutions of the {\bf 1D heat equation,}
% $u_t^{} = u_{xx}^{}$. Here
% is an {\bf initial boundary-value problem (IBVP)}
% for this equation with a particular choice of initial data,
% a pulse centered at $x=1$:
% $$ u_t^{} = u_{xx}^{}, \quad x\in [-\pi,\pi], ~ u(\pm\pi,t) = 0, ~
% u(x,0) = \exp(-50(x-1)^4). \eqno (22.2) $$
% Physically, the zero solution just mentioned represents the
% effect that eventually, all the heat
% described by (22.2) will flow out the ends of the
% interval, where the temperature is held at zero.
% A plot of the initial condition together with
% the solution at times $t = 0.01$, $0.1$, and $1$ shows the
% beginning of this process.\footnote{As the experiments of this
% chapter show, Chebfun can be very effective in solving
% certain PDE\kern .5pt s. PDE\kern .5pt s are not Chebfun's main focus, however,
% and its syntax is not as simple nor its capabilities as
% comprehensive for PDE\kern .5pt s as for ODE\kern .5pt s.}
% In this and other figures of this chapter displaying
% PDE solutions at various times~$t$, we plot the curves in
% orange, with the initial curve in brown.
%
u0 = chebfun('exp(-50*(x-1)^4)',[-pi,pi]);
pdefun = @(t,x,u) diff(u,2);
bc.left = @(t,u) u; bc.right = @(t,u) u;
opts = pdeset('plot','off');
t = [0 .01 .1 1]'; [t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0), hold off
title(['Fig.~22.1.~~Solution to (22.2) ' ...
'at $t=0, 0.01, 0.1, 1$'],FS,11), ylim([-.2 1.2])
%%
% \vskip 1.02em
%%
% It is apparent in this figure that until the
% last curve, the boundaries have had little effect. We can
% confirm this by integrating the curves to determine the
% total heat at each time, which doesn't diminish
% much until $t = 1$.
disp(' t heat'), disp([t sum(u)'])
%%
%
% \noindent
% As $t$ increases further, the left boundary begins to be
% important as well as
% the right one, as is clear from these images for
% $t = 1$ (again) and $2,4,8,16.$
%
u = u(:,end);
t = [1 2 4 8 16]'; [t,u] = pde15s(pdefun,t,u,bc,opts);
plot(u,CO,ibvp), hold on, plot(u(:,1),CO,ibvp0), hold off
title('Fig.~22.2.~~Solution continued to $t=1,2,4,8,16$',FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% Note how the pulse is now moving to the center, approaching its
% asymptotic form $C \exp(-t/4) \cos(x/2)$, as determined
% originally by Joseph Fourier around 1807.
% The integrals confirm that the total heat content is
% now diminishing speedily, with $98\%$ of the heat
% gone by $t=16$.
%
disp(' t heat 0.76exp(-t/4)')
disp([t sum(u)' 0.76*exp(-t/4)])
%%
%
% The Laplace and heat equations may be the archetypes, but
% these simple linear equations
% give only a hint of the rich relationship between
% time-dependent PDE\kern .5pt s and associated ODE BVPs.
% We will now look at some further
% examples illustrating the phenomena of stability,
% instability, and symmetry-breaking bifurcation.
%
%%
%
% \smallskip
% {\em Stability and instability.}\ \ Let us add
% an exponential term to (22.1) and (22.2) to make them
% into ODE and PDE versions of the nonlinear
% {\em Bratu equation}, considered
% in Chapters~16 and~18.\ \ Here is (16.4) again, with the constant
% on the nonlinear term
% adjusted to compensate for the change in domain from
% $[\kern .3pt 0,1]$ to $[-\pi,\pi]$.
% $$ y_{xx}^{} + {3e^y\over 4\pi^2} = 0 ,
% \quad x\in [-\pi,\pi], ~ y(\pm\pi) = 0 . \eqno (22.3) $$
% As shown in Chapter 16, there are two solutions to
% this ODE BVP.\ \ We plot one of them as a dashed line, the
% one we shall find is unstable.
%
N = chebop(-pi,pi); N.op = @(x,y) diff(y,2) + (.75/pi^2)*exp(y);
N.lbc = 0; N.rbc = 0; y1 = N\0; plot(y1,CO,bvpnl)
hold on, x = chebfun('x',[-pi pi]); N.init = .2*(pi^2-x^2);
y2 = N\0; plot(y2,'--',CO,bvpnl), ylim([-.5 2.5]), hold off
title(['Fig.~22.3.~~Two solutions of the ' ...
'Bratu equation (22.3)'],FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% Now let us embed the problem in a time-dependent PDE with
% initial condition $u(x,0) = u_0^{}(x)$,
% $$ u_t^{} = u_{xx}^{} + {3e^u\over 4\pi^2},
% \quad x\in [-\pi,\pi], ~ u(\pm\pi,t) = 0,~u(x,0) = u_0^{}(x) . \eqno (22.4) $$
% First of all, suppose we take $u_0^{}(x) = 0$. Then as $t\to\infty$,
% the solution increases and
% approaches a steady state, the lower solution of (22.3) just plotted.
%
u0 = 0*x; pdefun = @(t,x,u) diff(u,2)+(.75/pi^2)*exp(u);
t = (0:2:12)'; [t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0)
plot(y1,CO,bvpnl), plot(y2,'--',CO,bvpnl), hold off, ylim([-.5 2.5])
title('Fig.~22.4.~~Solution to (22.4) at $t=0, 2,\dots,12$',FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% Many other initial conditions also lead to the same steady
% state, which is a stable solution of this PDE.\ \ For example,
% in this plot we choose $u_0^{}$ to have
% a few wiggles and include extra values $t = 0.05, 0.2$
% to show how transient their influence is.
%
u0 = sin(x)+.3*sin(3*x)+.1*sin(13*x); t = [0 0.05 0.2 2:12]';
[t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0), ylim([-1.5 2.5])
plot(y1,CO,bvpnl), plot(y2,'--',CO,bvpnl), hold off
title(['Fig.22.5.~~Solution to (22.4) at ' ...
'$t=0,.05,.2$ and $2,4,\dots,20$'],FS,11)
%%
% \vskip 1.02em
%%
%
% As another example, next
% we take $u_0^{}$ to be $0.9$ times the {\em upper\/}
% solution of the BVP (22.3). Despite that starting point,
% the PDE converges to the other, lower solution. The
% upper solution of the BVP is an unstable solution of the PDE.
%
u0 = 0.9*y2; t = (0:4:32)'; [t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0), ylim([-.5 2.5])
plot(y1,CO,bvpnl), plot(y2,'--',CO,bvpnl), hold off
title(['Fig.~22.6.~~Solution to (22.4) at ' ...
'$t=0, 4,\dots,32$'],FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% If we take $u_0^{}$ to be a little larger than the upper
% solution of the BVP, the same instability is manifested as
% divergence to $\infty$ in finite time, at $t\approx 6.914$
% for this particular initial condition.
%
u0 = 1.1*y2; t = [0 4 6 6.5 6.8 6.9 6.913]';
[t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0)
plot(y1,CO,bvpnl), plot(y2,'--',CO,bvpnl), hold off
title(['Fig.~22.7.~~Solution to (22.4) at ' ...
'$t=0,4,6,6.5,6.8,6.9,6.913$'],FS,11)
%%
% \vskip 1.02em
%%
%
% As we have discussed at many points in this book, phenomena
% of stability and instability are of pervasive importance in
% the mathematical sciences. Generally speaking, a stable
% solution is likely to be observable in an experiment and
% an unstable one is not. Also, generally speaking, a stable
% solution is likely to emerge as a steady state of a
% time-dependent process and an unstable one is not.
% Still, a study of unstable solutions may be
% important to a full understanding of a problem, and they
% can often be determined by solving time-independent BVPs.
%
%%
%
% \smallskip
% {\em Symmetry-breaking bifurcation.}
% We have seen in Chapters 16--18 that the number of solutions
% of an ODE BVP may change as a parameter varies.
% Often when new solutions appear, this corresponds
% to a loss of stability of an old solution, and
% the new solutions are less symmetric
% than the original one.\footnote{Symmetry-breaking
% bifurcations are at the root of {\em phase transitions} in
% physics. Physicists believe that some particularly important
% phase transitions occurred soon after the Big Bang, when a
% single unified force separated into the different fundamental
% forces recognized today such as the electromagnetic and weak forces.}
%
%%
%
% To illustrate this effect, let us look at the
% time-dependent {\bf Allen--Cahn equation},
% $$ u_t^{} = \varepsilon\kern .3pt u_{xx}^{} + u - u^3, \quad
% x\in [-1,1], ~ u(\pm 1,t) = 0, ~u(x,0) = u_0^{}(x). \eqno (22.5) $$
% For $\varepsilon > \varepsilon_1^{} \approx 0.406$,
% the unique steady solution
% of this problem is $u\equiv 0$, as we illustrate here with
% $\varepsilon = 0.6$.
% The associated ODE BVP, which appeared in Exercise~18.1, is
% $$ \varepsilon\kern .3pt y_{xx}^{} + y - y^3 = 0, \quad
% x\in [-1,1], ~ y(\pm 1) = 0, \eqno (22.6) $$
% and we show its (zero) solution superimposed on the plot.
%
x = chebfun('x'); ep = 0.5; pdefun = @(t,x,u) ep*diff(u,2)+u-u.^3;
u0 = 0.5*cos(pi*x/2)+.7*sin(pi*x); t = (.5:.5:4)';
[t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0), ylim([-1.3 1.3])
title(['Fig.~22.8.~~Convergence to 0 for Allen--Cahn ' ...
'eq.~(22.5) with $\varepsilon=0.6$'],FS,11)
y1 = 0*x; plot(y1,CO,bvpnl), hold off
%%
% \vskip 1.02em
%%
%
% Reducing $\varepsilon$ to $0.15$, on the other hand, illustrates
% a different situation. Here, the zero solution is still
% valid, but it is unstable. With the same initial condition
% as before, the curve approaches a positive form as $t\to\infty$.
% By symmetry, if the initial condition were negated, the
% steady solution as $t\to\infty$ would be negated too.
%
ep = 0.15; pdefun = @(t,x,u) ep*diff(u,2)+u-u.^3;
[t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0), ylim([-1.3 1.3])
title(['Fig.~22.9.~~New solutions and broken symmetry ' ...
'with $\varepsilon=0.15$'],FS,11)
N = chebop(-1,1); N.lbc = 0; N.rbc = 0;
N.op = @(x,y) ep*diff(y,2) + y - y^3;
N.init = u0; y2 = N\0; y3 = -y2;
plot(y1,CO,bvpnl,'--'), plot([y2 y3],CO,bvpnl), hold off
%%
% \vskip 1.02em
%%
% What if one starts from an initial condition that is an odd function
% of $x$,
% hence equally close to the upper and lower solutions? For this
% value of $\varepsilon$, the solution quickly converges to zero.
u0 = .7*sin(pi*x); [t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0), ylim([-1.3 1.3])
title(['Fig.~22.10.~~Odd initial function: ' ...
'convergence to a saddle point'],FS,11)
plot(y1,CO,bvpnl,'--'), plot([y2 y3],CO,bvpnl), hold off
%%
% \vskip 1.02em
%%
%
% \noindent
% The zero solution is not truly stable, however. It is a saddle
% point, and if it is perturbed by a function that is not odd, the
% perturbations will generally grow.
%
%%
% If $\varepsilon$ is reduced further for this problem,
% below another bifurcation point $\varepsilon_2^{}\approx 0.1$, a further
% set of solutions appears.
ep = 0.02;
pdefun = @(t,x,u) ep*diff(u,2)+u-u.^3; u0 = .3*sin(pi*x);
[t,u] = pde15s(pdefun,t,u0,bc,opts);
plot(u,CO,ibvp), hold on, plot(u0,CO,ibvp0), ylim([-1.3 1.3])
title(['Fig.~22.11.~~Two more unstable steady states ' ...
'with $\varepsilon=0.02$'],FS,11)
N.op = @(x,y) ep*diff(y,2) + y - y^3;
N.init = cos(pi*x/2); y2 = N\0; y3 = -y2;
N.init = sin(pi*x); y4 = N\0; y5 = -y4;
plot(y1,CO,bvpnl,'--'), plot([y2 y3],CO,bvpnl)
plot([y4 y5],'--',CO,bvpnl), hold off
%%
% \vskip 1.02em
%%
%
% \noindent
% We have marked the new solutions with dashed curves, because
% in fact these again are saddle points, not stable points.
% If the initial condition were not an odd function, the solution
% would eventually converge to one of the stable steady states
% of a single sign.
%
%%
%
% \smallskip
% {\em Time-harmonic solutions and resonance.}
% Up to now we have looked at PDE\kern .5pt s of first order in $t$.\ \ An
% ODE BVP may also be relevant to time-oscillatory
% solutions of a PDE that is of second order in $t$.
% (We considered such problems in Exercise~6.3 and
% in the Application of Chapter~20,
% ``Why is New York hotter than San Francisco?'')
% Suppose we have a linear PDE
% $$ u_{tt}^{} = Lu + e^{i\kern .3pt\omega\kern .3pt t} f , \eqno (22.7) $$
% where $L$ is a linear operator acting on functions of $x$,
% $\omega$ is frequency,
% and $f$ is a function of $x$.
% A {\bf time-harmonic} solution of (22.7) is a solution
% of the form
% $$ u(x,t) = e^{i\kern .3pt\omega\kern .3pt t} y(x) \eqno (22.8) $$
% for some function $y$. Inserting
% (22.8) in (22.7) gives
% $$ -\omega^2 u = Lu + e^{i\kern .3pt\omega\kern .3pt t} f , $$
% or after dividing by $e^{i\kern .3pt\omega\kern .3pt t}$,
% $$ -\omega^2 y = Ly + f , \eqno (22.9) $$
% known as the {\bf reduced equation} associated with (22.7).
% One of the fundamental techniques of mathematical physics
% is to solve reduced equations to find time-harmonic solutions
% of time-dependent equations.\footnote{In quantum mechanics,
% the technique is applied to the Schr\"odinger equation,
% which differs from (22.7) in being of first order in $t$. We saw
% this in the Application of Chapter~6, ``Eigenstates of
% the Schr\"odinger equation.''}
%
%%
%
% \smallskip
% {\em Further links between ODE\kern .5pt s and PDE\kern .5pt s.}
% The difference between ODE\kern .5pt s and PDE\kern .5pt s is nothing more than
% the number of independent variables, so it is hardly surprising
% that connections between the two arise in many ways. In this
% chapter we have emphasized steady-state forms that a solution
% of a PDE may settle down to as $t\to\infty$.
% A related situation is that
% sometimes as $t\to\infty$ the solution of a PDE approaches not
% a fixed function but a {\em traveling wave,} a function
% $u(x,t) = y(x-c\kern .7pt t)$ for some constant wave velocity $c$
% that in most cases is not known a priori. Here again ODE\kern .5pt s play
% a natural role. Another way ODE\kern .5pt s arise from PDE\kern .5pt s is in
% {\em separation of variables,} and the time-harmonic solutions
% just looked at fall in this category.
% Separation of variables applies also in purely spatial
% problems, for example, if a radially-symmetric elliptic PDE is
% factored into an ODE in $r$ and an ODE in $\theta$\kern .5pt: this is the
% most familiar way in which one encounters the Bessel equation (Exercise~22.1).
% Another PDE-ODE link arises for problems that have {\em similarity
% solutions.} And if we move up one or more dimensions, it is worth
% noting that just as a PDE in one space dimension may reduce to
% an ODE in
% some limit, so a PDE in two or more space dimensions may reduce
% to a PDE in one dimension less. For example, steady solutions of
% the heat equation in a 2D spatial domain will be described by solutions
% of the Laplace equation, a PDE, in the same domain.
%
%%
%
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: solitons and the KdV equation} \\[-3pt]
% \hrulefill
% \end{center}
%
%%
%
% The simplest of all time-dependent PDE\kern .5pt s is the {\em first-order
% 1D linear wave equation} (or {\em advection equation}),
% $$ u_t^{} = -u_x^{}, \eqno (22.10) $$
% whose solutions consist of an arbitrary wave form $y(x)$ traveling
% rightward at speed~$1$ without changing shape,
% $$ u(x,t) = y(x-t). \eqno (22.11) $$
% For example, here we start with a Gaussian of amplitude
% $50$ centered at $x=-1$ and run to $t=2$ using
% the Chebfun {\tt spin} command (``stiff PDE integrator'').
% The wave slides along to $x=1$ without changing shape.
%
dom = [-5,5]; x = chebfun('x',dom);
u0 = 50*exp(-30*(x+1)^2); S = spinop(dom,[0 2]);
S.lin = @(u) -diff(u);
S.init = u0; u1 = spin(S,160,.01,'plot','off');
plot(u0,'r',u1,'k',LW,1), ylim([-20 130])
text(-1,57,'$t=0$',FS,9,HA,CT), text(1,57,'$t=2$',FS,9,HA,CT)
title(['Fig.~22.12.~~Advection with the linear ' ...
'wave equation (22.10)'],FS,11)
%%
% \vskip 1.02em
%%
%
% For a nonlinear equation, on the other
% hand, change of shape will be the
% rule. In this application we look at
% the most famous of all nonlinear wave
% equations, the {\em KdV} ({\em Korteweg--de Vries\/}) equation,
% which we write in the form
% $$ u_t^{} = -0.03(u^2)_x^{} - 0.01 u_{xxx}^{}. \eqno (22.12) $$
% The equation blends first-order nonlinear {\em advection}, involving
% $(u^2)_x^{}$ or equivalently $2\kern .3pt uu_x^{}$, with third-order
% linear {\em dispersion}, involving $u_{xxx}^{}$.
% The constants $0.03$ and $0.01$ have no great significance
% and have been chosen to make the behavior approximately match
% that of $u_t^{} = -u_x^{}$ for our initial condition.
% In the next figure we solve (22.12) to $t=2$ as before with {\tt spin},
% which is made for nonlinear problems like this.
%
S = spinop(dom,[0 2]);
S.lin = @(u) -0.01*diff(u,3); S.nonlin = @(u) -0.03*diff(u.^2);
S.init = u0; u1 = spin(S,160,.01,'plot','off');
plot(u0,'r',u1,'k',LW,1), ylim([-20 130])
text(-1,57,'$t=0$',FS,9,HA,CT), text(.78,50,'$t=2$',FS,9,HA,CT)
title(['Fig.~22.13.~~Nonlinear propagation with ' ...
'the KdV equation (22.12)'],FS,11)
%%
% \vskip 1.02em
%%
%
% To first approximation, Figure~22.13 is like Figure~22.12.
% However, the peak is moving a little more slowly, and it is losing
% amplitude. Meanwhile a tail of dispersive
% oscillations is appearing behind. (The small oscillations
% ahead of the main pulse are due to the use of periodic
% boundary conditions in all {\tt spin} simulations.) A complicated mix of effects
% like this is what one must expect to see, in general, with a nonlinear PDE.
%
%%
%
% But the KdV equation has special properties that
% this experiment does not reveal. It has special solutions,
% called {\em solitary waves\/} or {\em solitons},
% that travel at fixed speed with exactly
% uniform shape. For any constant $c>0$, the associated soliton
% profile is given by the formula
% $$ y(x) = 50\kern .7pt c\kern 1.3pt \hbox{sech}(\sqrt{100\kern .3pt c}
% \kern .8pt (x+1)/2)^2 , \eqno (22.13) $$
% where $\hbox{sech}(x) = 1/\cosh(x)$, the reciprocal of the hyperbolic
% cosine.
% Moreover, the speed at which the soliton moves
% is equal to $c$, which means that tall solitons move
% faster than short ones. Here we see this behavior with
% a simulation for $t\in [\kern .3pt 0,2\kern .2pt ]$ and $c=1$. Note that the
% $\hbox{sech}^2$ function is much like a Gaussian, the difference
% in shape being too little to notice by eye.
%
c = 1; soliton1 = 50*c*sech(sqrt(100*c)*(x+1)/2)^2;
S.init = soliton1; u1 = spin(S,160,.01,'plot','off');
plot(S.init,'r',u1,'k',LW,1), ylim([-20 130])
text(-1,57,'$t=0$',FS,9,HA,CT), text(1,57,'$t=2$',FS,9,HA,CT)
title(['Fig.~22.14.~~If the initial condition is a soliton, ' ...
'it translates unchanged'],FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% Since we took $c=1$, this soliton has traveled at speed $1$.
% If we double the height, we double the speed.
%
c = 2; soliton2 = 50*c*sech(sqrt(100*c)*(x+2)/2)^2;
S.init = soliton2; u1 = spin(S,160,.005,'plot','off');
plot(S.init,'r',u1,'k',LW,1), ylim([-20 130])
text(-2,107,'$t=0$',FS,9,HA,CT), text(2,107,'$t=2$',FS,9,HA,CT)
title('Fig.~22.15.~~A taller soliton moves faster',FS,11)
%%
% \vskip 1.02em
%%
%
% Another property of the KdV equation is truly extraordinary.
% {\em Solitons can pass through one another.} For example, if
% we start with a solution consisting of two well separated solitons,
% the rear one taller than the front one, then the rear
% one will catch up, pull
% ahead, and eventually return to the same shape and
% speed it had before. Here is an illustration.
%
S.init = soliton1 + soliton2; u1 = spin(S,160,.005,'plot','off');
plot(S.init,'r',u1,'k',LW,1), ylim([-20 130])
text(-1.3,64,'$t=0$',FS,9,HA,CT), text(1.4,64,'$t=2$',FS,9,HA,CT)
title('Fig.~22.16.~~A fast soliton passes through a slower one',FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% Note the locations of the peaks in this figure
% for the curve at $t=2$. The speeds have returned to
% essentially their initial values 2 and 1, but
% the tall soliton is a bit further
% along than $x=2$, and the short soliton is not as far along
% as $x=1$. This is because a shift of each soliton is
% introduced during the time of their interaction around $t=1$.
% The details of these nonlinear interactions are subtle and fascinating,
% as the reader can easily explore with further computations.
% A waterfall plot gives the idea.
%
spin(S,160,.005,'plot','waterfall'); set(gca,FS,9)
title('Fig.~22.17.~~Two solitons',FS,11), colormap([0 0 0])
xlabel('$x$',FS,10), ylabel('$t$',FS,10), zlabel(' ')
%%
% \vskip 1.02em
%%
%
% A soliton is an example of a traveling wave, which we
% defined earlier as a solution of
% a time-dependent PDE of the form
% $$ u(x,t) = y(x-c\kern .7pt t) \eqno (22.14) $$
% for some function $y(x)$.
% Of course we know that in the case of the KdV
% equation, the formula for $y$ must be (22.13), but
% let us suppose we do not know this.
% Inserting (22.14) in (22.12) gives
% $-c\kern .7pt y_x^{} = -0.03 \kern .5pt (y^2)_x^{} - 0.01 y_{xxx}^{},$
% that is,
% $$ y_{xxx}^{} = 100\kern .7pt c
% \kern .7pt y_x^{} -3\kern .5pt (y^2)_x^{}. $$
% This looks like a nonlinear 3rd-order ODE.\ \ However,
% all three terms are differentiated, so we
% can integrate it once to get the second-order equation
% $$ y_{xx}^{} = 100\kern .7pt c\kern .7pt y-3\kern .5pt y^2 + C $$
% for some constant $C$. Now let us use this ODE to find the
% shape of a traveling wave solution. If we assume that $y$ and
% its derivatives approach 0 as $|x|\to \infty$, then we have $C = 0$,
% and the equation becomes
% $$ y_{xx}^{} = 100\kern .7pt c\kern .7pt y-3\kern .5pt y^2.
% \eqno (22.15) $$
% This could be dealt with analytically, or one
% could try various things on the computer. For example, suppose
% we regard (22.15) as an IVP on the interval $x\in [-2,2\kern .2pt ]$, with
% $c=1$. With initial conditions $y(-2)=0$ and
% $y'(-2)= 10^{-4}$ we get a soliton of just the expected amplitude
% $50$.\ \ The horizontal location has no particular significance; it would be
% different with a different value of $y'(-2)$.
%
N = chebop(-2,2); N.op = @(y) diff(y,2) - 100*y + 3*y^2;
N.lbc = [0; 0.0001]; y = N\0;
plot(y,CO,ivpnl), axis([-5 5 -20 130])
title(['Fig.~22.18.~~Soliton travelling wave profile ' ...
'determined from an ODE'],FS,11)
%%
% \vskip 1.02em
%%
%
% We have not mentioned the scientific uses of the KdV equation.
% These began with water waves, and
% the story is often told of John Scott Russell's chase on horseback of a
% solitary wave in
% the Union Canal in Scotland in 1834. Two centuries later, the
% important applications are in optics. Light pulses travelling through
% optical fibers experience nonlinear effects combined with dispersion, and
% the KdV equation and its relatives are fundamental nowadays in
% optical technology.
%
%%
%
% \smallskip
% {\sc History.} There are few equations with more impact in
% the history of mathematics than the KdV equation.
% The understanding of solitons and their interactions came in the 1960s
% and led to sensational new discoveries related to so-called
% completely integrable dynamical systems; some of the key names
% were Zabusky, Kruskal, and Lax.
% The technical details are exhilarating, and so are the wider
% implications, for it was work around the KdV equation that led to
% some of our deepest understanding of the significance of nonlinearity.
% All this was sparked by computers. Though the KdV equation goes
% back to the 19th century, the great advances began
% with a 1955 paper by Fermi, Pasta, and Ulam,
% who worked on the Maniac I computer at Los Alamos.
% \smallskip
%
%%
%
% {\sc Our favorite reference.}
% {\em The PDE Coffee Table Book} was a group
% project at Oxford in the late 1990s coordinated by
% Trefethen together with Kristine Embree.
% The idea was to present 100 PDE\kern .5pt s in the simplest
% possible terms, each getting exactly a two-page spread with
% just the right information about mathematics, applications, and history,
% and always with a colorful computed illustration. Alas, the
% project halted when it was only $34\%$ complete, but that
% $34\%$ is worth checking out: see
% {\em The (Unfinished) PDE Coffee Table Book,}
% \verb|https://people.maths.ox.ac.uk/trefethen/pdectb.html|.
% \smallskip
%
%%
%
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 22.} Many ODE BVPs are
% related to associated time-dependent PDE\kern .5pt s in one space dimension.
% The ODE\kern .5pt s may arise as limits\/ $t\to\infty$ in which solutions
% approach a steady state, a steady oscillation, or a traveling wave.
% ODE\kern .5pt s may also arise from purely spatial PDE\kern .5pt s by the process
% of separation of variables.
% \vspace{2pt}}}
% \end{displaymath}
%
%%
% \smallskip\small\parskip=1pt\parindent=0pt
% {\em Exercise $22.1$. Bessel equation for a circular drum.}
% A 2D membrane with a fixed boundary $\Gamma$ oscillates according
% to the wave equation $u_{tt}^{} = \Delta u$, where $\Delta$ is the
% Laplacian operator, $\Delta u = u_{xx}^{} + u_{yy}^{}$. This is a PDE
% in two space variables and time.
% {\em (a)}~Assume that $u$ takes the form $u(x,y,t) =
% \exp(i\kern .5pt \omega\kern .3pt t) v(x,y)$ for a fixed frequency~$\omega$ and function
% $v$. Show that this assumption
% reduces the equation to a PDE eigenvalue equation in two space variables.
% {\em (\kern .7pt b)} Now assume that $\Gamma$ is the unit circle and
% that~$v$ takes the polar coordinates
% form $v(r,\theta) = \exp(i k \theta) y(r)$ for
% a fixed integer wave number~$k$. Using the polar coordinates
% representation of the Laplacian,
% $\Delta v = v_{rr}^{} + r^{-1} v_r^{} + r^{-2} v_{\theta\theta}^{}$,
% show that this reduces the problem
% to an ODE eigenvalue problem involving Bessel's equation.
% \par
% {\em \underline{Exercise $22.2$}.
% Asymptotic behavior of advection-diffusion problem.}
% {\em (a)} Execute {\tt chebgui} and run the demo labeled
% Advection-diffusion equation 2 under the tab Demos/PDE-scalar.
% The IBVP here is $u_t^{} = 0.3\kern .2pt u_{xx}^{} + 10\kern .3pt u_x^{}$ for
% $x\in[-1,1]$ with initial condition $u(x,0) = \exp(-10x^4/(1-x^2))$
% and boundary conditions $u(\pm 1) = 0$.
% Describe in words how the solution behaves as $t\to\infty$.
% {\em (\kern .7pt b)} This behavior can be analyzed by looking
% for solutions of this equation of the form $u(x,t) = e^{\lambda t}
% y(x)$, where~$y$ is a fixed function of $x$. Write down the
% formulas showing that $y$ and $\lambda$ are solutions of an
% eigenvalue problem. Determine $y$ and $\lambda$, either
% analytically or numerically, for the value of $\lambda$ corresponding to
% the slowest decay of $u(x,t)$ as $t\to\infty$. Show
% that~$y$ has a shape corresponding to the observations of part {\em (a)}.
% {\em (c)} By pressing the ``Export solution'' button in Chebgui and
% then typing \verb|max(u(:,end))|, find the value of
% $\max_{x\in[-1,1]} u(x,0.25)$ (since the final time of this simulation
% is $t=0.25$). Combine this figure with the result of part {\em (\kern .7pt b)}
% to estimate $\max_{x\in[-1,1]} u(x,5)$ to an accuracy of within
% $1\%$. {\em (d)} The asymptotic rate of decay is different if the advection
% term $10\kern .3pt u_x^{}$ is removed from the equation. (Physically,
% heat transfer is faster with convection than just conduction --- the
% oceanic effect that we found makes San Francisco cooler than New York in
% the summer.)
% Estimate $\max_{x\in[-1,1]} u(x,5)$ in this case.
% \par
% {\em \underline{Exercise $22.3$}. Metastability of Allen--Cahn equation.}
% {\em (a)} Execute {\tt chebgui} or {\tt pde15s}, as you
% prefer, for the Allen--Cahn equation
% $u_t^{} = 0.05\kern .3pt u_{xx}^{} + u = u^3$ for $x\in [-5,5]$ and
% $t\in [\kern .3pt 0,50\kern .3pt ]$ with boundary
% conditions $u(-5) = u(5) = -1$ and initial condition
% $u(x,0) = -1+2\exp(-x^2)$. Show how the initial function quickly changes
% to a ``stalagmite'' form that looks like a steady state.
% {\em (\kern .7pt b)} Now run the calculation for
% $t\in [\kern .3pt 0,1000\kern .3pt ]$ and show that the apparent steady
% state was not truly steady, only metastable. Estimate the time at
% which the stalagmite is extinguished --- say, the time $t$ at which
% $\max_{x\in[-5,5]} u(x,t)$ becomes negative. How does this critical time change
% if the diffusion coefficient $0.5$ is increased to $0.6$ or $0.7\kern .3pt$?
% (With values of the diffusion coefficient much less than $0.5$ the
% time soon becomes nearly infinite in practice --- too long to be
% measured by a Chebfun computation --- though it is always finite in
% principle.)
% {\em (c)}~Write down the ODE in the variable $y$ for a time-independent
% solution of the Allen--Cahn equation. One solution of the ODE is
% $y(x) = -1$. The waveform observed for $t=50$ in part {\em (a)}
% is a ``pseudo-solution'': it nearly satisfies the ODE, though not
% exactly. Insert this
% solution in the ODE to quantify ``nearly.''
%