Matlab Ignoring Half of Equation When Solving ODE
2 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
回答(1 个)
Akhil Gopinath
2020-1-16
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
2 个评论
Akhil Gopinath
2020-1-27
编辑:Akhil Gopinath
2020-1-27
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
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!