NDSolve not finding the right solution for a mass hanging of a thread

My NDSolve isn't predicting the right trajectory for one particular case of a problem.

System

Consider a mass vertically hanging by an almost inelastic thread. The thread's force acts acting to the potential

springV[x_] := If[x > l, k (x - l)^2/2, If[x < -l, k (-x - l)^2/2, 0]]


where x is the mass' position. Basically if the mass stretches the string beyond a length l, the string pulls back like a Hook spring$$^1$$ while if it goes limp, it exerts no force at all.

Let's say the string has an unstretched lengthl. When the mass starts to fall under gravity from initial positions greater than or less than l, NDSolve, gives the right trajectory. However for what is IMO the simplest case $$-$$ starting from rest at l, NDSolve fails. The solution it produces is one of pure free fall instead of the expected oscillations about equilibrium.

Code

(*aliases*)
echo = Echo;
cls = ClearAll;
first = First;
flatAll[x___] := Flatten[x, Infinity];

(*package - using lagrangian*)
Needs@"VariationalMethods";

Module[
{springV, eqns, initConds, x0, v0, x, v, L},

(*the spring potential*)
springV[x_] :=
If[x > l, k (x - l)^2/2, If[x < -l, k (-x - l)^2/2, 0]];

(*the lagrangian - T-V; V=0 @x1=0 +ve down*)
L = m1 x1'[t]^2/2 - springV[x1[t]] + m1 g x1[t];

(*the eqn to solve ma=mg-T*)
eqns = EulerEquations[L, {x1[t]}, t];
(*inital conditions*)
initConds = {x1[0] == x10, x1'[0] == v10};
(*combine them*)
eqns = flatAll@{eqns, initConds} // echo;
(*
sim params: everything is one except
gravity(=9.8)
inital velocity(=0)
k=10^4
*)
vals = {m1 -> 1, k -> 10^4, l -> 1, x10 -> 1., v10 -> 0, g -> 9.8};
sol = NDSolve[eqns /. vals // echo, {x1[t]}, {t, 0, 100}];
xsol = {x1[t]} /. first@sol
]


Getting $$\ldots$$

The yielded solution is a parabola instead of the expected$$^2$$ oscillations about some $$x_0^{[3]}$$

Unmagnified, the solution is 1.00098 - 0.00098 Cos[100. t]$$^4$$ (ignore the plotting artifact around the origin in the magnified expected solution)

Even though in the above I have taken an almost inelastic string by taking a large $$k$$ relative to $$m g$$, the same problem exists even with k=1.

Questions

1. Why is NDSolve not giving the right answer and how to correct it?

I have (blindly) tried basic NDSolve variants under "Possible Issues" from docs. to no avail.

Test cases (in case your results differ)

In case you have come up with your own code, do test to see it matches the behavior at other x10 values

Footer

$$^1$$ For a hook spring the restorative force is $$\propto$$ $$-$$extension$$^1$$ i.e. its linear.

$$^2$$using DSolve[{m1 x''[t] == m1 g - k (x[t] - l), x[0] == x10, x'[0] == v10} /. vals, x[t], t]

$$^3x_0=l+mg/k=1.00098$$ for the given sim

$$^4$$ i.e. it oscillates about the equilibrium position

final modification

potential springV[x] isn't unique at x==1. NDSolve has problems to solve with starting value x1[0]==1.

Try option Method -> {Automatic, "DiscontinuityProcessing" -> False}

sol = NDSolve[eqns /. vals // echo,{x1[t]}, {t, 0, 100} ,
Method -> {Automatic,"DiscontinuityProcessing" -> False}];
xsol = {x1[t]} /. first@sol]
Plot[xsol, {t, 0, 10}, PlotRange ->{-1.1, 1.1}]
`

• what was the issue with my code? I rechecked it .. seems to be working (wrongly) on my machine...btw running Mathematica 12.2....were the err codes any help? Jun 25 at 7:05
• strange.. ran your code ditto....its giving a parabola instead....yep re ran still parabola Jun 25 at 7:05
• btw the curve you show isn't the expected answer for the given params...the mass should stay $\approx 1$ Jun 25 at 7:08
• Some strange things happened. After correcting the potential the result is as expected. See my modified answer. Jun 25 at 7:24
• at least i am getting what you are getting...but the solution you got is still wrong (do carefully observe the sol I have provided)....I'll add some test cases Jun 25 at 7:33