Springs and Balls

AP1601
Columbia University

One ball; Two springs

Imagine a ball of mass, m, suspended between two fixed walls by two ideal springs.

As the ball moves, the springs exert restoring forces that accelerate the ball towards the origin, (x,y) = (0,0). When the ball is situated exactly at the origin, the ball is in "equilibrium" because there is no net force on the ball due to the springs. (And, for the first time, we will be ignoring gravity in our physics modeling.)

The force on the ball by an ideal spring is given by the simple formula
	F = - k (l – (Subscript[l, 0]))
where (Subscript[l, 0]) is the "unstreched" length of the spring and k is the ideal spring constant.

If we restrict the ball to move only the (x,y) plane, then the equations of motion for the ball are:
	(∂_t) x[t]  =  vx[t]
	(∂_t) y[t]  =  vy[t]
	(∂_t) vx[t]  =  ((Δ + x)/(Subscript[l, 1]))(Subscript[a, 1]) - ((Δ-x)/(Subscript[l, 2]))(Subscript[a, 2])
	(∂_t) vy[t]  =  (y/(Subscript[l, 1]))(Subscript[a, 1]) + (y/(Subscript[l, 2]))(Subscript[a, 2])
where the springs are anchored at either +Δ or -Δ, and the length of each spring is
	(Subscript[l, 1]) = (Sqrt[(((Δ + x))^2) + (y^2)])
	(Subscript[l, 2]) =  (Sqrt[(((Δ - x))^2) + (y^2)])
Notice that ratios of the length, y/l and (Δ+x)/l, multiply the accelerations, a, because the represent the ratios of the y and x directed components of the force to the total force.
The acceleration (i.e. the force divided by the ball's mass, m) applied by each spring is
	(Subscript[a, 1]) = - (k/m)((Subscript[l, 1]) - (Subscript[l, 0]))
	(Subscript[a, 2]) = - (k/m)((Subscript[l, 2]) - (Subscript[l, 0]))

We can determine the motion of the ball by using NDSolve and a particular initial condition.

Equations of Motion (One ball)

p = {x[t],y[t],vx[t],vy[t]}

l1 =(Sqrt[(((Δ + x[t]))^2)+ ((y[t])^2)]);
l2 =(Sqrt[(((Δ - x[t]))^2)+ ((y[t])^2)]);
a1 = - (k/m) (l1 - l0);
a2 = - (k/m) (l2 - l0);
dp = {vx[t], vy[t],(Δ+x[t]) a1/l1-(Δ-x[t]) a2/l2,y[t] a1/l1 + y[t] a2/l2}

p0 = {x0, y0,vx0, vy0};

Thread[{x[0], y[0], vx[0], vy[0]} == p0]

eqs = Thread[D[p,t] ==  dp]~Join ~ Thread[{x[0], y[0], vx[0], vy[0]} ==  p0]

ball1 = NDSolve[eqs /. k->1 /. m -> 1/.Δ -> 1 /. l0 -> 0.5/.x0 -> 0.1 /. y0 -> 0.05 /. vy0->0/. vx0 -> 0, p, {t, 0, 50}]

Basic Graphics

Plot[Evaluate[{x[t],y[t]} /. ball1], {t, 0, 50}, PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}}]

ParametricPlot[Evaluate[{x[t], y[t]} /. ball1], {t, 0, 50},AspectRatio -> Automatic]

Animation

ListAnimate[Table[ListPlot[Evaluate[{x[t],y[t]}/. ball1],
			PlotRange->{{-0.1,0.1},{-0.1,0.1}},
			PlotStyle->PointSize[0.04],DisplayFunction->Identity],{t,0,50,0.2}]]

Simplification of the One Ball Motion

For small displacements of the ball,

eqs1 = Series[dp,{x[t],0,1},{y[t],0,1}]

Linear motion of the ball is described by a very simple set of equations.

dp1 = {vx[t],vy[t],
		-2 (k/m) x[t],
		-2 (k/m) (1 - l0/Δ) y[t]}

eqs1a = Thread[D[p,t] == dp1]

eqs1a = eqs1a /. k -> m (w^2)/2 // Simplify

eqs1a = eqs1a /. (w^2)(l0 - Δ)/Δ -> -(wy^2)

dsol1 = DSolve[eqs1a~ Join ~
		Thread[{x[0], y[0], vx[0], vy[0]} == {x0,y0,0,0}],p,t] // ExpToTrig

Mathematica does not know what type of numbers are represented by w, wy, x0, and y0. If we assign particular real integers to these constants, we can see more easily the form of the linear solution.

dsol1 /. w -> 1 /. wy -> 1/5 /. x0 -> 1/20 /. y0 -> 1/10// Simplify

Notice that there are two independent motions: the x and the y-directions.

Two balls; Three springs

Let's add a ball and spring. To simplify our problem, we will restrict the motion of the balls to the y direction.

The equations of motion are:
	(∂_t) y[1][t]  =  vy[1][t]
	(∂_t) y[2][t]  =  vy[2][t]
	(∂_t) vy[1][t]  =  - (ω^2) y[1][t] -  (ω^2) (y[1][t] - y[2][t])
	(∂_t) vy[2][t]  =  - (ω^2) y[2][t] -  (ω^2) (y[2][t] - y[1][t])
where (ω^2)= (k/m) (1 - l0/Δ). If the position of one of the balls is held fixed, y = 0, these equations reduce to our previous equations.

Equations of Motion (Two balls)

p = {y[1][t],vy[1][t],y[2][t],vy[2][t]}

dp = {vy[1][t],-(ω^2) (2 y[1][t] - y[2][t]),vy[2][t], -(ω^2)(2 y[2][t] - y[1][t])}

p0 = {y0[1], vy0[1],y0[2], vy0[2]};

eqs = Thread[D[p,t] == dp]~ Join ~Thread[{y[1][0], vy[1][0], y[2][0], vy[2][0]} == p0]

ball2 = NDSolve[eqs /.ω->1 /.y0[1] -> 1 /.y0[2] -> 0 /. vy0[1]->0/. vy0[2] -> 0, p, {t, 0, 50}] // Flatten

Basic Graphics

Plot[Evaluate[{y[1][t],y[2][t]} /. ball2], {t, 0, 50},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}}]

Animation

ListAnimate[Table[ListPlot[Evaluate[{{-.7,y[1][t]},{.7,y[2][t]}}/. ball2],
			PlotRange->{{-1.5,1.5},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,DisplayFunction->Identity],{t,0,50,0.2}]]

Many Balls & Springs

Generalizing to N balls (and N + 1 springs):

The equations of motion are:
	(∂_t) y[i][t]  =  vy[i][t]  (* for all i *)
	(∂_t) vy[1][t]  =  - 2(ω^2) y[1][t] +  (ω^2)  y[2][t]
	(∂_t) vy[i][t]  =    (ω^2)  y[i-1][t] - 2(ω^2) y[i][t] +  (ω^2)  y[i+1][t]
	(∂_t) vy[N][t]  =  (ω^2)  y[N-1][t] - 2(ω^2) y[N][t] 
For simplicity, we can normalize time to 2π/ω which is equivalent to assigning ω-> 1.

Equations of Motion (N balls)

nBalls = 30;

p = Table[{y[i][t],vy[i][t]}, {i,1,nBalls}] // Flatten

dp = {vy[1][t],-2 y[1][t] + y[2][t]}~ Join ~Table[{vy[i][t],y[i-1][t]-2 y[i][t] + y[i+1][t]}, {i, 2, nBalls-1}] ~Join~{vy[nBalls][t], y[nBalls-1][t]-2 y[nBalls][t]} // Flatten;

aBall = 3;
p0 = Table[{y[i][0] == If[i == aBall,1,0], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

Animation

ListAnimate[Table[ListPlot[Evaluate[Table[y[i][t],{i,1,nBalls}]/. ballN],
			PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
			DisplayFunction->Identity],{t,0,100,0.5}]]

Practice

Different Initial Conditions

Examine wave motion with different initial conditions. For example, start with a "Gaussian" or a "plucked guitar string", where a Gaussian vary's like Exp[-(i-i0)^2/(di^2)] and a plucked string varies linearly with distance, If[i>i0, i/i0, 1 - (i-i0)/(nBalls-i0)].

Gaussian Example

The Gaussian example is:

i0 = 10;
di = 4;
initial[i_] := Exp[-(i-i0)^2/(di^2)]

p0 = Table[{y[i][0] == initial[i], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

Animation

ListAnimate[Table[ListPlot[Evaluate[Table[y[i][t],{i,1,nBalls}]/. ballN],
			PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
			DisplayFunction->Identity],{t,0,100,0.5}]]

Frequency Spectrum

The frequency spectrum of the string can be displayed using the Fourier built-in function. Read the online help about the Fourier function, and use it to determine the frequency spectrum of the sring for various initial conditions.

Example

The example below shows the spectrum for a gaussian initial condition.

aBallList = Table[Evaluate[y[10][t] /. ballN], {t, 0, 100, 0.2}];

ListPlot[aBallList]

ballFourier = Fourier[aBallList];

ListPlot[Abs[ballFourier],PlotRange->{{0,60},{0,3}},Joined->True]

The frequency spacing is equal to 2π/N, where N is the number os samples.

Short[aBallList]

aBallFreq = Table[{2 π (i-1)/500, Abs[ballFourier[[i]]]}, {i, 1, 100}];

ListPlot[aBallFreq,Joined->True,PlotRange->{0,3},AxesLabel->{"Normalized Frequency"," "}]

Notice, the natural frequency of the "oscillations" on the string is much lower than the natural frequency of a single ball, ω = 1.

Driven String

Consider one end of the string of balls driven by an oscillator or a pulsation.
How would you change the equations such that one of the ends changes as a function of time?

Solution

dp = {vy[1][t],-2 y[1][t] + y[2][t] + drive[t]}~ Join ~ Table[{vy[i][t],y[i-1][t]-2 y[i][t] + y[i+1][t]}, {i, 2, nBalls-1}] ~Join~
		{vy[nBalls][t], y[nBalls-1][t]-2 y[nBalls][t]} // Flatten;

drive1[t_] := If[t<π, Sin[t],0]

p0 = Table[{y[i][0] ==0, vy[i][0]==0}, 
		{i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs /. drive[t] -> drive1[t] , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

Animation

ListAnimate[Table[ListPlot[Evaluate[Table[y[i][t],{i,1,nBalls}]/. ballN],
			PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
			DisplayFunction->Identity],{t,0,100,0.5}]]

Friction

How does the string oscillations change when the oscillations of the balls are damped with friction?

The equations of motion are changed to:
	(∂_t) y[i][t]  =  vy[i][t]  (* for all i *)
	(∂_t) vy[1][t]  =  - 2(ω^2) y[1][t] +  (ω^2)  y[2][t] - k vy[1][t]
	(∂_t) vy[i][t]  =    (ω^2)  y[i-1][t] - 2(ω^2) y[i][t] +  (ω^2)  y[i+1][t] - k vy[i][t]
	(∂_t) vy[N][t]  =  (ω^2)  y[N-1][t] - 2(ω^2) y[N][t]  - k vy[N][t]
We can still normalize time to 2π/ω which is equivalent to assigning ω-> 1 and expressing k (usally much less than one) as the damping time relative to the natural frequency of a single ball.

Solution

dp = {vy[1][t],-2 y[1][t] + y[2][t] -k vy[1][t]}~ Join ~ Table[{vy[i][t],y[i-1][t]-2 y[i][t] + y[i+1][t] -k vy[i][t]}, {i, 2, nBalls-1}] ~Join~{vy[nBalls][t], y[nBalls-1][t]-2 y[nBalls][t] -k vy[nBalls][t]} // Flatten;

p0 = Table[{y[i][0] == initial[i], vy[i][0]==0}, {i,1,nBalls}] // Flatten;

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

ballN = NDSolve[eqs /.k ->0.05 , p, {t, 0, 100}] // Flatten;

Plot[Evaluate[y[3][t] /. ballN], {t, 0, 100}]

Animation

ListAnimate[Table[ListPlot[Evaluate[Table[y[i][t],{i,1,nBalls}]/. ballN],
			PlotRange->{{0,nBalls+1},{-1.2,1.2}},
			PlotStyle->PointSize[0.04],
			Frame -> True,
			Ticks -> None,
			FrameTicks -> None,
			DisplayFunction->Identity],{t,0,100,0.5}]]