Question regarding MATLAB code to solve for electron in constant electric and magnetic field

7 views (last 30 days)
Hello all the kind soul, I am trying to simulate the trajectory of an electron in constant electric and magnetic field. However, when I am running the code, I keep getting the error "Unable to perform assignment because the left and right sides have a different number of elements." Can anyone please help me, any kind of help will do. Thank you so much. Below is my code.
function dsdt=ExBfunc(~, s_vect)
%s_vect(1) = x
%s_vect(2) = y
%s_vect(3) = z
%s_vect(4) = dx/dt
%s_vect(5) = dy/dt
%s_vect(6) = dz/dt
global Bx; global By; global Bz;
global Ex; global Ey; global Ez;
global q; global m;
dsdt=zeros(6, 1);
dsdt(1) = s_vect(4);
dsdt(2) = s_vect(5);
dsdt(3) = s_vect(6);
dsdt(4) = (q/m)*(Ex + s_vect(5)*Bz-s_vect(6)*By);
dsdt(5) = (q/m)*(Ey + s_vect(6)*Bx-s_vect(4)*Bz);
dsdt(6) = (q/m)*(Ez + s_vect(4)*Bz-s_vect(5)*Bx);
end
clear
close all
clc
%Define constant variables q and m
q=1.602*10^(-19);
m=9.11*10^(-31);
%Define the value for constant electric field
Ex=0;
Ey=0;
Ez=1;
%Define the value for constant magnetic field
Bx=1;
By=0;
Bz=0;
%Define the initial position and initial velocity of the electron as a
%column vector
r0=[0 0 0]';
v0=[0 0 0]';
%Pass in the constants to the defined function E_3D and B_3D
%E=E_3D(Ex, Ey, Ez);
%B=B_3D(Bx, By, Bz);
s0=[r0; v0];
%Define the timespan
t0=0;
T=10;
tspan=[t0 T];
disp('Simulation begins');
%Solve the differential equation by numerical method using ode45
[t, s_vect]=ode45(@ExBfunc, tspan, s0');
disp('Simulation done');
%Plot the trajectory of electron of on each axis against time
plot(t,s_vect(:,1));
plot(t,s_vect(:,2));
plot(t,s_vect(:,3));
disp('Plotting done');

Answers (1)

Alan Stevens
Alan Stevens on 14 Oct 2021
  1. If you are going to use global variables they need to be declared outside the function as well as inside.
  2. Your values of q and m are such that the simulation takes forever, so scale S = s*m/q and then rescale s_vector*q/m at the end.
global Bx; global By; global Bz;
global Ex; global Ey; global Ez;
%Define constant variables q and m
q=1.602*10^(-19);
m=9.11*10^(-31);
%Define the value for constant electric field
Ex=0;
Ey=0;
Ez=1;
%Define the value for constant magnetic field
Bx=1;
By=0;
Bz=0;
%Define the initial position and initial velocity of the electron as a
%column vector
r0=[0 0 0]';
v0=[0 0 0]';
%Pass in the constants to the defined function E_3D and B_3D
%E=E_3D(Ex, Ey, Ez);
%B=B_3D(Bx, By, Bz);
s0=[r0; v0];
%Define the timespan
t0=0;
T=10;
tspan=[t0 T];
disp('Simulation begins');
Simulation begins
%Solve the differential equation by numerical method using ode45
[t, s_vect]=ode45(@ExBfunc, tspan, s0);
disp('Simulation done');
Simulation done
s_vect = s_vect*q/m;
%Plot the trajectory of electron of on each axis against time
subplot(3,1,1)
plot(t,s_vect(:,1));
subplot(3,1,2)
plot(t,s_vect(:,2));
subplot(3,1,3)
plot(t,s_vect(:,3));
disp('Plotting done');
Plotting done
function dsdt=ExBfunc(~, s_vect)
%s_vect(1) = x
%s_vect(2) = y
%s_vect(3) = z
%s_vect(4) = dx/dt
%s_vect(5) = dy/dt
%s_vect(6) = dz/dt
global Bx; global By; global Bz;
global Ex; global Ey; global Ez;
dsdt=zeros(6, 1);
dsdt(1) = s_vect(4);
dsdt(2) = s_vect(5);
dsdt(3) = s_vect(6);
dsdt(4) = (Ex + s_vect(5)*Bz-s_vect(6)*By);
dsdt(5) = (Ey + s_vect(6)*Bx-s_vect(4)*Bz);
dsdt(6) = (Ez + s_vect(4)*Bz-s_vect(5)*Bx);
end
  2 Comments
Song Hang Chai
Song Hang Chai on 14 Oct 2021
Hi Alan, I really appreciate your help. My code is now running fine!
However, I would like to ask one silly question, how did you know that my simulation takes forever to run with those values of q and m(That is really the case). Is it because that my electron is doing something unphysical, like it is moving exceeding or close the speed of light, but I didn't take into the account of realistivitic manner? Once again, I really appreciate your help
Alan Stevens
Alan Stevens on 15 Oct 2021
The ratio of q/m is around 10^11. I think this might give ode45 convergence problems, though I don't know for sure.

Sign in to comment.

Categories

Find more on MATLAB in Help Center and File Exchange

Tags

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!