I would like to simulate a point mass that is moving on a circular rail considering kinetic friction. The mass point starts at the bottom of the rail with an initial velocity v0. for the derivation of the algebraic differential equation system i used standard newtonian mechanics. fn is the normal force, that acts on the mass point and ϕ is the only degree of freedom. My computed solution behaves quite weird, e.g. the point mass never stops or it moves in a strange way (see Animation). Here is my code:

v0 := 14
m := 1
g := 9.81
μ := 0.05
R := 5
ϕ0 := 0
tmax := 50
eqn1[t_] := -fn[t] +m g Cos[ϕ[t]]^2 == -m R ϕ'[t]^2
eqn2[t_] := -μ Abs[fn[t]]*Sign[ϕ'[t]] -m g Sin[ϕ[t]] == m R ϕ''[t]
sol = NDSolve[{eqn1[t],eqn2[t], ϕ[0] == ϕ0, ϕ'[0] ==v0/R},
{ϕ,fn}, {t, 0, tmax}, MaxSteps -> Infinity, MaxStepSize -> 0.001,AccuracyGoal -> 6, PrecisionGoal -> 6,InterpolationOrder -> All][[1]];
Plot[ϕ[t] /. sol, {t, 0, tmax}]
Plot[fn[t] /. sol, {t, 0, tmax}, PlotRange -> All]
Plot[(fn[t] - g Cos[ϕ[t]]^2 - R ϕ'[t]^2) /. sol, {t,0,tmax},PlotRange -> All]
Plot[(-μ Abs[fn[t]]*Sign[ϕ'[t]] - g Sin[ϕ[t]]-R ϕ''[t]) /. sol, {t, 0, tmax}, PlotRange -> All]
phi[t_] := ϕ[t] /. sol
Animate[ParametricPlot[{R*Sin[phi[t]], -R*Cos[phi[t]]}, {t, 0, tmax},Epilog ->{Red, PointSize[Medium], Point[{R*Sin[phi[t]], -R*Cos[phi[t]]}]},PlotRange -> All, PerformanceGoal -> "Quality"], {t, 0,tmax}, AnimationRunning -> False]`

How can I improve the code, so the simulation gets appropriate?

  • $\begingroup$ Do you have a problem in writing the equations, or in solving an ODE with Mathematica? In the first case, your question does not concern Mathematica and is off-topic; in the second, could you provide the ODE you want to solve and explain why you believe Mathematica is wrong at solving it? $\endgroup$ – anderstood Dec 21 '16 at 21:56

The code can be improved in several ways. Here is my version, using ListAnimate instead of Animate to get better control of the individual frames. I also speed up the plotting by using Circle instead of ParametricPlot. This allows me to squeeze in $1000$ animation frames. The entire calculation is now done like this:

 {v0 = 14,
  m = 1,
  g = 9.81,
  μ = 0.05,
  R = 5,
  ϕ0 = 0,
  tmax = 50},
 eqn[t_] := -μ Abs[m g Cos[ϕ[t]] 
    + m R ϕ'[t]^2  ]*Sign[ϕ'[t]] - m g Sin[ϕ[t]] == m R ϕ''[t];
 phi = NDSolveValue[{eqn[t], ϕ[0] == ϕ0, 
                   ϕ'[0] == v0/R}, ϕ, {t, 0, tmax}, MaxSteps -> Infinity, 
   MaxStepSize -> 0.001, AccuracyGoal -> 6, PrecisionGoal -> 6, 
   InterpolationOrder -> All];
 Column[{Plot[phi[t], {t, 0, tmax}],
    Table[Graphics[{Circle[{0, 0}, R], {Red, PointSize[Medium], 
        Point[{R*Sin[phi[t]], -R*Cos[phi[t]]}]}}], {t, 0, tmax,tmax/1000}],
        AnimationRunning -> False]}]



In the differential equation, there was no need to have two equations because the normal force is really obtainable purely from $\phi$ alone. I also think you had a square in the normal force that doesn't belong there.

| improve this answer | |
  • $\begingroup$ oh… i didn't see the square in equation for the normal force. thank you for mentioning! $\endgroup$ – freddy90 Dec 21 '16 at 22:33
  • $\begingroup$ one more question. For some reason the "faster"-button in the animation doesn't work. is it just me or why is that? $\endgroup$ – freddy90 Dec 21 '16 at 23:15
  • $\begingroup$ Yes, I think you're right. I basically never use those buttons because their behavior isn't intuitive, and instead adjust the rate by changing the DefaultDuration option and the number of frames. $\endgroup$ – Jens Dec 21 '16 at 23:35
  • $\begingroup$ didn't know about DefaultDuration. Thanks again! $\endgroup$ – freddy90 Dec 21 '16 at 23:38

The solution is just fine, the problem here is Animate computing too few frames.

Monitor[frames = 
   Table[ParametricPlot[{R*Sin[phi[x]], -R*Cos[phi[x]]}, {x, 0, tmax},
      Epilog -> {Red, PointSize[Medium], 
       Point[{R*Sin[phi[t]], -R*Cos[phi[t]]}]}, PlotRange -> All, 
     PerformanceGoal -> "Quality"], {t, 0, tmax, 0.1}], t];

 Export["test.gif", frames]



both work.

enter image description here

You might also try Animate with the du specified in the iterator.

Note we can speed the whole thing up by not regenerating the parametric plot for each frame, now Animate works pretty well.

p = ParametricPlot[{R*Sin[phi[x]], -R*Cos[phi[x]]}, {x, 0, tmax}, 
  PlotRange -> All, PerformanceGoal -> "Quality"]
   Graphics[{Red, PointSize[Medium], 
     Point[{R*Sin[phi[t]], -R*Cos[phi[t]]}]}]}], {t, 0, tmax, 0.1}, 
 AnimationRate -> 1]
| improve this answer | |

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.