TransWikia.com

Piecewise Numerical results used in NDsolve

Mathematica Asked by Yuuu on April 5, 2021

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

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

One Answer

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

Answered by Yuuu on April 5, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP