# NDSolve cannot do a discontinuity replacement for event surfaces that depend only on time

I'm trying to solve a system of DAEs with discrete conditions:

muk = 0.2;
mus = 0.3;
m1 = 1;
m2 = 2;
m3 = 3;
Fn12 = 3;
Fn23 = 2;
minvel = 0.1;
F1[x_] := 2*Sin[5*x];
F2[x_] := 2*Sin[7*x];
Fs12[x_] := (m2*F1[x] - m1*F2[x])/(m1 + m2);
Off[NDSolve::pdord]
s = NDSolve[{m1*x1''[t] == F1[t] - Ff12[t], x1' == 0, x1 == 0,
m2*x2''[t] == F2[t] + Ff12[t], x2' == 0, x2 == 0,
Ff12[t] ==
If[Abs[x1''[t] - x2''[t]] < 0.1 && Abs[Fs12[t]] < mus*Fn12 ,
Fs12[t], muk*Fn12*Sign[x1''[t] - x2''[t]]]}, {x1, x2, Ff12}, {t,
10}]


However when solving I get two weird errors

NDSolve::tddisc

NDSolve cannot do a discontinuity replacement for event surfaces that depend only on time

NDSolve::nderr

Error test failure at t==...; unable to continue which searching on google does not offer that much of the results!

I have found some relevant posts and documentations (one, two, three). It seems that I have to use a helper function or something like that using Piecewise and/or WhenEvent functions but I can't get my head around it!

I would appreciate if you could explain why this error happens and how I can modify my code to avoid it?

• It evaluates without any errors after setting mus=10 and adding the following options to NDSolve: Method -> {Automatic, "DiscontinuityProcessing" -> False}, MaxStepSize -> 10^-2. Check if this produces desired results.
– K.J.
Nov 14 '17 at 12:54
• @K.J. mus can't be 10, it is the static friction coefficient :) I will try the other idea now and will let you know. thanks.
Nov 14 '17 at 13:37
• @K.J. nope, not helping :(
Nov 14 '17 at 13:39

By trial and error, integration as far as t == 1.3773405394574811 can be achieved with

s = NDSolveValue[{m1*x1''[t] == F1[t] - Ff12[t], x1' == 0, x1 == 0,
m2*x2''[t] == F2[t] + Ff12[t], x2' == 0, x2 == 0,
Ff12[t] == If[Abs[x1''[t] - x2''[t]] < .1 && Abs[Fs12[t]] < mus*Fn12,
Fs12[t], muk*Fn12*Sign[x1''[t] - x2''[t]]]}, {x1, x2, Ff12}, {t, 10},
MaxStepSize -> 2 10^-4];
Show[Plot[Fs12[t], {t, 0, 1.377}, PlotStyle -> Directive[Black, Dashed]],
Plot[Evaluate@Through[s[t]], {t, 0, 1.377}, PlotRange -> All],  ImageSize -> Large,
AxesLabel -> {t, "x1, x2, Ff12, Fs12"}, LabelStyle -> Directive[Bold, Black, Medium]] although with four messages.

NDSolveValue::pdord: Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations.

NDSolve::tddisc: NDSolve cannot do a discontinuity replacement for event surfaces that depend only on time. (twice)

NDSolveValue::reinitfail: Unable to reinitialize the system at t = 1.3773405394574811 within specified tolerances.

Comparison between the Ff12 curve (green) and Fs12 curve (black dashed) suggests that Ff12 can be approximated by

Ff12[t] == If[Abs[Fs12[t]] < mus*Fn12, Fs12[t], muk*Fn12*Sign[Fs12[t]]]


which allows integration to proceed as far as t == 10.

ss = NDSolveValue[{m1*x1''[t] == F1[t] - Ff12[t], x1' == 0, x1 == 0,
m2*x2''[t] == F2[t] + Ff12[t], x2' == 0, x2 == 0,
Ff12[t] ==  If[Abs[Fs12[t]] < mus*Fn12, Fs12[t], muk*Fn12*Sign[Fs12[t]]]},
{x1, x2, Ff12}, {t, 10}, MaxStepSize -> 10^-3]
`

A plot over the same domain as in the first plot suggests that this approximation is fairly good. The plot over the entire domain is This approach is by no means perfect but nonetheless produces reasonably credible results.