## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Computation of symmetric, time-periodic solutions of the vortex sheet with surface tension

Edited by Alexandre J. Chorin, University of California, Berkeley, CA, and approved December 7, 2009 (received for review September 21, 2009)

## Abstract

A numerical method is introduced for the computation of time-periodic vortex sheets with surface tension separating two immiscible, irrotational, two-dimensional ideal fluids of equal density. The approach is based on minimizing a nonlinear functional of the initial conditions and supposed period that is positive unless the solution is periodic, in which case it is zero. An adjoint-based optimal control technique is used to efficiently compute the gradient of this functional. Special care is required to handle singular integrals in the adjoint formulation. Starting with a solution of the linearized problem about the flat rest state, a family of smooth, symmetric breathers is found that, at quarter-period time intervals, alternately pass through a flat state of maximal kinetic energy, and a rest state in which all the energy is stored as potential energy in the interface. In some cases, the interface overturns before returning to the initial, flat configuration. It is found that the bifurcation diagram describing these solutions contains several disjoint curves separated by near-bifurcation events.

Many complex and rich phenomena in nature are controlled by coherent time-periodic or traveling structures. Although the computation of traveling waves is often straightforward, most numerical methods for computing time-periodic solutions were designed with ordinary differential equations in mind and are too expensive for partial differential equations (PDEs). We develop an adjoint-based optimal control algorithm for solving general two-point boundary value problems and use it to perform a computational study of the existence of time-periodic solutions of the vortex sheet with surface tension, which is the interface between two incompressible, irrotational, inviscid, immiscible fluids shearing past each other. This system poses many technical challenges for the method, and is of considerable interest in mathematical fluid mechanics.

We were drawn to several unique features of this problem. First, although the initial value problem is locally well posed (1–5), singularities can form in finite time due to self-intersection (6, 7) or, more speculatively, through the development curvature singularities (8, 9); hence, periodic solutions are special in that they are global solutions that remain smooth for all time. Second, asymptotic models of interface problems in fluid mechanics are often integrable; the KdV and Benjamin–Ono equations are two such examples. Craig and Worfolk have disproved a conjecture of Dyachenko and Zakharov on the integrability of free surface hydrodynamics (10). Nevertheless, studies of periodic solutions should help illuminate the connection between free surface flows for the full Euler equations and integrable model equations obtained in various asymptotic limits. Third, there are many interesting questions in fluid mechanics regarding ergodicity, recurrence (11–13), and the role of viscosity in fluid mixing. Along these lines, we note that periodicity is beginning to play an important role in the study of turbulence (14–16). Finally, interface problems in fluid mechanics generally suffer from small divisor problems that require variants of Nash–Moser and Kolmogorov–Arnold–Moser theory (17, 18) to study time periodicity. By developing such tools, Plotnikov and Toland (19) and Iooss et al. (20) have proved existence of time-periodic gravity-driven water waves. We aim to learn more about time-periodic interface problems by developing robust numerical methods capable of solving such problems whether or not small divisors are present.

Our numerical method involves two key ideas. First, by adapting adjoint-based optimal control methods (21–24) originally developed in the shape optimization community, we are able to use quasi-Newton line search algorithms such as the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method (25) to solve two-point boundary value problems rather than the standard methods of orthogonal collocation (26) or shooting (27). This leads to a tremendous reduction in computational cost, especially when approximate Hessian information from the previous solution is used in the continuation algorithm. The method is a variant of the one developed by the authors in (28, 29) for the Benjamin–Ono equation, but is necessarily more complex as the motion of the vortex sheet with surface tension is described by a coupled system of nonlinear equations rather than a single equation, and involves singular integrals. Second, to solve the forward and adjoint problems, we use a fourth-order additive Runge–Kutta method (30, 31) rather than an implicit–explicit multistep method (32) such as adopted by Hou et al. (6, 7). In either approach, a small-scale decomposition (developed in ref. 6) is employed in which the most singular terms in the evolution equations are treated implicitly to remove stiffness, whereas nonlinear terms are treated explicitly. The advantage of the additive Runge–Kutta framework is that the implicit part of the method is *L* stable. By contrast, high-order multistep methods lack *A* stability and must be filtered when used for dispersive problems.

## Equations of Motion

Following refs. 1, 6, 7, we consider two irrotational, ideal fluids of equal density separated by a sharp interface, which is a curve (*x*(*α*,*t*),*y*(*α*,*t*)) parametrized by *α*∈[0,2*π*) and time. We assume the curve is 2*π* periodic in the horizontal direction, i.e., *x*(*α* + 2*π*,*t*) = *x*(*α*,*t*) + 2*π*, *y*(*α* + 2*π*,*t*) = *y*(*α*,*t*). The jump in pressure across the interface is [*p*] = *τκ*, where *τ* > 0 is the (constant) coefficient of surface tension and *κ* is the curvature of the interface. We define the arclength element of the curve, *ds* = *σdα*, and tangent angle, *θ*, by (*σ* cos *θ*,*σ* sin *θ*) = (*x*_{α},*y*_{α}). We denote the tangent and normal vectors to the curve by and . We let *U* denote the normal velocity of the curve and *V* the tangential velocity of the curve. We denote by *L*(*t*) = 2*πσ*(*t*) the length of one period of the curve.

The evolution of *θ* and *σ* can be inferred from the evolution : [1]The curve is initially parametrized by arclength, normalized so that *α*∈[0,2*π*). We choose the tangential velocity to be a nonphysical velocity that maintains this normalized arclength parametrization at all positive times. Thus, we have [2]where [3]are orthogonal projections onto the mean and onto the space of zero-mean functions, respectively, and is the zero-mean antiderivative operator.

The normal velocity, , is determined by the fluid dynamics via the Birkhoff–Rott integral, **W** = (*W*_{1},*W*_{2}): [4]Here *z* = *x* + *iy*. The vortex sheet strength, *γ*, and true vortex sheet strength, , are related via or . The complex cotangent comes from summing over periodic images, . The reader could consult, for instance, the book of Saffman (33) for details on the Birkhoff–Rott integral. The evolution equation for *γ* is [5]where *V*_{1} is the average tangential velocity of the fluid across the interface, (*u* + *iv*)^{±} = [∓(*γ*/2*σ*) + *V*_{1} + *iU*]*e*^{iθ}.

We will need to reconstruct the curve from *θ* and *σ*. This is done by integrating the identity *z*_{α} = *σe*^{iθ} and using *z*_{t} = (*V* + *iU*)*e*^{iθ} to evolve the integration constant: [6]So we see that to reconstruct the curve, we only need to evolve one point in addition to *θ* and *γ*, namely, *z*_{0} = *P*_{0}*z* - *π*.

If *A* ≠ 0 and (*σ*(*t*),*θ*(*α*,*t*),*γ*(*α*,*t*),*z*_{0}(*t*)) is a solution with surface tension *τ*, then (*σ*(*At*),*θ*(*α*,*At*),*Aγ*(*α*,*At*),*z*_{0}(*At*)) is a solution with surface tension *A*^{2}*τ*. Thus, by rescaling time and vortex sheet strength, we may assume *τ* = 1. One may also show that is another solution.

If at any moment *γ*(*α*,*t*) and *θ*(*α*,*t*) are both even functions of *α*, then the real and imaginary parts of *e*^{iθ(α,t)} will be even and those of [*z*(*α*,*t*) - *z*_{0}(*t*)] will be odd. A change of variables in Eq. **4** then shows that **W**, *U*, and *V*_{1} are odd functions. Because *V* is also odd, *γ*_{t} and *θ*_{t} are even, whereas *z*_{0t} = 0. If, in addition to being even, *γ* and *θ* change sign upon translation by *π*, this will also remain true for all time. In this paper, we look for time-periodic solutions with initial conditions of the form [7]where {*c*_{k}∶*k* = 1,3,5,…} are real numbers and a hat denotes a Fourier coefficient, . By the above symmetry arguments, and remain real for all time, and remain zero if *k* is even. If at some time *T*/4 the solution with initial conditions [**7**] evolves to a state in which *γ* ≡ 0, a time-reversal argument (with *A* = -1 above) shows that the solution will evolve back to a flat state at *T*/2 with the sign of *γ* reversed. The evolution of *γ* and *θ* from *T*/2 to *T* will be identical to that from 0 to *T*/2, but with opposite signs, ending at the original state.

## Results

The standard approach to proving existence of time-periodic solutions of nonlinear PDE is to build periodicity into the solution space and use a Newton iteration (17, 18) to solve a system of lattice equations for the spatial and temporal Fourier modes. Newton’s method converges rapidly enough that “small denominators” can be dealt with via small numerators. Our numerical method is based instead on searching for *c*_{k} and *T* such that the Cauchy problem with initial conditions [**7**] satisfies *γ*(·,*T*/4) ≡ 0. We define and wish to solve *F* = 0. If the standard approach were turned into a numerical method, it would resemble spectral collocation (26), which is very expensive for PDE. If our approach were used for analytical purposes, it would also suffer from small divisor problems.

We begin our search for time-periodic solutions by linearizing the equations about the flat rest state: *z*_{0} = 0, *σ* = 1, , *γ*_{t} = *τθ*_{αα}. Because *H*[- sin(*kα*)] = cos(*kα*), the Fourier modes satisfy . . Thus, the solution of the linearized problem with initial conditions [**7**] is when *k* is even, and [8]when *k* is odd. Here and are the angular frequency and period of the *k*th Fourier mode.

If we linearize *F* about the flat rest state (*c*_{k} = 0) with any period *T* > 0, we find from Eq. **8** that the Jacobian of *F* with respect to the *c*_{k} is an infinite diagonal matrix *J* (indexed by positive odd integers) with entries *J*_{kk} = cos(*ω*_{k}*T*/4), while . A necessary condition for a bifurcation to occur is that *J* have a nontrivial kernel, i.e., there must exist *j*, *k* odd and positive so that *T* = *jT*_{k}. Fixing *T* (i.e., *k* and *j*), the other entries of *J* satisfy [9]Thus, the kernel of *J* is infinite dimensional (as *A*_{kn2} = 0 for odd *n*) and the range of *J* is not closed (as *A*_{m} accumulates at zero, being uniformly distributed (34) over [0, 1]).

Both of these properties prevent a rigorous bifurcation analysis of solutions of *F* = 0 via the Liapunov–Schmidt reduction (35). Nevertheless, in spite of zero and small divisors, our numerical method has no difficulty finding time-periodic solutions. We use a solution of the linearized problem (with *k* = 1) as a starting guess for our optimal control algorithm to find a solution of the nonlinear problem near the flat rest state. We then use numerical continuation to increase the amplitude beyond the realm of linear theory. The continuation algorithm consists of varying one of the Fourier modes *c*_{k0} in [**7**] of the initial vortex sheet strength, *γ*_{0}, and solving for the other *c*_{k} and *T* to minimize the deviation from time-periodicity, , defined in Eq. **14** below. For each new value of *c*_{k0}, we use linear extrapolation from two previously computed solutions as a starting guess for the remaining *c*_{k}.

In Fig. 1, we show the result of varying *c*_{1} from zero (the flat rest state) to a turning point at -1.08207, and then back up to -0.8321. The solution labeled A on the diagram remains qualitatively similar to the linearized solutions [**8**] with *k* = 1, but higher frequency Fourier modes become increasingly significant as we continue along the bifurcation curve. This diagram contains 1,704 time-periodic solutions, each computed down to *G* ≈ 10^{-24}, with the number of Fourier modes, *M*, ranging from 32 to 512. The simulations took 4 weeks running simultaneously on five machines with a total of 32 cores (running OpenMP, a shared memory parallel programming language, on each machine). Most of the running time was devoted to resolving the more complicated solutions beyond the turning point in the bifurcation curve and exploring near-bifurcation events (described below). The part of the curve connecting the flat rest state to the point labeled A contains 439 solutions with *G* ≈ 10^{-30}, but only took 4 h to compute on an eight core machine. A few of these solutions were recomputed in double–double precision arithmetic to *G* ≈ 10^{-63} to be sure the algorithm continues to converge when roundoff error is decreased.

We interpret the turning point as a transition from *c*_{1} being the dominant mode to *c*_{3} being the dominant mode. In fact, we used *c*_{3} as the bifurcation parameter to traverse this region of the curve. As shown in Fig. 2 (ignoring side branches), *c*_{3} decreases monotonically through the turning points in *T* and *c*_{1}. As we continue along the curve, the solutions develop visibly active secondary oscillations superimposed on the main carrier wave. In some cases, the vortex sheet briefly overturns before returning to its initial flat state.

We noticed small wobbles in some of the plots of *c*_{k} versus *T*. By refining the stepsize in the continuation algorithm near each wobble, we discovered that these curves actually consist of several disjoint branches. The 13th Fourier mode *c*_{13} of the initial vortex sheet strength gives a particularly nice representation of the “near-bifurcation” events that separate the various branches. As shown in Fig. 3, these near-bifurcations appear as perturbed pitchforks (35). To our surprise, numerical continuation of the side branches from one of the pitchforks led to reconnections with the side branches of another pitchfork. One of the branches appears to be a closed loop.

Bifurcation diagrams of still higher Fourier modes reveal additional near bifurcations not visible to the first, third, or 13th mode. We illustrate this with the 43rd mode in Fig. 4. A sudden jump in the curve indicates a transition to a new branch of solutions. Following side branches of similar anomalies in lower Fourier modes led to the four branches shown in Figs. 1–3. We stopped following the side-branches of the last two pitchforks in Fig. 3 (and did not follow the new side branches in Fig. 4) as the running time grew to more than a day per data point (running on eight cores).

We are confident that the disconnection of bifurcation branches is a true feature of solutions of the PDE rather than a numerical artifact; the curves remain identical (to 9–10 digits of accuracy) if we cut the mesh size in half. We also emphasize that the same simulations are shown in all four figures; the additional bifurcations visible in the 43rd mode are a result of looking at a higher-frequency mode, not a result of running the simulations with a smaller mesh size.

We suspect that the disconnection of the bifurcation curves is related to resonances and small divisors. In previous studies of nonlinear wave equations (17, 18), it was found that periodic solutions may not occur in smooth families—their existence could only be established for values of a parameter in a Cantor set. We seem to be observing exactly this phenomenon. The remarkable thing is that low-frequency modes are mostly determined by their interaction with each other; a sudden jump in a high-frequency mode has little effect. This is why it is possible to compute these solutions numerically.

## Numerical Method

We now describe our algorithm for computing time-periodic solutions of the vortex sheet with surface tension. For the symmetric solutions studied in this paper, *z*_{0}(*t*) remains zero for all time, so we drop it from the equations in the interest of brevity. Let *q* = (*σ*,*θ*,*γ*) and define the inner product [10]We adapt the small-scale decomposition (SSD) algorithm (6, 7) from the multistep framework to the additive Runge–Kutta framework and write the vortex sheet system in the form [11]where Here is the Hilbert transform, which has symbol . We have desingularized the Birkhoff–Rott integral by writing *U* = *U*_{1} + *U*_{2} with [12]We have suppressed the dependence of *γ* and *z* on time in the notation. The idea behind the decomposition [**11**] is to treat the nonlinear operator *f*_{2}(*q*) explicitly and the linear operator *f*_{1}(*q*), which is the source of stiffness, implicitly. This is done using two *s*-stage Butcher arrays (36), one for *f*_{1} and another for *f*_{2}. In the more general case that *f*_{1} and *f*_{2} depend on time (e.g., in the adjoint system described below), we define two sets of stage derivatives and an update step via [13]Here *h* = Δ*t* is the timestep, spatial derivatives and the Hilbert transform are computed via the fast Fourier transform (FFT), and multiplications are done in physical (as opposed to Fourier) space. The trapezoidal rule is used to evaluate the integral in Eq. **12**, using *K*(*α*,*α*) = *z*_{αα}(*α*)/2*z*_{α}(*α*). We do not simplify *z*_{αα}(*α*)/2*z*_{α}(*α*) = (*i*/2)*θ*_{α}(*α*) as this identity only holds to *O*(*h*^{2}) in internal Runge–Kutta stages. (The final Runge–Kutta update is nevertheless fourth order, i.e., *O*(*h*^{5}).) The Butcher array for *f*_{1} is diagonally implicit (*a*_{ij} = 0 for *i* < *j*), whereas that for *f*_{2} is explicit ( for *i* ≤ *j*). This allows the stage derivatives to be solved for in order: *k*_{1},*ℓ*_{1},…,*k*_{s},*ℓ*_{s}. In our code, we used the six-stage fourth-order scheme ARK4(3)6L[2]SA described in ref. 31. If *f*_{2} = 0, this scheme is stiffly accurate (36), and hence *L* stable.

Next we define a functional *G*(*q*_{0},*T*) of the initial conditions and supposed period that is zero if and only if the solution is time periodic. Following previous work on the Benjamin–Ono equation (28, 29), we could define , where *q* solves Eq. **11** with initial condition *q*_{0}. Instead, to achieve a factor of four improvement in speed and to emphasize that the method will work for any two-point boundary value problem (beyond the computation of time-periodic solutions), we define [14]where *γ* is the third component of *q*, which satisfies Eq. **11** with initial conditions *q*(0) = *q*_{0} to be determined. As in [**7**], we take *q*_{0} of the form *σ*_{0} = 1, *θ*_{0} ≡ 0, and , . We note that *T* is now one-quarter of the period, which is our convention in this section only.

We vary *T* and the *c*_{k} in [**7**] to minimize *G* using an arbitrary precision C++ version of the limited memory BFGS algorithm (25) we wrote for this project. BFGS is a quasi-Newton line search algorithm that builds an approximate (inverse) Hessian matrix from the sequence of gradient vectors it encounters during the course of the line searches. In our continuation algorithm, we initialize the approximate Hessian with that of the previous minimization step (rather than the identity matrix), which leads to a tremendous reduction in the number of iterations required to converge (by factors of 10–20 in many cases). We use the limited memory feature of the code for the opposite reason it was originally intended: We store twice as many Hessian updates as there are columns in the matrix before cyclically overwriting them, which gives the algorithm more time to achieve superlinear convergence in the final iterations. The cost of the linear algebra in the BFGS algorithm is dwarfed by the PDE solves required to compute *G* and ∇*G*, so there is no benefit to using fewer Hessian updates. On the other hand, using more than twice as many columns does not seem to improve convergence rates.

It remains to explain how to compute ∇*G*, which is needed by the BFGS algorithm. The *T* derivative is easily found by evaluating using the trapezoidal rule. Both quantities *γ*(·,*T*) and *γ*_{t}(·,*T*) are already known from solving Eq. **11**. One way to compute with *k* a positive, odd integer would be to define and solve the variational equation [15]with initial conditions to obtain in [16]To avoid the expense of solving Eq. **15** repeatedly (for each value of *k*), we solve a single adjoint PDE to find the function such that [17]Here are adjoint variables, whereas is the *k*th Fourier series coefficient of . The function is chosen so that [18]is independent of *t*. When *t* = *T*, we put so that [**18**] is equal to in Eq. **16**. When *t* = 0 in [**18**], we recover Eq. **17**. A sufficient condition for [**18**] to remain constant may be obtained by differentiation. This yields the adjoint equation [19]where *s* = *T* - *t* denotes “reversed” time. Like the variational equation, Eq. **15**, the adjoint equation is linear and nonautonomous due to the presence of the solution *q*(*t*) in the equation. Note that Eq. **19** only needs to be solved once to obtain all the derivatives ∂*G*/∂*c*_{k} simultaneously (after one additional FFT in Eq. **17**). Thus, ∇*G* can be computed in approximately the same amount of time as *G*.

To solve the adjoint equation numerically in the additive Runge–Kutta framework, the values of *q*(·,*T* - *s*) are needed between timesteps (due to *τ*_{i} and in Eq. **13**), and a small-scale decomposition must be chosen. We use cubic Hermite interpolation to compute *q* at these intermediate times, having stored *q* and *q*_{t} at each timestep when Eq. **11** was solved. This is enough to achieve fourth-order accuracy in the adjoint problem. Our SSD algorithm is described below.

Due to the presence of singular integrals in Eq. **11**, the variational and adjoint equations are rather complicated. To write down the adjoint equation, Eq. **19**, we must first find formulas for in Eq. **15**. This requires the intermediate quantities , , , and to be computed. As always, a dot indicates a directional derivative with respect to *q* in the direction. From Eq. **6**, we have [20]where all factors to the right of a projection are multiplied before applying the projection. Next, from Eq. **12**, we obtain [21]The last term is found by writing *K*(*α*,*β*) in Eq. **12** as an *α*-derivative and interchanging the order of differentiation when the dot is applied. As *β* → *α*, the derivative of the term in brackets approaches , so it is not a singular integral. Next, from *U* = *U*_{1} + *U*_{2}, , and , we obtain [22]It then follows from Eq. **11** that [23]where *V*_{2} = *V* - *V*_{1} and .

Our final task is to identify the adjoint operator *Df*(*q*)^{∗}. Eqs. **20**–**23** can be combined into a composition of linear operators, *Df*(*q*) = *ABC*, where [24]We then have *Df*(*q*)^{∗} = *C*^{∗}*B*^{∗}*A*^{∗}. When computing adjoints, the middle two spaces in [**24**] are treated as real inner product spaces with the imaginary component of the last entry acting as another real dimension, e.g., Multiplication of by *i* is interpreted as a rotation by 90° in this real vector space. From Eq. **20**, we obtain [25]Without parentheses, operators and multiplication are always resolved right to left in our formulas. Similarly, from Eq. **21**, we find that where in Eq. **21**. Note that as *P*_{0} is not enclosed in parentheses in the formula for , we multiply (*U*_{1} - *iV*_{1}) by *w* before applying *P*_{0} and then taking the real part. Next, we seek such that [26]for all sufficiently smooth test functions *w*(*α*) in *L*^{2}(0,2*π*). Here and Δ*z* are shorthand for and *z*(*α*) - *z*(*β*), respectively. As it stands, the singularity in cot(Δ*z*/2) as *β* → *α* is cancelled by . However, we must separate from to achieve the desired form on the left-hand side of Eq. **26**, which gives rise to singular integrals. One approach is to write with *K* as in Eq. **12** to convert the singular part of the integral [**26**] into a Hilbert transform before separating . Instead, we use the fact that remains constant if *α* and *β* are interchanged. Thus, Eq. **26** may be written We convert this to a principal value integral over the region with *ε* → 0, and then use the fact that (⋆) changes sign when *α* and *β* are interchanged to conclude that may be replaced by . Because *γ*(*α*) is real valued, we get the desired formula We evaluate this integral numerically using the trapezoidal rule [37]where a subscript *k* indicates evaluation at one of the grid points *α*_{k} = 2*πk*/*M* (0 ≤ *k* < *M*), and primes are *α*-derivatives computed via the FFT. The *j* = *k* term is the *ε* → 0 limit of the average of the two values of the integrand at *β* = *α*_{k} ± *ε*, weighted by 1/*σM*. The same formulas are obtained if the integrand is desingularized before applying the trapezoidal rule; hence, the method is spectrally accurate.

The operator *A* in [**24**] may be found by combining Eqs. **22** and **23**. Although tedious, the procedure of forming *A* and computing the adjoint *A*^{∗} term by term is routine. The result is given in Fig. 5. The terms in boxes are separated from the rest and treated implicitly in the Runge–Kutta method; these terms propagate through *B*^{∗} and *C*^{∗} unaltered. Note that, although *DF*(*q*)^{∗} is linear, a fully implicit approach is impractical as the full operator cannot be inverted via the FFT.

## Acknowledgments

This work was supported in part by the National Science Foundation through Grant DMS-0926378 (to D.M.A.) and by the Director, Office of Science, Computational and Technology Research, US Department of Energy under Contract DE-AC02-05CH11231 (to J.W.).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: wilken{at}math.berkeley.edu.

Author contributions: D.M.A. and J.W. designed research; J.W. performed research; D.M.A. and J.W. contributed new reagents/analytic tools; J.W. analyzed data; and D.M.A. and J.W. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/cgi/content/full/0910830107/DCSupplemental.

## References

- ↵
- ↵
- Ambrose D,
- Masmoudi N

- ↵
- Iguchi T,
- Tanaka N,
- Tani A

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Brush SG

- Poincaré H

- ↵
- Arnold VI

- ↵
- Li YC

- ↵
- ↵
- Kawahara G,
- Kida S

- ↵
- ↵
- ↵
- Bourgain J

- ↵
- ↵
- ↵
- Bristeau MO,
- Pironneau O,
- Glowinsky R,
- Periaux J,
- Perrier P

- ↵
- ↵
- ↵
- Mohammadi B,
- Pironneau O

- ↵
- Nocedal J,
- Wright SJ

- ↵
- ↵
- Stoer J,
- Bulirsch R

- ↵
- Ambrose DM,
- Wilkening J

- ↵
- ↵
- ↵
- ↵
- ↵
- Saffman PG

- ↵
- Kuipers L,
- Niederreiter H

- ↵
- Golubitsky M,
- Schaeffer DG

- ↵
- Hairer E,
- Norsett SP,
- Wanner G

## Citation Manager Formats

## Sign up for Article Alerts

## Jump to section

## You May Also be Interested in

### More Articles of This Classification

### Physical Sciences

### Related Content

- No related articles found.

### Cited by...

- No citing articles found.