# NDsolve alarms

Posted 8 years ago
4252 Views
|
0 Replies
|
0 Total Likes
|
 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
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments