Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I'm trying to make a simulation of an electron that goes through a magnetic and electric field perpendicular to each other and to the velocity of the electron but I'm running in some troubles with the Mathematica language. I'm a beginner at Mathematica language and I ended up probably using too much procedural programming and I would like suggestions on how I can improve it using real Mathematica language. For those who saw this question before I edited it, the second problem was fixed and it was a simple logic mistake, so I took it of the question.

vel = {0.05*300000000, 0, 0}; (* velocity *)
pos = {0, 0, 0}; (* position *)
acel = {0, 0, 0}; (* acceleration *)
ce = {0, -5000, 0}; (* electric field *)
cm = {0, 0, -0.0003}; (* magnetic field *)
dt = 0.000000000001;
l = 0.05; (* length of the region where there are magnetic and electric fields *)
me = 9.10938*10^-31; (* electrons mass *)
q = -1.60217657*10^-19; (* elementary charge *)
data = Table[0, {i, 1000000}, {j, 2}];
For[count = 1, pos[[1]] <= 1.5*l, count++,
     If[Abs[pos[[1]]] <= l, acel = q/me*(ce + (vel\[Cross]cm)), acel = {0, 0, 0}]; 
     pos = pos + vel*dt + (acel*dt^2)/2; vel = vel + acel*dt;
     data[[count, 1]] = pos[[1]]; data[[count, 2]] = pos[[2]]
For[r = 100000; s = 0.00001, r > 0.001, r = (r/10); s = (s/10),
     While[Abs[vel[[2]]] > r, cm[[3]] = cm[[3]] - s; 
             vel = {0.05*300000000, 0, 0}; 
             pos = {0, 0, 0}; 
             data = Table[0, {i, 1000000}, {j, 2}];
                 For[count = 1, pos[[1]] <= 1.5*l, count++,
                      If[Abs[pos[[1]]] <= l, acel = q/me*(ce + (vel\[Cross]cm)), acel = {0, 0, 0}]; 
                      pos = pos + vel*dt + (acel*dt^2)/2; vel = vel + acel*dt;
                      data[[count, 1]] = pos[[1]]; data[[count, 2]] = pos[[2]]

I know that there are duplicate variable declarations and redundant statements in the code, its because I pasted everything together to post here, but it was originally apart.

share|improve this question
Please post code, not images of code, and make sure that code evaluates properly in a fresh kernel session. – Yves Klett Oct 13 '13 at 8:00
I can't even read the tiny text in your code image. Please post your code with proper markdown formatting. – m_goldberg Oct 13 '13 at 8:40
So what if I do? Do you expect me, then, to manually retype all your code into Mathematica so I can work with it? No way. – m_goldberg Oct 13 '13 at 8:50
Markdown and formatting help: – Yves Klett Oct 13 '13 at 8:59
Most users will only bite if code is posted here, which is sensible, too. – Yves Klett Oct 13 '13 at 9:33
up vote 2 down vote accepted

I adapted a very primitive procedure I wrote to mimick the computation of the orbit of Saturn as seen on volume 1 of the Feynman's Lectures on Physics. Let's define the constants:

me = 9.10938*10^-31;
q = -1.60217657*10^-19;

The state of the system at a given time is described by the vectors position (r) and mom--- velocity (v). The state of the system a quid dt after that time is given by:

nextStep[dt_][{r0_, v0_}] :=
    Block[{a, v, r},
      a = q/me(Ef + Cross[v0, Bf]);
      v = v0 + a dt;
      r = r0 + v dt;
      {r, v}

(Two caveats, here: This is not the smartest way to compute a, v and r. I also have a version that follows Feynman's advice to make it a but smarter but I do not want to spoil your fun. Also, you'll have to add the condition in the computation of a. Something like a = If[ Abs[r0[[1]]]<1, that 'a' there, else 0].) EDIT: You might also want to impose similar conditions on the fields, in which case you will have to make them a function of position, so that a will use Ef[r0] and Bf[r0].

Let's write a procedure that computes all the steps you want. Since I totally lack fantasy, I will call it "compute".

compute[r0_, v0_, dt_, n_Integer] := NestList[nextStep[dt], {r0, v0}, n]

(NestList and Nest are function you want to learn). Let's compute a few values, shall we?

Ef = {0, -5000, 0}; Bf = {0, 0, -0.0003};
compute[{0, 0, 0}, {0.05*3 10^8, 0, 0}, 0.000000000001, 5]

you get a list in the form


From which you can extract the trajectory. For Feynman's example I had written a procedure to do just that (but I used only scalars, so since it is easy to write one, I won't bother you with that).

Oh, in case you are wondering what was that Feynman suggested to increase the precision of the computation, well... Look it up on the Lectures: Book I, Chapter 9, Psalm 7: Planetary Motion.

share|improve this answer

Your Answer


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.