# I'd like to display field lines for a point charge in 3 dimensions

I'd like to display field lines for a point charge in 3 dimensions. Not a force field (short arrows) but continuous field lines that start on the charge.

-
VectorPlot3D may be? reference.wolfram.com/mathematica/ref/VectorPlot3D.html – Nasser Jan 25 '12 at 21:54
I think he's looking for a 3D version of StreamPlot. – David Jan 25 '12 at 22:07
@David I was just coming to the same conclusion. That's a good question. (+1) – Mr.Wizard Jan 25 '12 at 22:09
@David seems so. – acl Jan 25 '12 at 22:19
Michael Trott gives some code for visualizing field lines of charges in his book. I can transcribe it if needed. – J. M. Jan 25 '12 at 22:44
show 1 more comment

This is something I have used for my classes. Over time, I've tried to make it more and more user friendly, but that's also made it a little longish. I'll post the complete set of functions, with apologies if it's a bit unwieldy...

As you'll see, I found it does indeed work better in my use cases if I normalize the field, so that we advance along the field lines in more balanced steps. The hardest part in applying these functions is to choose the appropriate seed points.

fieldSolve::usage =
"fieldSolve[f,x,x0,\!$$\*SubscriptBox[\(t$$, $$max$$]\)] \
symbolically takes a vector field f with respect to the vector \
variable x, and then finds a vector curve r[t] starting at the point \
x0 satisfying the equation dr/dt=\[Alpha] f[r[t]] for \
t=0...\!$$\*SubscriptBox[\(t$$, $$max$$]\). Here \[Alpha]=1/|f[r[t]]| \
for normalization. To get verbose output add debug=True to the \
parameter list.";

fieldSolve[field_, varlist_, xi0_, tmax_, debug_: False] := Module[
{xiVec, equationSet, t},
If[Length[varlist] != Length[xi0],
Print["Number of variables must equal number of initial conditions\
\nUSAGE:\n" <> fieldSolve::usage]; Abort[]];
xiVec = Through[varlist[t]];
(* Below, Simplify[equationSet] would cost extra time
and doesn't help with the numerical solution, so   don't try to simplify. *)

equationSet = Join[
Map[D[#, t] &, xiVec] ==
],
(xiVec /. t -> 0) == xi0
]
];
If[debug,
Print[Row[{"Numerically solving the system of equations\n\n",
TraditionalForm[(Simplify[equationSet] /. t -> "t") //
TableForm]}]]];
(* This is where the differential equation is solved.
The Quiet[] command suppresses warning messages because numerical precision isn't crucial for our plotting purposes: *)

Quiet[NDSolve[
equationSet,
xiVec,
{t, 0, tmax}
]]], 2]
]

fieldLinePlot::usage = "fieldLinePlot[field,varlist,seedList] plots \
3D field lines of a vector field (first argument) that depends on the \
symbolic variables in varlist. The starting points for these \
variables are provided in seedList. Each element of \
seedList={{\!$$\*SubscriptBox[\(p$$, $$1$$]\), \!$$\*SubscriptBox[\(T\$$, $$1$$]\)},{\!$$\*SubscriptBox[\(p$$, $$2$$]\), \
\!$$\*SubscriptBox[\(T$$, $$2$$]\)}...} is a tuple where \
\!$$\*SubscriptBox[\(p$$, $$i$$]\) is the starting point of the i-th \
field line and \!$$\*SubscriptBox[\(T$$, $$i$$]\) is the length of \
that field line starting from \!$$\*SubscriptBox[\(p$$, $$i\$$]\).";

fieldLinePlot[field_, varList_, seedList_, opts : OptionsPattern[]] :=
Module[{sols, localVars, var, localField},
If[Length[seedList[[1, 1]]] != Length[varList],
Print["Number of variables must equal number of initial conditions\
\nUSAGE:\n" <> fieldLinePlot::usage]; Abort[]];
localVars = Array[var, Length[varList]];
localField =
ReleaseHold[
Hold[field] /.
(* Assume that each element of seedList specifies a point AND the \
length of the field line: *)
Show[ParallelTable[
(
ParametricPlot3D[
Evaluate[Through[#[t]]], {t, #[[1, 1, 1, 1]], #[[1, 1, 1, 2]]},
Evaluate@
Apply[Sequence,
FilterRules[{opts}, Options[ParametricPlot3D]]]] &
)@
fieldSolve[localField, localVars, seedList[[i, 1]],
seedList[[i, 2]]]
, {i, Length[seedList]}]]
];

Options[fieldLinePlot] = Options[ParametricPlot3D];

SyntaxInformation[fieldLinePlot] = {"LocalVariables" -> {"Solve", {2, 2}},
"ArgumentsPattern" -> {_, _, _, OptionsPattern[]}};

SetAttributes[fieldSolve, HoldAll];


The main function is fieldLinePlot, but I split it into two functions to be more modular. Also, the problem of where to start drawing the field lines is treated separately because it depends a lot on the particular application.

fieldSolve[f,x,x0,Subscript[t, max]] symbolically takes a vector field f with respect to the vector variable x, and then finds a vector curve r[t] starting at the point x0 satisfying the equation dr/dt = α f[r[t]] for t=0...tmax. Here α = 1/|f[r[t]]| for normalization. To get verbose output add debug=True to the parameter list.

fieldLinePlot[field,varlist,seedList] plots 3D field lines of a vector field (first argument) that depends on the symbolic variables in varlist. The starting points for these variables are provided in seedList.

Each element of seedList={{p1, T1},{p2, T2}...} is a tuple where pi is the starting point of the $i^\mathrm{th}$ field line and Ti is the length of that field line in both directions from Pi.

Here are some examples:

### 1) Coulomb field of two opposite charges at $\vec{r} = \vec{0}$ and $\vec{r} = (1, 1, 1)$:

Look at the form of seedList to see how the field line starting points and lengths are specified.

seedList =
With[{vertices = .1 N[PolyhedronData["Icosahedron"][[1, 1]]]},
Join[Map[{#, 2} &, vertices],
Map[{# + {1, 1, 1}, -2} &, vertices]]];

Show[fieldLinePlot[{x, y, z}/
Norm[{x, y, z}]^3 - ({x, y, z} - {1, 1, 1})/
Norm[{x, y, z} - {1, 1, 1}]^3, {x, y, z}, seedList,
PlotStyle -> {Orange, Specularity[White, 16], Tube[.01]},
PlotRange -> All, Boxed -> False, Axes -> None],
Background -> Black]


### 2) Magnetic field of an infinite straight wire:

With[{seedList = Table[{{x, 0, 0}, 6.5}, {x, .1, 1, .1}]
},
Show[fieldLinePlot[{-y, x, 0}/(x^2 + y^2), {x, y, z},
seedList, PlotStyle -> {Orange, Specularity[White, 16], Tube[.01]},
PlotRange -> All, Boxed -> False, Axes -> None],
Graphics3D@Tube[{{0, 0, -.5}, {0, 0, .5}}], Background -> Black]]


-
 Welcome to Mathematica.SE and thank you for posting an answer! You can surround code or variable names with backticks (as in  f[x] ) to format them more readably, and you can also use LaTeX (surround it with $ signs: $\int x^2 \, dx$ renders as$\int x^2 \, dx$) – Szabolcs Jan 26 '12 at 18:22 Thanks for taking the time to tweak the formatting. I just tried it and it's nice to see the real-time update in the preview... – Jens Jan 26 '12 at 18:35 Yes, the real time preview is really nice. Unfortunately if the post is long, it may be a bit slow. If this happens to you, you can try right clicking a (rendered) math expression, and choose Settings -> Math Renderer -> MathML. (This will only work in browsers that support MathML, such as Firefox.) The formatting quality will be a little different (some would say worse), but the performance will be much better too. I have an old and slow machine, so I always use MathML. We also have one-click graphics uploading now :-) – Szabolcs Jan 26 '12 at 18:39 Assuming that the motion of the particle is governed by some force field, you could use NDSolve together with ParametricPlot3D to plot the individual field lines. For example, consider the force field force[p_] := ({1, 1, 1} - p)/Norm[{1, 1, 1} - p]^3 - p/Norm[p]^3  The equations of motion are given by p[t_] := {p1[t], p2[t], p3[t]} eqs := Thread[p''[t] == force[p[t]]]  We also need some initial conditions, for example ic := {Thread[p[0] == {1, 0, 0}], Thread[p'[0] == {0, 1, 0}]}  The system can then be solved according to sol = NDSolve[{eqs, ic}, {p1, p2, p3}, {t, 0, 20}]  And the plot of this solution looks like ParametricPlot3D[p[t] /. sol, {t, 0, 20}]  If you want several field lines, you would need to rerun NDSolve for a list of initial conditions. Edit As rcollyer pointed out, this doesn't actually plot the field lines. For that you would need to solve p'[t] == force[p[t]]/Norm[force[p[t]] for which you can still use the method above, e.g. eq := Thread[p'[t] == -force[p[t]]/Sqrt[force[p[t]].force[p[t]]]] (* seeds for field lines *) seeds = Flatten[{0.1 #, 0.1 # + {1, 1, 1}} & /@ N[PolyhedronData["Icosahedron"][[1, 1]]], 1]; sol = NDSolve[{eq, Thread[p[0] == #]}, {p1, p2, p3}, {t, 0, 40}][[1]] & /@ seeds; ParametricPlot3D[p[s] /. sol, {s, 0, 20}, PlotRange -> {{-3, 3}, {-3, 3}, {-3, 3}}]  -  It's cool anyway (mind you I just stopped myself from posting a very similar answer, realizing my mistake just before clicking "Post your answer"...) – acl Jan 25 '12 at 22:56 You can get the coordinates of the polyhedron easier (and more readable) by using the built-in parameter VertexCoordinates (... /@ PolyhedronData["Icosahedron", "VertexCoordinates"] ...). – David Jan 26 '12 at 4:29 One more thing, the field lines are the integral curves of the vector field, i.e. solutions of$p'(t) = X(p(t))\$, there is no division by a norm. (This equation holds on any differentiable manifold, and most of them don't even have a norm.) – David Jan 26 '12 at 4:43

The closest thing in the help files that I can see is this example. It might be of some use.

stackPlots[plots2D : {__Graphics},
dh_, o : OptionsPattern[Graphics3D]] := Graphics3D[
MapIndexed[
Function[{g, index}, g[[1]] /. Arrow[v2d_] :> Arrow[v2d /.
{x_?NumericQ, y_?NumericQ} :> {x, y, dh index[[1]]}]],
plots2D], o]

plotSpacing = 5;

values = {-1, 0, 1};
plots = MapIndexed[
Function[{\[Lambda], i},
StreamPlot[{y, \[Lambda] - x^2}, {x, -3, 3}, {y, -3, 3},
StreamPoints -> 16, StreamScale -> 0.07,
StreamStyle -> ColorData["SolarColors"][0.3 i[[1]]]]], values]

stackPlots[plots, plotSpacing, Axes -> True, Boxed -> False,
Ticks -> {Automatic, Automatic,
MapIndexed[{plotSpacing #2[[1]], Row[{"\[Lambda] = ", #1}]} &,
values]}]


-