I am attempting to change the sign of a constant when a certain condition is met during a numerical integration. Here is the code:

μk = 0.3; g = 9.81; R = 1; ti = 0; m = 0.3;

EOM =  -μk (R θ'[t]^2 + g Sin[θ[t]]) + 
    g Cos[θ[t]] == R θ''[t];
IC1 =  θ'[0] == 0;  IC2 =  θ[0] == 0;

sol = NDSolve[{EOM, IC1, IC2, 
    WhenEvent[θ'[t] == 0, {μk = -μk; tf = t; 
    "RestartIntegration"}]}, θ[t], {t, ti, ∞}];

θ = θ[t] /. sol;
θd = D[θ[t] /. sol, t];
Plot[{θ, θd}, {t, ti, tf}]

The problem lies in the fact that it is not dissipating as it should. θ should start to approach Pi/2. Instead it oscillates on to Infinity.

Note: I discovered that I cannot simply use Sign[θ'[t]] since the value is crossing zero. A zero value is undesirable in my case.

Slider System

Here is an image of the system being modeled.

  • $\begingroup$ If I understand your last sentence correctly, you do not wish μk == 0 when θ'[t] == 0, What value do you wish it to have? $\endgroup$ – bbgodfrey Sep 19 '18 at 14:13
  • $\begingroup$ It should simply have a value that is +μk or -μk. $\endgroup$ – MGoforth Sep 19 '18 at 18:02
  • $\begingroup$ Yes, but which? The issue is, what value should μk have, if a time-step lands exactly on θ'[t] == 0, which appears to be happening at t == 2.062520756856687? You will notice that μk flips sign there, but θ'[t] does not. I can change this behavior, so that uk does not flip signs there, in which case the answer is entirely different. $\endgroup$ – bbgodfrey Sep 19 '18 at 18:56
  • $\begingroup$ Ah I understand you question now. This is modeling a slider that oscillates in a vertical plane on a semicircular guide. Thus, when θ'[t] == 0 the friction force is changing its direction. Long story short, I am using when θ'[t] == 0 as a way to change the direction of the friction force. Moving forward with my model I need to add a check that indicates if the direction of friction 'should' change as well as check static friction to see if the slider slips after the slider comes to rest. I added a photo of the system to the original post. $\endgroup$ – MGoforth Sep 19 '18 at 19:54
  • $\begingroup$ In that case, why do you not wish to use Sign[θ'[t]], which gives no friction when the slider is stationary? Also, why is friction applied to the potential energy, g Sin[θ[t]]? $\endgroup$ – bbgodfrey Sep 19 '18 at 20:15

The code has three issues. First, μk cannot be changed during the computation unless it is designated by DiscreteVariables. Second, the list of actions to be taken by WhenEvent needs to be separated by commas, not semicolons. Third, the upper limit of integration cannot be infinity. With these changes, the code becomes

g = 9.81; R = 1; ti = 0; m = 0.3;

EOM =  -μk [t] (R θ'[t]^2 + g Sin[θ[t]]) + g Cos[θ[t]] == R θ''[t];
IC1 =  θ'[0] == 0;  IC2 =  θ[0] == 0;

sol = NDSolve[{EOM, IC1, IC2, μk [0] == 0.3, 
    WhenEvent[θ'[t] == 0, {μk[t] -> -μk[t], tf = t, "RestartIntegration"}]}, 
    {θ[t], θ'[t], μk[t]}, {t, ti, 10}, DiscreteVariables -> μk];

tmax = Flatten[θ[t] /. sol /. t -> "Domain"] // Last;
Plot[Evaluate[{θ[t], θ'[t], μk[t]} /. sol], {t, ti, tmax}, ImageSize -> Large]

enter image description here

Note that the computation terminates with a NDSolve::ndsz error. Basically, the computation has become unstable. Method -> "StiffnessSwitching" does not help. Incidentally, "RestartIntegration" has no substantive effect in this computation and can be omitted.

Improved Solution

Based on additional information provided by the OP in comments below and in the question itself, a better approach is

g = 981/100; R = 1; ti = 0; m = 3/10; μ = 3/10;

EOM =  -μ Sign[θ'[t]] (R θ'[t]^2 + g Sin[θ[t]]) + g Cos[θ[t]] == R θ''[t];
IC1 =  θ'[0] == 0;  IC2 =  θ[0] == 0;

sol = NDSolve[{EOM, IC1, IC2}, {θ[t], θ'[t]}, {t, ti, 10}, WorkingPrecision -> 30];

Plot[Evaluate[{θ[t], θ'[t], μ  Sign[θ'[t]]} /. sol], {t, ti, tmax}, ImageSize -> Large]

enter image description here

The two solutions are identical until θ'[t] == 0 at t == 2.062520756856687159564857914. In the earlier computation, WhenEvent causes μk to reverse sign, even though θ'[t] does not, thereby introducing negative friction. I have tried numerous WhenEvent options, as well as WorkingPrecision as high as 60 to resolve this issue, but without success. In conclusion, only the second solution is credible physically.

  • 3
    $\begingroup$ +1. Two remarks: (1) The upper limit of integration can be Infinity, but one cannot actually integrate to infinity. Usually one adds an event that will "StopIntegration", if the stopping time is unknown. (In this case, your solution seems the right approach.) (2) My guess, based on $y'=\pm y^2$, is that the solution blows up (singularity) once θ'[t]^2 grows too large, rather than numerical instability. $\endgroup$ – Michael E2 Sep 19 '18 at 10:36

It is possible to study the dependence of the model on the parameters, but even for a small value of the parameter $\mu k$ there is no convergence of the solution to the limiting value of $\pi/2$, since this is the point of unstable equilibrium. Below is a model and data showing the transition through a point of unstable equilibrium.

g = 981/100; R = 1; ti = 0; tf = 20;
f[x1_?NumberQ, x_?NumberQ] := -\[Mu]k*(R x1^2 + g Sin[x])
sol = Block[{\[Mu]k = 1/20}, 
   NDSolve[{R \[Theta]''[t] == 
      f[\[Theta]'[t], \[Theta][t]] + g Cos[\[Theta][t]], \[Theta]'[
       0] == 0, \[Theta][0] == 0, 
     WhenEvent[\[Theta]'[t] == 0, \[Mu]k = -\[Mu]k; 
      "RestartIntegration"]}, {\[Theta], \[Theta]'}, {t, ti, tf}, 
    WorkingPrecision -> 30]];

{Plot[Evaluate[\[Theta][t] /. sol], {t, ti, tf}], 
 Plot[Evaluate[\[Theta]'[t] /. sol], {t, ti, tf}]}



In case anyone is curious, here is the final product of the problem on which I was working. I found joy in the solution, so here you go:

enter image description here

muk = 0.3; mus = 0.4; g = 9.81; R = 1; ti = 0; thi = 0; thf = 0; sgn = 1; 
n=10^3; cnt = 0;
thList = {}; thdList = {}; TimeList = {};
th1 = ArcTan[1/mus];
th2 = Pi - th1;
 cnt = cnt + 1;
 EOM =  \[Theta]''[t] + sgn muk \[Theta]'[t]^2 + 
    g/R (sgn muk Sin[\[Theta][t]] - Cos[\[Theta][t]]) == 0;
 IC1 =  \[Theta]'[ti] == 0;  IC2 =  \[Theta][ti] == thi;
 sol = NDSolve[{EOM, IC1, IC2, 
    WhenEvent[\[Theta]'[t] == 0, tf = t; "StopIntegration"]}, \[Theta][t], 
     ti, \[Infinity]}];
 \[Theta] = \[Theta][t] /. sol;
 thd = D[\[Theta][t] /. sol, t];
 a = Flatten[Table[\[Theta], {t, ti, tf, (tf - ti)/n}]];
 b = Flatten[Table[thd, {t, ti, tf, (tf - ti)/n}]];
 time = Flatten[Table[t, {t, ti, tf, (tf - ti)/n}]];
 AppendTo[thList, a];
 AppendTo[thdList, b];
 AppendTo[TimeList, time];
 Clear[\[Theta], thd, a, b, time];
 ti = tf;
 thi = (\[Theta][t] /. sol /. t -> tf)[[1]];
 If[thi > Pi/2, sgn = -1, sgn = 1];
 If[th1 < thi < th2, Break[]];
 If[cnt == 20, {Break[], Print["count limit reached"]}];]
thList = Flatten[thList];
thdList = Flatten[thdList];
TimeList = Flatten[TimeList];
thData = Thread[{TimeList, thList}];
thdData = Thread[{TimeList, thdList}];
SetOptions[ListPlot, BaseStyle -> {FontFamily -> "Cambria", FontSize -> 12}, 
  GridLines -> Automatic, Frame -> True];
ListPlot[{thData, thdData}, Joined -> True, 
 PlotLegends -> {"\[Theta][t]  (rad)", "thd[t]  (rad/s)"}, 
 FrameLabel -> {"Time  (s)", ""}, PlotLabel -> "Slider Phenomenon", 
 PlotStyle -> {Lighter[Red, 0.3], Lighter[Blue, 0.3]}]
Print["The System Comes to Rest:  {{ thf = ", 
 thf = Last[thList], " rad   OR   thf = ", thf*180/Pi, "\[Degree] }}"]

You can change the values of "muk" and "mus" (Static and Kinetic Friction Coefficients) to observe different behavior of the slider.

This specifically is NOT an explicit answer to the question. I used the answer by "bbgodfrey" to come to the solution to this problem.

Also, this probably can be achieved in a much simpler way, which could be explored... Thanks.


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.