Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
1  
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
1  
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
    
Hello Rainforest Frog and welcome to the Mathematica StackExchange. Don't forget to upvote good answers (and other people's questions) using the triangle above the number next to the post, and use the checkmark to "accept" the answer to your question that you think best answers it. You now have enough "reputation" (points) to visit the chat room and chat if you would like. –  Verbeia Jan 25 '12 at 22:45

3 Answers 3

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[
    Thread[
     Map[D[#, t] &, xiVec] == 
      Normalize[field /. Thread[varlist -> xiVec]]
     ],
    Thread[
     (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: *)

  Map[Head, First[xiVec /.
     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] /. 
     Thread[Map[HoldPattern, Unevaluated[varList]] -> localVars]];
  (* 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]

enter image description here

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

enter image description here

share|improve this answer

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

Mathematica graphics

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

Mathematica graphics

share|improve this answer
    
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]}]

Mathematica graphics

Mathematica graphics

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.