differential equations – WhenEvent disabled in NDsolve

I try to simulate particle motion in electric field. Ex and Ey is the intensity of electric field in X axis and y axis resepectively, which is related to the position of particle.

``````g = 9.8; q = 1.08 10^-13;
s = NDSolve({m x''(t) + 6 (Pi) (Eta) R x'(t) == q Ex(x(t), y(t), t),
m y''(t) + 6 (Pi) (Eta) R y'(t) == q Ey(x(t), y(t), t) - m g,
x(0) == 0, y(0) == 100 b, x'(0) == y'(0) == 0,
WhenEvent(1 b < y(t) < 5 b, y'(t) -> -0.35 y'(t))}, {x, y}, {t, 0,
0.1}, Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 2,
PrecisionGoal -> 4, Method -> {"DiscontinuityProcessing" -> False},
MaxStepSize -> 10^-6);
StepDataPlot(s)
ParametricPlot(Evaluate({x(t), y(t)} /. s), {t, 0, 0.1})
``````

WhenEvent describes that when particle land on the plane, it will bounce. Initially, I wrote like this
`WhenEvent(y(t)==0, y'(t) -> -0.35 y'(t))}`

I thought Ndsovle might miss the point. So I enlarge the the section of WhenEvent and minimize the Max step size. But it still doesnot work.

I do not know Whether I set the parameters wrong or use the wrong method.

Sometimes, it will notes that Singularity or stiffness problesm, so I add “StiffnessSwitching”. And the function of Ex and Ey is piecewise periodic function, so I add `Method -> {"DiscontinuityProcessing" -> False}`

I am a rookie in using NDsolve. I am still learning it in dealing with particle trajectories. Thank you very much for any constructive ideas.

Here is complete code. Most of the previous ones are for calculating the electric field. If you have time to comb through the code, just run it. Thanks again.

``````Needs("DifferentialEquations`NDSolveProblems`");
Needs("DifferentialEquations`NDSolveUtilities`");
Needs("DifferentialEquations`InterpolatingFunctionAnatomy`");
(Epsilon) = 8.854187817 10^-12;
vc0 = 0.8 10^3;
a = 0.3 10^-3;
b = 18 10^-6;
dx = a;
l = 4 a + 4 dx;
n1 = 40;
n2 = 6;
k1 = 1/6;
k2 = 1/6;
TT = 0.2;
(Delta)1 = a/n1; (Delta)2 = b/n2; s1 = k1 (Delta)1; s2 =
k2 (Delta)2; c = a - 2 s2; d = b - 2 s1; (Delta)3 =
c/n1; (Delta)4 = d/n2; n = 2 (n1 + n2);
(*1*)xa1 =
Join(Table(i (Delta)1, {i, 0, n1}), Table(a, {i, 1, n2 - 1}),
Table(a - i (Delta)1, {i, 0, n1}), Table(0, {i, 1, n2 - 1}));
ya1 = Join(Table(b, {i, 0, n1}),
Table(b - i (Delta)2, {i, 1, n2 - 1}), Table(0, {i, 0, n1}),
Table(i (Delta)2, {i, 1, n2 - 1}));
xa2 = Join(Table(s2 + i (Delta)3, {i, 0, n1}),
Table(a - s2, {i, 1, n2 - 1}),
Table(a - s2 - i (Delta)3, {i, 0, n1}), Table(s2, {i, 1, n2 - 1}));
ya2 = Join(Table(b - s1, {i, 0, n1}),
Table(b - s1 - i (Delta)4, {i, 1, n2 - 1}), Table(s1, {i, 0, n1}),
Table(s1 + i (Delta)4, {i, 1, n2 - 1}));

(*2*)xb1 =
Join(Table(i (Delta)1, {i, 0, n1}), Table(a, {i, 1, n2 - 1}),
Table(a - i (Delta)1, {i, 0, n1}), Table(0, {i, 1, n2 - 1})) +
a + dx;
yb1 = Join(Table(b, {i, 0, n1}),
Table(b - i (Delta)2, {i, 1, n2 - 1}), Table(0, {i, 0, n1}),
Table(i (Delta)2, {i, 1, n2 - 1}));
xb2 = Join(Table(s2 + i (Delta)3, {i, 0, n1}),
Table(a - s2, {i, 1, n2 - 1}),
Table(a - s2 - i (Delta)3, {i, 0, n1}),
Table(s2, {i, 1, n2 - 1})) + a + dx;
yb2 = Join(Table(b - s1, {i, 0, n1}),
Table(b - s1 - i (Delta)4, {i, 1, n2 - 1}), Table(s1, {i, 0, n1}),
Table(s1 + i (Delta)4, {i, 1, n2 - 1}));

(*3*)xc1 =
Join(Table(i (Delta)1, {i, 0, n1}), Table(a, {i, 1, n2 - 1}),
Table(a - i (Delta)1, {i, 0, n1}), Table(0, {i, 1, n2 - 1})) +
2 (a + dx);
yc1 = Join(Table(b, {i, 0, n1}),
Table(b - i (Delta)2, {i, 1, n2 - 1}), Table(0, {i, 0, n1}),
Table(i (Delta)2, {i, 1, n2 - 1}));
xc2 = Join(Table(s2 + i (Delta)3, {i, 0, n1}),
Table(a - s2, {i, 1, n2 - 1}),
Table(a - s2 - i (Delta)3, {i, 0, n1}),
Table(s2, {i, 1, n2 - 1})) + 2 (a + dx);
yc2 = Join(Table(b - s1, {i, 0, n1}),
Table(b - s1 - i (Delta)4, {i, 1, n2 - 1}), Table(s1, {i, 0, n1}),
Table(s1 + i (Delta)4, {i, 1, n2 - 1}));

(*4*)xd1 =
Join(Table(i (Delta)1, {i, 0, n1}), Table(a, {i, 1, n2 - 1}),
Table(a - i (Delta)1, {i, 0, n1}), Table(0, {i, 1, n2 - 1})) +
3 (a + dx);
yd1 = Join(Table(b, {i, 0, n1}),
Table(b - i (Delta)2, {i, 1, n2 - 1}), Table(0, {i, 0, n1}),
Table(i (Delta)2, {i, 1, n2 - 1}));
xd2 = Join(Table(s2 + i (Delta)3, {i, 0, n1}),
Table(a - s2, {i, 1, n2 - 1}),
Table(a - s2 - i (Delta)3, {i, 0, n1}),
Table(s2, {i, 1, n2 - 1})) + 3 (a + dx);
yd2 = Join(Table(b - s1, {i, 0, n1}),
Table(b - s1 - i (Delta)4, {i, 1, n2 - 1}), Table(s1, {i, 0, n1}),
Table(s1 + i (Delta)4, {i, 1, n2 - 1}));

(*5*)xe1 =
Join(Table(i (Delta)1, {i, 0, n1}), Table(a, {i, 1, n2 - 1}),
Table(a - i (Delta)1, {i, 0, n1}), Table(0, {i, 1, n2 - 1})) +
4 (a + dx);
ye1 = Join(Table(b, {i, 0, n1}),
Table(b - i (Delta)2, {i, 1, n2 - 1}), Table(0, {i, 0, n1}),
Table(i (Delta)2, {i, 1, n2 - 1}));
xe2 = Join(Table(s2 + i (Delta)3, {i, 0, n1}),
Table(a - s2, {i, 1, n2 - 1}),
Table(a - s2 - i (Delta)3, {i, 0, n1}),
Table(s2, {i, 1, n2 - 1})) + 4 (a + dx);
ye2 = Join(Table(b - s1, {i, 0, n1}),
Table(b - s1 - i (Delta)4, {i, 1, n2 - 1}), Table(s1, {i, 0, n1}),
Table(s1 + i (Delta)4, {i, 1, n2 - 1}));

(*6*)xf1 =
Join(Table(i (Delta)1, {i, 0, n1}), Table(a, {i, 1, n2 - 1}),
Table(a - i (Delta)1, {i, 0, n1}), Table(0, {i, 1, n2 - 1})) +
5 (a + dx);
yf1 = Join(Table(b, {i, 0, n1}),
Table(b - i (Delta)2, {i, 1, n2 - 1}), Table(0, {i, 0, n1}),
Table(i (Delta)2, {i, 1, n2 - 1}));
xf2 = Join(Table(s2 + i (Delta)3, {i, 0, n1}),
Table(a - s2, {i, 1, n2 - 1}),
Table(a - s2 - i (Delta)3, {i, 0, n1}),
Table(s2, {i, 1, n2 - 1})) + 5 (a + dx);
yf2 = Join(Table(b - s1, {i, 0, n1}),
Table(b - s1 - i (Delta)4, {i, 1, n2 - 1}), Table(s1, {i, 0, n1}),
Table(s1 + i (Delta)4, {i, 1, n2 - 1}));

(*7*)xg1 =
Join(Table(i (Delta)1, {i, 0, n1}), Table(a, {i, 1, n2 - 1}),
Table(a - i (Delta)1, {i, 0, n1}), Table(0, {i, 1, n2 - 1})) +
6 (a + dx);
yg1 = Join(Table(b, {i, 0, n1}),
Table(b - i (Delta)2, {i, 1, n2 - 1}), Table(0, {i, 0, n1}),
Table(i (Delta)2, {i, 1, n2 - 1}));
xg2 = Join(Table(s2 + i (Delta)3, {i, 0, n1}),
Table(a - s2, {i, 1, n2 - 1}),
Table(a - s2 - i (Delta)3, {i, 0, n1}),
Table(s2, {i, 1, n2 - 1})) + 6 (a + dx);
yg2 = Join(Table(b - s1, {i, 0, n1}),
Table(b - s1 - i (Delta)4, {i, 1, n2 - 1}), Table(s1, {i, 0, n1}),
Table(s1 + i (Delta)4, {i, 1, n2 - 1}));

(*8*)xh1 =
Join(Table(i (Delta)1, {i, 0, n1}), Table(a, {i, 1, n2 - 1}),
Table(a - i (Delta)1, {i, 0, n1}), Table(0, {i, 1, n2 - 1})) +
7 (a + dx);
yh1 = Join(Table(b, {i, 0, n1}),
Table(b - i (Delta)2, {i, 1, n2 - 1}), Table(0, {i, 0, n1}),
Table(i (Delta)2, {i, 1, n2 - 1}));
xh2 = Join(Table(s2 + i (Delta)3, {i, 0, n1}),
Table(a - s2, {i, 1, n2 - 1}),
Table(a - s2 - i (Delta)3, {i, 0, n1}),
Table(s2, {i, 1, n2 - 1})) + 7 (a + dx);
yh2 = Join(Table(b - s1, {i, 0, n1}),
Table(b - s1 - i (Delta)4, {i, 1, n2 - 1}), Table(s1, {i, 0, n1}),
Table(s1 + i (Delta)4, {i, 1, n2 - 1}));

x1 = Join(xa1, xb1, xc1, xd1, xe1, xf1, xg1, xh1);
y1 = Join(ya1, yb1, yc1, yd1, ye1, yf1, yg1, yh1);
x2 = Join(xa2, xb2, xc2, xd2, xe2, xf2, xg2, xh2);
y2 = Join(ya2, yb2, yc2, yd2, ye2, yf2, yg2, yh2);
p0(i_, k_) := -(1/(4 (Pi) (Epsilon))) Log((x1((i)) -
x2((k)))^2 + (y1((i)) - y2((k)))^2);
p = Table(p0(i, k), {i, 1, 8 n}, {k, 1, 8 n});
vc1 = Join(Table(vc0, n), Table(vc0, n), Table(-vc0, n),
Table(-vc0, n), Table(vc0, n), Table(vc0, n), Table(-vc0, n),
Table(-vc0, n));
vc2 = Join(Table(-vc0, n), Table(vc0, n), Table(vc0, n),
Table(-vc0, n), Table(-vc0, n), Table(vc0, n), Table(vc0, n),
Table(-vc0, n));
vc3 = Join(Table(-vc0, n), Table(-vc0, n), Table(vc0, n),
Table(vc0, n), Table(-vc0, n), Table(-vc0, n), Table(vc0, n),
Table(vc0, n));
vc4 = Join(Table(vc0, n), Table(vc0, n), Table(-vc0, n),
Table(vc0, n), Table(vc0, n), Table(-vc0, n), Table(-vc0, n),
Table(vc0, n));(*Voltage of electrodes*)
(Lambda)1 = LinearSolve(p, vc1, Method -> "Multifrontal");
(Lambda)2 = LinearSolve(p, vc2, Method -> "Multifrontal");
(Lambda)3 = LinearSolve(p, vc3, Method -> "Multifrontal");
(Lambda)4 = LinearSolve(p, vc4, Method -> "Multifrontal");
ClearAll((Lambda));
(Lambda)(t_?NumericQ) :=
Which(t >= TT, (Lambda)(t - TT), t < 0, (Lambda)(t + TT),
0 <= t < TT/4, (Lambda)1, TT/4 <= t < TT/2, (Lambda)2,
TT/2 <= t < 3 TT/4, (Lambda)3, 3 TT/4 <= t < TT, (Lambda)4);
ex0(x_, y_, t_, k_) :=
Indexed((Lambda)(t), k)/(2 (Pi) (Epsilon)) (
x - x2((k)))/((x - x2((k)))^2 + (y - y2((k)))^2);
ey0(x_, y_, t_, k_) :=
Indexed((Lambda)(t), k)/(2 (Pi) (Epsilon)) (
y - y2((k)))/((x - x2((k)))^2 + (y - y2((k)))^2);
Ex(x_, y_, t_) := Sum(ex0(x, y, t, k), {k, 1, 8 n});
Ex(x_?NumericQ, y_, t_) :=
which(x >= 6 a + 11 dx/2, Ex(x, y, t) = Ex(x - 4 a - 4 dx, y, t),
x < 2 a + 3 dx/2, Ex(x, y, t) = Ex(x + 4 a + 4 dx, y, t),
2 a + 3 dx/2 <= x < 6 a + 11 dx/2, Ex(x, y, t));
Ey(x_, y_, t_) := Sum(ey0(x, y, t, k), {k, 1, 8 n});
Ey(x_?NumericQ, y_, t_) :=
which(x >= 6 a + 11 dx/2, Ey(x, y, t) = Ey(x - 4 a - 4 dx, y, t),
x < 2 a + 3 dx/2, Ex(x, y, t) = Ex(x + 4 a + 4 dx, y, t),
2 a + 3 dx/2 <= x < 6 a + 11 dx/2, Ey(x, y, t));
m = 1.8 10^-9; (Eta) = 1.8 10^-5; R = 5 10^-5;
g = 9.8; q = 1.08 10^-13;
s = NDSolve({m x''(t) + 6 (Pi) (Eta) R x'(t) == q Ex(x(t), y(t), t),
m y''(t) + 6 (Pi) (Eta) R y'(t) == q Ey(x(t), y(t), t) - m g,
x(0) == 0, y(0) == 100 b, x'(0) == y'(0) == 0,
WhenEvent(y(t) == 0, y'(t) -> -0.35 y'(t))}, {x, y}, {t, 0, 0.02},
AccuracyGoal -> 2, PrecisionGoal -> 4,
Method -> {"DiscontinuityProcessing" -> False});
ParametricPlot(Evaluate({x(t), y(t)} /. s), {t, 0, 0.02})
``````