A Supplementary Material: Objective function

# Optimizing snake locomotion in the plane. I. Computations

## Abstract

We develop a numerical scheme to determine which planar snake motions are optimal for locomotory efficiency, across a wide range of frictional parameter space. For a large coefficient of transverse friction, we find that retrograde traveling waves are optimal. The optimal snake deflection scales as the -1/4 power of the coefficient of transverse friction, in agreement with an asymptotic analysis. At the other extreme, zero coefficient of transverse friction, we propose a triangular direct wave which is optimal. Between these two extremes, a variety of complex, locally optimal motions are found. Some of these can be classified as standing waves (or ratcheting motions).

Keywords: snake; friction; sliding; locomotion; optimization.

## I Introduction

Snake locomotion has a long history of study by biologists, engineers, and applied mathematicians ([6]; [1]; [2]; [3]; [4]; [5]). A lack of appendages makes snake motions simpler in some respects than those of other locomoting animals ([7]). However, snakes can move through a wide range of terrestrial ([8]; [6]; [9]), aquatic ([10]), and aerial ([11]) environments, with different kinematics corresponding to the different dynamical laws which apply in each setting. Here we focus on terrestrial locomotion, with motions restricted to two dimensions, and a Coulomb frictional force acting on the snake. Even with these simplifying assumptions, there are many possible snake kinematics to consider. One way to organize our understanding of a diversity of locomotory kinematics is to propose a measure of locomotory performance (here, efficiency), and determine which motions are optimal, and how their performance varies with physical parameters. Well-known examples are optimization studies of organisms moving in low- ([12]; [13]; [14]; [15]; [16]; [17]) and high-Reynolds-number fluid flows ([18]; [19]; [20]; [21]; [22]; [23]).

In this work, we use a recently-proposed model for snake locomotion in the plane which has shown good agreement with biological snakes ([2]; [5]). The snake is a slender body whose curvature changes as a prescribed function of arc-length distance along the backbone, and time. Its net rotation and translation are then determined by Newton’s laws, with Coulomb friction acting between the snake and the ground. This model is perhaps one of the simplest which can represent arbitrary planar snake motions. Previous studies have used the model to find optimal snake motions when the motions are restricted to sinusoidal traveling waves ([5]), or the bodies are composed of two or three links ([24]). These motions are represented by a small number of parameters (2-10); at the upper end of this range, the optimization problem is difficult to solve using current methods. In ([1]), Guo and Mahadevan developed a model which included the internal elasticity and viscosity of the snake, and studied the effects of these and other physical parameters on locomotion for a prescribed sinusoidal shape and sinusoidal and square-wave internal bending moments. In the present work, we neglect the internal mechanics of the snake, and consider only the work done by the snake on its environment, the ground. However, we consider a more general class of snake motions. Hu and Shelley assumed a sinusoidal traveling-wave snake shape, and computed the snake speed and locomotory efficiency across the two-parameter space of traveling-wave amplitude and wavelength for two values of frictional coefficients ([5]). They also considered the ability of the snake to redistribute its weight during locomotion. Jing and Alben considered the locomotion of two- and three-link snakes, and found analytical and numerical results for the scaling of snake speed and efficiency ([24]). Among the results were the optimal temporal function for actuating the internal angles between the links, expressed as Fourier series with one and two frequencies.

Here we address the optimization problem for planar snakes more generally, by using a much larger number of parameters (45 in most cases) to represent the snake’s curvature as a function of arc length and time. Although our snake shape space is limited, it is large enough to represent a wide range of shapes. To solve the optimization problem in a 45-dimensional shape space, and simultaneously across a large region of the two-dimensional friction parameter space, we develop an efficient numerical method based on a quasi-Newton method. The simulations show clearly that two types of traveling-wave motions, retrograde and direct waves, are optimal when the coefficient of friction transverse to the snake is large and small, respectively. An asymptotic analysis in [25] shows how the snake’s traveling wave amplitude should scale with the transverse coefficient of friction when it is large. The analysis gives a power law for the amplitude which agrees well with the numerical results shown here. When the transverse coefficient of friction is zero, we give another optimally-efficient motion here which is a direct wave of deflection along the body. Between these extremes, the numerics show a more complicated set of optima which include standing wave (or ratcheting) motions, among a wide variety of other locally-optimal motions. Taken together, these results show which planar snake motions are optimal within a large space of possible motions and frictional coefficients.

## Ii Model

We use the same frictional snake model as ([2]; [5]; [24]), so we only summarize it here. The snake’s position is given by , a planar curve which is parametrized by arc length and varies with time . A schematic diagram is shown in figure 1.

The unit vectors tangent and normal to the curve are and respectively. The tangent angle and curvature are denoted and , and satisfy , , and . We consider the problem of prescribing the curvature of the snake as a function of time, , in order to obtain efficient locomotion. When is prescribed, the tangent angle and position are obtained by integration:

 θ(s,t) =θ0(t)+∫s0κ(s′,t)ds′, (1) x(s,t) =x0(t)+∫s0cosθ(s′,t)ds′, (2) y(s,t) =y0(t)+∫s0sinθ(s′,t)ds′. (3)

The trailing-edge position and tangent angle are determined by the force and torque balance for the snake, i.e. Newton’s second law:

 ∫L0ρ∂ttxds =∫L0fxds, (4) ∫L0ρ∂ttyds =∫L0fyds, (5) ∫L0ρX⊥⋅∂ttXds =∫L0X⊥⋅fds. (6)

Here is the snake’s mass per unit length and is the snake length. The snake is locally inextensible, and and are constant in time. is the force per unit length on the snake due to Coulomb friction with the ground [2]:

 f(s,t)=−ρgμt(ˆ∂tX⋅^n)^n−ρg(μfH(ˆ∂tX⋅^s)+μb(1−H(ˆ∂tX⋅^s)))(ˆ∂tX⋅^s)^s. (7)

Here is the Heaviside function and the hats denote normalized vectors. When we define to be . According to (7) the snake experiences friction with different coefficients for motions in different directions. The frictional coefficients are , , and for motions in the forward (), backward (), and transverse (i.e. normal) directions (), respectively. In general the snake velocity at a given point has both tangential and normal components, and the frictional force density has components acting in each direction. A similar decomposition of force into directional components occurs for viscous fluid forces on slender bodies ([26]).

We assume that the snake curvature is a prescribed function of and that is periodic in with period . Many of the motions commonly observed in real snakes are essentially periodic in time ([2]). We nondimensionalize equations (4)–(6) by dividing lengths by the snake length , time by , and mass by . Dividing both sides by we obtain:

 LμfgT2∫10∂ttxds =∫10fxds, (8) LμfgT2∫10∂ttyds =∫10fyds, (9) LμfgT2∫10X⊥⋅∂ttXds =∫10X⊥⋅fds. (10)

In (8)–(10) and from now on, all variables are dimensionless. For most of the snake motions observed in nature, ([2]), which means that the snake’s inertia is negligible. By setting this parameter to zero we simplify the problem considerably while maintaining a good representation of real snakes. (8)–(10) become:

 b=(b1,b2,b3)⊤=0;b1 ≡∫10fxds, (11) b2 ≡∫10fyds, (12) b3 ≡∫10X⊥⋅fds. (13)

In (11)–(13), the dimensionless force is

 f(s,t)=−μtμf(ˆ∂tX⋅^n)^n−(H(ˆ∂tX⋅^s)+μbμf(1−H(ˆ∂tX⋅^s)))(ˆ∂tX⋅^s)^s (14)

The equations (11)–(13) thus involve only two parameters, which are ratios of the friction coefficients. From now on, for simplicity, we refer to as and as . Without loss of generality, we assume . This amounts to defining the backward direction as that with the higher of the tangential frictional coefficients, when they are unequal. may assume any nonnegative value. The same model was used in ([2]; [5]; [24]), and was found to agree well with the motions of biological snakes in ([2]).

Given the curvature , we solve the three nonlinear equations (11)–(13) at each time for the three unknowns , and . Then we obtain the snake’s position as a function of time by (1)–(3). The distance traveled by the snake’s center of mass over one period is

 d= ⎷(∫10x(s,1)−x(s,0)ds)2+(∫10y(s,1)−y(s,0)ds)2. (15)

The work done by the snake against friction over one period is

 W =∫10∫10f(s,t)⋅∂tX(s,t)dsdt. (16)

We define the cost of locomotion as

 η=Wd (17)

and our objective is to find which minimize at values which range widely over . In our computations, instead of , we minimize a different function which is more convenient for physical and numerical reasons:

 F=−dWe2cos(θ(0,1)−θ(0,0)). (18)

To obtain from we have substituted for , to avoid infinities in the objective function for a common class of motions with and of order 1. The exponential term in (18) penalizes rotations over a period, to ensure that the snake moves in a straight path rather than a circular path over many periods. In the Supplementary Material section A we give a more detailed explanation for the choice of .

## Iii Numerical Minimization Method

Our goal is to determine the snake shape, , which minimizes for a given . Since has -period one, we represent it via a double-series expansion:

 κ(s,t)=m1−1∑j=0n1−1∑k=0(αjkcos(2πjt)+βjksin(2πjt))Tk(s). (19)

Here is a Chebyshev polynomial in :

 Tk(s)=cos(karccos(s)). (20)

The expansion (19) converges as for a large class of which includes differentiable functions, and with a convergence rate that increases rapidly with the number of bounded derivatives of ([27]). It is reasonable to hope that minimizing over a class of functions (19) with small and gives a good approximation to the minimizers with infinite and . Since the terms in (19) with coefficients are automatically zero, the functions given by (19) are a -dimensional space. In most of this work we use , so we are minimizing in a 45-dimensional space. We give a few results at large with for comparison, and find similar minimizers in this region of parameter space. The minimization algorithm converges in a smaller portion of space as and increase, and convergence is considerably slowed by the increased dimensionality of the search space. We find that is small enough to allow convergence to minima in a large portion of space, but large enough to approximate the gross features of a wide range of . Another reason to expect that good approximations to minimizers can be obtained with small is that the dynamical equations (11)–(13) depend only on spatial integrals of body position and velocity.

Using the dynamical equations (11)–(13), we show in Supplementary Material section B that , , and are the same for any reparametrization of time. Such a reparametrization can be used to reduce the high-frequency components of a motion while keeping the efficiency the same. Therefore, it is reasonable to expect that a good approximation to any minimizer of can be obtained with low temporal frequencies, that is, with an expansion in the form of (19) with small .

We use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm ([28]) to minimize numerically, by computing a sequence of iterates which converge to a local minimizer of . The algorithm requires routines for computing the values and gradients of as well as a line search method for finding the next iterate at each step, using the search direction provided within the BFGS algorithm. The BFGS algorithm is a quasi-Newton algorithm, and therefore forms an approximation to the Hessian matrix of second derivatives, without the expense and difficulty of computing the Hessian matrix explicitly.

We have developed an efficient method for evaluating and its gradient with respect to the shape parameters with just a single simulation of the snake trajectory over one period. Computing the value of requires computing the trajectory of the snake over one period. We discretize the period interval uniformly with time points: . At each time point we write (11)–(13) as a nonlinear algebraic system of three equations in three unknowns: . To formulate the system, we set to at each time point, which means that we set the snake in an artificial frame rotated from the lab frame. is the same in the two frames, so after solving for in the artificial frame, we integrate in time to obtain in the lab frame. This function is the angle between the artificial frame and the lab frame at each time, and we use it to rotate to the lab frame. We then integrate in the lab frame to obtain , and integrate using to obtain the body trajectory a posteriori. Using as the unknowns instead of has two advantages: we avoid cancellation error associated with computing discrete time derivatives of , and we solve decoupled systems of three equations in three unknowns instead of a large coupled system of equations in unknowns ( at all times points), which is more expensive to solve.

We solve the system (11)–(13) using Newton’s method with a finite-difference Jacobian matrix. We solve it at each time sequentially from , using a random initial guess for at and an extrapolation from previous time points as an initial guess at other . An important element of our method is the quadrature used to compute the integrals in (11)–(13) using (14). Here we introduce notation for the tangential and normal components of the snake velocity:

 us=∂tX⋅^s;un=∂tX⋅^n. (21)

In (14) we can write

 ˆ∂tX⋅^s=us√u2s+u2n;ˆ∂tX⋅^n=un√u2s+u2n. (22)

These terms have unbounded derivatives with respect to and when . The local error when integrating (14) using the trapezoidal rule (for example) with uniform mesh size is . To obtain convergence in this case it is necessary to have , so the trapezoidal rule (and other classical quadrature rules) needs to be locally adaptive when becomes small.

We use a different approach which allows for a uniform mesh even when . We perform the quadrature analytically on subintervals using locally linear approximations for and . Using and , we compute , , and then and on a uniform mesh. On each subinterval, we approximate and as linear functions of the form using their values at the interval endpoints. We approximate , , and as constants using their values at the midpoints of the intervals. Then the integrals in (11)–(13), using (14), are of the form

 C0∫baAs+B√(As+B)2+(Cs+D)2ds (23)

on each subinterval. We evaluate such integrals analytically using

 I1 Missing or unrecognized delimiter for \left (24) I2 ≡∫10sds√s2+a2s+a3=−√a3+√1+a2+a3−a22I1. (25)

When used to evaluate (23), the logarithms in (24)–(25) have canceling singularities when the linear approximations to and are proportional (including the case in which one is zero), i.e. when . When this occurs, we use different formulae because the integrands are constants or sign functions, which are easy to integrate analytically. This quadrature method is for mesh size , due to the use of linear and midpoint approximations of the functions in the integrands.

With this method for evaluating (11)–(13), we compute by Newton’s method as described above, and obtain the snake trajectory and , , and for the prescribed curvature function.

The BFGS method also requires the gradient of with respect to the curvature coefficients . For notation, let be a vector whose entries are the curvature coefficients. We compute a finite-difference gradient using

 Extra open brace or missing close brace (26)

where is the unit basis vector in the th coordinate direction. (26) requires evaluations of , or computations of the snake motion over a period. We employ a very accurate approximation which requires only a single computation of the snake motion over a period.

Computing requires solving , i.e. (11)–(13), for at each time using the vector of curvature coefficients . is a function of and the curvature coefficients, so we write it as and solve for at each time . The value of which satisifies these equations at each may be written . Having simulated the snake motion and evaluated when the curvature coefficient vector is , we have for each . solves , a set of equations whose terms are extremely close to (within about of) those of at each . Thus is an extremely good initial guess for , and we need only perform a single Newton iteration with this initial guess to obtain to machine accuracy (about ). The one-step Newton iteration is as follows. We define for our description. We start with the usual Taylor series expansion from which Newton’s method is obtained:

 0=b(ck(a+ϵj),a+ϵj)=b(ck(a),a+ϵj)+J––––(ck(a),a+ϵj)(ck(a+ϵj)−ck(a))+O(∥ϵj∥2). (27)

Here is the 3-by-3 Jacobian matrix of partial derivatives of with respect to its leading 3-by-1 vector argument:

 J––––(x,y)ij≡∂xjbi(x,y),i=1,2,3,j=1,2,3. (28)

Having previously simulated the snake motion with the curvature coefficient vector to evaluate , we have computed via a separate instance of Newton’s method, and we have also calculated and stored at each . We therefore modify in (27) to use this information, while retaining the same order of accuracy in :

 0=b(ck(a+ϵj),a+ϵj)=b(ck(a),a+ϵj)+J––––(ck(a),a)(ck(a+ϵj)−ck(a))+O(∥ϵj∥2). (29)

We solve (29) to obtain , our approximation to using a single Newton iteration:

 cjk≡ck(a)−J––––(ck(a),a)−1b(ck(a),a+ϵj)=ck(a+ϵj)+O(∥ϵj∥2). (30)

(30) requires only a single computation of the snake motion over a period, which gives and at each . It also requires evaluations of the function . These are much less expensive than even a single computation of the snake motion over a period. The error in is of the same order of magnitude as machine precision. The value of the cost function computed from is approximately within machine precision of (which corresponds to ), and therefore contributes an error of the same order as the finite-difference error in (26) when it is used in place of in (26). This level of error leads to essentially no deterioriation in the convergence of BFGS up to machine precision, relative to the version with the exact gradient ([29]). Therefore we obtain the gradient of for approximately the cost of computing . Each computation of a search direction in BFGS requires only one gradient computation. These ingredients give an efficient optimization method, which is critical to finding optima in the high-dimensional space of curvature coefficients.

## Iv Large-μt results

In Figure 2 we show the result of one run of the optimization routine with and , starting with a randomly-chosen initial vector of curvature coefficients. In panel a we plot a map of curvature values over one period, and find that the curvature approximates a sinusoidal traveling wave with wavelength approximately equal to the snake body length. In this case the waves appear to travel with nearly constant speed, shown by the straightness of the curvature bands. The bands deviate from straightness near the ends, where the largest curvatures occur. In panel b we plot the objective function and of the 2-norm of its gradient over the iteration count of the optimization routine. We find that the gradient norm decreases by four orders of magnitude within the first 100 iterations, and the objective function reaches a plateau at the same time. In panel c we give snapshots of the position of the snake at equal instants in time, moving from left to right, and see that the snake approximately follows a sinusoidal path.

Figure 3 shows three different converged results for the same parameters and but with different randomly-chosen initial iterates. Again we have traveling waves of curvature shown by bands, but unlike in Figure 2a, the bands are not straight. Also, the bands travel from right to left instead of left to right. The bands are not straight because the traveling wave speed changes over time. However, at any given time, the distances in from the maximum to the minimum of curvature are nearly the same in Figure 3a–c and Figure 2a, and so are the shapes. In all four cases, the snake passes through the same sequence of shapes, but at different speeds, and different (integer) numbers of times in one period. Thus the four motions are essentially the same after a reparametrization of time. Correspondingly, the values of for all four are in the range . The distances traveled, , are within 0.2% of integer multiples of 0.597, with the multiples 4,3,2, and 1 for the motions in Figure 2a and Figures 3a–c, respectively, which correspond to the number of times they move through the sequence of shapes represented by Figure 3c. The values of are the same whether the waves move from left to right or right to left in , because (or , so there is a symmetry between forward and backward motions.

Over an ensemble of 20 random initializations, this sinusoidal motion yields the lowest value of among the local optima found. We have repeated the search with twice as many spatial and temporal modes ( = 10, = 10) and 48 random initializations. About three-quarters of the searches converge, and the majority of the converged states are very close to those shown in Figures 2 and 3. The values of for these states lie in the range . The mean is systematically lower than the mean for = 5, = 5 by about 2.5%. The decrease in the objective function is slight, and the optimal shape is fairly smooth, which indicates that the = 5, = 5 minimizer is representative of the minimizer for much larger and , at least for and . With = 10 and = 10, the curvature bands of Figures 2a and 3 become straighter near the ends, showing the effect of the truncation in the Chebyshev series. However, it is more difficult to obtain convergence with four times as many modes, because of the higher dimensionality of the space. Since a primary goal here is to obtain the optima across as broad a range of () space as possible, we proceed with = 5 and = 5.

Performing searches with and a wide range of , we find that similar sinusoidal traveling-wave optima occur for . In Figure 4a we plot the values of corresponding to the optima over a range of on a logarithmic scale, and find that the optima are more efficient at higher . In panel b, the curvatures of the optimal shapes are plotted for the same , at the times when the curvatures have a local maximum at . The corresponding snake body positions are shown in panel c, and are vertically displaced to make different bodies easier to distinguish. In all cases the bodies approximate sinusoidal shapes with slightly more than one wavelength of curvature on the body. The curvature amplitudes decrease monotonically with increasing .

Performing searches with a wide range of and now, we find that similar sinusoidal traveling-wave optima are found for and all used (). For these motions, the curvature traveling waves move backward in , and the tangential motion is purely forward, so they have the same at all . However, for , a separate class of (local, not global) optima are also found. These have curvature traveling waves which move forward in , and for which the tangential motion is purely backward. Since friction is higher in the backward direction, these local optima have higher (are less efficient) than the oppositely-moving optima, and the values of for these optima increase with increasing , for the same reason.

In figure 5a we plot the values of for these “Reverse” sinusoidal motions at various and , with leftward-pointing triangles. The triangles connected by a given solid line are data for various at a particular value of , labeled to the right of the line. When , the reverse sinusoidal motions have the same values as the forward sinusoidal motions shown in figure 4. These have the lowest values across all for , and are thus the most efficient motions we have found.

Three other sets of symbols are used to denote other classes of local optima found by our method. All are less efficient than the forward sinusoidal motions, but we present examples of them to show some of the types of local optima that exist for this problem. For these three sets, shown by circles, downward pointing triangles, and plus signs, we do not join the optima with a given by lines, to reduce visual clutter. The circles, labeled “Curl,” give values for shapes which are like sinusoidal traveling waves, but they have a curl—the snake body winds around itself once in a closed loop. Snapshots of the snake motion for such a case are shown in figure 5b. Here the snake moves horizontally, but the snapshots are displaced vertically to make them distinguishable. The arrow shows the direction of motion. The tail of the snake is marked by a circle. The snake body moves rightward while the curvature wave and curl both move leftward. Because the snake overlaps itself, this motion is not physically valid in 2D. However, a large variety of such optima are found numerically, with various small-integer numbers of curls (and with both directions of rotation). The efficiency of such motions decreases as the number of curls increases. The downward-pointing triangles, labeled “Curl, Reverse,” are for traveling waves with curls which move in the forward direction along the snake. These have a yet smaller efficiency for the same reason that the reverse sinusoidal motions are less efficient than the forward sinusoidal motions, in the absence of curls. The plusses, labeled “Other,” are for other local optima, which have a variety of motions. Panel c shows snapshots of a curling (but not self-intersecting) motion in which the snake moves horizontally leftward. The horizontal displacement is increased by 8/9 of a snake length between snapshots to make them distinguishable. Panel d shows a different curling motion in which the snake moves horizontally rightward, and the snapshots are displaced vertically to make them distinguishable. These alternative minima increase greatly in number as decreases below 6, and some of them outperform the forward sinusoidal motions at these transitional values. Before discussing the regime of below 6, we describe some scaling properties for the traveling-wave solution to the problem in the large- regime. In a separate paper [25], we derive analytically the form of the optimal prescribed curvature in the limit of large . In that paper, we calculate that the optimal curvature function is a periodic traveling wave with any functional form, but in the limit that the wavelength of the traveling wave tends to 0, and such that the snake deflection has a root-mean-square amplitude of . The numerically-determined optima in figure 4 do not have small wavelengths, however, due to the finite number of modes used. This is explained further in Supplementary Material section C.

We compare the numerically-found optima to the theoretical results of [25] in figure 6. Here we take the plots of figure 4 and transform them according to the theory. Figure 6a transforms from figure 4a into by the formula , since rotations over one period are negligible for the traveling-wave solutions. In figure 6a we plot versus for the numerical solutions (squares) from figure 4a, along with the small-wavelength theoretical result (solid line), which is

 η→1+√2μ−1/2t+O(μ−1t). (31)

The numerics follow the same scaling as the theory, but with a consistent upward shift. Some degree of upward shift is expected from the finite-mode truncation, which makes the numerical result underperform the analytical result. The uniformity of the shift may be due to the shape consistently assumed by the numerical optima for various . For we have also computed via the dynamical simulation for a much shorter-wavelength shape given by with wavelength and chosen so that the snake deflection has a root-mean-square amplitude of . The result is given by the open circle, which outperforms the optimum identified in the truncated-mode space, as expected. The circle is 0.013 above the solid line. This deviation is due in part to the neglect of the next term in the asymptotic expansion, proportional to or 1/300 here. Numerically, it is difficult to simulate motions with much smaller than 1/8. Figure 6b shows the curvature, rescaled by to give a collapse according to [25]. The data collapse well compared to the unscaled data in 4b. Figure 6c shows the body shapes from 4c with the vertical coordinate rescaled, and the shapes plotted with all centers of mass located at the origin. We again find a good collapse, consistent with the curvature collapse.

## V Moderate to small μt: from retrograde waves to standing waves to direct waves

Near there is a qualitative change in the numerically-computed optima. At values of , the global minimizers of all have backward-moving traveling waves which propel the snake forward. In addition, there are local minimizers in which the waves move forward along the snake (and the snake moves backward), and overlapping solutions with integer numbers of curls appearing in the traveling waves. The values of at the local minima, shown in figure 5, are generally well-separated. In addition, there are a small number of local minima with ratcheting types of motions. For in figure 5, there are an increasing number of local minimizers, denoted “Other,” quite close to the global minimizer. These motions are perturbed traveling waves. The number of local minima increases greatly as decreases below 6, and their values are no longer well-separated.

Figure 7 shows the values of for all computed local optima at four values of near the transition: 2, 3, 6, and 10. For each , the values are plotted against . Also, symbols are used to distinguish whether the local maximum in curvature travels forward or backward on average, and hence gives the direction of the curvature traveling wave (or perturbed traveling wave). At we find that the smallest is nearly constant with respect to . The same is true for except at the largest , for which there are several closely spaced values of which are slightly below the minimal at smaller . These values mark the first appearance of the perturbed traveling waves as decreases. For , the distribution of values is quite different, with a much stronger decrease in as increases. There are many clusters of closely-spaced but distinct values of , showing a large increase in the number of local equilibria. At each the smallest values of are spread over a band of about 0.5 in width. At , the same qualitative trend holds, except that the decrease in with is much greater. There is more scatter of values overall and of those representing forward-moving traveling waves. There is randomness in the sample of values, particularly evident for , because they are values for a finite sample of local minimizers starting from random initial coefficients.

Intuitively, the minima for represent perturbations of traveling wave motions towards “ratcheting”-type motions which involve small amounts of backward motion in a portion of the snake to “push” the rest of the snake forward. The appearance of backward motion explains the variation of with . As decreases from 2000 to 6, for the traveling waves increase, and for , other motions outperform the traveling wave motions at sufficiently large . For , the motions tend more toward ratcheting-type motions, and we will investigate this further subsequently.

Another way of viewing the transition is to consider the large- solution as decreases. In the snake motion of figure 2c, the snake “slips” tranversely to itself as it slithers forward (rightward). This transverse motion has a component in the backward (leftward) direction, and the resulting friction with the ground propels the snake forward as a thrust force from transverse friction, which balances the drag force from the forward component of tangential friction. The thrust force is proportional to , and increases as the slope of the snake increases, which increases the component of transverse friction in the horizontal direction. If is decreased, the thrust force can be maintained by increasing the slope of the snake. This is indeed observed in the sinusoidal optima as decreases, in figure 4c. When decreases to 3, the slope of the snake becomes nearly vertical at its maximum, and no further decreases in can be accommodated by increasing the slope of the snake. Thus, near this value of , motions different from the sinusoidal traveling wave become preferred.

In figure 8a we plot a function of the shape dynamics which quantifies the transition away from traveling wave motions near . The quantity is

 ψ[κ]≡mintmax0.2≤s≤0.8|κ(s,t)|maxtmax0.2≤s≤0.8|κ(s,t)|. (32)

For which is a traveling wave with a wavelength sufficiently small that a maximum or minimum in curvature occurs in at all times, and for which the maximum and minimum are equal in magnitude, equals unity. For , the optimal are approximately such traveling waves; we restrict away from the ends to avoid slight intensifications in curvatures near the ends. In figure 8a we plot the maximum of over all equilibria located at a given pair.

We find for . increases slightly from about 0.8 to 0.9 as increases from 6 to 200. For = 6, decreases slightly but noticeably, from 0.8 to 0.7, as exceeds 200, showing that at transitional , the perturbations away from the traveling wave solutions appear first at the largest . At = 3, the values of have decreased to around 0.6 for most ; is as large as 0.9 at smaller and as small as 0.4 at larger . For , the values of are generally much smaller. There is considerable randomness in the values of plotted for because they are the maxima over samples of local equilibria starting from 10-20 random initial coefficient vectors.

In figure 8b we plot the values of for the optima whose -values are shown in panel a. We find that varies smoothly over the region of transitional , , attaining its highest values there. The retrograde traveling-wave motions and other motions are similarly efficient at a transitional , and the other motions become preferable as decreases.

A global picture of the optimal states is presented in the “phase diagram” of figure 9. At large , the retrograde traveling waves are optimal. At small , all of the optima have the form of direct traveling waves. Unlike for the large- waves, these direct waves all have large amplitudes, and their specific form varies with . Near there is a region where standing waves with non-zero mean deflection (or ratcheting motions) are seen. Examples will be given shortly. These optima co-exist with retrograde-wave optima at larger . We note that there are other optima which do not fit easily into one of these categories, throughout the phase plane. These other optima are generally found less often by the computational search, and are generally suboptimal to the traveling waves, but there are exceptions. For the sake of brevity we do not discuss such motions here.

To obtain a clear picture of the snake kinematics across the transitional region shown in figure 9, we give figures showing snapshots of snake kinematics, together with spatiotemporal curvature plots, at = 1 and 3 in this section and at = 2, 5, and 20 in Supplementary Material section D. Each of these figures represents one horizontal “slice” through the phase diagram in figure 9.

In figure 10, we show the optimal snake kinematics for and eight spanning the transition. Each set of snapshots is arranged in a single column, moving downward as time moves forward over one period, from top to bottom (shown by arrows). Similarly to figure 5b and d, the snapshots are shown in a frame in which the net motion is purely horizontal, from left to right. However, the snapshots are given an extra uniform vertical spacing to make them visable (i.e. to prevent overlapping). Below the snapshots, the corresponding curvature plot is shown, in the same format as in figures 2 and 3, with increasing from 0 to 1 on the horizontal axis, from left to right, and increasing from 0 to 1, from bottom to top. We have omitted the axes, and the specific curvature values for each plot, instead using a single color bar (shown at bottom) with curvature values which range from the maximum to the minimum for each plot. By omitting the axes we are able to show several of the optimal motions side-by-side. By comparing a large number of motions in a single figure, it is easier to see the general trends.

We begin by discussing the rightmost column in figure 10, with . A retrograde traveling wave is clearly seen, so this optimum is a continuation of the large- solutions. At this , just above the transition, the wave amplitude is large, and the wave clearly moves backwards in the lab frame. Thus the snake slips backward considerably even though it moves forward, because is not very large. Moving one column to the left, the optimum at shows a clear difference. At certain instants, the curvature becomes locally intensified near the snake midpoint. This occurs when the local curvature maximum reaches the snake midpoint. The motion is still that of a retrograde traveling wave, but now the wave shape changes over the period. Hence the displacement is not a traveling wave of the form for a fixed periodic function . Comparing the curvature plots for and 3, we see that the bands have nearly constant width across time for . At , by contrast, the curvature bands become modulated. The bands show bulges near the curvature intensifications. Moving leftward again, to and then to 1.5, the curvature intensifications increase further, showing a general trend. In all cases, however, the curvature peaks move backward along the snake, even as the peak values change, so these motions are retrograde waves. At , a novel feature is seen: the endpoints of the snake curl up (but do not self-intersect), forming “anchor points.” Between these points, the remainder of the snake moves forward. In the curvature plot, bands can still be seen corresponding to retrograde waves. The three remaining columns, moving leftward, all show direct waves. That is, the wave moves in the same direction as the snake [30]. No optimum is shown for . The point , in which friction is equal in all directions, is in some sense a singular point for our computations. Near this point, it is difficult to obtain convergence to optimal motions, and optima are difficult to classify clearly as retrograde or direct waves (or another simple motion). Nonetheless, we have computed one optimum which can be classifed as a direct wave, but it has complicated and special features which are not easy to classify. In the three left columns, the direct waves can be recognized by the bands in the curvature plots. The bands have the direction of direct waves. However, the bands are not nearly as simple in form as those for the retrograde waves, at large .

Figure 11 shows optima for now increased to 3. The same transition is seen, but now in two middle columns, with and 1.5, standing waves are seen. The curvature plots are nearly left-right symmetric. In the corresponding motions, the snake bends and unbends itself. Because , there is a preference for motion in the forward direction even though the curvature is nearly symmetric with respect to . These optima did not occur for (figure 10) because symmetric motions with equal forward and backward friction would yield zero net distance traveled. It is somewhat surprising that such motions can be completely ineffective for and yet represent optima for as small as 2, as shown in figure 9.

We have presented only a limited selection of optima for moderate , and largely with pictures, because it is difficult to describe and classify many of the motions we have observed. Future work may consider this interesting region of parameter space further.

We now move on to the limit of zero . In this limit, we find a direct-wave motion with zero cost of locomotion. This explains the preference for direct waves at smaller in this section.

## Vi Zero μt: direct-wave locomotion at zero cost

When , it is possible to define a simple shape dynamics with zero cost of locomotion. The shape moves with a nonzero speed while doing zero work. It is a traveling wave with the shape of a triangle wave which has zero mean vertical deflection:

 y(s,t)=A0−∫sgn(sin(2π(s−t)))ds. (33)

The vertical speed corresponding to (33) is a square wave:

 ∂ty(s,t)=−A0sgn(sin(2π(s−t))). (34)

The wave of vertical deflection moves forward with constant speed one, and the horizontal motion is also forward, with a constant speed :

 ∂tx(s,t)=U (35)

We determine by setting the net horizontal force on the snake to zero:

 ∫ˆ∂tX⋅^ssxds=0. (36)

with

 ^s=⎛⎝√1−A20A0sgn(sin(2π(s−t)))⎞⎠ (37)

Solving (36) for we obtain

 U=A20√1−A20. (38)

With this choice of , is identically zero on the snake, so the motion is purely in the normal direction. Therefore, there are no forces and torques on the snake, and no work done against friction. However, the snake moves forward at a nonzero speed . This solution has infinite curvature at the peaks and troughs of the triangle wave, and we may ask if zero cost of locomotion () can be obtained in the limit of a sequence of smooth curvature functions which tend to this singular curvature function. To answer this question, we have set to the Fourier series approximation to (33) with wave numbers, for ranging from 2 to 1000. At each time , we have then calculated a net horizontal speed , net vertical speed and a net rotational speed such that , and . We find very clear asymptotic behaviors for , , , and :

 U(t)=A20√1−A20+O(n−1),V(t)=O(n−1),dRdt=O(n−1),η=O(n−1). (39)

Therefore zero cost of locomotion does in fact arise in the limit of a sequence of smooth shape dynamics. This type of traveling wave shape dynamics, which yields locomotion in the direction of the wave propogation, has been called a “direct” wave in the context of crawling by snails ([30]) and worms ([31]).

The triangular traveling wave motion is particularly interesting because it can represent an optimal motion at both zero and infinite , and its motion can be computed analytically for all . We show how this motion embodies many aspects of the general problem in Supplementary Material section E.

## Vii Conclusion

We have studied the optimization of planar snake motions for efficiency, using a fairly simple model for the motion of snakes using friction. Our numerical optimization is performed in a space of limited dimension (45), but which is large enough to allow for a wide range of motions. The optimal motions have a clear pattern when the coefficient of transverse friction is large. Our numerical results indicate that a traveling-wave motion of small amplitude is optimal, which agrees with an asymptotic analysis in [25]. When this coefficient is small, another traveling-wave motion, of large amplitude, is optimal. When the coefficient of transverse friction is neither large nor small, the optimal motions are far from simple, and we have only begun to address this regime.

A natural extension of this work is to include three-dimensional motions of snakes. In other words, one would give the snake the ability to lift itself off of the plane ([1]; [2]; [5]), as occurs in many familiar motions such as side-winding ([9]). Another interesting direction is to consider the motions of snakes in the presence of confining walls or barriers ([4]; [32]).

###### Acknowledgements.
We would like to acknowledge helpful discussions on snake physiology and mechanics with David Hu and Hamidreza Marvi, helpful discussions with Fangxu Jing during our previous study of two- and three-link snakes, and the support of NSF-DMS Mathematical Biology Grant 1022619 and a Sloan Research Fellowship.

Supplementary material for:
Optimizing snake locomotion in the plane. I. Computations

## Appendix A Supplementary Material: Objective function

To understand our choice of , we first describe the class of motions which result from periodic-in-time . Let and be given for the moment. Let us rotate by an angle and translate by a vector , uniformly over , so we have a rigid body motion. Then and are also uniformly rotated by . Let us also rotate uniformly by . Then by (14), is uniformly rotated by . Hence, if and are such that the force and torque balances (11)–(13) hold, they will continue to hold if we apply a rigid-body rotation and translation to and rotate by the same angle. Now, since , is obtained from by a rigid-body rotation by an angle and translation by a vector . In fact, the trajectory of the body during each period is the same as during the previous period but with a rigid translation by and rotation by . Over many periods such a body follows a circular path when , and because of the rotation, the distance traveled over periods may be much smaller than times the distance traveled over one period. To obtain a motion which yields a large distance over many periods, we restrict to body motions for which is obtained from by a translation only, with rotation angle . We implement this constraint approximately in our optimization problem, using a penalization term. We multiply in (17) by a factor which penalizes rotations over a single period:

 d→de2cos(θ(0,1)−θ(0,0)) (40)

(40) involves the change in angle at the snake’s trailing edge over a period, which is the same as . The exponential factor has a maximum of when , and a minimum of when . Using the factor of 2 in the exponent of (40), all of the optima found by our computations involve very little rotation—less than 0.5 degrees over all of space. When the factor is decreased from 2, some of the optima begin to involve nontrivial amounts of rotation. In short, the factor of 2 represents a choice of how much to weight lack of rotation relative to distance and work in the objective function. As described in the results section, for and all , we find optima which have an analytical traveling wave form, and in the limit of large , rotations tend to zero for these solutions. These optima are found even with no rotation penalty, so the rotation penalty plays a negligible role in this region of parameter space. For smaller , however, in the absence of the rotation penalty, many of the optima have significant rotations.

A second modification to is made to avoid large values of . There are snake motions (i.e. with the symmetry and ) which have and , giving . The minimization search commonly encounters motions which are nearly as ineffective throughout shape space, and has a large gradient in the vicinity of such motions. It is challenging to minimize such a function numerically. We therefore substitute for . The minima are the same for both because each is an increasing function of the other, but varies much more smoothly in most of shape parameter space. There are cases where tends to , but these only occur in a small region of parameter space which is rarely approached by our computations.

Combining these ideas, the function we minimize in place of is

 F=−dWe2cos(θ(0,1)−θ(0,0)). (41)

## Appendix B Supplementary Material: Invariance under reparametrization of time

Here we show that , , and are the same for any reparametrization of time—in other words, for with any nondecreasing differentiable mapping from to with and . Such a reparametrization can be used to reduce the high-frequency components of a motion while keeping the efficiency the same. Therefore, it is reasonable to expect that a good approximation to any minimizer of can be obtained with low temporal frequencies, that is, with an expansion in the form of (19) with small . To show that , , and are the same for a reparametrization of time, let be given and be the aforementioned parametrization function. Let be the motion obtained by integrating from as in (1)–(3) with the integration constants and chosen to satisfy the dynamical equations (11)–(13). Now we show that also satisfies (11)–(13) at each time . By the chain rule of differentiation,

 ˆ∂t~X(s,t)=dX(s,γ(t))/dt∥dX(s,γ(t))/dt∥=γ′(t)∂γX(s,γ)|γ=γ(t)∥γ′(t)∂γX(s,γ)|γ=γ(t)∥=ˆ∂γX(s,γ)|γ=γ(t) (42)

using . Let be the force distribution corresponding to . It is given by (14) with