Motion of a Massive Spring
Fixed At One End
and with a
Rigid Mass At the Other
Bill Watts

This solves the equation of motion of a spring of mass m/unit length, fixed at one end and containing a mass M at its free end.
The mass M is initially at rest and is displaced by an amount U0 such that the spring is in tension.  The spring has spring constant k, natural length L.
Because the spring has mass, this system is considerably more complicated than the usual mass - spring system of high school physics.
However, Mathematica solves the problem rather easily.

ClearAll["Global`*"];
Off[General::spell]
Off[General::spell1]
SetDirectory[ToFileName[("FileName"/.NotebookInformation[EvaluationNotebook[]])[[1]]]];

$Assumptions=n∈ Integers

We will make use of the differential equation originally applied to longitudinal vibrations of a thin rod.

pdeRod=(∂_x) (Ε A[x](∂_x) u[x,t])==m[x] (∂_(t,t)) u[x,t]

Where
u[x,t] is the longitudinal displacement of the rod from equilibrium as a function of position on the bar and time.
E is Young's modulus
A[x] is the cross-sectional area of the rod
m[x] is the mass/ unit length of the rod.

We convert the equation to that of a spring with the relationship E A[x] -> k L.  Also use contant mass/ unit m length along the spring.

pde=pdeRod//.{A'[x]->0,Ε->((k L)/(A[x])),m[x]->m}

pde=Thread[pde/m,Equal]

Look for separable solutions for u.

u[x_,t_]=X[x]T[t]

eq=Thread[(pde/(X[x]T[t])),Equal]

Since the left side of the above equation depends on x alone and the right side depends on t alone, each side
must be equal to the same constant.  That constant may be positive, negative, or 0, each with different forms of
the solution.

0 constant first
the T side

t0eq=eq[[2]]==0

DSolve[t0eq,T[t],t]//Flatten

T0=T[t]/.%/.{C[1]->c1,C[2]->c2}

The X side

x0eq=eq[[1]]==0

DSolve[x0eq,X[x],x]//Flatten

X0=X[x]/.%/.{C[1]->c3,C[2]->c4}

Now do the same with a negative constant

t1eq=eq[[2]]==-(ω^2)

DSolve[t1eq,T[t],t]//Flatten

T1=T[t]/.%/.{C[1]->c5,C[2]->c6}

x1eq=eq[[1]]==-(ω^2)

DSolve[x1eq,X[x],x]//Flatten

X1=X[x]/.%/.{C[1]->c7,C[2]->c8}

Using a positive constant will give us exponential solutions instead of sinusoidal, and will obviously not
satisfy our boundary and initial conditions, so we do not need to consider such solutions.
Now combine the pieces of the solutions we have found.

u[x_,t_]=X0 T0+X1 T1

Apply the boundary conditions

u[0,t] = 0				The spring is fixed at x = 0
(∂_t) u[x, 0] = 0				The system is initially at rest
kL (∂_x) u[x,L]=-M (∂_(t,t)) u[x,L] 		The force at x = L.  F = k Δx = M × acceleration. The negative sign indicates that a positive force (tension) pulls M in the negative x direction.
u[x,0] = U0 x/L				Initial displacement.  Linear such that the displacement at x = 0 is 0, and the displacement at x = L is U0.
u[x,L] = U0				Seemingly redundant to the above, but necessary to take into account orthogonality of our resulting Fourier series.

The spring is fixed at x = 0

u[0,t] = 0

(u[0,t]/.x->0)==0

To satisfy the above set:

c3=0;
c7=0;

u[x,t]

Too many constants, so combine them by setting:

c4=1;
c8=1;

u[x,t]

The initial velocity is 0

(∂_t) u[x,0]=0

((∂_t) u[x,t]/.t->0)==0

Set:

c2=0;
c6=0;

u[x,t]

At the end with mass M, F = M a.

k L (∂_x) u[L,t]=-M (∂_(t,t)) u[x,L]

(k L(∂_x) u[x,t]/.x->L)==-(M (∂_(t,t))u[x,t]/.x->L)//FullSimplify

Set c1 = 0 and set the expression inside the parentheses to be 0

c1=0;

ωeq=(Sqrt[k]) (Sqrt[L]) (Sqrt[m]) Cos[(((Sqrt[L]) (Sqrt[m]) ω)/(Sqrt[k]))]-M ω Sin[(((Sqrt[L]) (Sqrt[m]) ω)/(Sqrt[k]))]==0

There are an infinite number of ω's and the above is the equation that determines their value.  We will solve numerically
for ω when we have assigned values to the input parameters.  For future use set:

F1[ω_]=%[[1]]

For later simplification.

ωrule= Cos[(Sqrt[((L m)/k)])  ω]->((M ω)/((Sqrt[k]) (Sqrt[L]) (Sqrt[m]) )) Sin[(((Sqrt[L]) (Sqrt[m]) ω)/(Sqrt[k]))]

We need to solve the ωeq numerically when we have assigned k, L, m, and M and sum the solutions over all ω's.

u[x,t]

u[x,0]=U0 x/L

((U0 x)/L)==u[x,0]

Multiply and integrate in the manner of Fourier series.

Thread[m Sin[(((Sqrt[m]) x ω)/((Sqrt[k]) (Sqrt[L])))]%,Equal]

$Assumptions=$Assumptions&&L>0&&m>0&&ω>0&&k>0

eq1=((m U0)/L)(∫_0^L) x Sin[(((Sqrt[m]) x ω)/((Sqrt[k]) (Sqrt[L])))]\[DifferentialD]x==c5 m(∫_0^L) Sin[(((Sqrt[m]) x ω)/((Sqrt[k]) (Sqrt[L])))]Sin[(((Sqrt[m]) x ω)/((Sqrt[k]) (Sqrt[L])))]\[DifferentialD]x

Solving the above equation for c5 will not get us what we want because the right hand side of eq1 does
not satisfy the orthogonality condition for the ω's.
We need to add the seemingly redundant condition

u[L,0] = U0 

because the orthogonality relationship for this system is

(∫_0^L) m u[x, ω1]u[x,ω2]dx + M u[L, ω1] u[L, ω2] = 0

 when ω1 ≠ ω2

U0==u[L,0]

Multiply and add to the original equation.

eq2=Thread[M Sin[(((Sqrt[L]) (Sqrt[m]) ω)/(Sqrt[k]))]%,Equal]

eq=Thread[eq1+eq2,Equal]

Now the right side is c5 × the orthogonality condition.

Solve[eq,c5]//Flatten

c5=c5/.%//Simplify

And do some further simplification

c5=c5/.ωrule//ExpandAll//FullSimplify

c5=(c5//TrigExpand)/.ωrule//Simplify

u[x_,t_]=u[x,t]//Simplify

The velocity of any point on the spring

v[x_,t_]=(∂_t)u[x,t]//FullSimplify

The spring force on any point on the spring.

F[x_,t_]=k L (∂_x)u[x,t]//FullSimplify

Assign indices to the ω's since our solution will be a sum over the ω's.

u[x_,t_]=u[x,t]/.ω->ω[[n]]

v[x_,t_]=v[x,t]/.ω->ω[[n]]

F[x_,t_]=F[x,t]/.ω->ω[[n]]

Since there are an infinite number of ω's our solution needs to be an infinite sum.  We will use a sum with a finite number of terms, however
in order to compute numerical values

u[x_,t_,mm_]:=(∑_(n=1)^mm)u[x,t]

v[x_,t_,mm_]:=(∑_(n=1)^mm)v[x,t]

F[x_,t_,mm_]:=(∑_(n=1)^mm)F[x,t]

Plug in some numbers, so we can make plots.  First consider the case where the mass of the
spring mL is the same as the mass M at the end of the spring.  Recall that m = mass/unit length
which we are taking to be constant.

L=10;
m=1;
M=10;
U0=1;
k=100;

Clear[ω]

Plot[Evaluate[F1[ω]],{ω,0,15}]

The ω's are the roots of F1[ω].  Examine the curve and manually input the first 2 approximate roots into the
first two FindRoots in the routine below.  The routine will compute the number of ω's specifed by the 
value of n.

n=500;
ω=Table[0,{i,n}];
a=FindRoot[F1[α]==0,{α,3},WorkingPrecision->500]
ω[[1]]=α/.a;
a=FindRoot[F1[α]==0,{α,11},WorkingPrecision->500]
ω[[2]]=α/.a;
For[i=3,i<=n,
	a=FindRoot[ F1[α]==0,{α,2ω[[i-1]]-ω[[i-2]]},WorkingPrecision->500 ];
  ω[[i]]=α/.a;	
i++]
Clear[n]

Spot check the orthogonality condition with 2 different ω's

m(∫_0^L) Sin[(((Sqrt[m]) x ω[[1]])/((Sqrt[k]) (Sqrt[L])))]Sin[(((Sqrt[m]) x ω[[2]])/((Sqrt[k]) (Sqrt[L])))]\[DifferentialD]x+M Sin[(((Sqrt[m]) L  ω[[1]])/((Sqrt[k]) (Sqrt[L])))]Sin[(((Sqrt[m]) L  ω[[2]])/((Sqrt[k]) (Sqrt[L])))]

Pretty close to zero.

Compute the period of oscillation.  Since we are stretching the end of the spring a certain distance and releasing it,
the displacement should be zero at (1/4) of a period.

FindRoot[u[L,t,100]==0,{t,.5}]

And the period becomes.

T=4t/.%

Check the initial displacement at x = L.  Should be close to U0.

u[L,0,50]

The displacement of the end of the spring as a function of time.

Plot[Evaluate[u[L,t,100]],{t,0,T},PlotRange->All,AxesLabel->{"t","u"},PlotLabel->"Displacement of mass M vs t"]

Check solution against the given initial condition.

Plot[{((U0 x)/L),Evaluate[u[x,0,100]]},{x,0,L},PlotRange->All,PlotLabel->"Displacement at t = 0",AxesLabel->{"x","u"}]

Check that F = k U0 initially.

Plot[Evaluate[F[x,0,200]],{x,0,L},PlotRange->{0,110},AxesLabel->{"x","Force"},PlotLabel->"Initial spring force"]

The velocity of the end of the spring.

Plot[Evaluate[v[L,t,400]],{t,0,T},PlotRange->All,AxesLabel->{"t","Velocity"},PlotLabel->"Velocity of mass M"]

The tension at the fixed end of the spring. x = 0

Plot[Evaluate[F[0,t,400]],{t,0,2T},AxesLabel->{"t","Force"},PlotLabel->"Force at fixed end"]

Plot[Evaluate[F[L,t,400]],{t,0,2T},AxesLabel->{"t","Force"},PlotLabel->"Spring Force at x = L"]

The above shows the effects of compression/tension waves traveling along the spring.

The tension (spring force) as a function of position and time.

Table[Plot[Evaluate[F[x,t,200]],{x,0,L},PlotRange->{{0,L},{-150,150}},ImageSize->{600,400},AxesLabel->{"x","Force"},PlotLabel->"Spring Force vs x and t"],{t,0,2T,T/48}];
ListAnimate[%]

The above shows that the spring is not in uniform tension or compression along its length other than initially.

The velocity along the length of the spring as a function of time.

Table[Plot[Evaluate[v[x,t,100]],{x,0,L},PlotRange->{{0,L},{-4,4}},ImageSize->{600,400},AxesLabel->{"x","Velocity"},PlotLabel->"Velocity vs x and t"],{t,0,2T,T/24}];
ListAnimate[%]

And the displacement

Table[Plot[Evaluate[u[x,t,200]],{x,0,L},PlotRange->{{0,L},{-1,1}},ImageSize->{600,400},AxesLabel->{"x","Displacement"},PlotLabel->"Displacement vs x and t"],{t,0,2T,T/24}];
ListAnimate[%]

Let's try the case where the mass of the spring is very small compared to the mass at the end.
This case will approach the simple harmonic motion of the mass-spring system of high school physics.

L=10;
m=.0001`500;
M=10;
U0=1;
k=100;

Clear[ω]

Plot[Evaluate[F1[ω]],{ω,0,1000}]

n=500;
ω=Table[0,{i,n}];
a=FindRoot[F1[α]==0,{α,3},WorkingPrecision->500]
ω[[1]]=α/.a;
a=FindRoot[F1[α]==0,{α,1000},WorkingPrecision->500]
ω[[2]]=α/.a;
For[i=3,i<=n,
	a=FindRoot[ F1[α]==0,{α,2ω[[i-1]]-ω[[i-2]]},WorkingPrecision->500 ];
  ω[[i]]=α/.a;	
i++]
Clear[n]

In this case the motion is so close to simple harmonic motion(which has only one fundamental frequency)
that only ω[[1]] has much effect on the solution.  However, we will use more ω's to show this.

Spot check the orthogonality condition with 2 different ω's

m(∫_0^L) Sin[(((Sqrt[m]) x ω[[1]])/((Sqrt[k]) (Sqrt[L])))]Sin[(((Sqrt[m]) x ω[[2]])/((Sqrt[k]) (Sqrt[L])))]\[DifferentialD]x+M Sin[(((Sqrt[m]) L  ω[[1]])/((Sqrt[k]) (Sqrt[L])))]Sin[(((Sqrt[m]) L  ω[[2]])/((Sqrt[k]) (Sqrt[L])))]

Pretty close to zero.

Compute the period of oscillation.  Since we are stretching the end of the spring a certain distance and releasing it,
the displacement should be zero at (1/4) of a period.

One term

FindRoot[u[L,t,1]==0,{t,.5}]

Fifty terms

FindRoot[u[L,t,50]==0,{t,.5}]

And the period becomes.

T=4t/.%

Check the initial displacement at x = L.  Should be close to U0.

One term

u[L,0,1]//N[#,10]&

Fifty terms

u[L,0,50]//N[#,10]&

The displacement of the end of the spring as a function of time.

Plot[Evaluate[u[L,t,50]],{t,0,T},PlotRange->All,AxesLabel->{"t","u"},PlotLabel->"Displacement of mass M vs t"]

Check solution against the given initial condition.

Plot[{((U0 x)/L),Evaluate[u[x,0,50]]},{x,0,L},PlotRange->All,PlotLabel->"Displacement vs x at t = 0",AxesLabel->{"x","u"}]

Check that F = k U0 initially.

Plot[Evaluate[F[x,0,50]],{x,0,L},PlotRange->{0,110},AxesLabel->{"x","Force"},PlotLabel->"Initial spring force"]

The velocity of the end of the spring.

Plot[Evaluate[v[L,t,400]],{t,0,T},PlotRange->All,AxesLabel->{"t","Velocity"},PlotLabel->"Velocity of mass M vs t"]

The tension at the fixed end of the spring.

Plot[Evaluate[F[0,t,400]],{t,0,2T},AxesLabel->{"t","Force"},PlotLabel->"Spring Force at fixed end vs t"]

The force on mass M at the end of the spring

Plot[Evaluate[F[L,t,400]],{t,0,2T},AxesLabel->{"t","Force"},PlotLabel->"Spring Force at free end vs t"]

The tension (spring force) as a function of position and time.

Table[Plot[Evaluate[F[x,t,200]],{x,0,L},PlotRange->{{0,L},{-150,150}},ImageSize->{600,400},AxesLabel->{"x","Force"},PlotLabel->"Spring Force"],{t,0,2T,T/48}];
ListAnimate[%]

From the above plot, the spring maintains nearly uniform tension or compression along its length.

The velocity along the length of the spring as a function of time.

Table[Plot[Evaluate[v[x,t,100]],{x,0,L},PlotRange->{{0,L},{-4,4}},ImageSize->{600,400},AxesLabel->{"x","Velocity"},PlotLabel->"Velocity"],{t,0,2T,T/24}];
ListAnimate[%]

And the displacement

Table[Plot[Evaluate[u[x,t,200]],{x,0,L},PlotRange->{{0,L},{-1,1}},ImageSize->{600,400},AxesLabel->{"x","Displacement"},PlotLabel->"Displacement"],{t,0,2T,T/24}];
ListAnimate[%]

Now try the case where the mass on the end is very small compared to the mass of the spring.

L=10;
m=1;
M=0.0001`500;
U0=1;
k=100;

Clear[ω]

Plot[Evaluate[F1[ω]],{ω,0,15}]

n=500;
ω=Table[0,{i,n}];
a=FindRoot[F1[α]==0,{α,5},WorkingPrecision->500]
ω[[1]]=α/.a;
a=FindRoot[F1[α]==0,{α,15},WorkingPrecision->500]
ω[[2]]=α/.a;
For[i=3,i<=n,
	a=FindRoot[ F1[α]==0,{α,2ω[[i-1]]-ω[[i-2]]},WorkingPrecision->500 ];
  ω[[i]]=α/.a;	
i++]
Clear[n]

Spot check the orthogonality condition with 2 different ω's

m(∫_0^L) Sin[(((Sqrt[m]) x ω[[1]])/((Sqrt[k]) (Sqrt[L])))]Sin[(((Sqrt[m]) x ω[[2]])/((Sqrt[k]) (Sqrt[L])))]\[DifferentialD]x+M Sin[(((Sqrt[m]) L  ω[[1]])/((Sqrt[k]) (Sqrt[L])))]Sin[(((Sqrt[m]) L  ω[[2]])/((Sqrt[k]) (Sqrt[L])))]

Pretty close to zero.

Compute the period of oscillation.  Since we are stretching the end of the spring a certain distance and releasing it,
the displacement should be zero at (1/4) of a period.

FindRoot[u[L,t,100]==0,{t,.5}]

And the period becomes.

T=4t/.%

Check the initial displacement at x = L.  Should be close to U0.

u[L,0,500]

The displacement of the end of the spring as a function of time.

Plot[Evaluate[u[L,t,100]],{t,0,T},PlotRange->All,AxesLabel->{"t","u"},PlotLabel->"Displacement of the free end vs t"]

Interesting sawtooth shape.

Check solution against the given initial condition.

Plot[{((U0 x)/L),Evaluate[u[x,0,100]]},{x,0,L},PlotRange->All,PlotLabel->"Displacement at t = 0",AxesLabel->{"x","u"}]

Check that F = k U0 initially.

Plot[Evaluate[F[x,0,200]],{x,0,L},PlotRange->{0,110},AxesLabel->{"x","Force"},PlotLabel->"Initial spring force"]

The velocity of the end of the spring.

Plot[Evaluate[v[L,t,400]],{t,0,T},PlotRange->All,AxesLabel->{"t","Velocity"},PlotLabel->"Velocity of free end vs t"]

The tension at the fixed end of the spring.

Plot[Evaluate[F[0,t,400]],{t,0,2T},AxesLabel->{"t","Force"},PlotLabel->"Spring Force of fixed end vs t"]

Plot[Evaluate[F[L,t,400]],{t,0,2T},AxesLabel->{"t","Force"},PlotLabel->"Spring Force of free end vs t"]

With very small mass at the free end, the force is very nearly zero.

The tension (spring force) as a function of position and time.

Table[Plot[Evaluate[F[x,t,200]],{x,0,L},PlotRange->{{0,L},{-150,150}},ImageSize->{600,400},AxesLabel->{"x","Force"},PlotLabel->"Spring Force"],{t,0,2T,T/48}];
ListAnimate[%]

The velocity along the length of the spring as a function of time.

TP1=Table[Plot[Evaluate[v[x,t,100]],{x,0,L},PlotRange->{{0,L},{-4,4}},ImageSize->{600,400},AxesLabel->{"x","Velocity"},PlotLabel->"Velocity"],{t,0,2T,T/24}];
ListAnimate[TP1]

And the displacement

TP1=Table[Plot[Evaluate[u[x,t,200]],{x,0,L},PlotRange->{{0,L},{-1,1}},ImageSize->{600,400},AxesLabel->{"x","Displacement"},PlotLabel->"Displacement"],{t,0,2T,T/24}];
ListAnimate[TP1]