13
$\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.

Improve this question
$\endgroup$
1
  • $\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
    Commented Feb 7, 2016 at 18:49

4 Answers 4

Reset to default
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

Improve this answer
$\endgroup$
2
  • $\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
    Commented Feb 7, 2016 at 19:03
  • $\begingroup$ @drN Good point. I elaborated the answer somewhat. $\endgroup$
    – Sjoerd C. de Vries
    Commented Feb 7, 2016 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]
Improve this answer
$\endgroup$
3
  • $\begingroup$ EquationSimplification->"Residual" is magic? $\endgroup$
    – dearN
    Commented Feb 7, 2016 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
    Commented Feb 7, 2016 at 20:56
  • 1
    $\begingroup$ This method no longer works at least since version 12.3. $\endgroup$
    – xzczd ♦
    Commented Jul 31 at 12:06
7
$\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

Improve this answer
$\endgroup$
1
$\begingroup$

At least since version 12.3, OP's code gives the desired output even if the equation has actually been ruined by Listable attribute:

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

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, tend}];
ParametricPlot3D[Evaluate[pos[t] /. sol], {t, 0, tend}]

enter image description here

And the equation is no longer ruined since version 14.1, because the Cross function now owns a NonThreadable attribute:

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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