A Small-wavelength integrals

# Optimizing snake locomotion in the plane. II. Large transverse friction

## Abstract

We determine analytically the form of optimal snake locomotion when the coefficient of transverse friction is large, the typical regime for biological and robotic snakes. We find that the optimal snake motion is a retrograde traveling wave, with a wave amplitude that decays as the -1/4 power of the coefficient of transverse friction. This result agrees well with our numerical computations.

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

## I Introduction

In [1] we considered the problem of optimizing snake locomotion in the plane numerically. The mechanics of snake locomotion has been described previously by biologists and engineers, and modeled by applied mathematicians ([7]; [2]; [3]; [4]; [5]; [6]). As with other terrestrial locomoting animals ([8]), a terrestrial snake pushes against the ground with its body, and obtains a reaction force from the ground which propels it forward. However, a wide range of possible kinematics can be employed, and determining which are most effective (i.e. efficient) in different environments, and why, has been a major theme of locomotion studies ([9]; [10]; [11]; [12]; [13]; [14]; [15]; [16]; [17]; [18]; [19]; [20]).

A first approximation to snake locomotion is to consider motions confined to two dimensions [2]; [3]; [6]. Here the reaction force in the plane of snake motion is due to friction, and Coulomb friction provides a simple model. Hu and Shelley [3]; [6] emphasized the importance of frictional anisotropy in snake locomotion. In particular, they found that when the curvature of the snake backbone is prescribed as a sinusoidal traveling wave, high speed and efficiency is obtained when the coefficient of transverse friction is large compared to the coefficient of forward friction. Jing and Alben [21] found that the same holds for two- and three-link snakes. In biological snakes, the ratio of transverse to forward friction is thought to be greater than one [3]; [6], although how much greater is not clear. In [3]; [6], a ratio of only 1.7 was measured for anaesthetized snakes. These works noted that active snakes use their scales to increase frictional anisotropy, so the ratio in locomoting snakes is potentially much larger, and indeed, a ratio of 10 was found to give better agreement between the Coulomb friction model and a biological snake [6]. Wheeled snake robots have been employed very successfully for locomotion [22]; [23], and for most wheels the frictional coefficient ratio (with forward friction defined by the rolling resistance coefficient of the wheel) is much greater than 10 [24].

In this paper, we develop an analytical solution for the optimal locomotion of a snake in the plane, when the ratio of transverse to forward friction is large. We make few assumptions on the snake kinematics at the outset, but use numerical computations from [1] to provide guiding intuition. We find first that the optimal form of the prescribed backbone curvature is a traveling wave. We then find that the root-mean-square (RMS) amplitude of the curvature should decay as the power of the ratio of transverse to forward friction coefficients. Surprisingly, any periodic traveling wave motion can achieve the optimal efficiency, subject to the aforementioned restriction on its RMS amplitude, and to the requirement that its wavelength, normalized by the snake length, tend to zero. In the limit of large transverse friction, the power required to move a snake optimally is simply that needed to tow a straight snake forward.

## Ii Model

We use the same frictional snake model as ([3]; [6]; [21]; [1]), 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 [3]:

 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 ([25]).

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 ([3]). 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, ([3]), 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 ([3]; [6]; [21]), and was found to agree well with the motions of biological snakes in ([3]).

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 here is to find which minimizes as .

## Iii Large-μt analysis

We now analytically determine the optimal snake motion in the limit of large . Figure 2 shows a few numerical results from [1], which provide some intuition as we develop the analytical solution. In [1] we show that for , the numerically-computed optimal motions are always (retrograde) traveling waves, an example of which is shown in panel a for . A sequence of snapshots of the snake is shown over a period of motion. The snake moves from left to right, and appears to follow a sinusoidal path, very similar to what has been observed biologically and studied numerically [6]. We have computed similar traveling-wave optima at a range of . In panel b we compare the curvatures from these optimal motions at a particular instant when a curvature maximum occurs at the snake midpoint, for ten different values of : 5, 6, 7, 10, 20, 30, 60, 100, 200, and 300. The profiles are similar, but the amplitudes decay monotonically as increases. Panel c shows the ten snake shapes corresponding to the curvatures in panel b. Consistently, we see the decay of the traveling-wave amplitude as increases. With this picture from the numerics, we now proceed to derive the analyical solution in the asymptotic limit of large .

We begin with the position of the snake in terms of its components:

 X(s,t)=(x(s,t),y(s,t)). (18)

If is large, the snake moves more easily in the tangential direction than in the transverse direction. Furthermore, the most efficient curvature function will give motion mainly in the tangential direction, to avoid large work done against friction. We assume that the mean direction of motion is aligned with the -axis, and that the deflections from the -axis are small—that is, and higher derivatives are for some negative . This has been observed in the numerical solutions, and makes intuitive sense. Small deflections allow the snake to move along a nearly straight path, which is more efficient than a more curved path, which requires more tangential motion (and work done against tangential friction) for a given forward motion, i.e. .

We expand each of the terms in the force and torque balance equations in powers of and retain only the terms which are dominant at large . We expand the cost of locomotion similarly, and find the dynamics which minimize it, in terms of .

We first expand . We decompose into its -average and a zero--average remainder:

 x(s,t)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯x(s,t)+−∫scosθ(s′,t)ds′ (19)

where means the constant of integration is chosen so that the integrated function has zero -average. Hence

 ∂tx(s,t)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂tx(s,t)+−∫s−∂tθ(s′,t)sinθ(s′,t)ds′. (20)

We denote the -averaged horizontal velocity (the horizontal velocity of the snake center of mass) by

 U(t)≡¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂tx(s,t). (21)

We expand the integrand in (20) in powers of and its derivatives:

 −∂tθ(s,t)sinθ(s,t)=−∂sty(s,t)∂sy(s,t)+O(|∂sy|4). (22)

The quartic remainder term in (22), , actually includes other quartic terms involving time derivatives of , but for brevity we use the assumption (which will be correct for our solution) that and all its derivatives are of the same order of smallness with respect to . We define

 h(s,t)=−∫s−∂s′ty(s′,t)∂s′y(s′,t)ds′ (23)

and with the small-deflection assumption, (20) becomes

 ∂tx(s,t)=U(t)+h(s,t)+O(|∂sy|4), (24)

and

 ∂tX(s,t)=(U(t)+h(s,t),∂ty(s,t))+O(|∂sy|4). (25)

We will see that for a general class of shape dynamics near the optimum, is (1), i.e. of the same order as one snake length per actuation period, even in the limit of large . is therefore the dominant term in (24) and (25). We use this assumption on to expand the velocity 2-norm:

 ∥∂tX(s,t)∥ =√∂tx2+∂ty2=U(1+hU+12∂ty2U2)+O(|y|4). (26)

In writing we have again lumped terms with four or more powers of and/or its - and -derivatives. Dividing (25) by (26), we obtain the expansion for the normalized velocity:

 ˆ∂tX(s,t)=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝1−12∂ty2U2+O(|y|4)∂tyU(1−hU−12∂ty2U2)+O(|y|5)⎞⎟ ⎟ ⎟ ⎟ ⎟⎠. (27)

We next expand the unit tangent and normal vectors:

 ^s=(∂sx∂sy)=(√1−∂sy2∂sy) =(1−12∂sy2∂sy)+O(|y|4). (28) ^n =(−∂sy1−12∂sy2)+O(|y|4). (29)

Using (27) – (29),

 ˆ∂tX⋅^s =1−12(∂sy−∂tyU)2+O(|y|4), (30) ˆ∂tX⋅^n =(∂tyU−∂sy)(1−12∂ty2U2)+∂tyU(−hU−12∂sy2)+O(|y|5). (31)

Using (30) and (31) we can write the force and torque balance equations (11)–(13) and the cost of locomotion in terms of , , and derivatives of .

### iii.1 Traveling-wave optimum

We first expand and use it to argue that the optimal shape dynamics is approximately a traveling wave. This allows us to neglect certain terms, which simplifies the formulae for the rest of the analysis. If the power expended at time is , the cost of locomotion is

 η =∫10P(t)dt/∫10U(t)dt. (32) =∫10U(t)∫10[(ˆ∂tX⋅^s)2+μt(ˆ∂tX⋅^n)2]∥∂tX∥U(t)dsdt/∫10U(t)dt. (33)

In (33) we have assumed that the tangential motion of the snake is entirely forward, with no portion moving backward, so the coefficient in front of the tangential term is unity (i.e. , rather than a term involving both and ). This assumption is not required but it simplifies the equations from this point onward and will be seen to be correct shortly. The numerator of (33) can be separated into a part involving tangential velocity and a part involving normal velocity with a factor of . Using (26), (30), and (31), we can list the terms with the lowest powers of in each part:

 (ˆ∂tX⋅^s)2∥∂tX∥U(t) =1+O(|y|2) (34) μt(ˆ∂tX⋅^n)2∥∂tX∥U(t) =μt(∂sy−∂tyU)2+O(μt|y|4). (35)

The optimal shape dynamics is a that minimizes subject to the constraints of force and torque balance. The leading-order term in (34) is 1, which is independent of the shape dynamics. At the next (quadratic) order in , terms appear in both (34) and (35); only that in (35) is given explicitly. It is multiplied by , so it appears to be dominant in the large limit. Thus, to a first approximation, minimizing is achieved by setting

 (∂sy−∂tyU)=0, (36)

so a first guess for the solution is a traveling wave, i.e. a function of with constant. However, such a function does not satisfy the -component of the force balance equation (11), which is:

 ∫[(ˆ∂tX⋅^s)sx+μt(ˆ∂tX⋅^n)nx]ds=0. (37)

Inserting (30) and (31) and retaining only the lowest powers in from each expression,

 ∫101+μt(∂sy−∂tyU)∂syds=0. (38)

which cannot be satistied by a function of . Physically, the first term (unity) in (38) represents the leading-order drag force on the body due to forward friction. For a function of , in which the deflection wave speed equals the forward speed, the snake body moves purely tangentially (up to this order of expansion in ) along a fixed path on the ground, with no transverse forces to balance the drag from tangential frictional forces. Also, when a function of is inserted into (38) all factors of cancel out and is left undetermined. This is related to a more fundamental reason why posing the shape dynamics as a function of is invalid: the problem of solving for the snake dynamics given in section II is not then well-posed. In the well-posed problem we provide the snake curvature and curvature velocity versus time as inputs, and solve the force and torque balance equations to obtain the average translational and rotational velocities as outputs. Thus is one of the three outputs (or unknowns) that is needed to solve the three equations.

Instead we pose the shape dynamics as

 y(s,t)=g(s+Uwt), (39)

which is a traveling wave with a prescribed wave speed , different from in general. Since is periodic in with period 1, is periodic with a period of . Inserting (39) into (38), we obtain an equation for in terms of and :

 1+μt(1−UwU)∫10g′(s+Uwt)2ds=0. (40)

In (40), as , . In the solution given by (39) and (40), the shape wave moves backwards along the snake at speed , which propels the snake forward at a speed , with slightly less than . Therefore the snake slips transversely to itself, in the backward direction, which provides a forward thrust to balance the backward drag due to forward tangential friction. We refer to in (40), a measure of the amount of backward slipping of the snake, as the “slip.” In figure 2a, the snake slips transversely to itself, and backward (leftward), even as it moves rightward. Here the slip is fairly small because is 30, fairly large.

We now solve (40) for and insert the result into (33) to find and which minimize . Using the lowest order expansions (34) and (35) with (40) we obtain:

 η =∫10P(t)dt/∫10U(t)dt. (41) =∫10U(t)∫101+μt(∂sy−∂tyU)2+O(|y|2,μt|y|4)dsdt/∫10U(t)dt. (42) =1/∫1011+1μt⟨g′(s+Uwt)2⟩+O(|g|2,μt|g|4)dt. (43)

where

 ⟨g′(s+Uwt)2⟩≡∫10g′(s+Uwt)2ds. (44)

If we minimize (43) over possible for a given , using only the terms in the expansion that are given explicitly, we find that tends to a minimum of 1 as the amplitude of diverges. However, in this limit our small- power series expansions are no longer valid. We therefore add the next-order terms in our expansions of and the equation and search again for an -minimizer.

### iii.2 Optimal wave amplitude

We now expand to higher order the components of , (34) and (35),

 (ˆ∂tX⋅^s)2∥∂tX∥U(t) =1−(∂sy−∂tyU)2+hU+12∂ty2U2+O(|y|4) (45) μt(ˆ∂tX⋅^n)2∥∂tX∥U(t) =μt[(∂tyU−∂sy)+∂tyU(−hU−12∂sy2)]2(1+O(|y|2)), (46)

and equation (11),

 ∫1+μt(∂sy−∂tyU+∂tyU(hU+12∂sy2))∂syds=0. (47)

We again assume the traveling wave form of (39), and use this to evaluate given by (23) in terms of :

 h(s,t)=−Uw2g′(s+Uwt)2+Uw2⟨g′(s+Uwt)2⟩. (48)

Inserting the traveling-wave forms of (39) and (48) into (47) we obtain

 1+μt[(1−UwU)⟨g′(s+Uwt)2⟩+12⟨g′(s+Uwt)2⟩2]=0. (49)

To obtain (49) we’ve used the fact that the slip as , which will be verified subsequently. One can also proceed without this assumption, at the expense of a lengthier version of (49). We insert the traveling-wave forms of and into (45) and (46) to obtain , again with the assumption that as to simplify the result:

 η=1+∫1012⟨g′(s+Uwt)2⟩+μt(1−UwU+12⟨g′(s+Uwt)2⟩)2⟨g′(s+Uwt)2⟩dt (50)

Solving for the slip from (49) and inserting into (50), we obtain an improved version of (43):

 η=1/∫1011+12⟨g′(s+Uwt)2⟩+1μt⟨g′(s+Uwt)2⟩+O(|g|4,μt|g|8)dt. (51)

(51) is a more accurate version of (43), with an additional term in the denominator of the integrand which penalizes large amplitudes for . If we approximate as constant in time, we obtain

 η=1+12⟨g′2⟩+1μt⟨g′2⟩+O(|g|4,μt|g|8). (52)

which is minimized for

 ⟨g′2⟩1/2=21/4μ−1/4t. (53)

Thus the RMS amplitude of the optimal should scale as . (53) represents a balance between competing effects. At smaller amplitudes, the forward component of normal friction is not large enough to balance the forward component of tangential friction, so the snake slips more, which reduces its forward motion and does more work in the normal direction. At larger amplitudes, as already noted, the snake’s path is more curved, which requires more work done against tangential friction for a given forward distance traveled. (53) is only a criterion for the amplitude of the motion. To obtain information about the shape of the optimal , we now consider the balances of the -component of the force and the torque.

### iii.3 Small-wavelength shapes

So far we have searched for the optimal shape dynamics in terms of , but in fact it is the curvature which is prescribed. We obtain from the curvature by integrating twice in :

 y(s,t) =y(0,t)+∫s0sinθ(s′,t)ds′ (54) =y(0,t)+∫s0θ(s′,t)ds′+O(y3) (55) =y(0,t)+∫s0[θ(0,t)+∫s′0κ(s′′,t)ds′′]ds′+O(y3). (56) =y(0,t)+θ(0,t)s+∫s0∫s′0κ(s′′,t)ds′′ds′+O(y3). (57) ≡Y(t)+sR(t)+k(s,t)+O(y3). (58)

where and are defined for notational convenience. Prescribing the curvature is equivalent to prescribing . We set

 k(s,t)=g(s+Uwt) (59)

so we have the same form for as before, with an additional translation and rotation . and are determined by the -force and torque balance equations:

 ∫10[(ˆ∂tX⋅^s)sy+μt(ˆ∂tX⋅^n)ny]ds=0. (60) ∫10[(ˆ∂tX⋅^s)(xsy−ysx)+μt(ˆ∂tX⋅^n)(xny−ynx)]ds=0. (61)

We expand these to leading order in and obtain

 ∫10∂sy−∂tyU+∂tyU(hU+12∂sy2)ds=0. (62) ∫10s(∂sy−∂tyU+∂tyU(hU+12∂sy2))ds=0. (63)

We insert (58) with (59) into the three equations (47), (62), and (63) to solve for the three unknowns, , , and in terms of and . We obtain:

 ∫101+μt[(−UwU+1+12⟨g′2⟩)g′2+(−Y′U−R′sU+R)g′]ds=0. (64) ∫10(−UwU+1+12⟨g′2⟩)g′−Y′U−R′sU+Rds=0. (65) ∫10s[(−UwU+1+12⟨g′2⟩)g′−Y′U−R′sU+R]ds=0. (66)

We solve (65) and (66) for and in terms of :

 Y′U−R =4B−6A (67) R′U =12A−6B (68) A ≡(−UwU+1+12⟨g′2⟩)⟨sg′⟩ (69) B ≡(−UwU+1+12⟨g′2⟩)⟨g′⟩ (70)

We then solve (64) for in terms of :

 Missing or unrecognized delimiter for \left (71)

We now form , (33) with terms given by (45) and (46), using the updated from of ((58) with (59)) and (71). We obtain a more correct version of (51):

 η=1/∫1011+12⟨g′2⟩+1μt(⟨g′2⟩−⟨g′⟩2−3(⟨g′⟩−2⟨sg′⟩)2)+O(|g|4,μt|g|8)dt. (72)

We can find -minimizing in a few steps. First, let be the family of orthonormal polynomials with unit weight on (essentially the Legendre polynomials):

 ∫10LiLjds=δij;L0≡1,L1=√12(s−1/2), …. (73)

At a fixed time we expand in the basis of the :

 g′(s+Uwt)=∞∑k=0ck(t)Lk(s). (74)

Then we have:

 ⟨g′2⟩ =∞∑k=0ck(t)2, (75) ⟨g′⟩2 =c0(t)2, (76) 3(⟨g′⟩−2⟨sg′⟩)2 =⟨L1g′⟩2=c1(t)2. (77)

Inserting into (72), becomes

 η=1/∫1011+12(c0(t)2+c1(t)2+∑∞k=2ck(t)2)+1μt(∑∞k=2ck(t)2)+O(|g|4,μt|g|8)dt. (78)

The denominator of the integrand in (78) is minimized when

 c0(t)=0;c1(t)=0;∞∑k=2ck(t)2=√2μ−1/2t. (79)

The function that minimizes is that for which (79) holds for all . Then the integral in (78) is maximized, so is minimized. Recall that is a periodic function with period . The relations (79) hold for any such in the limit that , as long as is normalized appropriately. Define

 A≡(1Uw∫Uw0g′(x)2dx)1/2. (80)

Then

 ∫10g′(s+Uwt)2ds =A2+O(Uw) (81) ∫10g′(s+Uwt)ds =O(Uw) (82) ∫10sg′(s+Uwt)ds =O(Uw). (83)

(81) – (83) are fairly straightforward to show, and we show them in Appendix A for completeness. If (81) – (83) hold, then in the limit that , (79) holds with , by (75) – (77). A physical interpretation of the small-wavelength limit is that the net vertical force and torque on the snake due to the traveling wave alone () become zero in this limit, so no additional heaving motion () or rotation () are needed, and thus the additional work associated with these motions is avoided. Similar “recoil” conditions were proposed as kinematic constraints in the context of fish swimming [15]. We also note two other scaling laws. For a given , by (71) the slip and by (67)–(70), .

### iii.4 Optimal cost of locomotion and comparison with numerics

To obtain an -minimizing function , let be any periodic function with period . Define

 g∗(s+Uwt)≡21/4μ−1/4tAg(s+Uwt) (84)

where is given in (80). Then in the limit ,

 η(g∗)→1+√2μ−1/2t+O(μ−1t), (85)

by the estimate in (72). Thus the optimal traveling wave motions become more efficient as , which agrees with the numerical solutions in [1]. The numerically-determined optima in figure 2 do not have small wavelengths, due to the finite number of modes used. This is explained further in the Supplementary Material of [1].

We compare the numerically-found optima to the theoretical results in figure 3. Here we take the plots of figure 2 and transform them according to the theory. In figure 3a we plot versus for the numerical solutions (squares) from [1], along with the small- theoretical result (solid line), which is given by (85). 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 in the numerical computation, which makes the numerical result underperform the analytical result, as explained in [1]. The circle shows the efficiency for a curvature function with a higher wavelength than those which could be represented in the numerical optimization, and its distance from the solid line is , of the order of the next term in the asymptotic expansion of .

Figure 3b shows the curvature, rescaled by to give a collapse according to (53). The data collapse well compared to the unscaled data in figure 2b. Figure 3c shows the body shapes from figure 2c 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.

## Iv Conclusion

We have studied the optimization of planar snake motions for efficiency, using a fairly simple model for the motion of snakes using friction. When the coefficient of transverse friction is large, our analysis shows that a traveling-wave motion of small amplitude is optimal. The amplitude tends to zero as the transverse friction coefficient tends to infinity, scaling as the transverse friction coefficient to the -1/4 power, and the efficiency tends to unity in this limit. The corresponding power is that for a straight snake towed forward.

In [1] we were able to compute many optimal motions at moderate and small also. Those found at small also corresponded to traveling waves (direct now, not retrograde), and at zero , a simple triangular-wave solution was found with optimal efficiency. It remains to determine optimal motions with small but nonzero analytically. At moderate , a region of standing-wave or ratcheting motions was found. It may be possible to compute some of these motions analytically using a set of assumptions specialized to the moderate- regime.

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

## Appendix A Small-wavelength integrals

To show (81) – (83), we first note that since is periodic, is periodic with zero average over a period. To show (81), we decompose the integral in (81) into two intervals:

 ∫10g′(s+Uwt)2ds=∫Uw⌊1/Uw⌋0g′(s+Uwt)2ds+∫1Uw⌊1/Uw⌋g′(s+Uwt)2ds. (86)

The first integral on the right of (86) is over an integral number of periods of , and is thus