# Solving a Nonlinear Complementary Problem (plasticity)

I would like to solve the following: given $$t\mapsto\sigma(t)$$ and $$E>0$$, $$\sigma_y>0$$, find $$\epsilon$$ such that:

$$\left\lbrace\begin{array}{l}g(t,\epsilon)\geq 0,\\ \phi(t,\epsilon,\epsilon')\geq0, \\g(t,\epsilon)\phi(t,\epsilon,\epsilon')=0,\\ \epsilon(0)=0. \end{array}\right.$$

with $$g(t,\epsilon(t))=\sigma_y - |\sigma(t) - E\epsilon(t)|$$ and $$\phi(t,\epsilon(t),\epsilon'(t)) = (\sigma(t) - E\epsilon(t))\epsilon'(t)$$.

Example:

sigma[t_] := Sin[t];
sigmay = 0.5;
E0 = 1;
g[t_?NumericQ, epsi_] := sigmay - Abs[sigma[t] - E0*epsi]
phi[t_?NumericQ, epsi_, dotepsi_] := (sigma[t] - E0*epsi)*dotepsi
epsisol = NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, epsi[0] == 0}, epsi, {t, 0, 100}]

(* NDSolveValue::ntdvdae: Cannot solve to find an explicit formula for the derivatives.
NDSolve will try solving the system as differential-algebraic equations.

NDSolveValue::ndcf: Repeated convergence test failure at t == 1.5717016618338167; unable to continue. *)


Attempt with WhenEvent:

events = {WhenEvent[g[t, epsi[t]] == 0, coef[t] -> 1], WhenEvent[phi[t, epsi[t], epsi'[t]] == 0, coef[t] -> 0]}

epsisol = First@NDSolveValue[{g[t, epsi[t]]*coef[t] + phi[t, epsi[t], epsi'[t]]*(1 - coef[t]) == 0, epsi[0] == 0,
coef[0] == 0, events}, {epsi, coef}, {t, 0, 10}, DiscreteVariables -> {coef}]

(* NDSolveValue::ntdvdae: Cannot solve to find an explicit formula for the derivatives.
NDSolve will try solving the system as differential-algebraic equations.

... then integration stops at t = 1.54 with no further error *)


Context and details

I am interested in plasticity, in particular the computation of the displacement $$\varepsilon$$ induced by an applied force $$\sigma$$ to a system composed of a slider and a spring in parallel, such as:

This problem is addressed in Solving a discontinuous differential-algebraic equation system for plasticity behaviour ($$C_2$$ is $$H$$ and $$C_1$$ is $$\infty$$) --- I believe there is a slight mistake in the equations but it still works after correction.

However, with both answers, I did not manage to adapt the code so that it would work with two such devices in series. Also, I wanted to derive the equation from "standard" plasticity theory:

• a plasticity criterion $$f =|\sigma -X| - \sigma_y \leq 0$$ ($$g=-f$$ above, and $$X=E \epsilon$$ is the stress in the spring);
• the positiveness of dissipation: $$\phi = \epsilon' (\sigma_y -E\epsilon)\geq 0$$
• an orthogonality condition implying that energy is dissipated iff there is plasticity ($$f=0$$): $$f\times \phi = 0$$.

This is often written altogether: $$0\leq (-f) \perp \phi \geq 0$$

Such kinds of formulation are also found in intermittent contact dynamics: the reaction force is always non-negative, and can be non-zero only if there is contact, i.e. when the distance is zero. Conversely, if the distance is non-negative, the reaction for can only be zero.

Such problems are numerically challenging, even though there are dedicated numerical methods. Even the formulation involving the derivative $$\epsilon'$$ is wobbly, because $$\epsilon$$ is not differentiable everywhere (just like the velocity of a bouncing ball is not defined at impact times).

Anyway, WhenEvent works very fine for bouncing balls with few contacts, so I would have expected WhenEvent to be efficient here.

I think the behavior of WhenEvent encountered by OP is a bug. Anyway, here's a working WhenEvent based solution:

sigma[t_] = Sin[t];
sigmay = 0.5;
E0 = 1;
g[t_, epsi_] = sigmay - Abs[sigma[t] - E0*epsi];
phi[t_, epsi_, dotepsi_] = (sigma[t] - E0 epsi) dotepsi;

events = {WhenEvent[g[t, epsi[t]] < phi[t, epsi[t], epsi'[t]], coef[t] -> 1],
WhenEvent[phi[t, epsi[t], epsi'[t]] < g[t, epsi[t]], coef[t] -> 0]};

epsisol = First@
NDSolveValue[{g[t, epsi[t]] coef[t] + phi[t, epsi[t], epsi'[t]] (1 - coef[t]) == 0,
epsi[0] == 0, coef[0] == 0, events}, {epsi, coef}, {t, 0, 100},
DiscreteVariables -> coef, SolveDelayed -> True]

Plot[epsisol[t], {t, 0, 15}]


• It looks like your solution not limited by $\pm 0.5$ as solution by Steffen Jaeschke. – Alex Trounev Oct 27 '20 at 13:31
• @Alex Seems that the DAE solver has turned to i.c. epsi[0] == -0.038287, this might be the reason, so this method isn't that good… – xzczd Oct 27 '20 at 15:02
• @AlexTrounev I find another method that gives the correct solution, have a look. – xzczd Oct 28 '20 at 3:23
• @AlexTrounev Actually I did manage to make it work with regularization by reformulating the problem with sign (sign[x_] = 2/Pi*ArcTan[10000*x]; epsipsol = NDSolveValue[{E1*epsip[t] - sigma[t] == -sigmay*sign[epsip'[t]], epsip[0] == 0}, epsip, {t, 0, 100}]). But I disagree that it is supposed to be better (see my previous comment). It's supposed to be simpler, for sure, but a high cost! – anderstood Oct 28 '20 at 13:45
• @anderstood You need to increase points for plotting, try e.g. PlotPoints -> 1000 or MaxRecursion -> 15 :) – xzczd Oct 30 '20 at 4:57

A trick to obtain the full result.

sigma[t_] := Sin[t];
sigmay = 0.5;
E0 = 1;
tmax = Pi;
g[t_?NumericQ, epsi_] := sigmay - Abs[sigma[t] - E0*epsi]
phi[t_?NumericQ, epsi_, dotepsi_] := (sigma[t] - E0*epsi)*dotepsi
tmax = Pi;
tmin = 0;
epsisolant = sigma[tmin];
GR = {};

While[tmax < 100,
epsisol = NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, epsi[tmin] == epsisolant}, epsi, {t, tmin, tmax}, Method -> {"EquationSimplification" -> "Residual"}];
AppendTo[GR, Plot[epsisol[t], {t, tmin, tmax}]];
epsisolant = epsisol[tmax];
tmin = tmax;
tmax += Pi/4
]

Show[GR, PlotRange -> All]

• Mmh, that works because for this specific simple sigma the singular points are known. But its not very robust. Still the idea of breaking manually the integration might do the trick. – anderstood Oct 23 '20 at 9:20

You input is to my knowledge applied to it correctly. Well done.

But this is a discretized attempt to solve the problem.

sigma[t_] := Sin[t];
sigmay = 0.5;
E0 = 1;
g[t_?NumericQ, epsi_] := sigmay - Abs[sigma[t] - E0*epsi]
phi[t_?NumericQ, epsi_, dotepsi_] := (sigma[t] - E0*epsi)*dotepsi
epsisol =
NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0,
epsi[0] == 0}, epsi, {t, 10^-13, 100}]


The second message opens a page ndsolve::ndcf with the direct invitation to contact the technical support of Wolfram Inc.

I found that the domain depends with a rapid jump on the starting time at a little bit more than 10^-13 for example a quarter I reproduce Your results and around that less again. Might be this is a match for the domain length 4.71. That can even be got again at higher starting times as 0.0001 or so.

My output is:

Plot[epsisol[t], {t, 0.005, 4.71}, PlotRange -> Full]


From that on I agree with [@cesareo]5 it might go on delayed quasi-periodic. This might already be chaotic not only in the starting time but in the period. The rise and fall may be characteristic. Somehow this is similar to a sawtooth. Therefore and because the switch function suggests it, I make the solution ideation that this might be solved with Fourier or Laplace methodologies for more domain. This will only work in approximation.

But the curious idea changed my plans: make the domain smaller arbitrarly:

epsisol =
NDSolveValue[{Min[g[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0,
epsi[0] == 0}, epsi, {t, 10^-13, 10}]


Plot[epsisol[t], {t, 0.005, 10}, PlotRange -> Full]


Hope that helps. This is done with V12.0.0 on iMac Catalina.

This can be solved up to 10.99639 if the Method -> {"EquationSimplification" -> "Residual"} is used. The message remains: ndcf. The repeated convergence test does not accept the rapid stagnation of the growth of the solution at -0.5. But it suffices for the full period of the graph. Perhaps the treatment as a differential-algebraic equation.

Seems that a better match for sigmay and sigma give a longer domain in the capabilities for of-the-shelf differential-algebraic methods. Maybe this is on the other-hand a question designed for th fail of the adaptibility of the repeated convergence test.

I made a comparison between both solution, mine and from @xzczd.

Plot[{epsisol[t], epsisolu[t],
1.25 10^7 (epsisol[t] - epsisolu[t])}, {t, 0.00001, 10.99},
PlotRange -> Full]


Despite both solutions look at first sight very similiar they are different.

Plot[{epsisol[t] - .5, epsisolu[t] - .5}, {t, 1.8, 3.14},
PlotRange -> Full, PlotLegends -> "Expressions"]


Mine stays a little bit, an order of one magnitude further away from the limiting 0.5. This is even bigger for the negative border and bigger at the second constant interval. Then my solution goes over to fail. The even very small error oscillates up and the end is the test fails.

With InterpolationOrder->All the oscillation get much smaller and more repetitive:

But the domain is not larger.

To each Accuracy 9,10,11,... there is interval near zero for which the integration is successful.

epsisol = NDSolveValue[{Min[gi[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] == 0, epsi[0] == 0}, epsi, {t, 10^-10.1295, 11}, Method -> {"EquationSimplification" -> "Residual"}, InterpolationOrder -> All, AccuracyGoal -> 10]

Plot[{epsisol[t], epsisolu[t],
1.25 10^7 (epsisol[t] - epsisolu[t])}, {t, 0.00001, 7.85},
PlotRange -> Full]


For Accuracy 11 the domain has a much large interval for which my solution gets much closer to the reference solution and the oscillation get tamed. At -0.5 mine is better than that of the competitor. But the oscillations remain still order 10^-7.

Fast and dirty as Mathematica built-ins are these days. The behavior is a clear hint that Mathematica uses StiffnessSwitching internally for the calculation of the solution.

ListLinePlot@
Quiet@Table[(epsisol =
NDSolveValue[{Min[gi[t, epsi[t]], phi[t, epsi[t], epsi'[t]]] ==
0, epsi[0] == 0}, epsi, {t, 10^expon, 11},
Method -> {"EquationSimplification" -> "Residual"},
InterpolationOrder -> All, AccuracyGoal -> 13])[[1, 1,
2]], {expon, -5, -16, -.01}]
`

There are many possible start values for Accuracy 12. The result is still stiffness switching wildly but the exactness increases strongly.

• Fail to reproduce the last result in v11.3 and v12.1.1, which version are you in? – xzczd Oct 27 '20 at 14:43