Hi everyone,
I read an interesting tuto on mathematica. Pretty interesting for electrostatic: On the importance of being edgy in electrostatic
I'm trying to find the trajectories of charged particles with random initial conditions in the vicinity of a uniformly charged disk with sharp edge.
I defined the electrostatic field, then the gradient. After that I want to solve numerically the equations of the motion for several random conditions.
I define 1st the potential as a sum of elliptic functions:
Q = 1; \[Epsilon] = 1; R = 1; r[x_, y_, z_] := Sqrt[(x)^2 + y^2];
p[x_, y_, z_] := Sqrt[(r[x, y, z] + R)^2 + (z)^2];
n[x_, y_, z_] := (4 r[x, y, z] R)/(r[x, y, z] + R)^2;
m[x_, y_, z_] := (4 r[x, y, z] R)/p[x, y, z]^2;
Ud[x_, y_, z_] :=
2/p[x, y,
z] (p[x, y, z]^2*EllipticE[m[x, y, z]] - (r[x, y, z]^2 - R^2)*
EllipticK[m[x, y, z]] - (z)^2*(r[x, y, z] - R)/(r[x, y, z] + R)
EllipticPi[n[x, y, z], m[x, y, z]] );
Ui[x_, y_, z_] :=
Piecewise[{{-2 \[Pi]*Sqrt[(z)*(z)],
r[x, y, z] \[LessSlantEqual] R}, {0, r[x, y, z] > R}}]; \[Phi][x_,
y_, z_] :=
1/(4 \[Pi] \[Epsilon])*Q/(\[Pi] R^2) (Ud[x, y, z] + Ui[x, y, z]);
Then I define the gradient:
Linedisk[{x_, y_,
z_}] := {-D[\[Phi]t[x, y, z], x], -D[\[Phi][x, y, z],
y], -D[\[Phi][x, y, z], z]};
After that I want to solve the equation of motion for several random conditions (here 5):
SeedRandom[111]; Show[{(*the charged disk *)
Graphics3D[{Directive[GrayLevel[0.3], Specularity[Yellow, 25]],
Cylinder[{{0, 0, 0.01}, {0, 0, -0.01}}, 1]}],(*sample orbits *)
Table[nds =
NDSolve[Join[
Thread[{x''[t], y''[t],
z''[t]} == -Linedisk[{x[t], y[t], z[t]}]],
Thread[{x[0], y[0], z[0]} == 2 RandomReal[{-0.5, 0.5}/2, {3}]],
Thread[{x'[0], y'[0], z'[0]} ==
RandomReal[{-0.5, 0.5}, {3}]]], {x, y, z}, {t, 0, 5},
MaxSteps -> 50000];
ParametricPlot3D[
Evaluate[{x[t], y[t], z[t]} /. nds[[1]]], {t, 0.1,
nds[[1, 1, 2, 1, 1, 2]]},
PlotStyle -> {ColorData["RedBlueTones"][RandomReal[]],
Tube[0.012]}, PlotRange -> All], {5}]}, PlotRange -> All]
After that I have many alarms like: Power::infy: "Infinite expression 1/0. encountered." or Infinity::indet: "Indeterminate expression ComplexInfinity+ComplexInfinity encountered."
At the end, I have a 3D graph but I don't trust the results because of the alarms.
There is a singularity at r=1 but normally it should be ok.
I guess the problems are from the elliptic functions that are maybe too difficult and heavy for NDsolve.
Any help will be appreciated to better understand the alarms from Mathematica.
Thanks