# Message Boards

0
|
4803 Views
|
0 Replies
|
0 Total Likes
View groups...
Share
 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