%% 18. Continuation and path-following
%
% \setcounter{page}{229}
%
%%
% As we have seen in the last chapter, a central
% idea in the study of ODE\kern .5pt s is the investigation of the dependence
% of solutions on parameters. This is important both
% scientifically, since varying parameters may reveal the range
% of behaviors of a problem,
% and algorithmically, since solutions
% can often be calculated most reliably
% by varying a parameter incrementally
% from one range to another. In this chapter we pursue
% these ideas.
%%
%
% Dependence on parameters can be of interest for linear problems,
% as we shall see in our investigation of boundary and interior
% layers in Chapter 20. The heart of this subject, however,
% is nonlinear problems, where the varying of parameters may be
% essential to the computation of some solutions and to the elucidation
% of the structure of a problem.
%
%%
%
% As a starting example, suppose we are faced with the problem
% $$ \varepsilon\kern .3pt y'' +y- y^2 = 0, \quad x\in [\kern .3pt 0,1], ~
% y(0) = 1, ~y(1) = 0 \eqno (18.1) $$
% and seek a solution with $\varepsilon = 2^{-10}$.\ \ If
% we try calling Chebfun directly,
% it fails, as Newton iterations often do when started from an initial
% guess not close enough to a solution. Larger values
% of $\varepsilon$ lead to ready convergence, however.
% Taking advantage of this
% effect, we can creep up on the solution we want,
% decreasing~$\kern 1pt\varepsilon\kern 1pt$ steadily in
% the process known as {\bf continuation}.\footnote{Other terms
% include {\bf embedding} and {\bf homotopy} methods. Such ideas are
% important not just for ODE\kern .5pt s (and PDE\kern .5pt s), but also for finding solutions of
% nonlinear algebraic systems of equations.}
% The number
% $\varepsilon$ is called the {\bf continuation parameter}.
% The following code starts with
% $\varepsilon = 2^{-3}$ and reduces the value successively
% by factors of $\sqrt 2$ until
% $2^{-10}$ is reached.
%
ODEformats
N = chebop(0,1); N.lbc = 1; N.rbc = 0; epsvec = 2.^-(3:.5:10);
for ep = epsvec
N.op = @(x,y) ep*diff(y,2) + y - y^2;
y = N\0; plot(y,CO,bvpnl,LW,0.6), hold on, ylim([-1.5 1.5])
N.init = y;
end
title(['Fig.~18.1.~~A family of solutions of (18.1) found' ...
' by continuation in $\varepsilon$'],FS,11), hold off
%%
% \vskip 1.02em
%%
%
% The key feature of this code segment is that for each value
% of $\kern 1pt\varepsilon\kern 1pt$ after the first,
% an initial guess is specified
% in the \verb|N.init| field corresponding to the converged
% solution from the previous value.
% Note that the solution obtained shows a
% boundary layer at $x=1$; its width
% is $O(\varepsilon^{1/2})$ (see Exercise 18.2).
%
%%
%
% The family of solutions of (18.1) we have just plotted is
% not the only one. If we start from the initial
% guess $y(x) = 1 - x - 3e^x (x-x^2)$, another family emerges,
% which we plot in a new color. Again there is a boundary
% layer at $x=1$.
% Images like these show that continuation may offer a powerful
% tool for investigating problems with multiple solutions.
%
x = chebfun('x',[0,1]); N.init = 1-x-3*exp(x)*(x-x^2);
for ep = epsvec
N.op = @(x,y) ep*diff(y,2) - y^2 + y;
y = N\0; plot(y,CO,purple,LW,0.6), hold on, ylim([-1.5 1.5])
N.init = y;
end
hold off
title('Fig.~18.2.~~Another family of solutions of (18.1)',FS,11)
%%
% \vskip 1.02em
%%
% Next, let us turn to equation (16.5),
% $$ \varepsilon\kern .3pt y'' + y + y^2 = 1, \quad x\in [-1,1], ~
% y(\pm 1) = 0. \eqno (18.2) $$
% Suppose we want a solution with $\varepsilon = 2^{-10}$. This time, a call
% to Chebfun is successful:
N = chebop(-1,1); N.lbc = 0; N.rbc = 0; ep = 2^-10;
N.op = @(x,y) ep*diff(y,2) + y + y^2;
y = N\1; plot(y,CO,bvpnl), ylim([-2.5 2.5])
title(['Fig.~18.3.~~A solution of (18.2) ' ...
'with $\varepsilon=2^{-10}$'],FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% This solution is valid, but it is far from the
% simplest solution for this choice of $\varepsilon$.
% To find another, the following code segment uses continuation
% with $\varepsilon = 1/2, 1/4, \dots, 1/1024.$
% Here we obtain solutions with boundary layers
% of width $O(\varepsilon^{1/2})$ at both ends.
%
N = chebop(-1,1); N.lbc = 0; N.rbc = 0; epsvec = 2.^-(1:10);
for ep = epsvec
N.op = @(x,y) ep*diff(y,2) + y + y^2;
y = N\1; plot(y,CO,bvpnl,LW,0.6), hold on, ylim([-2.5 2.5])
N.init = y;
end
title(['Fig.~18.4.~~A family of solutions of (18.2) found by' ...
' continuation in $\varepsilon$'],FS,11), hold off
%%
% \vskip 1.02em
%%
%
% Using a continuation parameter may have a number of benefits, including these:
% \par\vskip 9pt
% 1. The computation may be faster because fewer Newton steps are needed.
% \par\vskip 3pt
% 2. Convergence may be achieved in cases where
% a cold start would fail.\footnote{When we
% use an initial condition from a nearby problem
% to improve our results, this is called a {\em warm start}.}
% \par\vskip 3pt
% 3. In problems with multiple solutions,
% continuation may help pick out the solution of interest.
% \par\vskip 9pt
%
%%
%
% \noindent
% In the experiment we have just shown, benefit (3) has been the crucial one.
%
%%
% For a third example, let us look at the Bratu equation,
% $$ y'' + \lambda \kern -.7pt \exp(y) = 0, \quad
% x\in [\kern .3pt 0,1], ~y(0) = y(1) = 0 . \eqno (18.3) $$
% In (16.4), $\lambda$ took the value
% $3$. Now, suppose we want it to range over the values
% $0, 0.25, \dots, 3.5$. The continuation approach gives
% a simple family of curves.
N = chebop(0,1); N.lbc = 0; N.rbc = 0;
lamvec1 = 0:.25:3.5; mvec1 = [];
for lam = lamvec1
N.op = @(x,y) diff(y,2) + lam*exp(y);
y1 = N\0; plot(y1,CO,bvpnl,LW,0.6), hold on
N.init = y1; mvec1 = [mvec1 y1(0.5)];
end
title(['Fig.~18.5.~~Solutions to the Bratu eq.\ (18.3)'...
' for various $\lambda$'],FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% Starting from another initial guess, we can find the
% other branch of solutions, which we superimpose on
% the same plot with a new vertical scale.
%
N = chebop(0,1); N.lbc = 0; N.rbc = 0;
lamvec2 = 0.5:.25:3.5; mvec2 = [];
x = chebfun('x',[0 1]); N.init = 6*sin(pi*x);
for lam = lamvec2
N.op = @(x,y) diff(y,2) + lam*exp(y);
y2 = N\0; plot(y2,CO,purple,LW,0.6)
N.init = y2; mvec2 = [mvec2 y2(0.5)];
end
title('Fig.~18.6.~~Another set of solutions to (18.3)',FS,11), hold off
%%
% \vskip 1.02em
%%
%
% At this point we have shown five examples of families of solutions
% to parameter-dependent nonlinear ODE\kern .5pt s.
% In each case we superimposed
% collections of solutions on a single plot. There is another approach
% to such problems that is often fruitful, however, and that is to focus
% on a single scalar quantity that is characteristic of the
% solution of interest, an idea we have already explored
% in the last chapter on bifurcation.
%
%%
% For the Bratu problem, a natural scalar to measure is
% the maximum of $y(x)$, that is, $y(0.5)$. This quantity was stored at
% each step of the iterations above
% in the vectors |mvec1| and |mvec2|. Thus we can immediately plot
% $y(0.5)$ as a function of $\lambda$ for the two families
% of solutions just computed, giving a bifurcation diagram for (18.3).
plot(lamvec1,mvec1,'.-',CO,bvpnl,MS,10), axis([0 4 0 7]), hold on
plot(lamvec2,mvec2,'.-',CO,purple,MS,10), hold off
xlabel('$\lambda$',FS,10,IN,LT), ylabel('$y(0.5)$',FS,9,IN,LT)
title(['Fig.~18.7.~~Dependence of amplitude ' ...
'on $\lambda$'],FS,11), ylim([0,6])
%%
% \vskip 1.02em
%%
%
% The point of interest is around $\lambda = 3.5$, where we have
% what is a called a {\bf fold} or a {\bf saddle-node bifurcation},
% looking like a pitchfork without the central tine.
% To investigate the structure near such points,
% it is time to advance from {\em continuation}\/ to
% {\bf path-following}. The idea of path-following is that
% we will not just vary a parameter such as $\lambda$, but
% we will follow a path of solutions.
% For example, what happens to one of the solution branches in this plot
% as $\lambda$ approaches $3.5$? Clearly in some
% sense it bends around to turn back in the
% other direction, making $y(0.5)$ a double-valued function of $\lambda$.
% Let us show another plot, then explain what is going on.
%
H = chebop(0,1); H.op = @(x,y,lam) diff(y,2) + lam*exp(y);
H.lbc = @(y,lam) y; H.rbc = @(y,lam) y; lam0 = 0;
[y,lamvec,mvec,lamfun,mfun] = ...
followpath(H,lam0,'measure',@(y) y(0.5),'maxstepno',17);
plot(lamfun,mfun,'-k'), axis([0 4 0 8])
hold on, plot(lamvec, mvec,'.k',MS,12), hold off
xlabel('$\lambda$',FS,10,IN,LT), ylabel('$y(0.5)$',FS,9,IN,LT)
title(['Fig.~18.8.~~The same curve traced ' ...
'by path-following'],FS,11), ylim([0,6])
%%
% \vskip 1.02em
%%
%
% The plot shows a trajectory of solutions that passes around
% the turning point.
% Obviously more was involved in generating this plot than just
% varying~$\lambda$. The data are the dots (not the curve, which
% is just a smooth interpolant), and they have been generated
% by a process called {\em pseudo-arclength continuation,} implemented
% in the Chebfun code {\tt followpath}, whose
% details we shall not describe. The idea is to
% extrapolate from one dot to the next not just with respect
% to~$\lambda$, but with respect to $\lambda$ and also the solution
% $y$. A trajectory is in fact being followed in a high-dimensional
% space, though it is convenient to imagine that it is being followed
% in the two-dimensional space described by $\lambda$ and
% the scalar measure $y(0.5)$.
%
%%
%
% The Bratu equation has just two solutions, but (18.2) has
% four or more, for $\varepsilon = 0.2$ at least, as shown in Fig.~16.13.
% To get a sense of a more complicated bifurcation structure
% let us apply path-following to (18.2).
% We will use $\varepsilon^{-1}$ as a parameter to track them,
% with $y'(1)$ as the scalar measure.
%
%%
% First, we compute the smoothest family of solutions, which
% arises from Chebfun's default initial guess.
H = chebop(-1,1); H.op = @(x,y,epi) diff(y,2)/epi +y+y^2-1;
H.lbc = @(y,ep) y; H.rbc = @(y,lam) y; epi0 = 0.01;
meas = @(y) feval(diff(y),1);
MSN = 'maxstepno'; DI = 'direction';
SM = 'stepmax'; ME = 'measure';
[y1,epi1,m1,epif1,mf1] = followpath(H,epi0,ME,meas,MSN,12,SM,1);
plot(y1,CO,green,LW,0.6)
title('Fig.~18.9.~~First family of solutions to (18.2)',FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% Here is the second symmetric family of solutions,
%
x = chebfun('x'); H.init = 5*(1-x^2); epi0 = 0.5;
[y2,epi2,m2,epif2,mf2] = followpath(H,epi0,ME,meas,MSN,27);
plot(y2,CO,blue,LW,0.6)
title('Fig.~18.10.~~Second family of solutions to (18.2)',FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% Here is the first nonsymmetric family of solutions.
%
epi0 = 10; N.op = @(x,y) diff(y,2)/epi0 + y+y^2-1;
H.init = sin(pi*x);
[y3,epi3,m3,epif3,mf3] = followpath(H,epi0,ME,meas,MSN,13,DI,-1);
plot(y3,'r',LW,0.6)
title(['Fig.~18.11.~~Third family of solutions: ' ...
'broken symmetry'],FS,11)
%%
% \vskip 1.02em
%%
%
% \noindent
% The second nonsymmetric family of solutions is the same, but with
% $x$ replaced by $-x$.
%
H.init = -sin(pi*x);
[y4,epi4,m4,epif4,mf4] = followpath(H,epi0,ME,meas,MSN,13,DI,-1);
plot(y4,CO,purple,LW,0.6)
title(['Fig.~18.12.~~Fourth family of solutions: ' ...
'reflection of the third'],FS,11)
%%
% \vskip 1.02em
%%
%
% Let us now superimpose the $\varepsilon\,$-$\,y'(1)$ data from all
% four of these computations, retaining the same colors.
% We see an approximation to quite
% an interesting bifurcation diagram.
%
plot(epif1,mf1,CO,green), hold on, plot(epi1,m1,'.',CO,green,MS,10)
plot(epif2,mf2,CO,blue), plot(epi2,m2,'.',CO,blue,MS,10)
plot(epif3,mf3,'-r'), plot(epi3,m3,'.r',MS,10)
plot(epif4,mf4,CO,purple), plot(epi4,m4,'.',CO,purple,MS,10)
hold off
xlabel('$\varepsilon^{-1}$',FS,10,IN,LT), ylabel('$y''(1)$',FS,9,IN,LT)
title(['Fig.~18.13.~~Part of the bifurcation ' ...
'diagram for (18.2)'],FS,11)
axis([0 10 -7 7])
%%
% \vskip 1.02em
%%
%
% \noindent
% Each curve in the diagram has something to tell us.
% For $\varepsilon^{-1}= 0$, (18.2) reduces to the linear equation
% $y''=0$ with Dirichlet boundary conditions, with the unique solution
% $y = 0$. For any $\varepsilon>0$, this solution evolves into
% the smooth mode captured by the green curves, always an even
% function of $x$.
% This solution satisfies $y(x)< 0$ for all $x\in (-1,1)$.
%
%%
% The blue curve corresponds to a different even solution that exists for
% any $\varepsilon >0$ but not for $\varepsilon =0\kern 1pt$; as
% $\varepsilon \to 0,$ its amplitude diverges to infinity. This mode
% persists too for all values of $\varepsilon$.
%%
%
% The bifurcation diagram shows that something new happens when
% $\varepsilon^{-1}$ rises above a bifurcation point at
% about $4.5$, that is, when $\varepsilon$ falls below
% about $0.2$ (more precisely, about
% $0.2139$). In a symmetry-breaking pitchfork bifurcation, just like
% what we saw in the final pages of the last chapter, the
% blue solution splits into three, two of which
% are no longer even functions of $x$.
%
%%
%
% As $\kern 1pt\varepsilon\kern 1pt$ shrinks further,
% more and more solution branches appear.
% We shall not attempt to track them, but from Fig.~18.3
% we know they must exist.
%
%%
%
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: Arrhenius chemical reaction}\\[-3pt]
% \hrulefill
% \end{center}
%
%%
% The phenomenon of spontaneous combustion or thermal runaway
% was mentioned in Chapter 3, and it can be modeled by the
% Bratu equation (18.3). Let us rewrite that equation
% with $y$ replaced by $T$, indicating temperature,
% $$ T'' + \lambda \kern -.7pt \exp(T) = 0, \quad
% x\in [\kern .3pt 0,1], ~T(0) = T(1) = 0 . \eqno (18.4) $$
% The function $T(x)$ represents
% the temperature distribution in a one-dimensional
% body with coordinate $x$. In (18.4),
% the body is undergoing a heat-generating
% chemical reaction whose reaction rate
% increases exponentially with $T$. The
% equation balances heat diffusion,
% the second derivative, against the reaction term, and a solution
% to the BVP corresponds to a steady-state temperature distribution.
% If $\lambda$ is small enough, such a solution is possible, with
% heat leaving at the boundaries at $x=0$ and $1$
% fast enough to balance the heat generation in the interior.
% For larger $\lambda$, the heat transfer at the boundaries is
% not fast enough. The medium
% keeps heating up until the temperature
% explodes to $\infty$ in a finite time, and there is no
% steady state. The physics of such
% processes was worked out by the Soviet physicist
% Frank-Kamenetskii in 1939.
%%
%
% Of course, this must be an idealized model if it predicts
% an infinite temperature. As a step towards more
% realistic chemistry,
% one may replace the exponential law in (18.4) by what is
% known as the {\em Arrhenius reaction rate,} which grows
% exponentially for smaller temperatures but then levels off.
% We will pick an explicit constant and write the equation in the
% form
% $$ T'' + \lambda \kern -.7pt \exp(T/(1+0.2\kern .3pt T)) = 0, \quad
% x\in [\kern .3pt 0,1], ~T(0) = T(1) = 0 . \eqno (18.5) $$
% Note that for small $T$, the exponential term grows exponentially
% as before, but as $T\to \infty$, this term is limited by the value $\exp(5)$.
%
%%
%
% An analysis of equations like (18.5) was presented by
% J. R. Parks of the Monsanto Chemical Company
% in ``Criticality criteria
% for various configurations of a self-heating chemical as
% functions of activation energy and temperature of
% assembly,'' {\em Journal of Chemical Physics,} 1961.
% More details of the relevant mathematics can be found in Brown,
% Ibrahim, and Shivaji, ``S-shaped bifurcation
% curves,'' {\em Nonlinear Analysis, Theory, Methods \&
% Applications,} 1981.
% What we find for this equation is that the bounded reaction rate
% shuts off the explosion, allowing solutions for all $\lambda$.
% Moreover, the dependence on $\lambda$ reveals an elegant
% {\em S-shaped bifurcation curve,} which we can track
% with {\tt followpath}.
%
H = chebop(0,1); H.op = @(x,T,lam) diff(T,2) + lam*exp(T/(1+.2*T));
H.lbc = @(T,lam) T; H.rbc = @(T,lam) T; lam0 = 1;
[T,lamvec,mvec,lamfun,mfun] = ...
followpath(H,lam0,'measure',@(T) T(0.5),'stepmax',2,'maxstepno',35);
semilogy(lamfun,mfun,'-k',LW,1), axis([1 7 0.1 300])
xlabel('$\lambda$',FS,9,IN,LT), ylabel('$T(0.5)$',FS,9,IN,LT)
title('Fig.~18.14.~~S-shaped bifurcation curve for (18.5)',FS,11)
%%
% \vskip 1.02em
%%
% The curve shows that for most values of $\lambda$,
% there is a single solution to (18.4). For $\lambda$ approximately
% between $3.5$ and $4.5$, however, there are three solutions.
% For $\lambda = 4$, the amplitudes $T(0.5)$ will be
% approximately 1, 7, and 30, and we can use this information to
% find the solutions and plot them.
N = chebop(0,1); N.lbc = 0; N.rbc = 0;
lam = 4; N.op = @(x,T) diff(T,2) + lam*exp(T/(1+.2*T));
x = chebfun('x',[0 1]);
N.init = sin(pi*x); T = N\0; plot(T,LW,1), hold on
N.init = 7*sin(pi*x); T = N\0; plot(T,LW,1)
N.init = 30*sin(pi*x); T = N\0; plot(T,LW,1), hold off
title('Fig.~18.15.~~Three solutions for $\lambda = 4$',FS,11)
axis([0 1 0 40]), set(gca,XT,0:.2:1,YT,0:10:40)
%%
% \vskip 1.02em
%%
%
% {\sc History.} Pseudo-arclength continuation is an idea
% closely associated with a software package, AUTO, written by
% Eusebius Doedel and growing out of his work with Herbert Keller
% at Caltech in the late 1970s. AUTO is written in Fortran and it
% has evolved and been used for research and applications now
% for four decades. A standard reference on continuation methods
% is Allgower and Georg, {\em Numerical Continuation
% Methods: An Introduction}, Springer, 2012.
% \smallskip
%
%%
%
% {\sc Our favorite reference.}
% A complicated bifurcation diagram is traced
% in extraordinary detail in S. J. Chapman and P. E. Farrell,
% Analysis of Carrier's problem, {\em SIAM Journal on Applied
% Mathematics} 77 (2017), pp.~924--950.
% Gradually reducing the value of $\kern 1pt \varepsilon \kern 1pt$
% in the Carrier equation (16.6), while also using the technique
% of deflation to discover new branches (see Exercise~16.4), the authors follow
% the curves connecting a
% pitchfork bifurcation at $\varepsilon\approx 0.2198$,
% a fold bifurcation at $\varepsilon\approx 0.0814$,
% another pitchfork bifurcation at $\varepsilon\approx 0.0551$,
% another fold bifurcation at $\varepsilon\approx 0.0295$,
% and so on in an infinite sequence of ever increasing
% multiplicity and complexity.
% \smallskip
%
%%
%
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 18.} A powerful technique for the study of
% nonlinear problems is continuation, in which one
% moves incrementally from one solution to
% nearby solutions by varying a parameter. A refinement of this
% idea is path-following, in which a trajectory is followed
% involving not just a parameter but also the solution itself.
% These are fundamental tools of bifurcation analysis.
% \vspace{2pt}}}
% \end{displaymath}
%
%%
% \smallskip\small\parskip=1pt\parindent = 0pt
% {\em \underline{Exercise $18.1$}. Allen--Cahn equation.}
% One solution of the BVP $\varepsilon\kern .3pt y'' + y - y^3=0$, $y(\pm 1) = 0$
% is $0$, but there is another one you will find by setting $\varepsilon=0.1$
% and starting from the initial guess $y(x) = 1-x^2$.
% {\em (a)} Plot this
% solution and report its slope $y'(-1)$.
% {\em (\kern .7pt b)} Plot the solution of the same structure and report $y'(-1)$
% with $\varepsilon = 0.01$ and $\varepsilon = 0.001$.
% {\em (c)} In the other direction, what happens if you increase
% $\varepsilon$ above $0.1$? How far can you go and still find a
% nonzero solution? Based on your explorations, make a plot of
% $y'(-1)$ as a function of $\varepsilon$.
% \par
% {\em \underline{Exercise $18.2$}. Nonlinear boundary layer.}
% In the text it is stated that the width of the boundary layer in
% Figure 18.1 is $O(\varepsilon^{1/2})$. Verify this numerically by
% determining, for each $\varepsilon$ in the plot, the value $\delta$
% such that $y(1-\delta) = 0.5$. Make a table of $\delta(\varepsilon)$
% and also the quotient $\delta(\varepsilon)/\varepsilon^{1/2}$, and
% plot $\delta(\varepsilon)$ against $\varepsilon$
% on a log-log scale.
% \par
% {\em \underline{Exercise $18.3$}. Deforming the S-shaped curve.}
% {\em (a)} Rerun the example of Figure 18.14 with the coefficient $0.2$ of
% (18.4) changed to $0.15$ and $0.35$.
% {\em (\kern .7 pt b)} To two digits of accuracy at least, what is
% the largest coefficient choice that gives a curve that is triple-valued
% for some values of $\lambda\kern .7pt ?$
% \par
% {\em \underline{Exercise $18.4$}. Hysteresis with the S-shaped curve.}
% The subcritical pitchfork bifurcation of Fig.~17.16 led to a jump
% transition and then hysteresis in Fig.~17.18. Devise a similar experiment
% to generate a demonstration of hysteresis for the S-shaped
% bifurcation curve of Fig.~18.14.
%