Re: Newton e Hooke

In article <vnttl2$20mb1$1@dont-email.me>, Luigi Fortunati wrote

If I push the end A of a spring, do I compress it or accelerate it?
 
Obviously I compress it and accelerate it at the same time, because I
am not pushing the elastic body of the spring but only its point A.
 
I ask myself: how much of my force is dedicated to compression and how
much to acceleration?

In general the spring's center of mass accelerates AND the spring is
compressed.

But, it's not correct to say that only part of the applied force is
dedicated to acceleration -- actually, the spring's center-of-mass
acceleration is determined by *all* of the applied force, while at the
same time the spring compresses.

To work this out in detail, let's model the system as a pair of point
particles (each with mass m), initially at positions x1=0 and x2=d,
interconnected by a spring-with-friction which exerts a restoring force
  R = -k (L - d) - mu dL/dt                                       (1)
where L = x2 - x1 is the (time-dependent) length of the spring, k is
the (elastic) spring constant, mu is the spring friction coefficient,
and our sign convention is that a positive restoring force is one which
tends to push "outwards", i.e., one that tends to increase L.

Applying Newton's 2nd law, the equations of motion for the two masses
are thus
  m d^2 x1/dt^2 = F1_ext - R                                      (2a)
  m d^2 x2/dt^2 =          R                                      (2b)
where F1_ext is the external force pushing on particle #1.

Now let's work out the motion of the spring's center of mass.  The
center of mass is located at
  xc = (x1+x2)/2  .                                               (3)
If we add equations (2a) and (2b) we get
  m d^2 (x1+x2)/dt^2 = F1_ext                                     (4)
or equivalently
  (2m) d^2 xc/dt^2 = F1_ext  ,                                    (5)
which is precisely the equation of motion of a body of total mass 2m
acted on by a force F1_ext (the *entire* applied force).  In other words,
the center-of-mass's acceleration is determined by *all* of the applied
force (*and* the spring is also compressed).

(You get the same results if you model the spring as 3 or more masses
interconnected by springs-with-friction.  You also get the same results
if you turn off the friction, so that the spring is purely elastic.)

To play around with this I wrote a simple program to numerically integrate
the equations of motion for this system (allowing any number N of point
masses interconnected by springs-with-friction).  The program is short
enough that I'll give it below.

For the example above, with the program's default parameters m=1, d=1,
k=2.5, mu=0.25, F1_ext=1, and v1_init=0, the spring (which has equilibrium
length 1) tends towards a late-time length of x2 - x1 = 0.8, i.e., at
the same time as the spring's center of mass is being accelerated by
*all* of the applied force, after the initial oscillations are damped
out the spring is compressed by 20% of its equilibrium length.

Here is the simulation program.  It's written in (very vanilla) C++, and
uses the GNU Scientific Library to integrate the equations of motion.

--- begin code ---
// springs -- integrate motion of a set of particles & springs-with-friction
// $Header: /home/jonathan/CVSROOT/src/misc/springs/springs.cc,v 1.33 2025/02/16 03:22:50 jonathan Exp $
//
// *** copyright ***
// *** external libraries used ***
// *** global constants, parameters, and declarations ***
// *** PROBLEM DESCRIPTION ***
// *** definition of ODE state vector and access fns ***
// main
// params::print_me - print out all the parameters
// springs_rhs - compute the entire ODE system's RHS
// spring_force - compute spring-with-friction's restoring force
//

#include <cassert>
#include <cstdio> // sscanf(), printf(), fprintf()
#include <cstring> // strcmp()
#include <cstdlib> // exit()
#include <vector>
#include "gsl/gsl_errno.h"
#include "gsl/gsl_odeiv2.h"

using std::sscanf;
using std::printf;
using std::fprintf;
using std::strcmp;

//******************************************************************************

//
// Copyright 2025, Jonathan Thornburg
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//

//******************************************************************************

//
// *** external libraries used ***
//

//
// This program uses the GNU Scientific Library (GSL) for the ODE integration.
// This program compiles and runs successfully using Clang 16.0.6 and GSL
// version 2.7.1, but would probably work using any reasonably-modern C++
// compiler and GSL version.
//

//******************************************************************************

//
// *** global constants, parameters, and declarations ***
//

// printf() format for printing numbers
const char* const number_format = "%f";

// "fuzz" for fuzzy floating-point arithmetic
const double fp_fuzz = 1.0e-10;

// all the system parameters
struct params
{
int N; // number of particles
double m; // mass of each individual particle
double d; // equilibrium distance between particles
double k; // (elastic) spring constant
double mu; // spring friction coefficient
double F1_ext; // external force applied to leftmost particle
double v1_init; // initial velocity of leftmost particle
bool right_wall_flag; // is there a fixed "wall" to the
// right of the rightmost particle?
double dt; // time interval at which to print output
double t_max; // time at which to end the simulation

// print all the parameters to /stdout/
// in a (vaguely) human-readable format
// ... each line begins with '##' (& hence is a gnuplot comment)
void print_me() const;

// construct with default parameters
explicit
  params()
: N      (2  ),
  m      (1.0),
  d      (1.0),
  k      (2.5),
  mu     (0.25),
  F1_ext (1.0),
  v1_init(0.0),
  right_wall_flag(false),
  dt     ( 0.1),
  t_max  (20.0) // no comma
{ /* empty constructor body */ }
};

// declarations of functions local to this file
extern "C"
  int springs_rhs(const double /* t */,
  const double u[],
  double udot[],
  void* void_params_ptr);
double spring_force(const params& p,
    const double u[], const int i);

//******************************************************************************

//
// *** PROBLEM DESCRIPTION ***
//

const char* const help_msg =
"Usage:\n"
"  springs [ N=%d ]\n"
"          [ m=%g ]   [ d=%g ]\n"
"          [ k=%g ]   [ mu=%g ]\n"
"          [ F1_ext=%g ]   [ v_init=%g ]\n"
"          [ right_wall_flag={false|true} ]\n"
"          [ dt=%g ] [ t_max=%g ]\n"
"or:\n"
"   springs --print-default-params\n"
"or:\n"
"   springs { --help | --help-long | --help-examples }\n"
"\n"
"This program integrates the Newtonian equations of motion for a system\n"
"of $N$ point particles connected by springs-with-friction.  The command-line\n"
"arguments specify the system parameters and the output times, and may be\n"
"specified in any order.\n"
"\n"
"The output begins with a descriptive header, then gives the position of\n"
"each particle at each specified output time.  The output format is directly\n"
"readable by the /gnuplot/ graphics package.\n"
;

const char* const help_long_msg =
"*** Physical System ***\n"
"\n"
"In more detail, the physical system consists of $N$ point particles (each\n"
"with mass $m$); we label the particles by an integer index $i=1$, $i=2$,\n"
"$i=3$, ..., $i=N$.  Initially (i.e., at $t=0$) each particles is at the\n"
"position $x_i=(i-1)d$, the leftmost ($i=1$) particle has a velocity\n"
"$v_{1,init}$, and all the other particles are at rest.\n"
"\n"
"Each adjacent pair of particles is connected by a spring-with-friction,\n"
"which exerts a restoring force\n"
"  $R_{i,i+1} = -k (L_{i,i+1} - d) - \mu \dot{L}_{i,i+1}$,\n"
"where $L_{i,i+1} = x_{i+1} - x_i$ is the length of the spring connecting\n"
"the adjacent particles $i$ and $i+1$, i.e., the distance between particles\n"
"$i$ and $i+1$, $k \ge 0$ is the (elastic) spring constant, $\mu \ge 0$\n"
"is the spring friction coefficient, and our sign convention is that a\n"
"positive restoring force is one which tends to push \"outwards\", i.e.,\n"
"tends to increase $L_{i,i+1}$.\n"
"\n"
"Optionally, there may be a fixed \"wall\" at the position $x=Nd$; in this\n"
"case the rightmost particle is connected to the wall by a spring-with-friction\n"
"just like all the other springs.\n"
"\n"
"For $t \ge 0$ there is an external force $F{1,ext}$ applied to the\n"
"leftmost ($i=1$) particle.\n"
"\n"
"For example, for N=5 the system looks like this:\n"
"\n"
"                                                             ||\n"
"      * vvvvvvvv * vvvvvvvv * vvvvvvvv * vvvvvvvv * vvvvvvvv ||\n"
"                                                             ||\n"
"     i=1        i=2        i=3        i=4        i=5\n"
"\n"
"With the sign convention that positive forces push to the right, the\n"
"equations of motion are thus\n"
"  m\ddot{x}_1 = F1_ext - R_{1,2}\n"
"  m\ddot{x}_i = R_{i-1,i} - R_{i,i+1} for $1 < i \le N$\n"
"where we define $R_{N,N+1}$ as the spring-with-friction restoring force\n"
"between the rightmost particle and the wall if the wall is present, or 0\n"
"if there's no wall.\n"
;

const char* const help_examples_msg =
"Example:\n"
"  springs\n"
"This simulates a system with default parameters.\n"
"\n"
"Example:\n"
"  springs N=5 mu=2.5 right_wall_flag=true t_max=50\n"
"This simulates a system of N=5 point masses with a larger-than-default\n"
"spring friction coefficient mu=2.5 and a \"wall\" to the right of the\n"
"rightmost mass.  The simulation will run until t=t_max=50.\n"
;

//******************************************************************************

//
// *** definition of ODE state vector and access fns ***
//

//
// For the ODE system (which must be 1st order in time), we define the
// state vector (which is of size $2N$) as follows:
//   $u = (x_1, x_2, x_3, ..., x_N,
//         \dot{x}_1, \dot{x}_2, \dot{x}_3, ..., \dot{x}_N)$
//

// number of 1st-order-in-time equations
int N_eqns(const params& p)
{
return 2*p.N;
}

// 0-origin integer index in state vector corresponding to $x_i$
int posn_of_xi(const params& p, const int i)
{
assert(i >= 1); // sanity check
assert(i <= p.N); // ditto
return i-1;
}

// 0-origin integer index in state vector corresponding to $\dot{x}_i$
int posn_of_xdoti(const params& p, const int i)
{
assert(i >= 1); // sanity check
assert(i <= p.N); // ditto
return p.N + (i-1);
}

//******************************************************************************

int main(const int argc, const char* const argv[])
{
//
// parse the command-line arguments
//
params p;
if ((argc == 2) && (strcmp(argv[1], "--help") == 0))
{
printf("%s", help_msg);
exit(0); // *** --help exit
}
if ((argc == 2) && (strcmp(argv[1], "--help-examples") == 0))
{
printf("%s", help_examples_msg);
exit(0); // *** --help exit
}
if ((argc == 2) && (strcmp(argv[1], "--help-long") == 0))
{
printf("%s", help_msg);
printf("\n\n");
printf("%s", help_long_msg);
exit(0); // *** --help exit
}
if ((argc == 2) && (strcmp(argv[1], "--print-default-params") == 0))
{
p.print_me();
exit(0); // *** --help exit
}

for (int ap = 1 ; ap < argc ; ++ap)
{
if (sscanf(argv[ap], "N=%d", &p.N) == 1)
{ } // argument parsed ok ==> no-op here
else if (sscanf(argv[ap], "m=%lf", &p.m) == 1)
{ } // argument parsed ok ==> no-op here
else if (sscanf(argv[ap], "d=%lf", &p.d) == 1)
{ } // argument parsed ok ==> no-op here
else if (sscanf(argv[ap], "k=%lf", &p.k) == 1)
{ } // argument parsed ok ==> no-op here
else if (sscanf(argv[ap], "mu=%lf", &p.mu) == 1)
{ } // argument parsed ok ==> no-op here
else if (sscanf(argv[ap], "F1_ext=%lf", &p.F1_ext) == 1)
{ } // argument parsed ok ==> no-op here
else if (sscanf(argv[ap], "v1_init=%lf", &p.v1_init) == 1)
{ } // argument parsed ok ==> no-op here
else if (strcmp(argv[ap], "right_wall_flag=false") == 0)
{ p.right_wall_flag = false; }
else if (strcmp(argv[ap], "right_wall_flag=true") == 0)
{ p.right_wall_flag = true; }
else if (sscanf(argv[ap], "dt=%lf", &p.dt) == 1)
{ } // argument parsed ok ==> no-op here
else if (sscanf(argv[ap], "t_max=%lf", &p.t_max) == 1)
{ } // argument parsed ok ==> no-op here
else {
fprintf(stderr,
"***** springs: can't parse argument \"%s\"!\n",
argv[ap]);
fprintf(stderr, "%s", help_msg);
exit(1); // *** error exit
}
}

//
// setup state vector & initial conditions
//
// ... note state vector is a C++ /std::vector/ which GSL doesn't grok,
//     so we must explicitly pass a pointer to the beginning element when
//     passing the state vector to GSL
//
const int N = p.N;
std::vector<double> uvec;
for (int i = 1 ; i <= N ; ++i)
{
// initial position of particle /i/
const double xi_init = (i-1)*p.d;
uvec.push_back(xi_init);
}
for (int i = 1 ; i <= N ; ++i)
{
// initial velocity of particle /i/
const double xdoti_init = (i == 1) ? p.v1_init : 0.0;
uvec.push_back(xdoti_init);
}

//
// setup GSL ODE integration
//
gsl_odeiv2_system springs_ODE_system
= {
  springs_rhs, // RHS fn
  nullptr, // Jacobian fn (unused)
  static_cast<size_t>(N_eqns(p)), // number of ODEs to integrate
  static_cast<void*>(&p) // parameters to RHS fn
  };
const gsl_odeiv2_step_type* springs_step_type
= gsl_odeiv2_step_msadams;
const double dt_init = 0.01*p.dt; // guess for initial time step
const double eps_abs = 1.0e-8; // error tolerance (absolute)
const double eps_rel = 1.0e-8; // error tolerance (relative)
gsl_odeiv2_driver* driver
= gsl_odeiv2_driver_alloc_y_new(&springs_ODE_system,
springs_step_type,
dt_init,
eps_abs, eps_rel);

//
// print an output header
//
p.print_me();
printf("# column 1 = t\n");
for (int i = 1, col = 2 ; i <= N ; ++i, ++col)
{
printf("# column %d = x%d\n", col, i);
}

//
// main time evolution
//
double t = 0.0; // starting time
for (int k = 1 ; ; ++k) // k = time step number
{
// print the main output (current time & particle positions)
printf(number_format, t);
for (int i = 1 ; i <= N ; ++i)
{
printf("\t");
printf(number_format, uvec.at(posn_of_xi(p,i)));
}
printf("\n");

if (t >= p.t_max-fp_fuzz)
break; // *** loop control ***

// step to the next output time
const double t_goal = k*p.dt;
const int status = gsl_odeiv2_driver_apply(driver,
   &t, t_goal,
   &uvec.at(0));
if (status != GSL_SUCCESS)
{
fprintf(stderr,
"***** springs: gsl_odeiv2_driver_apply() returned error code %d\n"
"               trying to take time step k=%d t=%.10f t_goal=%.10f\n",
status, k, t, t_goal);
std::exit(1); // *** ERROR EXIT ***
}
}

gsl_odeiv2_driver_free(driver);
return 0;
}

//******************************************************************************

//
// This function prints all the parameters in a /params/ object to /stdout/
// in a (vaguely) human-readable format.  Each output line begins with '##'
// (& hence is a gnuplot comment).
//
void params::print_me()
  const
{
printf("## N=%d\t\t\t\t# number of particles \n", N);

printf("## m=");
printf(number_format, m);
printf("\t\t\t# mass of each individual particle\n");

printf("## d=");
printf(number_format, d);
printf("\t\t\t# equilibrium distance between particles\n");

printf("## k=");
printf(number_format, k);
printf("\t\t\t# (elastic) spring constant\n");

printf("## mu=");
printf(number_format, mu);
printf("\t\t\t# spring friction coefficient\n");

printf("## F1_ext=");
printf(number_format, F1_ext);
printf("\t\t# external force applied to leftmost particle\n");

printf("## v1_init=");
printf(number_format, v1_init);
printf("\t\t# initial velocity of leftmost particle\n");

printf("## right_wall_flag=");
printf(right_wall_flag ? "true" : "false");
printf("\t# is there a fixed \"wall\" to the\n");
printf("##\t\t\t\t# right of the rightmost particle?\n");

printf("## dt=");
printf(number_format, dt);
printf("\t\t\t# time interval at which to print output\n");

printf("## t_max=");
printf(number_format, t_max);
printf("\t\t# time interval at which to end the simulation\n");
}

//******************************************************************************

//
// This function computes the ODE system's RHS.
// This function's API is specified by GSL and must not be changed.
//
// Arguments:
// t = (unused) the time at which to evaluate the RHS
// u = the ODE state vector (a C-style array) at which to evaluate the RHS
// udot = a C-style array into which this function should store the RHS
// void_params_ptr = a pointer to our /params/ object, cast to /void */
//
// Results:
// This function returns /GSL_SUCCESS/ if the computation succeeds.
// If any errors occur, this function prints an error message and calls
// /exit()/ (which doesn't return).
//
extern "C"
  int springs_rhs(const double /* t */,
  const double u[],
  double udot[],
  void* void_params_ptr)
{
const params& p = *static_cast<const params*>(void_params_ptr);
const int    N = p.N;
const double m = p.m;
const double F1_ext = p.F1_ext;

for (int i = 1 ; i <= N ; ++i)
{
udot[posn_of_xi(p,i)] = u[posn_of_xdoti(p,i)];
}

for (int i = 1 ; i <= N ; ++i)
{
const double F_net_i
= (i == 1) ? F1_ext                - spring_force(p,u,i)
   : spring_force(p,u,i-1) - spring_force(p,u,i);
udot[posn_of_xdoti(p,i)] = F_net_i / m;
}

return GSL_SUCCESS;
}

//******************************************************************************

//
// This function computes the restoring force $R_{i,i+1}$ exerted by
// the spring between particles $i$ and $i+1$, or, if /i = N/ and
// /p.right_wall_flag/ is true, between particle N and the "wall".
// If /i = N/ and there's no wall, this function returns 0.0.
//
// Arguments:
// p = a reference to the /params/ object describing the system
// u = a pointer to the state-vector array (a C-style array)
// i = specifies which force is wanted
//
// Results:
// This function returns the restoring force $R_{i,+1}$.
// If any of the inter-particle spacings are (fuzzily) zero or negative,
// this function prints an error message and calls /exit()/ (which doesn't
// return).
//
//
double spring_force(const params& p,
    const double u[], const int i)
{
const int    N  = p.N;
const double d  = p.d;
const double k  = p.k;
const double mu = p.mu;
const bool right_wall_flag = p.right_wall_flag;

assert(i >= 1); // sanity check
assert(i <= N); // ditto

if ((i < N) || right_wall_flag)
{
const double    x_iplus1 = (i < N) ? u[posn_of_xi   (p,i+1)] : N*d;
const double xdot_iplus1 = (i < N) ? u[posn_of_xdoti(p,i+1)] : 0.0;

const double L    =    x_iplus1 - u[posn_of_xi   (p,i)];
const double Ldot = xdot_iplus1 - u[posn_of_xdoti(p,i)];
if (L <= fp_fuzz)
{
fprintf(stderr,
"***** springs: spring from i=%d to %d has length %.10f\n"
"               which is (fuzzily) <= 0.0!\n",
i, i+1, L);
exit(1); // *** ERROR EXIT ***
}

const double R = -k*(L-d) - mu*Ldot;
return R; // *** normal return
}
   else return 0.0; // *** normal return
}
--- end code ---

ciao,
--
-- "Jonathan Thornburg [remove -color to reply]" <dr.j.thornburg@gmail-pink.com>
   (he/him; currently on the west coast of Canada)
       The Three Laws of Thermodynamics:
       1) You can't win, only lose or break even.
       2) You can only break even at absolute zero.
       3) You can't get to absolute zero.