1
$\begingroup$

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

enter image description here

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

It seems the event is identified at the first red circle point. However, after that point, the event will be ignored.

Electric field Distribution. Red line reprensents Ey and blue line reprensents Ex.

enter image description here

enter image description here

Improve this question
$\endgroup$
8
  • $\begingroup$ Can you, please, post a MWE? Thanks! $\endgroup$
    – CA Trevillian
    Apr 3, 2021 at 7:21
  • $\begingroup$ @Yuuu There are several typos in your code, for example, which should be Which. It looks like you try to define piecewise function without Piecewise. $\endgroup$
    – Alex Trounev
    Apr 3, 2021 at 9:16
  • $\begingroup$ @ CA Trevillian, sorry, what's MWE? Do you mean megawatt electrical? I will post the eletric field distribution at a certian height. @Alex Trounev, Thanks for pointing out that error. I corrected it, but Ndsovle still seems to ignore the WhenEvent. $\endgroup$
    – Yue Yu
    Apr 3, 2021 at 10:09
  • $\begingroup$ @Yuuu What kind of electric field do you try to calculate? $\endgroup$
    – Alex Trounev
    Apr 3, 2021 at 10:18
  • $\begingroup$ @ Alex Trounev, I'd like calculate electric field of parallel arrayed electrodes. I precalculate four electrodes and use it as a cycle domain. The electric field distribution is added in the question, black line in the X axis are electrodes. Hope it can make the problem clearer. $\endgroup$
    – Yue Yu
    Apr 3, 2021 at 10:22

1 Answer 1

Reset to default
5
$\begingroup$

After several typos correction in Ex, Ey definition we run code up to final result at t=0.1 with Method -> "ExplicitEuler"

\[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}];
Ex1[x_?NumericQ, y_, t_] := 
  Which[x >= 6 a + 11 dx/2, Ex[x - 4 a - 4 dx, y, t], 
   x < 2 a + 3 dx/2, 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}];
Ey1[x_?NumericQ, y_, t_] := 
  Which[x >= 6 a + 11 dx/2, Ey[x - 4 a - 4 dx, y, t], 
   x < 2 a + 3 dx/2, Ex[x + 4 a + 4 dx, y, t], 
   2 a + 3 dx/2 <= x < 6 a + 11 dx/2, Ey[x, y, t]];

It looks like WhenEvent works good in this evaluation, but particle moving in region x<0

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 Ex1[x[t], y[t], t], 
   m y''[t] + 6 \[Pi] \[Eta] R y'[t] == q Ey1[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.1}, StartingStepSize -> 10^-4, Method -> "ExplicitEuler"
  ];

ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 0.1}, 
 Frame -> True, PlotRange -> All, AspectRatio -> 1/2]

Figure 1

If particle bouncing not only from the bottom, but also from the top, then we can use combined rule as follows

WhenEvent[{y[t] == 0, y[t] == 101 b}, y'[t] -> -0.35 y'[t]]

Figure 2

Modified code for long time computation

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) (Tanh[100 (y[t]/R - 1)] + 1)/2, x[0] == 0, y[0] == 100 b, 
   x'[0] == y'[0] == 0, 
   WhenEvent[y[t] == R && Abs[y'[t]] > R g, y'[t] -> -0.35 y'[t]], 
   WhenEvent[y[t] == R/2, y'[t] -> 0]}, {x, y}, {t, 0, 1}, 
  StartingStepSize -> 10^-4, 
  Method -> 
   "ExplicitEuler"](*Method\[Rule]{"DoubleStep", \
Method\[Rule]"ExplicitEuler"}*)


ParametricPlot[Evaluate[{x[t], y[t]} /. s[[1]]], {t, 0, 1}, 
 Frame -> True, PlotRange -> All, AspectRatio -> 1/2]

Figure 3

Improve this answer
$\endgroup$
23
  • $\begingroup$ Thank you very much for going through the complete code. It works well now. I am wondering the reason Ndsolve ignores the WhenEvent is the use of Method -> {"StiffnessSwitching",}. Particle move direction is effected by electric field distribution. Under some parameters, it does travels in X<0. $\endgroup$
    – Yue Yu
    Apr 3, 2021 at 11:57
  • $\begingroup$ In addition, when I set the solution time into 0.1 s.s = NDSolve[{m x''[t] + 6 \[Pi] \[Eta] R x'[t] == q Ex1[x[t], y[t], t], m y''[t] + 6 \[Pi] \[Eta] R y'[t] == q Ey1[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.1}, Method -> {"DiscontinuityProcessing" -> False}]; It will notes that NDSolve::ndsz: At t == 0.06778142044938072`, step size is effectively zero; singularity or stiff system suspected.. Do you know how I can solve this problem. Thank you very much. $\endgroup$
    – Yue Yu
    Apr 3, 2021 at 12:04
  • $\begingroup$ That is why I add Method -> {"StiffnessSwitching",} $\endgroup$
    – Yue Yu
    Apr 3, 2021 at 12:05
  • $\begingroup$ @Yuuu See update to my answer for integration time up to t=0.1. $\endgroup$
    – Alex Trounev
    Apr 3, 2021 at 16:33
  • $\begingroup$ Cool and thanks. A different method can make a great change to the solution. I also tried several different methods for several days, none of which worked. I still know too little about NDsolve.Thank you very much! $\endgroup$
    – Yue Yu
    Apr 3, 2021 at 17:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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