Problem plotting the trajectory of electron in static electromagnetic field

I. Content :

When the electron moves in an uniform electromagnetic field, it will be acted by Lorentz force:

We can thus determine the acceleration of the electron. If we know the initial position and velocity, we can determine the kinetic motion equations of electron . Eliminating t from mentioned motion equations, we can derive , which is the orbital equation of electron.

II.Task :

Write Matlab program to:

1. Input the initial position (), velocity of the electron (); constant static electric () and magnetic field () parallel to z-axis.
2. Use symbolic operations to calculate the electromagnetic force acting on the electron, thereby deriving the electron's acceleration, velocity, and motion equation.
3. Plot the trajectory of the electron

IIT.My take on this :

EXPECTED OUTPUT:

MY OUTPUT:

MY CODE:

function pjoject6_alpha
clc
clear

syms t vx(t) vy(t) vz(t) x(t) y(t) z(t) C

%[INPUT] FOR INITIAL POSITION: x0 y0 z0
xyz = [0 0 0];

x0= xyz(1); 
y0= xyz(2); 
z0= xyz(3);

%[INPUT] FOR INITIAL VELOCITY VECTOR
v0 = [1.75 1 1];

v0x= v0(1); 
v0y= v0(2); 
v0z= v0(3);

%[INPUT] FOR
% E : ELECTRIC FIELD VECTOR
% B : MAGNETIC FIEL VECTOR 
E = [0 0 1.5];
B = [0 0 1];

%[PREDEFINED CONSTANT]
% me : MASS OF ELECTRON
% q :  PARTICLE CHARGE OF ELECTRON
me= 9.10938*10^-31; q= -1.60276*10^-19;

%[CACULATING] ACCELERATION VECTOR
% USING F = me.accel = q.(E) + q.(v.B)
accel = (E + cross(v0, B))*q/me; % cross(v0,B) : multiply vector v0 with vector B 


%[SOLVE] FOR VELOCITY EQUATION vx(t), vy(t), vz(t)
vx  = dsolve([diff(vx(t),t) == accel(1),vx(0) == v0x]);
vy  = dsolve([diff(vy(t),t) == accel(2),vy(0) == v0y]);


vz  = dsolve([diff(vz(t),t) == accel(3),vz(0) == v0z]);


x= dsolve([diff(x) == vx, x(0) == x0 ]);
y= dsolve([diff(y) == vy, y(0) == y0 ]);
z= dsolve([diff(z) == vx, z(0) == z0 ]);


if vz== 0
    fplot(x, y); title('Trajectory ');
else
    fplot3(x, y, z); title('Trajectory ');
end
end

MY THOUGHT PROCESS:

Followed the instruction, I reckon that I need to find the Vector of acceleration (accel).

After that I dsolve them to get the equation for velocity , respectively

I do the same for another time and get .

The last step is to fplot3(x,y,z) (The plot function will eliminate t and return an equation for the relation of

MY SPECULATION REGARDING THE ERROR:

I think there was something wrong with the plotting.

When I tried to eliminate t manually using eliminate():

equa = [x;y;z]
sol = eliminate(equa,t)
disp(sol)

output:

sol =
 
Empty sym: 1-by-0

And I think this is why i got a straight line

Please help me, I'm a complete novice

Answer 1

[Bjorn Gustavsson]

When you let the symbolic solver handle your equations of motion it might be safer to simply give dsolve the second-order equation straight-off and see if that handles it all in one go. Also it might be very good if you sit down an solve these equations by hand first. They separate nicely into two equations perpendicular to B and one || to B. That at least gives you some understanding of the problem.

syms m_e q_e positive % remember that q_e now is the positive elementary charge
syms E [3,1] real
syms B [3,1] real
syms t
syms r(t) [3,1]    
D2rDt2 = diff(r,2)
syms r0 [3,1] real
syms v0 [3,1] real
DrDt = diff(r);
qwe = dsolve(D2rDt2 == -q_e/m_e*(E + cross(DrDt,B)),r(0)==r0,Dr(0) == v0);
% Check those general solutions and appreciate the power of orienting your
% coordinates!
qwe.r1
qwe.r2
qwe.r3
%% A slightly more sensible setup:
syms Bz real        
syms Ex Ez real
asd = dsolve(D2rDt2 == -q_e/m_e*([Ex;0;Ez] + cross(DrDt,[0;0;Bz])),r(0)==r0,Dr(0) == v0);
% this solution you will have to look into more carefully, preferably with
% pen and paper, or with repeated use of the simplify and expand functions.

Then for plotting I strongly suggest you don't bother "eliminating t", simply calculate the electron-gyro-frequency and the gyro-period, and plot the trajectories for a couple of gyro-periods (5-8 or something in that range).

HTH

[Khoa]

Thank you so much!! I think I will need some time to digest your code but you provided me with so many useful info.

[Khoa]

The "" is used to assign initial velocity vector into the (I think).

However, looks like there is a minor problem.

sol = dsolve(D2rDt2 == -q_e/m_e*(E + cross(DrDt,B)), r(0) == r0, Dr(0) == v0 );

output:

Unrecognized function or variable 'Dr'.

Error in pro_trajectory (line 34)
sol = dsolve(D2rDt2 == -q_e/m_e*(E + cross(DrDt,B)), r(0) == r0, Dr(0) == v0 );

I got the idea but don't know the right syntax to work this around. Please help me

Thank you in advance

[Bjorn Gustavsson]

My bad. I changed the notation for clarity, but didn't do it throughout. it should be:

sol = dsolve(D2rDt2 == -q_e/m_e*(E + cross(DrDt,B)), r(0) == r0, DrDt(0) == v0 );

For a second-order ODE you need to specify the initial conditions for both the solution and its derivative at the start (or two boundary-conditions on two points).

HTH

[Khoa]

Thanks!! I also figured that out.

However, I'm sorry to bother you again.

When I run the code, the terminal said that :

Warning: Unable to find explicit solution. 

Then I type into the terminal:

>> sol
 
sol =
 
[ empty sym ]

There is nothing, no function, no variable. :(

[Bjorn Gustavsson]

That's a bit peculiar. When I run the code I get this output (for asd.r1):

>> asd.r1
 
ans =
 
(2*m_e*(v01*tan((Bz*q_e*t)/(2*m_e)) + (Ex + Bz*v02)/Bz))/(Bz*q_e*(tan((Bz*q_e*t)/(2*m_e))^2 + 1)) - (- q_e*r01*Bz^2 + 2*m_e*v02*Bz + 2*Ex*m_e)/(Bz^2*q_e)

This is run on matlab R2020a, what version are you using, is it a student-license?

[Khoa]

I am currently using the MATLAB R2019b version

When I put this into the terminal:

sol.r1

i got:

Dot indexing is not supported for variables of this type.
Error in pro_trajectory (line 36)
sol.r1

And btw what is this .r1 does. I'm curious

[Bjorn Gustavsson]

The capabilities of the symbolic tools evolve through the releases. I'll check if I have an older release somewhere and will try to see if it produces the same result as you see.

In principle dsolve integrates the equations of motion for the three components of r(t), and stores them as

sol.r1, sol.r2 and sol.r3 - it is up to us to re-label them as x(t), y(t) and z(t) as we see fit. If your dsolve cannot handle this set of equations sol will be empty.

Perhaps you can try to dumb-it-down and solve the equations of motion for a charged particle in a constant and homogenous E-field (1-D problem so should be simplest). Then perhaps also an electron gyrating in a constant and homgenous B-field, since this is effectively a 2-D problem it should be slightly more complicated for dsolve to handle, you should already know the solution so you can easily check.

HTH

[Khoa]

So I test the code again but on R2022b (Web version).

And looks like it worked !!!!

>> qwe.r1
 
ans =
 
(m_e*(cos((q_e*t)/m_e) + sin((q_e*t)/m_e)))/q_e - (10*m_e - q_e)/(10*q_e)

>> qwe.r2
 
ans =
 
m_e/q_e - (m_e*(cos((q_e*t)/m_e) - sin((q_e*t)/m_e)))/q_e

>> qwe.r3
 
ans =
 
t - (q_e*t^2)/(2*m_e)
 

But now i have problem plotting it .-.

Error using fplot3
Input must be a function or functions of a single variable.

Error in PRO_6 (line 22)
fplot3(qwe.r1,qwe.r2,qwe.r3)

It confused me. Obviously the input was a single variable function. Am I missing something (the interval of t? maybe or maybe this is not the plotting function I am looking for)

Hope you can help. TY

[Bjorn Gustavsson]

The symbolic tools retrurn symbolic expressions. You can convert those convert those to ordinary matlab-functions. Either you do this "by hand" - which might be a good exercise for you as a student for this reasonably simple case, or you can create functions with the matlabFunction function - check the corresponding help and documentation. It could work something like this:

r1_fcn = matlabFunction(asd.r1); % this creates a function-handle for x(t). Cycl.

From there you'd have to define the inputs to the functions, the electron mass, magnetic field-strength etc:

m_e = 9.11e-31; % Store these physical constants in one file to have at 
q_e = 1.62e-19; % hand at all times! These are plucked from memory to act as a 
% detering example!
B_z = 5e-3;
E_x = 1;
T_gyro = pi;
t = linspace(0,7*T_gyro,701);
plot3(r1_fcn(q_e,m_e,B_z,E_x,t),...
      r2_fcn(q_e,m_e,B_z,E_x,t),...
      r3_fcn(q_e,m_e,B_z,E_x,t),...
     ',r.-')

You'll have to check the order of the argument to each function and such on yourself, just as the values of eveything else. This way you'll get the trajectory plotted with 100 points per gyration, and equal time between each timestep.

HTH

[Khoa]

OMG! THANK YOU FOR EVERYTHING SIR I FINALLY GET IT TO WORK

[Bjorn Gustavsson]

Good!

(If someone tasks you with solving the equations of motion for charged particles in E and B-fields have a look at the Boris-mover, it is a way better integration-method for these eq-s of motion than the standard ODE-solvers. The wikipedia-entry is sufficient.)