Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I found the behavoiur of Sign function weird in the code below.

when $T=10$

T = 10;
A = 1;
Fk0 = 2;
m1 = 1;
x0[t_] = A Sin[2 Pi t];
Fk[t_] = -Fk0  Sign[x1'[t]];
sol = NDSolve[{
  m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0}, 
 {x0, x1}, {t, 0, T}, MaxSteps -> Infinity]
Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True]

but when $T=40$ , I got

which is how it is supposed to be. I wonder what I did wrong here.

share|improve this question
    
no errors/warnings? –  george2079 Feb 2 at 13:53
    
It returns no error. –  Dadan Ari Wibowo Feb 2 at 13:57

3 Answers 3

Using WhenEvent[] to get a typical kinetic friction behavior:

A = 1;
m1 = 1;
x0[t_] := A Sin[2 Pi t];
Fk0 = m1;
Fk[t_] := -Fk0 a[t];
sol[T_?NumericQ] := 
 NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0, a[0] == 1,
         WhenEvent[x1'[t] == 0, a[t] -> Sign[x1'[t]]]}, {a, x1}, {t, 0, T}, 
         MaxSteps -> Infinity, DiscreteVariables -> a, MaxStepSize -> 10^-3]
GraphicsRow@{
       Plot[{x0''[t], x1[t] /. sol[10], x1'[t] /. sol[10]}, {t, 0, 10}, Frame -> True,
             Evaluated -> True,  PlotStyle -> {Red, Blue, Green}],
       Plot[{x0''[t], x1[t] /. sol[40], x1'[t] /. sol[40]}, {t, 0, 40}, 
            Frame -> True, Evaluated -> True, PlotStyle -> {Red, Blue, Green}]}

Mathematica graphics

share|improve this answer
    
It could be interesting to add a static friction component at those instants when x1'[t]==0 –  belisarius Feb 2 at 16:55
    
does WhenEvent detect a sign change, or only trigger if the numerical value is precisely zero? –  george2079 Feb 2 at 20:02
    
@george2079 You can adjust the detection method by using options, but NDSolve being a numerical method implementation can't rely in exact values –  belisarius Feb 2 at 20:04
    
@george2079 I believe it detects a sign change. –  Daniel Lichtblau Feb 2 at 21:43
    
@george2079: it is documented to detect sign changes: "f==0: the function f of the solution variables crosses zero" (from ref/WhenEvent). –  Albert Retey Feb 3 at 9:34

As it turns out the only thing wrong with the original solution is a stroke of bad luck.. Once you hit the zero point at 5, the default time step is one and the solution marches out values at 6,7,8,9,10.. where the forcing function is zero..

Simply limiting the step size fixes things..

 T = 10;
 A = 1;
 Fk0 = 2;
 m1 = 1;
 x0[t_] = A Sin[2 Pi t];
 Fk[t_] = -Fk0 Sign[x1'[t]];
 sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, 
                x1[0] == 0}, {x0, x1}, {t, 0, T}, MaxSteps -> Infinity, 
                MaxStepSize -> .01]
 Plot[{x0''[t], x1'[t], 10 Sign[x1'[t]]} /. sol, {t, 0, T}, 
          Evaluated -> True, Frame -> True]

enter image description here

Here is the original showing the actual solution points..

Show[{Plot[x1'[t] /. sol, {t, 0, T}, Evaluated -> True, Frame -> True],
    Graphics[{PointSize[.01], 
    Point@{#, x1'[t] /. First@sol /. t -> # } & /@ 
         First@Cases[sol, f_InterpolatingFunction, Infinity][[1, 3]]}]}]

enter image description here

By the way there is a cleaner way to display the solution points using EvaluationMonitor described in the documentation.

share|improve this answer
    
In v8.0.4, MaxStepSize -> .01 doesn't work, but MaxStepSize -> .03 does. I think that a more strict control for the step is needed, though I don't know how to achieve that… –  xzczd Feb 4 at 10:30
    
yes 9.0.1. I find it frankly a bit disturbing that the default time step algorithm has made such a poor choice here abruptly jumping by a factor of ~1000. There should be a way to limit the change in time step instead of just setting an absolute limit but I can't readily see a way. –  george2079 Feb 4 at 19:20
    
I found it, see the 3rd approach in my answer. –  xzczd Feb 5 at 8:28

Dadan Ari Wibowo, george2079 and belisarius are all in v9, right? With this question, I noticed some slight changes are made for NDSolve in v9.

The syntax of NDSolve in OP's code is incorrect for v8.0.4. x0, as a known function, isn't allowed to be the 2nd argument of NDSolve.

Then the behavior of NDSolve at $t=5$ is different:

T = 6;
A = 1;
Fk0 = 2;
m1 = 1;
x0[t_] = A Sin[2 Pi t];
Fk[t_] = -Fk0 Sign[x1'[t]];
sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0}, x1, {t, 0, T}]
Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True]

enter image description here

Notice that the option MaxSteps -> Infinity has been removed or it'll eat up all my computer memory.

Finally, I've found 3 solutions for this problem. The 1st one is to use option WorkingPrecision -> 16:

sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0}, 
              x1, {t, 0, T}, WorkingPrecision -> 16]
Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True]

enter image description here

I have to mention that only WorkingPrecision -> 16 will work in this case i.e. even WorkingPrecision -> 15 or WorkingPrecision -> 17 will generate the same error. This solution is probably a coincidence caused by numerical error so it's not that instructive.

The 2nd approach is to use Method -> "StiffnessSwitching":

sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0}, 
              x1, {t, 0, T}, Method -> "StiffnessSwitching"]
Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True]
(* The plot is same as the former so I'd like not to "repost" it. *)

I found this approach mainly by… er… intuition so I can't give a theoretical explanation for it. (From the failing plot we see that something seems to happen at zero point so I guess it might be a stiffness issue and it seems to be! ) Maybe someone more experienced in NDSolve can give a detailed instruction?

Inspired by this post, I've found the 3rd solution:

sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0}, x1, {t, 0, T}, 
              Method -> {"FixedStep", Method -> "ExplicitRungeKutta"}, 
              StartingStepSize -> 1, MaxStepFraction -> Infinity]
(* Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True] *)
<< "DifferentialEquations`NDSolveUtilities`"
StepDataPlot[sol]

enter image description here

Here I deliberately chose 1 as the step size so all of the zero point should have been hit, but correct result is still generated, so it's hard (at least for me) to tell what actually happened inside NDSolve, after all, it's one of the most complex functions in Mathematica.

Not sure if my 3 approaches will work for v9.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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