key: cord-0630705-zuv8mnz4
authors: Clark, William A.; Gomes, Mario W.; Rodriguez-Gonzalez, Arnaldo; Stein, Leo C.; Strogatz, Steven H.
title: Surprises in a classic boundary-layer problem
date: 2021-07-24
journal: nan
DOI: nan
sha: dcb05b762d7a6a22d6b2bf0c7090ed1351d4a2e6
doc_id: 630705
cord_uid: zuv8mnz4

We revisit a textbook example of a singularly perturbed nonlinear boundary-value problem. Unexpectedly, it shows a wealth of phenomena that seem to have been overlooked previously, including a pitchfork bifurcation in the number of solutions as one varies the small parameter, and transcendentally small terms in the initial conditions that can be calculated by elementary means. Based on our own classroom experience, we believe this problem could provide an enjoyable workout for students in courses on perturbation methods, applied dynamical systems, or numerical analysis.

1. Introduction. In many parts of mathematics, physics, engineering, and the life sciences, researchers have developed ingenious techniques for gaining insight into difficult problems by exploiting the presence of a small parameter in them. These techniques, known as perturbation methods, have shed light on all sorts of fascinating phenomena in fluid dynamics, mathematical biology, optics, chemical engineering, quantum mechanics, plasma physics, climate science, and many other disciplines. See [2, 12, 13, 17, 19, 20, 22] for just a few of the textbook introductions to perturbation methods and their applications. This is the story of a textbook problem that has surprised us over and over again, and for which we have come to feel genuine affection. The problem is to solve the following nonlinear differential equation, subject to the given boundary conditions: y = yy − y , y(0) = 1 , y(1) = −1 . (1.1) Here we use the notation y = dy/dx, and we want to solve for y(x) on the domain 0 ≤ x ≤ 1 under the assumption that the parameter is small: 0 < 1. We first met problem (1.1) in the classic textbook by Mark Holmes [13] . Soon after that book appeared in 1995, one of us (Strogatz) decided to discuss (1.1) in an introductory course on asymptotics and perturbation methods. It quickly became clear the problem contains unexpected subtleties and riches.

In the years since then, whenever Strogatz has had a chance to teach that course, he has revisited problem (1.1) and always learned something new about it, thanks to the questions and insights of the students, postdocs, and colleagues in attendance, most notably the co-authors on this paper. Together we think we may have finally gotten to the bottom of it.

But because this is meant to be an education paper, we will also be pointing out some of the false turns we took along the way. Confusion is a natural part of doing mathematics. We did not land on the right way to think about (1.1) initially; it required lots of trial and error and a willingness to be open and vulnerable about what puzzled us at any given time. Experienced researchers know this is all part of the process, but we mention it for the sake of students who may be misled by the confident presentations and pristine appearance of the mathematics they see in most textbooks and journal articles. We want to be a little more honest here about how the sausage is actually made.

Our analyses have relied on three parts of applied mathematics: perturbation theory, nonlinear dynamics, and numerical analysis. As such, we think problem (1.1) could be useful to students or teachers in any of those subjects. We assume that the reader is comfortable with asymptotics and perturbation methods at the level of the books by Holmes [13] or Bender and Orszag [2] , as well as nonlinear dynamics at the level of the book by Strogatz [24] . Not much exposure to numerical analysis is required; a basic knowledge of what it means to solve an ordinary differential equation numerically should be sufficient, say at the level of someone who knows how to use Mathematica or Matlab to solve an initial-value problem. For readers who want to immerse themselves in the details, we also provide supplementary notebooks: a jupyter notebook for the numerical methods, and a Mathematica notebook for the analytical calculations [1] .

Although we have grown enamored of problem (1.1) for its pedagogical value, it would be even more appealing if it also had some real-world applications. Alas, we have not found any so far. But it does have some close relatives of scientific interest. For example, if we change the sign of the right-hand side of the differential equation in problem (1.1), we get an equation known as the Lagerstrom-Cole equation [17, 20] , which has been studied in aerodynamics and fluid dynamics as a model problem in connection with shock layers [20] . The differential equation in (1.1) also pops up in the study of delicate nonlinear phenomena known as canards [3, 7, 9, 10, 15, 18, 19] . As defined by Krupa and Szmolyan [18] , "a canard solution is a solution of a singularly perturbed system which is contained in the intersection of an attracting slow manifold and a repelling slow manifold." In applications, canards often arise in the analysis of nonlinear oscillations having both fast and slow time scales, as seen in chemistry, neurobiology, electronics, and many other fields [19] .

So what are some of the surprises in problem (1.1)? The first is how many solutions it has. Holmes [13] argued that it has a unique solution, but it turns out that it actually has three. Figure 1 shows their graphs. We call these solutions B0, M, and B1, with the names chosen to indicate that they have a boundary layer at x = 0, an interior layer in the middle, or a boundary layer at x = 1.

The graph of the B0 solution has a very negative initial slope y (0) at the left endpoint of the domain, while the other two solutions have y (0) extremely close to 1. In fact, those slopes differ from 1 by a "transcendentally small" term, meaning a term that goes to zero faster than any positive power of as → 0 + . A common example is an exponentially small term of the form exp(−c/ ) where c > 0 is a constant. Normally, such transcendentally small terms are beyond the reach of basic perturbation theory (although sometimes they can be handled by more sophisticated techniques known as exponential asymptotics, superasymptotics, hyperasymptotics, or asymptotics beyond all orders [4, 5] ). What especially surprised us about problem (1.1) is that the leadingorder asymptotics of the initial slope y (0) could be found by elementary means, even when the difference between y (0) and 1 is transcendentally small. As we show in section 6,

for the M solution, while for B1 it is given by

These transcendentally small terms can be calculated using nothing more than higherorder matching, phase plane analysis, and a constant of motion for the associated flow. Not only were we surprised by this good fortune; we were also delighted and relieved by it, given that none of us knew anything about the more sophisticated techniques mentioned above! Yet another surprise is the behavior of the solutions with respect to . In section 7 we'll see that as we increase from zero, the three solutions persist until a critical value c ≈ 0.216, at which point they collide in a pitchfork bifurcation. For larger values of , problem (1.1) has a unique solution that resembles the M solution above.

summarize his analysis here. We are not going to show the mathematical details yet, because they will appear in later sections. For now we just want to emphasize a question that arises whenever one confronts a differential equation like (1.1) with a small parameter multiplying its highest derivative:

Question 1: Does problem (1.1) have any layers (meaning regions of rapid variation in y or its derivatives), and if so, where are they located?

This issue comes up midway through Holmes's chapter on matched asymptotic expansions (chapter 2 in Ref. [13] ). Earlier in the chapter he has already introduced the basic ideas of inner and outer solutions, boundary layers, matching, and uniformly valid composite solutions. There is also a discussion of higher-order matching and examples with multiple boundary layers. In all these prior examples, the layers occur at the endpoints of the domain; in other words, they are genuine "boundary layers." In contrast, problem (1.1) is offered as the first instance of a problem having an "interior layer." Holmes writes:

Generally, when one first begins trying to solve a problem it is not known where the layer(s) is. If we began this problem as we did the previous two and assumed there is a boundary layer at either one of the endpoints, we would find that the expansions do not match. This is a lot of effort for no results, but fortunately there is a simpler way to come to the same conclusion. He then offers a convexity argument to rule out boundary layers at either x = 0 or x = 1, but is careful to note that "these are only plausibility arguments and they do not prove anything. What they do is guide the analysis and hopefully reduce the work necessary to obtain the solution." 3.1. Holmes's plausibility argument. Holmes first considers a possible boundary layer at x = 0. He looks at the governing equation y = yy − y along with its left boundary condition, y(0) = 1. Then he sketches the graph of a candidate solution that looks like our B0 solution in Figure 1 , except that he also assumes that y > 0 in the boundary layer; in other words, y(x) is assumed to be concave up for all small x > 0. On the other hand, although y (hypothetically) has one sign in the layer, y itself clearly changes sign from positive to negative as it drops from y(0) = 1 at the boundary to a value y ≈ −2 where it matches the outer solution. With those sign considerations in mind, Holmes [13] continues:

If the solution behaves like the other example problems we have examined, then near x = 0 it would be expected that y < 0 and y > 0. This means that y > 0 but y(y − 1) is both positive and negative in the boundary layer. This is impossible . . . . It is possible to rule out a boundary layer at x = 1 in the same way. Having ruled out (convex or concave) boundary layers at either end, Holmes [13] next considers the possibility of an interior layer centered at a point 0 < x 0 < 1. He works out three possible forms for the inner solution (depending on whether a certain constant of integration is positive, negative, or zero) and shows that only one of these can be matched to the outer solution. Then he uses a symmetry argument to conclude that the layer must be centered at the midpoint of the domain, x 0 = 1/2, and demonstrates that the inner solution must have odd symmetry about that point. Finally he matches the inner and outer solutions, constructs a uniformly valid composite solution, and shows it agrees with a numerical solution obtained for = 0.01. In this way, Holmes convincingly demonstrates the existence and properties of the solution we called M above, a solution of (1.1) with an interior layer in the middle of the domain.

But what about the possibility of non-convex or non-concave boundary-layer solutions? Recall that the plausibility argument only rules out solutions with y having strictly one sign in the layer. As Holmes's careful wording suggests, that loophole could potentially allow for sneakier solutions where y changes sign within a boundary layer. As we will see in the next section, that little finesse is precisely what allows the solutions B0 and B1 to exist.

4. Another approach to locating the layers: Phase plane analysis. As we worked through Holmes's analytical approach to problem (1.1), we began to wonder if it might be helpful to supplement it with a more geometric style of reasoning known as phase plane analysis. After all, the equation appearing in (1.1), y = yy − y, is a nonlinear, second-order, autonomous differential equation, and phase plane analysis is a powerful tool for illuminating how the solutions to such equations behave. Plus, we have to admit, we have more experience with nonlinear dynamics than perturbation theory, so it felt like a more secure way to approach an unfamiliar problem.

To recast the problem into the language of dynamics, we replace the independent variable x with t, and we think of it as time. Then the dependent variable y becomes a function of time t. The advantage of this approach is that it allows us to use our physical intuition about time and motion and the difference between fast versus slow. Abstract solutions to the differential equation turn into easily pictured trajectories of (imaginary) particles moving around in a two-dimensional space known as the phase plane.

To construct the phase plane, we convert the second-order equation y = yy − y into a pair of first-order equations, and then view those as defining a vector field on the plane. We perform the first step by introducing a new dependent variable z, defined as z = y , and then we rewrite y = yy − y in terms of z to get z = yz − y. By solving for y and z and placing them on the left hand side of a pair of first-order differential equations, we obtain the following vector field:

Next we interpret (4.1) as a dynamical system. From this perspective, the vector (y , z ) then tells us the instantaneous "velocity" that an imaginary particle at (y, z) would have at time t. As the imaginary particle moves around in the (y, z) phase plane, it traces out a trajectory (y(t), z(t)), which is the geometric counterpart of a solution (y(x), y (x)) to the original problem (1.1).

This construction allows us to visualize how the solutions to (4.1) behave by imagining how particles move around in the phase plane. There is no need to be quantitatively precise just yet; a qualitatively correct picture is enough at this stage. By looking at the signs of y and z in (4.1) and sketching a few vectors in various parts of the (y, z) plane, we are led to the picture shown in the left panel of Figure 2 . This picture is called the phase portrait for the system. It shows that there are three qualitatively different types of trajectories for (4.1): parabolic-looking trajectories that flow from left to right above the horizontal line z = 1; a straight trajectory that flows from left to right along the invariant line z = 1; and periodic trajectories below z = 1 that form closed loops, and on which a particle would circulate round and round, always moving clockwise. The right panel of Figure 2 shows the same information quantitatively.

Incidentally, we can prove that the periodic-looking trajectories in Figure 2 are truly periodic and are not merely slowly-winding spirals in disguise. There are two standard ways of proving this, so we will not dwell on the details; see section 6.5 and 6.6 in Ref. [24] for an introduction. Briefly, one way is to note that (4.1) is a "conservative" system. To see this, rewrite it as y /z = dy/dz = z/[y(z − 1)] and then separate variables and integrate to obtain 2 [z + log(1 − z)] − y 2 = constant. This implicit equation can be shown to define closed curves for all z < 1. Another way is to observe that that (4.1) is also a "reversible" system: the vector field (4.1) is unaltered by the change of variables (x, y, z) → (−x, −y, z), which corresponds to a time reversal combined with a mirror reflection across the z-axis in Figure 2 . From this symmetry we can conclude the trajectories lying below z = 1 are composed of two left/right mirror-image halves that together form a bilaterally symmetric loop.

A stronger reversibility symmetry of (4.1) is worth noting: Both the differential equation y = yy − y and its boundary conditions

Hence the possible solutions of the original boundary-value problem (1.1) either come in symmetrical pairs y(x) andỹ(x), whereỹ(x) = −y(1 − x), or the pair degenerates to a single solution with this symmetry, y(x) = −y(1 − x). We already saw a visual manifestation of this symmetry in Figure 1 , where B0 and B1 form a symmetric pair and M is self-symmetric. Likewise, when those same solutions are plotted in the phase plane shown in Figure 3 , the red and blue curves are paired under the symmetry, and the green curve is self-symmetric. The vertical coordinate z denotes y . The red, green, and blue, trajectories correspond to the B0, M, and B1 solutions having layers at x = 0, x = 1/2,and x = 1 respectively. The beginning (x = 0) of each trajectory is denoted with a circle, and the end (x = 1) is denoted with a square; the flow is clockwise on each trajectory. All trajectories start on the dashed vertical line y = 1 and end on the dashed vertical line y = −1, corresponding to the original boundary conditions. Note that the B0 and B1 trajectories form a symmetric pair. As such, both of them lie on the same periodic orbit in the phase plane. The slow part of each trajectory (the "outer solution") occurs at the top, close to z = 1, while the much faster part (the "inner solution" in the layer) occurs everywhere else below that. y(0) = 1 means that at time t = 0 our imaginary particle must start somewhere on the vertical line y = 1 in the phase plane ( Figure 3 ). Its z coordinate on that line, however, is unspecified and remains to be determined; indeed, the key to solving (1.1) is to figure out the initial value of z that will enable the moving particle to satisfy the other boundary condition, y(1) = −1. In dynamical terms, this other boundary condition y(1) = −1 is a final condition, not an initial condition. It says that the particle must reach the vertical line y = −1 in the (y, z) plane after exactly one unit of travel time.

Thus, we see what a difficult challenge our imaginary particle is facing. It must find exactly the right place to start on the line y = 1, such that after it gets carried along by the flow determined by the vector field, somewhat like a tiny speck of leaf being carried downstream by a gentle brook, it manages to land somewhere on the line y = −1 precisely when the clock strikes time t = 1. Figure 2 immediately implies that the trajectories on the line z = 1 or above it are disqualified as candidate solutions because there's no way they can satisfy the boundary conditions: A particle starting on any of them would move monotonically to the right, and so would have no chance of making the leftward journey required to get from y = 1 to y = −1. That leaves the closed loops below z = 1 as our only hope. And indeed, we can imagine a particle starting somewhere on the line y = 1 and then flowing along an arc on one of the closed loops such that if everything is chosen just right (i.e., we pick the right loop to start on), the particle will reach y = −1 after exactly one unit of travel time. Figure 3 shows three arcs that do the trick.

4.4. Slow-fast structure. So far we have not used the assumption that is small, but now we will. For 1, we see from (4.1) that the vector field has large regions where the flow is very fast in the vertical direction, with vertical velocities z of order O(1/ ) occurring at all points in the (y, z) plane where z − 1 = O(1). This region of fast variation, as we will soon see, corresponds to the "inner region" in a perturbation treatment via boundary-layer theory. In the phase portrait, it consists of all points outside the thin gray strip shown schematically in Figure 2 . The strip does not have well-defined edges, but its blurriness does not matter; the key thing is that its thickness is O( ), however we define it. (To check its thickness, observe that if y is O(1) and z − 1 = O( ), then y and z are both O(1) in (4.1), indicating the flow is slow compared to the O(1/ ) speeds achieved everywhere outside the strip.) This strip in the phase plane where the motion is comparatively slow corresponds to the "outer region" in a boundary-layer treatment.

It was by contemplating the slow-fast structure of the flow that we originally came to suspect that there might be more than one solution to problem (1.1). As we reconsidered the solution discussed by Holmes, with its interior layer centered at x = 1/2, we pictured it as a particle moving with a slow-fast-slow trajectory in the phase phase, as shown in the middle panel of Figure 3 . Our imaginary particle spends about half of its travel time dawdling through the initial slow region at the top of the green arc in Figure 3 , then rockets down and around and up again through the fast region in almost no time at all, and then dawdles through the remaining slow region for the remaining half of its travel time. Why, we wondered, couldn't a particle spend nearly all its time in a slow region at the beginning? Or at the end? If the particle started sufficiently close to the invariant line z = 1, or ended up near there, it seemed like these sorts of solutions should also be possible.

This intuition turned out to be correct: such solutions do exist. The blue and red curves in Figure 3 show what they look like as trajectories.

4.5. More backstory: Puzzling over the initial slope. It took us considerable trial and error to find these boundary-layer solutions numerically the first time we looked for them in the computer, more than a decade ago. They eluded us completely when was very small. Fortunately, for only moderately small it was not difficult to find them. For = 0.1, for example, we found that y (0) ≈ 0.9999 yielded a trajectory that was slow at the beginning and fast at the end (the B1 solution), whereas y (0) ≈ −10.6942 gave a trajectory that was fast and then slow (the B0 solution).

The strikingly small difference between 0.9999 and 1 made us wonder what the formula for y (0) as a function of might be for the B1 solution. Likewise, given the size of −10.6942, we were curious how negative y (0) might get for smaller values of as we continued tracking the B0 solution. Our numerics couldn't answer these questions at the time, since some of us were naive about computational methods, so we did not know how to solve the boundary-value problem reliably for 1. Perhaps the initial slope could be found by asymptotic analysis? We felt sure it could but did not immediately see how to do it. We pose that as our next big question and call it the puzzle of the initial slope:

Question 2: For the three solutions of problem (1.1), how do their initial slopes y (0) depend on for 0 < 1?

We will work toward answering that question in the next few sections. But before we leave the phase plane, we should notice that it tells us one more thing of interest. It reveals that the B0 and B1 solutions sneak through the loophole in Holmes's plausibility argument by having layers that are non-convex. To see how that conclusion follows from the phase portrait, observe that on any of the trajectories shown in Figure 3 , the value of z (the vertical velocity) evidently changes sign as the particle goes down and then back up on its journey through the fast region. Since z = y , that change of sign means the concavity of y(x) changes sign in the layer! 5. Perturbation theory. In this section we solve problem (1.1) for 1, both inside and outside the boundary layers or interior layers. Then we match the inner and outer solutions and find composite solutions that give uniformly valid asymptotic approximations of y(x) over the whole domain 0 ≤ x ≤ 1. We perform the match to first order in (i.e., we go beyond the leading order of perturbation theory) because it turns out we need this higher-order information to solve the puzzle of the initial slope (Question 2). In that sense, the following analysis provides a motivational case study of why one would ever want to do higher-order matching. The details of this analysis are included in the supplementary Mathematica notebook [1] .

First we consider the outer region, where regular perturbation theory applies. In this region, we expand y(x, ) in the regular perturbation series

insert this into the differential equation y = yy − y appearing in (1.1), and collect terms having like powers of . At leading order we find

This equation has two possible solutions: y 0 = 0 (which cannot satisfy the boundary conditions), or y 0 = 1, yielding y 0 (x) = x + a for some real constant a. In fact, a bit of study shows that the higher corrections all satisfy y n = 0 for n > 0, and thus y n = a n for some constants a n . But the leading-order solution y 0 (x) already fixes the constant by satisfying the boundary condition, so all the higher constants vanish: a n = 0 for all n > 0. Therefore y 0 (x) = x + a is not merely the zeroth-order approximation to the outer solution; it is the outer solution at all orders of . We can also reach this conclusion by noting that y 0 (x) = x + a satisfies the original differential equation y = yy − y exactly for all values of .

For an outer solution that includes x = 0 in its domain, we can determine the constant a by applying the boundary condition at that endpoint, and similarly for an outer solution that includes x = 1. These two potential outer solutions that satisfy either the left or right boundary condition are

Now we move on to the inner solutions. Suppose there is a layer at x = x 0 . Holmes [13] shows that x 0 = 1/2 is the only possible location for an interior layer, and we are going to show that boundary layers can occur at x 0 = 0 or x 0 = 1 as well. As before, we refer to these three cases as M (layer in the middle), B0 (layer on the boundary at x = 0), and B1 (layer on the boundary at x = 1).

Introduce a layer thickness δ, which is a function of to be determined. Let

be a scaled independent variable that describes positions in the layer, and let Y (X) ≡ y(x) = y(x 0 + δX) (5.5) be a new dependent variable that describes how y varies in the layer. Derivatives of the new variable are

Now our original differential equation y = yy − y becomes

If δ is chosen correctly, then Y and all its derivatives should be O(1) as → 0. Thus we find a distinguished limit when /δ 2 = 1/δ, or simply δ = . Therefore the inner equation in the layer is given by

To solve the inner equation (5.8) asymptotically, we expand Y as

Inserting this series into (5.8) yields, at leading order,

This can be integrated to obtain

where A is an integration constant.

This result has a nice geometrical interpretation. If we recall that y = z in the phase plane, then (5.11) shows that, to leading order, the trajectories in the (y, z) plane follow parabolic arcs as they move through the inner region where the motion is fast. From the phase plane pictures shown earlier, we know that the only parabolic arcs of interest are those with a negative z-intercept, as these are the arcs that lie on the closed loops. Hence we see that A should be negative, say A = − 1 2 b 2 . Then separating the variables in (5.11) and integrating gives

with an additional integration constant c. The two constants b, c are determined by matching to the solution in the outer region.

Zeroth-order matching for the symmetric solution M. Let's see how the matching goes for the interior layer at x 0 = 1/2. We know we have a symmetric solution that satisfies y M (1/2) = 0 = Y M 0 (0) (we could also learn this from the matching alone). This symmetry condition tells us that Y M 0 is an odd function of X and hence c = 0. To get the value of b, we need to look at the large-X asymptotic behavior of the inner solution and match it to the asymptotic behavior of the outer solution as x approaches x 0 = 1/2 from either side. Thus we need to take the following limit, lim X→±∞ Y M (X), and match it to y L 0 (1/2) = 3/2 and y R 0 (1/2) = −3/2. If we recall that tanh z is odd and lim z→+∞ tanh z = 1, we see that both limits agree if and only if b = ±3/2. For either choice of b, our inner solution at zeroth order becomes

Finally we can construct a composite solution y c by the usual recipe: y c = y outer + Y inner − y match . Carrying out those steps for the zeroth-order approximation to M yields

5.5. Zeroth-order matching for the asymmetric solutions B0 and B1. We can proceed similarly for the B0 and B1 cases. In fact we only need to do the work for one of them, since we can get the other one from the symmetry transformation (x, y) → (1 − x, −y). Let's focus on the B0 case, which has a layer at x 0 = 0. Its zeroth-order inner solution (5.13) needs to be matched to y R 0 (0) = −2, which tells us that b = −2. We also need to satisfy the boundary condition y(0) = 1, which is now inside the layer, so Y B0 0 (0) = 1 determines the value of c as c = − tanh −1 1 2 . Hence

Similarly for B1 we get

Now constructing the leading-order composite solutions, we get

5.6. Remarks about the zeroth-order composite solutions. The leadingorder composite solutions for M, B0, and B1 look extremely similar to the numerical solutions plotted in Figure 1 . They are also uniformly valid, as can be checked by looking at the difference between the numerical and analytical solutions, as plotted in Figure 4 for B0. Moreover, the layers for all of the zeroth-order solutions are non-convex, as we expected from our earlier phase plane analysis. For example, the B0 solution (5.18) has an inflection point at x = tanh −1 1

2 . Yet informative as these leading-order solutions are, they are not accurate enough to allow us to calculate the initial slope y (0) correctly. So, let's proceed to the next order.

First-order matching for the symmetric solution M. As we noted earlier, the outer solutions are y 0 (x) = x + a to all orders in ; only the inner solution and the matching change as we proceed to higher orders. To study the first-order The error is defined as the difference between the asymptotic solution (5.18) and a very careful numerical solution (taken as a surrogate for the unknown exact solution). The scaled error is defined here as the error divided by ; this scaling is appropriate because we expect to incur errors of size O( ) in a leading-order solution. We examine the error within the boundary layer, x = O( ), i.e., the layer has a width linear in (the error outside of the layer is transcendentally small). Notice that the composite solution is uniformly valid, and the error is proportional to , as expected from a zeroth-order solution. The error behavior is similar for the other two asymptotic solutions with layers at x = 1/2 and x = 1.

correction Y 1 in the inner solution, we insert the series (5.9) into our inner equation (5.8) and collect the first-order terms. We find that Y 1 (X) satisfies

Compared to the asymmetric B0 and B1 solutions, the symmetric solution M yields the simplest expressions for Y 1 , so we focus on that calculation now and relegate the others to Appendix A.

Using the zeroth-order solution for Y M 0 from (5.14), we find that Y M 1 satisfies the following second-order inhomogeneous linear differential equation with variable coefficients:

We must solve this beast subject to the side condition that Y M 1 (0) = 0 (a condition that follows from symmetry, as Y M is an odd function). Impressively, Mathematica obliges and produces a long expression that we give in (A.1) of Appendix A. Of the two free constants of integration appearing in (A.1), we fix the constant c 2 by imposing the side condition Y M 1 (0) = 0, giving c 2 = π 2 /18. To determine c 1 , we need to perform a first-order match to the outer solution y R 0 on the right (there's no need to worry about additionally matching to y L 0 , the outer solution on the left; matching on the right automatically takes care of matching on the left, by the odd symmetry noted above).

To perform the matching on the right, we need to know the asymptotic behavior of Y M 1 as X → +∞. The relevant asymptotics are:

where Li s (z) is a special function called a "polylogarithm" of order s [8]. The polylogarithms can be defined by their power series or recursively from an integral,

For s < 2, they are elementary functions, for example Li 1 (z) = − log(1 − z).

With further help from Mathematica we eventually find

Miraculously-and yet not miraculously at all, if one believes in perturbation theorythis function has the exactly the right large-X behavior, Y M 1 ∼ X, needed to match onto the outer solution, if the constant is correct! Recall that at zeroth order, lim X→∞ Y M 0 (X) = −3/2 already matches the value of lim x→1/2 + y R 0 (x). Therefore we want the additive constant above to vanish (or else we would make an O( ) error in the matching). Hence

Finally, by again invoking the recipe y c = y outer +Y inner −y match for forming a composite solution, we obtain the composite solution up to terms of order 2 . Remarkably, the y match that we need to subtract off here is simply the exact outer solution y R 0 . So the whole composite solution boils down to the inner solution Y M 0 + Y M 1 . Thus, the first-order composite solution for M is:

where X = (x − 1 2 )/ . We proceed similarly for the B0 and B1 solutions, collecting the results in Appendix A. Figure 5 confirms that the error between the first-order asymptotic solution for B0 and a numerical solution (presumed to be close to exact) truly does shrink in proportion to 2 , as it should for a first-order match. This sort of test provides a reassuring check when doing complicated numerics and asymptotics. 6. Solving for the initial slope. Now that we have constructed asymptotic approximations to the three solutions M, B0, and B1 of our original problem (1.1), we can use those approximations to estimate the initial slope y (0) in each case. Knowing this initial slope is theoretically interesting since (as we'll see) it depends on in an intriguing way. But it's also practically useful information: having a good approximation of the initial slope helps us solve the boundary-value problem (1.1) numerically.

One computational approach to solving a boundary-value problem is called the "shooting method." A more thorough discussion can be found in several textbooks on numerical methods, for example [23] . We also provide a supplementary jupyter notebook which implements the numerics described in this section [1] . Figure 6 shows an example of how shooting applies to our problem. We start a trial solution at y(0) = 1 and launch it with some initial slope y (0), somewhat like shooting an artillery shell at an intended target. In our case, the target is the point at the other boundary condition: x = 1, y = −1. Incorrect choices of the initial slopes at x = 0 will produce solutions that fail to hit the boundary condition at x = 1. The three curves shown on the left in Figure 6 indicate what happens if we aim too high or too low or just right.

If we compute where we hit for many possible y (0), and plot the resulting y(1)'s versus y (0), we get the graph of the target function shown in the right panel of Figure 6 . To find a numerical solution of our problem, then, we just have to figure out where the graph crosses the dashed horizontal line y = −1. This is a standard numerical task; it amounts to a root-finding problem in a small neighborhood of the solution. Notice that for the value = 0.1 used to make Figure 6 , the graph of the target function crosses the dashed line in three places. Those are the desired initial slopes of our three solutions. To find analytical estimates of these slopes, we will see next why we needed to go to the trouble of doing higher-order matching. The leading behavior here is y B0 ∼ −3/(2 ), and this is indeed correct, as we will soon see. However, the slope in (6.1) has an error that is O(1). This can be seen graphically in Figure 4 : Within the layer of width O( ), there is an O( ) error in the value of y, leading to an O(1) error in the slope. Or, reiterating the point in another way: The derivative of a function's asymptotic series (at some order) is not necessarily the asymptotic series of the derivative of the function (to the same order). Since we have the composite solution to higher order, we can easily check what the correction is by using the results for B0 in Appendix A. Evaluation and differentiation of the first-order composite solution at x = 0 yield

This is now the full result for y B0 (0) up to errors of order O( ). As shown by the top curve of Figure 7 , this asymptotic result agrees nicely with the value of y (0) produced by the shooting method.

6.2. Initial slope: The tempting but wrong way. For the B0 case, we saw that the O(1) term in the initial slope changed from 1 to 1 + log 16 when we went to first order. This change in the O(1) term is a rather small effect as → 0, since the −3/(2 ) term dominates in this limit anyway. But to our surprise (probably because of our inexperience in these matters), we soon discovered that the B1 and M cases are much more subtle: differentiating the first-order composite solution does not give the correct initial slope, not even at leading order in .

It's an instructive trap to fall into, so let's take the plunge. For simplicity, let's work with M. If we calculate the initial value and initial slope of the zeroth-order composite solution y M c,0 given in (5.15), we find

where TST stands for "transcendentally small terms," i.e., terms that are smaller than any power of times the smallest reported term.

Should we trust these results? Let's check by going to the next order of perturbation theory. Our first-order composite solution yields

This doesn't look good at all! Going to the next order in y has resulted in a change at a lower order to the values of y(0) and y (0). Apparently we can't trust this. So for this problem at least, naively differentiating the composite solution does not give us the correct initial slope.

6.3. Initial slope: The right way. Instead of the approach above, let's turn to a more global analysis, using our knowledge of the structure of solutions in the phase plane. The key insight is that we can use a conserved quantity (also known as a constant of motion, or a first integral) to transfer trustworthy information from inside a layer to a distant point outside the layer where we want to calculate an initial slope. For example, we can transfer information from the M layer at x 0 = 1/2 all the way over to x = 0; this trick is how we are going to extract the leading (but still minuscule!) transcendentally small term in M's initial slope. The advantage of using a conserved quantity is that it is exact; it allows us to shuttle information around in the phase plane with perfect fidelity.

6.3.1. Conserved quantity for the trajectories. To set the stage to perform the desired transfers of information, we need to gather a few facts about the conserved quantity.

We have already mentioned that our vector field (4.1) is conservative. Recall that at each point in the (y, z) phase plane,

which can be separated and integrated to yield

where C is a constant that labels the trajectories. Since C is constant, the quantity 2 [z + log(1 − z)] − y 2 remains unchanged as y(t) and z(t) flow along a trajectory and hence is a "conserved quantity." For a given value of C, we can generate two explicit formulas for the trajectories as curves in the (y, z) plane. Either we can write y(z) in the right or left half plane by solving for y and using the relevant branch of the square root. Or we can solve for z(y) if we allow ourselves to use the implicitly-defined Lambert W function, which satisfies W (z) exp(W (z)) = z . (6.9) This equation has multiple solutions, giving the multiple branches W n (z). For more on the Lambert W function, see [6] . From (6.8) the explicit solution for z is

For 0 ≤ z < 1, we want the branch W 0 . For z ≤ 0, we want the branch W −1 .

The existence of a conserved quantity for Eq. (4.1) suggests that the dynamical system might actually be a Hamiltonian system in appropriate coordinates. Indeed, Eq. (4.1) is Hamiltonian for z < 1: Canonical variables are Q = y and P = log(1 − z), in terms of which the vector field corresponding to (4.1) becomes Q = 1−e P , P = Q/ , with a corresponding Hamiltonian H(P, Q) = P − e P − Q 2 /(2 ).

Having set the stage, we're now ready to explain how to transfer slope information reliably with the help of the conserved quantity 2 [z + log(1 − z)] − y 2 . Here's the idea: Recall that, by definition, z = y , so z represents a slope y on a graph of y versus x. Suppose we have one point (y 1 , z 1 ) on a solution where we trust the slope z 1 . Then, by using the constancy of C, we can use this information to get the slope z 0 at some other y 0 on the same solution. For the M and B1 solutions, we will take (y 1 , z 1 ) to be a convenient point inside the layer where the inner solution is trustworthy and we can compute z 1 = O( −1 ) accurately. Then we will transfer this information over to the initial condition at x = 0, where y 0 = 1 and we seek z 0 .

Taking two copies of (6.8) on the same curve of constant C and eliminating C, we find that our two points are related by

Let y 0 = 1. Solve for the desired initial slope z 0 by using the appropriate branch of the Lambert W function:

Equation (6.12) is the key formula for transferring information from one point to another. Now if we can find trustworthy values of y 1 and z 1 to plug in, we'll be in business. Let's illustrate the idea with the B1 solution. Remember that B1 is related to the B0 solution by the symmetry transformation (x, y, z) → (1 − x, −y, z). That means that the final slope in B1's layer at x = 1 is the same as the initial slope in B0's layer at x = 0. But we have already calculated B0's initial slope reliably! It's given in (6.2). (It's reliable because it was calculated inside a layer, namely, the layer at x = 0 for the B0 solution.)

Hence the trustworthy point (y 1 , z 1 ) to pick for the B1 solution is

When we plug this point into (6.12), we soon discover something fascinating about the initial slope z B1 (0): it deviates from 1 by a transcendentally small quantity. To see this, however, takes a few more steps. Upon performing the substitution we first obtain

which may look opaque to anyone unfamiliar with the Lambert W function. Fortunately this expression can be simplified by using the fact that W 0 (z) is analytic at z = 0 and its power series is convergent within a radius |z| < 1/e (see [6] for a thorough treatment). The first few terms are

as can be verified by substituting back into the defining equation (6.9). We only need to keep the first term since the second term is already transcendentally small in relative to the first term, in light of the form of the argument of W 0 in (6.14).

Thus, we finally arrive at the correct asymptotic behavior for B1's initial slope:

which can be further simplified to

To obtain this last equality, we simplified (6.16) by replacing exp(log 16) with 16 and exp(O( )) with 1 + O( ). That relative error at O( ) in (6.16) times the leading −1 term results in our ignorance of the subdominant O(1) term in the prefactor multiplying the controlling exponential in (6.17) . Nevertheless, we still have enough information to nail down the leading-order -dependence of the initial slope, given by the term [24/ ] exp(−3/2 ).

In retrospect, these error considerations clarify why we needed go to the bother of approximating the inner solution with higher-order perturbation theory: It's because we needed the error in the argument of the exponential in (6.16) to be O( ), which necessitated an error of O( 2 ) in the inner solution.

For the M solution, we redo the calculation above, except now we use the slope in the layer at x = 1/2 as our trustworthy value of z 1 . We get that slope with sufficient precision by simply evaluating the derivative of the first-order composite solution (5.28) at x = 1/2. We also recall that y vanishes at x = 1/2, by the odd symmetry of M about its midpoint. Hence the trustworthy values of y and its slope are Here, as in (6.16), we have enough precision to nail the leading-order term after expanding exp(O( )), but not enough to determine the subleading O(1) correction. Figure 7 compares these analytical initial slopes against numerical results obtained from the shooting method. In all three cases, the agreement between asymptotics and numerics is excellent. And for both M and B1, the agreement extends over more than twenty orders of magnitude. Wow! 7. Pitchfork bifurcation. As we've seen, the nonlinear boundary-value problem (1.1) has three solutions when is sufficiently small. The final surprise in this problem comes when we examine what happens to the three corresponding initial slopes, y (0), as we increase away from 0. We can mull over the three equations (6.2), (6.16), and (6.21), or examine Figure 7 that summarizes all three. For sufficiently small , we have the ordering y B0 (0) < y M (0) < y B1 (0). As increases, y B0 (0) increases (becomes less negative) and moves toward zero, while both y M (0) and y B1 (0) decrease away from 1 and also move toward zero. This behavior suggests that at sufficiently large , there is a possibility that two or even three solutions might approach each other and merge. It also suggests that the initial slopes of the merging solutions might lie somewhere close to y (0) = 0. We first discovered experimentally that this actually happens; a three-way merger occurs through a pitchfork bifurcation precisely when y (0) = 0. After the fact, we were able to establish what analytical conditions describe the bifurcation and identify the critical value c where it occurs. 7.1. Visualizing and explaining the pitchfork bifurcation. We visualize the merger of solutions by plotting contours of 1 + y(1) in the ( , y (0)) plane in Figure 8 . Satisfying the boundary-value problem means finding the contours where y(1) + 1 = 0, which occurs at the boundary between reds and blues in the figure.

Below the critical value < c , there are three solutions; these bifurcate at = c and y (0) = 0. For large values of , there is only one solution. In the broad view shown in the top panel, there is a rather obscure shape for the solution set. In the bottom panel, we zoom in to the sixth decimal place in , and see a classic pitchfork shape, which we will explain below.

The first condition for the bifurcation comes from understanding the relationship between the initial and final slopes for the B0 and B1 solutions, which are paired under the symmetry (x, y, z) → (1 − x, −y, z). Let's return to (6.11) , the relationship between two points (y 0 , z 0 ) and (y 1 , z 1 ) on the same trajectory. Take these two points to be the endpoints of the trajectory, where y 0 = 1 and y 1 = −1. Then the y 2 terms cancel, leaving a simpler condition relating the two endpoint slopes:

For convenience, let us define a function f (z) = z + log(1 − z) , (7.2) whose domain is z < 1. Equation (7.1) says that f (z 0 ) = f (z 1 ). Now there are two possibilities. For the self-symmetric solution M, the slopes at the endpoints agree automatically, so z 0 = z 1 and (7.1) is vacuously satisfied. This case tells us nothing new. However, for the asymmetric solutions B0 and B1, there are two distinct slopes, z 0 = z 1 , and we can show that they must have opposite signs.

This result follows from the shape of the graph of f . Note first that f (z) = 1 − 1/(1 − z), so f has only one local extremum, at z = 0, where f (0) = 0. The second derivative f (z) = −(1 − z) −2 is strictly negative everywhere, so z = 0 is in fact a global maximum. For z < 0, f (z) is monotonically increasing, while for 0 < z < 1, f (z) is monotonically decreasing. Therefore, in the case where we have two distinct slopes f (z 0 ) = f (z 1 ) but z 0 = z 1 , we can deduce that they must lie in the intervals z 0 < 0 and 0 < z 1 < 1, or vice versa, and hence have opposite signs, as claimed.

Graphically, every horizontal slice through the graph of f (z) below its global maximum cuts it in two places: one value z 0 < 0, and the other value at 0 < z 1 < 1. As we increase our horizontal slice towards the maximum, these two roots both approach and ultimately confluence at z = 0.

Finally let's apply this knowledge to the slopes z B0 and z B1 , when both solutions exist. The B0 solution has endpoint slopes y (0) = z B0 and y (1) = z B1 , and the B1 solution exchanges these two. By the previous analysis, f (z B0 ) = f (z B1 ) with z B0 < 0 while 0 < z B1 < 1. We know that as increases, these two slopes approach each other. The only place they can merge is at z B0 = z B1 = 0, which happens at some = c .

What about the intermediate M solution? How does its slope compare to the other two? We know that the three slopes are the three roots of the target function plotted in the right panel of Figure 6 . Is it possible for z M to confluence with one of z B0 or z B1 before the outer two roots meet each other? No: the M solution has the same slope at each endpoint, so if it confluences with either of the other two roots, all three must confluence simultaneously. Generically, this leads us to expect that in a small neighborhood of the bifurcation point, the error function y(1) + 1 plotted in Figure 6 should be approximately a cubic of the form

We further expect A > 0 and B > 0 by examining the shape of the curve in Figure 6 . For > c , there is only one real solution; for < c , there are three real solutions; and as → c from below, these three solutions degenerate into a triple root. All of these expectations are confirmed by the classic pitchfork scenario seen in the lower panel of Figure 8 . 7.2. Calculating the pitchfork bifurcation value. The preceding analysis told us that the bifurcation happens when (1.1) admits a solution with y (0) = 0. We can now use that result to determine c .

To do so, we first observe that the solution of (1.1) when = c must satisfy the boundary condition y(0) = 1 as well as the bifurcation condition y (0) = 0. These conditions then uniquely determine the corresponding critical trajectory: it starts at (y(0), z(0)) = (y(0), y (0)) = (1, 0) and also satisfies 2 [z + log(1 − z)] = y 2 − C 2 by (6.8). So by plugging in y(0) = 1 and z(0) = 0 we see that C = 1 on the critical trajectory.

Next, we use C = 1 to find two conditions on c . Both conditions relate c to z c , defined as the minimum (i.e. most negative) value of z on the critical trajectory. By solving those two conditions simultaneously, we find z c and c , as follows.

To obtain the first condition on c , we note that z = 0 when z(t) reaches its minimum along the trajectory, which implies (from (4.1)) that y = 0 there. Like all other points on the critical trajectory, this point (y, z) = (0, z c ) must satisfy (6.8) with C = 1. This gives us our first condition on c :

The second condition is that the time required for the critical trajectory to go from y = 1 to y = −1 is t = 1, just as it is for every solution of our original boundary-value problem. But for the critical trajectory, we can say more. The critical trajectory is self-symmetric, which implies that the time to go halfway is simply t = 1/2. And at that halfway point, the trajectory is at the point we have just been discussing, (y, z) = (0, z c ). To translate these observations into the condition we seek, we need to find a formula for the travel time. The trick is to write dt in terms of z on the trajectory and then integrate. Recall that z = y(z − 1)/ , from (4.1), so dt = dz/z = dz/[y(z − 1)]. Thus where we are thinking of y as a function of z on the critical trajectory. That function can be written explicitly by using the fact that C = 1 on the critical trajectory; solving 2 c [z + log(1 − z)] = y 2 − 1 for y and noting that y ≥ 0 on the first half of the trajectory, we take the positive square root and obtain y = 1 + 2 c [z + log(1 − z)]. Now we have two conditions, (7.4) and (7.5) , in terms of z c and c . We use (7.4) to eliminate c , leaving an integral equation for z c alone to satisfy. After a bit of algebra, this combined condition is 0 = g(z c ) = z c + log(1 − z c ) + This is now a root-finding problem for the function g(z c ). We find z c ≈ −3.9052637703. Plugging this root back into (7.4) yields c ≈ 0.2159869288903.

A deeper look at the pitchfork bifurcation. In retrospect, there were at least two reasons to expect that a pitchfork bifurcation could occur in our problem. The reasons have to do with the symmetry of the problem itself and the special structure of the associated vector field.

First, recall that the differential equation y = yy −y and its boundary conditions y(0) = 1, y(1) = −1 are left unchanged by the transformation (x, y) → (1 − x, −y). If we apply this transformation a second time, we get back to (x, y), as if toggling a switch or reflecting an image in a mirror twice. Such a transformation is known as a Z 2 symmetry. Now it turns out that the occurrence of a pitchfork bifurcation in a boundary-value problem with Z 2 symmetry is a codimension-1 phenomenon [21] , which means that we should expect to see it in a generic one-parameter family of such problems. So we should not be surprised to find a pitchfork occurring in our problem as we vary its single parameter, .

Second, recall from Sec. 6.3.1 that the dynamical system (4.1) is Hamiltonian for z < 1. McLachlan and Offen [21] give the generic bifurcations for Hamiltonian boundary-value problems and find that a pitchfork can be codimension-1 in planar problems even without the kind of symmetry seen in our system.

So if we had known what we know now, we should have expected a pitchfork all along. As always, everything becomes clearer in hindsight. Fortunately, before clarity comes, we have the pleasure of being surprised. where, as above, X = x − tanh −1 1 2 . For the B1 solution, we can use the symmetry to write y B1 (x) = −y B0 (1 − x).

BVP for a supplementary jupyter notebook demonstrating the numerical approach, and a supplementary Mathematica notebook demonstrating the analytical calculations

Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory

Weakly Nonlocal Solitary Waves and Beyond-All-Orders Asymptotics: Generalized Solitons and Hyperasymptotic Perturbation Theory

The devil's invention: asymptotic, superasymptotic and hyperasymptotic series

On the Lambert W function

The canard unchained or how fast/slow dynamical systems bifurcate

Canard Cycles and Center Manifolds

Relaxation oscillations including a standard chase on french ducks

Array programming with NumPy

Perturbation Methods, Cambridge Texts in Applied Mathematics

Introduction to Perturbation Methods

Matplotlib: A 2d graphics environment

On fast-slow consensus networks with a dynamic weight

mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 1.2.1)

Multiple Scale and Singular Perturbation Methods

Extending geometric singular perturbation theory to nonhyperbolic points-fold and canard points in two dimensions

Multiple Time Scale Dynamics

Matched Asymptotic Expansions: Ideas and Techniques

Bifurcation of solutions to hamiltonian boundary value problems

Singular Perturbation Methods for Ordinary Differential Equations

Numerical Recipes: The Art of Scientific Computing

Nonlinear Dynamics and Chaos: With Applications to Physics

Lecture 16: A tricky nonlinear boundary-value problem

SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python

Acknowledgments. This manuscript was written during the global COVID-19 pandemic, which has affected different members of society inequitably. The authors were quite privileged to have had the opportunity to collaborate and continue to do research-opportunities not afforded to everyone. At the same time, one impetus for this work was that S.H.S. taught his course on asymptotics and perturbation methods online, and L.C.S. became aware of Strogatz's lectures because of Twitter, and was able to follow along on YouTube (the problem at hand was discussed in [25] ). L.C.S. would therefore like to thank Twitter and YouTube for facilitating his initial involvement in this research. We would also like to thank John Boyd, Mathieu Desroches, John Guckenheimer, Martin Krupa, Christian Kuehn, Ian Lizarraga, and Richard Rand for their helpful advice. In our numerical work we made use of the python package mpmath [16] for arbitrary precision calculation, numpy [11] , scipy [26] , and matplotlib [14] .

Appendix A. Higher-order asymptotics.Here we collect the long expressions obtained at higher orders of perturbation theory. For the middle layer solution M, at first order, we findwhere X ≡ (x − 1/2)/ . As discussed in the main text below (5.21), after matching, the integration constants take the values c 1 = π 2 /18, and c 2 = 1 + log 4.For the B0 case, we condense the notation a bit by defining X ≡ x/ and