Slinky Take Two

mass = 212.0/1000;
pluslength = 55.4/1000;
density = mass/pluslength;
gravity=9.8;
coils = 84;
width = pluslength/coils;
length=pluslength-width;

mcdata={0,2+3/4,5+1/4,7+5/8,10,12+3/8,14+5/8,17,19+1/8, 21+1/8,23+3/16,25+1/8, 27,28+7/8,30+5/8,32+1/4,34,35.5,36+15/16,38+3/8,39+3/4,41, 42+1/8,43+1/4,44+3/8,45+3/8,46+1/4,47+1/8,48,48+3/4,49.5,50,50.5,51,51+3/8,51+3/4,52,52+1/4,52+3/8, 52.5,52+9/16,52+5/8}

Length[mcdata]

mkdata=-2.54 mcdata/100

step[a_,b_,m_,j_]=a+j (b-a)/m;

plotdata=Table[{-mkdata[[j+1]],step[0,length,41,j]},{j,0,41}];

ListPlot[plotdata]

constant elasticity model

Lagrangian :  eps (w’-1)^2/2 - rho g w  
with eps = elasiticy assumed constant

Euler-Lagrange ODE:  eps w’’ = - rho g

general solution with w’(L) = 1:   

                                             w(x) = - rho g x^2/(2 eps) + (rho g L/eps + 1) x + b.

fitdata=Table[{step[0,length,41,j],-mkdata[[j+1]]},{j,0,41}];

conone[x_]=Fit[fitdata,{1,x,x^2},x]

Show[ListPlot[plotdata],ParametricPlot[{conone[x],x},{x,0,length}]]

This suggests we can fit the data quite well with a quadratic function.  The coefficients of this function, however, 
do not fit the form of our model:

elasticconstantquadratic=NSolve[-density gravity/(2 findeps)==Coefficient[conone[x],x^2]][[1]][[1,2]]

elasticconstantlinear=NSolve[density gravity/findeps+1==Coefficient[conone[x],x]][[1]][[1,2]]

To see this another way, consider choosing the coefficient b = 0 in the model so that w(0) = 0.  Then there is precisely one value of the elasticity constant eps for which w(length) = the measured value = 52+5/8 inches = 1.337 m.  That is,

elasticconstantends=NSolve[density gravity length^2/( 2 findeps)+length==-mkdata[[42]]][[1]][[1,2]]

Notice that this value is between the other two we got from the best quadratic fit.

quadraticends[x_]=-density gravity x^2/(2 elasticconstantends)+(density gravity length/elasticconstantends+1) x;

Show[ListPlot[plotdata],ParametricPlot[{quadraticends[x],x},{x,0,length}]]

And this fit is not so great.  

Another possibility is to assume again that we have w(0) = 0, so that b = 0, and then write 

w(x) = (1/eps) (-rho g x^2/2 + rho g length x) + x .

We can then use our data set to consider w(x) as a linear function of (-rho g x^2/2 + rho g x) and solve for the 
best eps:

secondbestfitdata=Table[{-mkdata[[j+1]]-step[0,length,41,j],density gravity (-step[0,length,41,j]^2/2+length step[0,length,41,j])},{j,0,41}];

ListPlot[secondbestfitdata]

second[x_]=Fit[secondbestfitdata,{x},x]

secondeps=Coefficient[second[x],x]

secondquadratic[x_]=-density gravity x^2/(2 secondeps)+(density gravity length/secondeps+1) x;

Show[ListPlot[plotdata],ParametricPlot[{secondquadratic[x],x},{x,0,length}],PlotRange->All]

This doesn’t fit so great either especially near x = length where the extension is small, but the elasticity constant is rather close to the last one; around 0.043.  

Finally, we might use elasticity = 0.043 and attempt to fit a quadratic with nonzero constant term b:

finalquadratic[x_]=Fit[fitdata,{-density gravity x^2/(2 (0.043))+(density gravity length /0.043+1) x,1},x]

Show[ListPlot[plotdata],ParametricPlot[{finalquadratic[x],x},{x,0,length}],PlotRange->All]

That is perhaps about the best we can do with the constant elasticity model.  Some modelers have, at this point, looked at the nice quadratic fit and simply adopted that as a “model” without any really good explanation, except that there is a good fit of the data.  We will proceed somewhat differently.

Preliminary discussion of improvement

The basic assertion is that the elasticity should depend on the extension rate w’.  Here are some calculations intended to give some idea about how to find that dependence:
Given two data values (x_j,w_j) and (x_k,w_k) give 

                                            w_j = - rho g x_j^2/(2 eps) + (rho g L/eps + 1) x_j + b
                                            w_k = - rho g x_k^2/(2 eps) + (rho g L/eps + 1) x_k + b
                                            
          and
                                   w_k - w_j = - rho g [ x_k^2 - x_j^2 ]/(2 eps) + (rho g L/eps + 1)(x_k - x_j)
                                
           or
           
                                  w_k - w_j - (x_k-x_j) = (rho g/eps) [ - (x_k^2 - x_j^2/2 + L(x_k - x_j) ].
                                  
           which gives 
           
                                 eps = rho g (x_k - x_j)[ L - (x_k + x_j)/2] / [ w_k - w_j - (x_k-x_j) ]
                               
 Thus, we define

elas[xone_,xtwo_,wone_,wtwo_]=density gravity (xtwo-xone)( length - (xtwo+xone)/2 )/(wtwo-wone-(xtwo-xone));

elasticitydata=Table[{(-mkdata[[j+2]]+mkdata[[j+1]])/(step[0,length,41,j+1]-step[0,length,41,j]),elas[step[0,length,41,j],step[0,length,41,j+1],-mkdata[[j+1]],-mkdata[[j+2]]]},{j,0,40}]

The first coordinate here is approximating the extension rate w’ and the second is using the two points to determine a value of the elasticity constant.  The problem with using only two points (at a time) is that there is some instability in the resulting data, especially for small extensions when w’ is close to 1 at the end.  This is somewhat expected because, first of all, the coils are all bunched up there and it is difficult to get accurate height measurements.  In fact, it is probably the case that the accuracy of the measurement of the height of any given coil is significantly outside the width of the equilibrium coils, so one can expect this data to be rather bad (from point to point for the last few points).  In addition, the coil tends to “bounce” a bit, so reading the height is difficult for that reason too.  
At any rate, this is what you see in the data.  For relatively large extension, the elasticity is around our best constant value 0.043.  It gets fuzzy for small extension, i.e., w’ close to 1, but is uniformly higher.

Show[ListPlot[elasticitydata,AxesOrigin->{0,0},PlotRange->All],Plot[0.043,{x,0,55}]]

elaone[x_]=Fit[elasticitydata,{1,x},x]

Show[ListPlot[elasticitydata,AxesOrigin->{0,0},PlotRange->All],Plot[elaone[x],{x,0,55}]]

eladata=Table[{(-mkdata[[j+2]]+mkdata[[j+1]])/(step[0,length,41,j+1]-step[0,length,41,j]),elas[step[0,length,41,j],step[0,length,41,j+1],-mkdata[[j+1]],-mkdata[[j+2]]]},{j,0,35}];

elatwo[x_]=Fit[eladata,{1,x-1},x]

Show[ListPlot[elasticitydata,AxesOrigin->{0,0},PlotRange->All],Plot[elaone[x],{x,0,55}],Plot[elatwo[x],{x,0,55},PlotStyle->Red],ListPlot[eladata,PlotStyle->Green]]

Let us remember the coefficients of elatwo[w’] = alpha(w’-1) + beta

alphatwo = -0.00038192900950803854`;
betatwo=0.05618417202447363` ;

Variable elasticity model

As suggested above, we will consider a variable elasticity of the form eps = a (w’-1) + b  where presumably the coefficient a is small and negative and the coefficient b is somewhere around 0.043.  

Here is a derivation:

Lagrangian :      eps (w’-1)^2/2 - rho g w

Euler-Lagrange ODE:  -(d/dx)[ eps’ (w’-1)^2/2 + eps (w’-1) ] = rho g

Note that here eps’ = (d eps / dp) = a is constant where p denotes the p = w’ dependence.  

Thus, we get 

[ a(w’-1) + a (w’-1) + eps ] w’’ = - rho g.

That is,

[ 3 a (w’-1) + b ] w’’ = - rho g.

Integrating once from x to L (and remembering the condition w’(L) = 1) we get 

-[ 3 a w’^2/2 + (b - 3a) w’ ] + b-3a/2 = -rho g (L-x).

This can be written as 

3 a w’^2 + 2(b-3a) w’ - 2 rho g(L-x) + 3a-2b = 0.

Solving the quadratic, 

w’ = (1/(6a)) [ 2(3a-b) + 2 sqrt( (b-3a)^2 - 3a(-2 rho g(L-x) + 3a-2b) ) ]  
      = (1/(3a)) [  3a-b + sqrt( b^2 +6a rho g(L-x) ) ]
     
 We take the “+“ sign in the quadratic formula because 3a-b is assumed to be negative, b being rather larger than a and positive, and we need an opposite sign to get the quadratic model (constant elasticity) in the limit as a tends to 0, i.e., eps = a w’ + b tends to b.
 
 Integrating this time from 0 to x and assuming w(0) = 0 we get 
 
 w(x) = (1-b/(3a)) x - (1/(3a)) (2/3)  [ b^2 +6a rho g(L-x) ]^(3/2) (1/( 6 a rho g )) 
                                     +  (1/(3a)) (2/3)  [ b^2 +6a rho g L ]^(3/2) (1/( 6 a rho g )). 
                                     
 Substituting our choice elatwo(w’) for eps here we get

wtwo[x_] = (1 - betatwo/(3 alphatwo)) x-(1/(27 alphatwo^2 density gravity))  ((betatwo^2 + 6 alphatwo density gravity (length-x))^(3/2)-(betatwo^2 + 6 alphatwo density gravity length)^(3/2))

Show[ListPlot[plotdata],ParametricPlot[{wtwo[x],x},{x,0.02,length}],PlotRange->All]

alphatwo

Plot[betatwo^2 + 6 alphatwo density gravity (length-x),{x,0,length}]