%% 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. %