3
$\begingroup$

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? enter image description here

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"}]

The energy decrease after the impacts

$\endgroup$

closed as off-topic by xzczd, b.gatessucks, rhermans, MarcoB, bbgodfrey Sep 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.gatessucks, rhermans, MarcoB, bbgodfrey
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 6
    $\begingroup$ 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)}] $\endgroup$ – xzczd Sep 25 '18 at 9:24
2
$\begingroup$

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}$ fig1

On the second fig. the solution is shown

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

fig2

$\endgroup$

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