# Piecewise Numerical results used in NDsolve

I have a solution of LinearSolve. And there are four conditions in single period. Then I used Which to expand this solution.

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"];
\[Lambda][t_] :=
Which[t >= T, \[Lambda][t - T], t < 0, \[Lambda][t + T],
0 <= t < T/4, \[Lambda]1, T/4 <= t < T/2, \[Lambda]2,
T/2 <= t < 3 T/4, \[Lambda]3, 3 T/4 <= t < T, \[Lambda]4];


I can use this to calculate other functions and plot them.

ex0[x_, y_, t_, k_] := \[Lambda][t][[k]]/(2 \[Pi] \[Epsilon]) (
x - x2[[k]])/((x - x2[[k]])^2 + (y - y2[[k]])^2);
ey0[x_, y_, t_, k_] := \[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}];
Ey[x_, y_, t_] := Sum[ey0[x, y, t, k], {k, 1, 8 n}];
Plot[Ex[2 a, 2 b, t], {t, 0, 2}]


But when I use this function in NDsovlve,

m = 1.8 10^-9; \[Eta] = 1.8 10^-5; R = 5 10^-5;
g = 9.8; q = -5.4 10^-14;
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] == -m g + q Ey[x[t], y[t], t],
x[0] == 0, y[0] == 5 b, x'[0] == y'[0] == 0,
WhenEvent[y[t] == 0, y'[t] -> -0.95 y'[t]]}, {x, y}, {t, 0.1},
AccuracyGoal -> 10, Method -> {"DiscontinuityProcessing" -> False}];
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 0.1},
AspectRatio -> 1/2]


it always notes that

Part::partw: Part 13 of Which[t>=0.01,\[Lambda][t-T],t<0,\[Lambda][t+T],0<=t<T/4,\[Lambda]1,T/4<=t<T/2,\[Lambda]2,T/2<=t<(3 T)/4,\[Lambda]3,(3 T)/4<=t<T,\[Lambda]4] does not exist.

How can I solve this? Thank you very much!

Complete code

\[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;
T = 0.01;
T1 = 0.1;
\[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 >= T, \[Lambda][t - T], t < 0, \[Lambda][t + T],
0 <= t < T/4, \[Lambda]1, T/4 <= t < T/2, \[Lambda]2,
T/2 <= t < 3 T/4, \[Lambda]3, 3 T/4 <= t < T, \[Lambda]4];
ex0[x_, y_, t_, k_] := \[Lambda][t][[k]]/(2 \[Pi] \[Epsilon]) (
x - x2[[k]])/((x - x2[[k]])^2 + (y - y2[[k]])^2);
ey0[x_, y_, t_, k_] := \[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}];
Ey[x_, y_, t_] := Sum[ey0[x, y, t, k], {k, 1, 8 n}];
m = 1.8 10^-9; \[Eta] = 1.8 10^-5; R = 5 10^-5;
g = 9.8; q = -5.4 10^-14;
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] == -m g + q Ey[x[t], y[t], t],
x[0] == 0, y[0] == 0, x'[0] == y'[0] == 0,
WhenEvent[y[t] == 0, y'[t] -> -0.95 y'[t]]}, {x, y}, {t, 0, 0.1},
AccuracyGoal -> 10, Method -> {"DiscontinuityProcessing" -> False}];
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 0.1},
AspectRatio -> 1/2]

• What happens if you change the def. of lambda to ClearAll[\[Lambda]]; \[Lambda][t_?NumericQ] :=...<same Which[] statement>,,, ? – Michael E2 Mar 30 at 23:26
• vc0, n, p, T , b are not defined, so your code can not be tested. – Daniel Huber Mar 31 at 6:57
• I defined it before. It is a huge table, so I didnot coyp it to here. – Yuuu Mar 31 at 8:55
• ClearAll[\[Lambda]]; \[Lambda][t_?NumericQ] := Which[t >= T, \[Lambda][t - T], t < 0, \[Lambda][t + T], 0 <= t < T/4, \[Lambda]1, T/4 <= t < T/2, \[Lambda]2, T/2 <= t < 3 T/4, \[Lambda]3, 3 T/4 <= t < T, \[Lambda]4];Still notes that: Part::partw: Part 2 of [Lambda][t] does not exist. – Yuuu Mar 31 at 8:57
• Until someone has time to comb through your code, look up Indexed and try replacing the call(s) to Part that cause the error with Indexed. – Michael E2 Apr 1 at 1:39

As @ Michael E2 suggested, turn the call of Lambda from the form of \[Lambda][t][[k]] to Indexed[\[Lambda][t], k]