0
$\begingroup$

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]}$

plots showing a parabola over a 400x cosine

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

enter image description here


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]

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

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

$\endgroup$
0
$\begingroup$

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}] 

enter image description here

$\endgroup$
13
  • $\begingroup$ 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? $\endgroup$
    – lineage
    Jun 25 at 7:05
  • $\begingroup$ strange.. ran your code ditto....its giving a parabola instead....yep re ran still parabola $\endgroup$
    – lineage
    Jun 25 at 7:05
  • $\begingroup$ btw the curve you show isn't the expected answer for the given params...the mass should stay $\approx 1$ $\endgroup$
    – lineage
    Jun 25 at 7:08
  • $\begingroup$ Some strange things happened. After correcting the potential the result is as expected. See my modified answer. $\endgroup$ Jun 25 at 7:24
  • $\begingroup$ 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 $\endgroup$
    – lineage
    Jun 25 at 7:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.