Matlab Ignoring Half of Equation When Solving ODE

Hello,

I am trying to solve the Lorentz Force differential equation set in 3 space, however I am currently omitting the magnetic field entirely so ideally, the program would only look at the electric field. However, it is not and I know this because I have changed the electric field to be very different magnitudes while keeping the time step the same and the electron still travels the same distance when I look at the plot which it shouldn't. For some reason Matlab is completely ignoring the electric field entirely. Why is this the case? (All code displayed is also attached in seperate files)

PS - I have tried all the ODE solvers that I could and the same issue kept happening. I am also aware that in this specific case, the ODE solver is not exactly needed but it will be in the more complex version of the code

Main script:

%User inputs
%TSTEP MUST BE ABLE TO GO INTO TFIN AN ODD AMOUNT OF TIMES (Because simulation includes tstep @ 0 so it ends up being even)
tstep = 0.0625E-9; %Defining time step
tfin = 5.25E-8; %Defining final time
particlecount = 1;
%# of particles to generate energies for energy generation
s1 = 0.0200001; %Specify minimum radius particles may be generated
s2 = 0.0200002; %Specify maximum radius particles may be generated
phi1 = 90;
phi2 = -90;
%cs = 1.5E-20; %Cross-sectional area input in (m^2). (Held constant for sake of simplicity)
cs = 0;
energyaddition = 0; %Absolute value of the E-field strength of experiment in (eV)
inelasticenergyloss = 0; %Amount of energy an electron loses per collision in (eV)
scaleenergy = 0; %If need to scale all energies by a certain value (constant/unitless)
mTorr = 750; %Define pressure in mTorr (From metter reading)
ict = 290; %Temperature in Kelvin of neutral gas. Converted from C and @ room temp

%Natural Constants
m_e = 9.11E-31; %Mass of electron
jtoev = 6.24E+18; %Joules to eV conversion
epnaut = 8.854187E-12;
k = 1/(4*pi*epnaut);

%Create zeros matrices to populate later
energymat = zeros(1,1);
vmagmat = zeros(1,1);
intspan = [0:tstep:tfin]; %Total time span
intspancol = intspan(:); %inverts matrix's rows and columns
[introw,intcol] = size(intspancol);
icvelocity(1).matrix = zeros(1,3);
icposition(1).matrix = zeros(1,3);

%Convert additional energy to Joules
entrans = energyaddition*(1.6022E-19);
%Energy conversion to Joules for inelastic collisions
inelasticenergy = inelasticenergyloss*(1.6022E-19);

%Generates energy for particle
numin = 1E-19;

%Generate random positions for each particle between s1 and s2
%In Spherical: Generates random position with boundaries of: s1<r<s2
rrand = (s2 - s1)*rand() + s1; %Random # generation for component x between s1 and s2
thetarand = (2*pi)*rand(); %Random # generation for component y between s1 and s2
phirand = asin((sin(phi2rad) - sin(phi1rad))*rand() + sin(phi1rad)); %Random # generation for component z between s1 and s2
[xrand,yrand,zrand] = sph2cart(thetarand,phirand,rrand);
icposition = [xrand yrand zrand];

%Generates velocity components for each particle
vmagin = sqrt(2*energymat(1,1)/m_e);
vmagmat(1,1) = vmagin;
%vx = (icposition(1)/norm(icposition(1:3)))*vmagin;
%vy = (icposition(2)/norm(icposition(1:3)))*vmagin;
%vz = (icposition(3)/norm(icposition(1:3)))*vmagin;
vx = 1E+1;
vy = vx;
vz = vy;
icvelocity = [vx vy vz];

%Generate matrix that will be populated by positions and velocities (Each matrix within structure "Wposandvel" should be of size 2*introw-1 by 6)
Wposandvel = zeros(2*introw-1,6);

%This section is purely for generating matrices the program will populate later
tindex = 0:1:introw-2;
tindexcol = tindex(:); %inverts matrix's rows and columns
[tinrow,tincol] = size(tindex);
PCPmat = zeros(tincol,1); %Probability of collision solutions
choicemat = zeros(tincol,1); %Set of RNG choices system makes
colnewv = zeros(intcol,3); %New velocity components after collision
velocitytransitionmat = zeros(intcol,3);
energytransitionmat = zeros(intcol,1);

spart(1) = icposition(1,1);
spart(2) = icposition(1,2);
spart(3) = icposition(1,3);
vpart(1) = icvelocity(1,1);
vpart(2) = icvelocity(1,2);
vpart(3) = icvelocity(1,3);

for t = 0:1:introw-2 %complete time interval and time step
        
        [trow,tcol] = size(t);
		
		%Defining relative position and velocity
		x = spart(1);
		y = spart(2);
		z = spart(3);
		vx = vpart(1);
		vy = vpart(2);
		vz = vpart(3);
        
      r = sqrt(x.^2 + y.^2 + z.^2);
	vmag = sqrt(vx.^2 + vy.^2 + vz.^2);
	vmagtimemat(1).matrix(t+1,1) = vmag; %Populate vmagtimemat structure
      entrack = ((.5)*m_e*(vmag^2)); %Define energy at current time step
        
	if r <= 0.02 %Defining particle collision with sphere surface. If collision occurs, vmag = 0
          vx = 0;
          vy = 0;
          vz = 0;
	    break
      end
        
		%Coupled differential equation solver for trajectory within current time step
	icv = [x; y; z; vx; vy; vz]; %Initial conditions
      tspan = [intspan(t+1) ((intspan(t+2)-intspan(t+1))/2)+intspan(t+1) intspan(t+2)]; %Time span
        
        %Trajectory solver
    	[T,S] = ode15s(@bdipuniodefun, tspan, icv);
    	[rownum,colnum] = size(S);
       Wposandvel((1+2*t):(3+2*t),(1:6)) = S;
        
	%Redfine the velocity and position components to reference on next if-loop run
      spart(1) = Wposandvel((3+2*t),1);
	spart(2) = Wposandvel((3+2*t),2);
	spart(3) = Wposandvel((3+2*t),3);
	vpart(1) = Wposandvel((3+2*t),4);
	vpart(2) = Wposandvel((3+2*t),5);
	vpart(3) = Wposandvel((3+2*t),6);
        
end


%Plotting Stuff
[X,Y,Z] = meshgrid(-0.1:0.005:0.1,-0.1:0.005:0.1,-0.1:0.005:0.1);

%{
%Plotting dipole B-field
[Bx, By, Bz] = B_test();
Bfieldx = arrayfun(Bx,X,Y,Z);
Bfieldy = arrayfun(By,X,Y,Z);
Bfieldz = arrayfun(Bz,X,Y,Z);
%}

%{
%Plotting Efield
[Ex, Ey, Ez] = E_test();
Efieldx = arrayfun(Ex,X,Y,Z);
Efieldy = arrayfun(Ey,X,Y,Z);
Efieldz = arrayfun(Ez,X,Y,Z);
%}

%{
%Plotting B-field by itself
figure
%subplot(3,3,2)
%Graph labels
axis equal
xlabel 'x in m'
ylabel 'y in m'
zlabel 'z in m'
title('B-Field')
%Changes figure background color
set(gca,'color','k');
quiver3(X,Y,Z,Bfieldx,Bfieldy,Bfieldz,'color','b') %Plots B-field
%}

%{
%Plotting E-field by itself
figure
%subplot(3,3,3)
%Graph labels
axis equal
xlabel 'x in m'
ylabel 'y in m'
zlabel 'z in m'
title('E-field')
%Changes figure background color
set(gca,'color','k');
quiver3(X,Y,Z,Efieldx,Efieldy,Efieldz,'color','g') %Plots E-field
%}


%Plotting trajectories only
figure
%Graph labels
axis equal
xlabel 'x in m'
ylabel 'y in m'
zlabel 'z in m'
title('Electron Trajectory Solutions')
%Changes figure background color
set(gca,'color','k');
hold on
Wposandvel(~any(Wposandvel(1:2*introw-1,1:6),2),:) = []; %Clears all rows in collpos that contain only zeros
plot3(Wposandvel(:,1),Wposandvel(:,2), Wposandvel(:,3), '-r', 'LineWidth',2,'color','w') %Plot trajectories

%Generate a sphere consisting of 20 by 20 faces
[x,y,z]=sphere;
% use surf function to plot
SSurface = surf(sphere1radius*x+sphere1centerx,sphere1radius*y+sphere1centery,sphere1radius*z+sphere1centerz);
set(SSurface,'FaceColor',[0.5 0.5 0.5],'FaceLighting','gouraud','EdgeColor','none')
camlight

scatter3(icposition(:,1),icposition(:,2),icposition(:,3),1,'filled','y') %Plot initial position points
    
hold off
view(90,90)

"bdipuniodefun" script:

function bdip = bdipuniodefun(t,s)

%Using SI units
q = -1.60217662E-19;
m_e = 9.11E-31;

persistent Bx By Bz Ex Ey Ez

if isempty(Bx)
    [Bx, By, Bz] = B_test();
end

if isempty(Ex)
    [Ex, Ey, Ez] = E_test();
end

%bdip = [s(4); s(5); s(6); (q/m_e)*(Ex(s(1),s(2),s(3)) + s(5)*Bz(s(1),s(2),s(3)) - s(6)*By(s(1),s(2),s(3))); (q/m_e)*(Ey(s(1),s(2),s(3)) + s(6)*Bx(s(1),s(2),s(3)) - s(4)*Bz(s(1),s(2),s(3))); (q/m_e)*(Ez(s(1),s(2),s(3)) + s(4)*By(s(1),s(2),s(3)) - s(5)*Bx(s(1),s(2),s(3)))];
bdip = [s(4); s(5); s(6); (q/m_e)*(Ex(s(1),s(2),s(3))); (q/m_e)*(Ey(s(1),s(2),s(3))); (q/m_e)*(Ez(s(1),s(2),s(3)))];

end

"E_test" script:

function [Ex, Ey, Ez] = E_test()

syms x y z
R_s = 0.02;
V = 2000;
epnaut = 8.854187E-12;
k = 1/(4*pi*epnaut);
Q = (V*R_s)/k;
r = [x, y, z];
E = (k*Q/norm(r)^3)*r;
Ex = matlabFunction(E(1));
Ey = matlabFunction(E(2));
Ez = matlabFunction(E(3));

end

Answer 1

[Akhil Gopinath]

Hi Tom

I pressume the issue is with persistant variable you have created in the 'bdipuniodefun' these variables will cling to the value stored in them untill you clear the function 'bdipuniodefun'. you can find more about it here: https://www.mathworks.com/help/matlab/ref/persistent.html

[Tom Keaton]

Hi Akhil,

Thank you for responding. I am using the persistent variable definitions because it saves computation time and I don't want the for loop to have to calcualte the B and E fields in all space in every time step. Without the definition, the solver pretty much froze. What would be your proposed solution for this (that perhaps doesn't invovle persistent variables) so that I can both save computation time and avoid this issue (Or is that in the link you sent, because I have already looked at that documentation and have not found something related to this issue there)? I am new to the coding world in general so I have only thought up the method I am currently using.

Tom

[Akhil Gopinath]

Hi Tom,

I believe I was not clear enough, I think the issue is that persistant variables are deleted only when you do a

clear function_name

Else, even when you start a new simulation with different parameters, It will start with the last stored value in the memory.

I could not observe any clear function in the code, hence the suggestion