# WhenEvent - energy conservation not valid [closed]

I would like to simulate the following conservative system, with purely elastic impacts (coefficient of restitution is 1). Applying the WhenEvents the energy of the system is not conserved. How can I correct this?

T = 10;
a = 0.1;
b = 0.1;
c = 1
{solu[t_], {data}} =
Reap[NDSolveValue[{a*u1''[t] + b*(u1[t] - u2[t]) == 0,
u2''[t] + u2[t] + b*(u2[t] - u1[t]) == 0, u1[0] == 0,
u1'[0] == 3, u2[0] == 0, u2'[0] == 0,
WhenEvent[Abs[(u2[t] - u1[t])] == c, Sow[{t, u2[t] - u1[t]}]],
WhenEvent[
Abs[(u2[t] - u1[t])] ==
c, {u1'[t] -> ((a - 1)*u1'[t] + 2*u2'[t])/(1 + a),
u2'[t] -> (2*a*u1'[t] + (1 - a)*u2'[t])/(1 + a)} ]},
{u1[t], u2[t], t}, {t, 0, T}, MaxStepSize -> 0.001]];

Plot[a*(solu'[t][[1]])^2 + (solu'[t][[2]])^2 +
b*((solu[t][[1]]) - (solu[t][[2]]))^2 + (solu[t][[2]])^2, {t, 0, T},
PlotRange -> {{0, T}, All}, Frame -> True,
FrameLabel -> {"t", "Energy"}]


## closed as off-topic by xzczd, b.gates.you.know.what, rhermans, MarcoB, bbgodfreySep 29 '18 at 17:23

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – xzczd, b.gates.you.know.what, rhermans, MarcoB, bbgodfrey
If this question can be reworded to fit the rules in the help center, please edit the question.

• This has been mentioned in the Possible Issues of document of WhenEvent: "With sequential event actions, the variables are modified in turn, …, To swap the variable values, use simultaneous events." So the correct way to set the event is: WhenEvent[Abs[u2[t] - u1[t]] == c, {u1'[t], u2'[t]} -> {((a - 1) u1'[t] + 2 u2'[t])/(1 + a), ( 2 a u1'[t] + (1 - a) u2'[t])/(1 + a)}] – xzczd Sep 25 '18 at 9:24

Apparently, we need to show the corrected code and the result

T = 10;
a = 1/10;
b = 1/10;
c = 1;
{U1, U2, U11, U21} =
NDSolveValue[{a*u1''[t] + b*(u1[t] - u2[t]) == 0,
u2''[t] + u2[t] + b*(u2[t] - u1[t]) == 0, u1[0] == 0, u1'[0] == 3,
u2[0] == 0, u2'[0] == 0,
WhenEvent[
Abs[u2[t] - u1[t]] ==
c, {u1'[t],
u2'[t]} -> {((a - 1) u1'[t] + 2 u2'[t])/(1 +
a), (2 a u1'[t] + (1 - a) u2'[t])/(1 + a)}]}, {u1, u2, u1',
u2'}, {t, 0, T}, WorkingPrecision -> 30];
energy[t_] :=
a*(U11[t])^2 + (U21[t])^2 + b*((U1[t]) - (U2[t]))^2 + (U2[t])^2
Plot[energy[t] - energy[0], {t, 0, T}, PlotRange -> {{0, T}, All},
Frame -> True, FrameLabel -> {"t", "Energy"}]


Here, the deviation of energy from the initial value of about $$10^{-11}$$

On the second fig. the solution is shown

Plot[{U1[t], U2[t]}, {t, 0, 10}]