Problem with NDSolve and Sign

Question

If c1 cnd c2 are zero, the code responds very well; but if one of c1 or c2 is nonzero, the code doesn't work correctly. I get the error

NDSolve::ndinid: Initial condition {0} is not in the range specified by the discrete variable NDSolve`s$147246. >>

I think it's because of the Sign function. Without the Sign function, the code works properly. I would like it to work correctly, and also I would like to know why this code didn't work properly.

ClearAll[θ1, θ2];
m1 = 1;  m2 = 1; lc1 = 1;  lc2 = 1; l1 = 2; l2 = 2; Ix2 = 1; Iy2 = 1;  Iz1 = 1;
Iz2 = 1;  g = 9.81;  Im1 = 1; Im2 = 1; η1 = 28;  η2 = 28; b1 = 1; b2 = 2;
c1 = 0; c2 = 1;

τ1[t_] := 10;
τ2[t_] := 10;

u1[t_] := τ1[t] - b1*θ1'[t] - c1*Sign[θ1'[t]];
u2[t_] := τ2[t] - b2*θ2'[t] - c2*Sign[θ2'[t]];


m11[t_] := m1*(lc1)^2 + m2*(l1 + lc2*Cos[θ2[t]])^2 + 
           Ix2*(Sin[θ2[t]])^2 + Iy2*(Cos[θ2[t]])^2 + Iz1;
m22[t_] := m2*(lc2)^2 + Iz2;
V11[t_] := 0;
V12[t_] := (-2*m2*lc2*(l1 + lc2*Cos[θ2[t]])* Sin[θ2[t]] + 
           2 (Ix2 - Iy2)*Sin[θ2[t]]*Cos[θ2[t]])*θ1'[t];
V21[t_] := (m2*lc2*(l1 + lc2*Cos[θ2[t]])* Sin[θ2[t]] - 
           (Ix2 - Iy2)*Sin[θ2[t]]* Cos[θ2[t]])*θ1'[t];
V22[t_] := 0;
M = {{m11[t], 0}, {0, m22[t]}} + {{Im1*(η1)^2, 0}, {0, Im2*(η2)^2}};
V = {{V11[t], V12[t]}, {V21[t], V22[t]}};
G = {{0}, {m2*lc2*Cos[θ2[t]]*g}};
initu = q[0] == {{0}, {.1}};
initv = q'[0] == {{0}, {0}};
q[t_] := {θ1[t], θ2[t]};
sol = First@  NDSolve[{Thread[M.q''[t] + V.q'[t] + G == {u1[t], u2[t]}], initu, 
     initv}, q[t], {t, 0, 100}];

Plot[{θ1[t] /. sol, θ2[t] /. sol}, {t, 0, 100}, PlotTheme -> "Detailed"]

Answer 1

It's a good idea to include the error you got in your question. It gives a clue and might prompt someone to investigate:

NDSolve::ndinid: Initial condition {0} is not in the range specified by the discrete variable NDSolve`s$147246. >>

Now Sign is discontinuous and will cause NDSolve to invoke special processing of the ODE. I suspect that the strange internal variable NDSolve`s$147246 is added to the system along with an event to implement the change in Sign. Note that the message is helpful in calling this variable a discrete variable, which is just what one would need to implement the discontinuous changes in Sign. (Normally one wouldn't expect to see an internal local variable, and I'm still undecided whether you really should. It would have been nicer if the basename were not the single letter s but something like Sign, say.)

In any case, the problem is with the initial condition for q, that is to say, for {θ1, θ2}. It defines θ1 and θ2 to be 1D vectors instead of scalars, so that Sign[θ1'[t]] will be a vector. But apparently NDSolve thought Sign was going to be a scalar. The mismatch causes the error.

If it is acceptable for θ1 and θ2 to be scalars, then here is a fix:

initu = q[0] == Flatten@{{0}, {.1}};  (* <--- N.B. *)
initv = q'[0] == Flatten@{{0}, {0}}; (* <--- N.B. *)
sol = First@NDSolve[{Thread[M.q''[t] + V.q'[t] + G == {u1[t], u2[t]}], 
    initu, initv}, q[t], {t, 0, 100}];

Plot[{θ1[t] /. sol, θ2[t] /. sol}, {t, 0, 100}, PlotTheme -> "Detailed"]

Note: It is not essential, but what I often do, especially when debugging code on this site, is to examine the pieces. The error message complains about the initial conditions. Since I didn't write the code, one of the first things I'm curious about in NDSolve problems is what the ODE system actually looks like. So I evaluate it:

 {Thread[M.q''[t] + V.q'[t] + G == {u1[t], u2[t]}], initu, initv}

The first thing that struck me were those strange (to me) initial conditions for θ1 and θ2. I often do this sort of examination when I've written the code for my own system and it does not work as expected.