12
$\begingroup$

I'm trying to simulate a particle in an electric and magnetic fields, but numerically instead of analytically. This is basically solving the equation

$$q \cdot \left(p'\times B\right) + q\cdot E = m p'',$$

where $p(t)$ is the position in $(x,y,z)$ coordinates.

After viewing a few topics on this site, I've got a good idea on how to get the solution using NDSolve, but my program gets stuck, and doesn't come up with anything.

b = {1, 0, 0};
e = {0, 0, 1};
q = 1;
m = 1;

sol = NDSolve[ {q*e + q*Cross[D[pos[t], t], b] == m D[pos[t], {t, 2}],
     pos[0] == {0, 0, 0}, (D[pos[t], t] /. t -> 0) == {0, 0, 0}}, 
   pos, {t, 0, 1}];
ParametricPlot3D[Evaluate[pos[t] /. sol], {t, 0, 1}];

It is also worth mentioning that if you remove the $q\cdot E$ term, the calculation is finished, but nothing shows up in the plot.

$\endgroup$
  • $\begingroup$ I haven't worked with vectors in Mathematica so I may not be able to help with that. However, your ParametricPlot3D has incorrect syntax. Assuming that pos is the position, is your initial acceleration zero? $\endgroup$ – dearN Feb 7 '16 at 18:49
13
$\begingroup$

The main problem is that your pos is not seen as a 3D vector.

The cross product is therefore interpreted as a scalar:

q*Cross[D[pos[t], t], b]

Mathematica graphics

when adding this to the vector q.e this 'scalar' term is added to each of the vector components:

q*e + q*Cross[D[pos[t], t], b]

Mathematica graphics

This won't work, instead do:

b = {1, 0, 0};
e = {0, 0, 1};
q = 1;
m = 1;

Define pos as a 3D vector. Also take more time than a single second:

ClearAll[pos]
pos[t_] = {px[t], py[t], pz[t]};
sol = NDSolve[
  {
   q*e + q*Cross[D[pos[t], t], b] == m D[pos[t], {t, 2}],
   pos[0] == {0, 0, 0},
   (D[pos[t], t] /. t -> 0) == {0, 0, 0}
   }, pos[t], {t, 0, 20}]

ParametricPlot3D[Evaluate[pos[t] /. sol], {t, 0, 20}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ what is the difference with this. In the latter answer, the author did not have to make m a 3D vector. In the current problem, is it the addition of q*E that leads to two matrices that may not be "added" together? $\endgroup$ – dearN Feb 7 '16 at 19:03
  • $\begingroup$ @drN Good point. I elaborated the answer somewhat. $\endgroup$ – Sjoerd C. de Vries Feb 7 '16 at 19:27
8
$\begingroup$

Alternative method:

b = {1, 0, 0};
e = {0, 0, 1};
q = 1;
m = 1;

sol = NDSolve[{e + Cross[pos'[t], b] == m/q pos''[t], 
    pos[0] == {0, 0, 0}, pos'[0] == {0, 0, 0}}, pos, {t, 0, 10}, 
   Method -> {"EquationSimplification" -> "Residual"}];
ParametricPlot3D[pos[t] /. sol, {t, 0, 10}, PlotRange -> All]
$\endgroup$
  • $\begingroup$ EquationSimplification->"Residual" is magic? $\endgroup$ – dearN Feb 7 '16 at 18:57
  • 2
    $\begingroup$ @drN It's a bit of a mystery, isn't it, since the LHS of the pos''[t] evaluates to 3x3 matrix and the RHS is a 3-vector. But it works, I guess: Look at the diagonal of the matrix in Inactivate[{e + Cross[pos'[t], b] == m/q pos''[t], pos[0] == {0, 0, 0}, pos'[0] == {0, 0, 0}}, Equal] /. sol /. t -> 1. It probably minimizes the norm of the residual (the matrix minus the vector), which in this case would lead to the correct solution, but if so, you'd think the large norm would cause an error/warning. $\endgroup$ – Michael E2 Feb 7 '16 at 20:56
6
$\begingroup$

Another alternative is to package your constant vector parameters as DiscreteVariables. In the OP's case, it necessary only to chnage e since b occurs inside Cross, which will not evaluate until all its arguments are vectors. Note that in the equation we changed e to e[t] and set its value with e[0] == {0, 0, 1}.

b = {1, 0, 0};
(*e={0,0,1};*)
q = 1;
m = 1;

sol = NDSolve[{q*e[t] + q*Cross[D[pos[t], t], b] == m D[pos[t], {t, 2}], 
    pos[0] == {0, 0, 0}, (D[pos[t], t] /. t -> 0) == {0, 0, 0}, 
    e[0] == {0, 0, 1}}, pos, {t, 0, 10}, DiscreteVariables -> {e}];

ParametricPlot3D[pos[t] /. sol, {t, 0, 10}, PlotRange -> All]

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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