# Using NDSolve to find particle trajectory

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.

• 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? Commented Feb 7, 2016 at 18:49

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]


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]


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


• 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? Commented Feb 7, 2016 at 19:03
• @drN Good point. I elaborated the answer somewhat. Commented Feb 7, 2016 at 19:27

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]

• EquationSimplification->"Residual" is magic? Commented Feb 7, 2016 at 18:57
• @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. Commented Feb 7, 2016 at 20:56
• This method no longer works at least since version 12.3. Commented 2 days ago

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]


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


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