Morin 4.22 (Projectile on Spring)

A projectile of mass m is fired from the origin with speed v0 and angle θ0 in the presence of gravity g.
It is also attached to a zero-length massless spring with constant k=m ω^2 (with the other end at the origin).
(a) Find x(t) and y(t)
(b) Show that the motion is the expected parabola when the spring is weak, and the expected harmonic motion when the spring is strong.  (And note that "weak"and "strong" refer to the strength relative to gravity.)
(c) Find the ω such that the projectile hits the ground going straight down.

Begin as always with the forces: there is a -mg in the y direction, and a spring force of magnitude k*r
with x-component -kr*Cos[θ] = -kx and y component -ky:

eqnx = m x''[t] == -m ω^2 x[t]
eqny = m y''[t] == -m ω^2 y[t] - m g

We launch with given conditions:

s = DSolve[{eqnx,eqny,x[0]==0,y[0]==0,x'[0]== v0 Cos[θ0], y'[0]== v0 Sin[θ0]},{x[t],y[t]},t]

Let's package up the solution in the form of a vector-valued function of time:

pt[t_] = {x[t],y[t]}/.s[[1]]

When ω is small we get:

Limit[pt[t],ω->0]

which sure does look like regular old projectile motion.
And when ω is large, or more precisely when gravity is weak, we get:

Limit[pt[t],g->0]

which is simple harmonic oscillation along an axis tilted at θ0.

To be a little more precise about "large" and "small", this problem has a scale defined by
the time it takes a projectile to go up and down, i.e. T=2 v0 Sin[θ0] /g.
The spring parameter ω should be compared to 1/T.
So let us define the dimensionless parameter α=ω T and express ω as:

αrule = {ω-> α g/(2 v0 Sin[θ0])}

With this replacement, our solution is:

sα = DSolve[{eqnx,eqny,x[0]==0,y[0]==0,x'[0]== v0 Cos[θ0], y'[0]== v0 Sin[θ0]}/.αrule,{x[t],y[t]},t]//Simplify

ptα[t_] = {x[t],y[t]}/.sα[[1]]

And now we can take the more meaningful limits:

Limit[ptα[t],α->0]

This one gets stuck, despite all the assumptions

Limit[ptα[t],α->Infinity]

$Assumptions = {t>0,v>0,θ0>0,θ0<π/2,g>0,α>0}

$Assumptions=.

OK Mathematica got confused by the Sin[] and Cos[], but we can see that for large α the αSin[] term dominates the first two terms in y[t], yielding the harmonic motion above.

To determine the magical value of ω which makes the mass hit the ground straight down, we first find the time when y=0:

sT = Solve[pt[t][[2]]==0,t]

We want the positve time:

T = t/.sT[[3]]

And at that time we want the x-velocity:

pt'[T][[1]]

which we demand to be zero:

Solve[%==0,ω]

In other words we want our dimensionless α to be 2.
Finally, make some plots:

p = {g->1,v0->1,θ0-> π/3,α->2}

Find the final time:

Tf = T/.αrule /.p

Draw the curve, checking that it does crash to the ground vertically:

curve = ParametricPlot[ptα[t]/.p,{t,0,Tf}]

Here is a little function which creates a unit vector perpendicular to the input vector:

e2[{x_,y_}] = {-y,x}/Sqrt[x^2+y^2]

which we use to draw the wiggles in the spring:

draw[t_] := Show[
curve,
ListPlot[{ptα[t]/.p},PlotStyle->{PointSize[.03],Red}],
ParametricPlot[(s*ptα[t]+.003*Sin[20 π s]*e2[ptα[t]])/.p,{s,0,1},PlotStyle->Green]
]

Animate[draw[t],{t,0,Tf}]

By the way, if we ignore the ground and keep going, we get an ellipse centered at the true equilibrium point,
which is a distance mg/k below the ground (that being the place where a mass could hang motionless).
In other words, the exact motion is always just 2D simple harmonic motion; the effect of gravity is just
to shift the equilibrium point.

ParametricPlot[ptα[t]/.p,{t,0,10Tf},AspectRatio->2.5,PlotRange->{{-.5,.5},{-2,.5}}]