1
$\begingroup$

I was playing around with NDSolve and the following ODE

  s = NDSolve[{y''[x] == -Sin[y[x]], y[0] == \[Pi], y'[0] == 0 }, 
  y, {x, 0, 120}]

The equation and the boundary conditions have an immediate physical interpretation, that is a pendulum with $0$ velocity and rotated $\pi$ degrees from the equilibrium position. This is an unstable equilibrium and I understand that it might not be entirely trivial for Mathematica to "hold" it.

Indeed the solution loses stability pendulum orbit.

Solving instead

  s = NDSolve[{y''[x] == -Sin[y[x]], y[0] == \[Pi]-$MachineEpsilon, y'[0] == 0 }, 
  y, {x, 0, 120}]

I would expect to get the same solution as before, and indeed I do.

I then perturb the initial position by $+$ \$MachineEpsilon, it still computes the same solution. I would have thought it would have gone in the opposite direction, as the little "push" was toward increasing angles. I guess the reasoning does not apply, as we are still within machine precision range . If I perturb by say, $+1.1$ \$MachineEpsilon to get the expected behaviour (swing from $\pi$ to $ 3\pi$).

It all could make sense, but I would like to know what is the proper way to tighten tolerances to the max and hold the unstable position as much as possible. Further, what is exactly 'as much as possible'? I would love to know why the instability loss occurs around 60 seconds. Is there some probabilistic explanation, e.g. every time step the some expression is checked and after some iterations noise gets in the works?

Thanks

I checked somehow similar question, e.g. unstable BVP, but could not come to a conclusion,

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ You can increase the AccuracyGoal. With AccuracyGoal -> 24 the system remains at rest until time 120. Basically, round off errors cause the system to deviate from the unstable equilibrium. Look at LogPlot[Abs[y[x] - \[Pi]] /. s, {x, 0, 60}] $\endgroup$ – Chris K Mar 31 '20 at 21:09
  • 1
    $\begingroup$ as Chris says, at $\pi$ it is unstable equilibrium. This is in theory only. When using a computer, any tiny deviation from this point will result is system moving away from this point. Since you using numerical solver, this is expected to happen. $\endgroup$ – Nasser Apr 1 '20 at 2:08
  • $\begingroup$ Thanks for your replies, @Nasser. I understand that, but I was wondering, why 60 seconds? What relationship to the chosen time step? The AccuracyGoal is exactly what I was looking for, but are there other ways? Thanks $\endgroup$ – Smerdjakov Apr 1 '20 at 5:21
  • $\begingroup$ I can show you analytical solution if you like, but the analytical solution is not defined at exactly $\pi$. It is accurate for very close to $\pi$, as close as you want, but not exactly $\pi$ and hence the analytical solution will also not stay at $\pi$ for all time. I do not know why numerical solution deviates at 60 seconds or any other time. Must be again related to numerical issues. $\endgroup$ – Nasser Apr 1 '20 at 5:23
  • $\begingroup$ Obviously, pi cannot be represented exactly in floating point $\endgroup$ – mikado Apr 1 '20 at 5:30

Your Answer

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

Browse other questions tagged or ask your own question.