I fixed it, thankfully, using rootshuffle.m https://www.mathworks.com/matlabcentral/fileexchange/29464-rootshuffle-m?s_tid=srchtitle !!! Apparently, it is because matlab finds the roots, but doesn't necessarily put them in order, thereby causing sign flipping of the roots!!!
Getting strange vertical lines when plotting roots of biquadratic equation/ dispersion relation
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Katey Stevenson
el 21 de Sept. de 2022
Editada: Katey Stevenson
el 21 de Sept. de 2022
Hello,
I'm simply trying to plot the real and imaginary parts of the two positive roots of a biquadratic equation/dispersion relation, p^4 + C3*p^2 + C5 = 0, the plot looks correct apart from these weird vertical lines which seems to be from the roots changing sign...I've racked my brain and I am not sure of what I'm doing wrong. I've put in a divide by zero check, but it doesn't appear that dividing by zeros is the problem.
Here is the plot I'm getting:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1132190/image.jpeg)
and here is what it's supposed to look like:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1132195/image.jpeg)
Thank you in advance for any input!!
%% constants %%
%absolute value of electron charge
e = 1.60217662E-19; %[C]
%speed of light
c = 299792458; %[m/s]
boltzmann = 1.380649E-23; %[J/K]
m_e = 9.0194E-31; %[kg]
m_i = 6.6335209E-26; %[kg] %argon
% permittivity of free space
eps_0 = 8.8542E-12; %[F/m]
steps = 1e+3;
B = linspace(0.0002,0.1,steps);%[TESLA]
e_coll_freq = 9e+4;
ion_coll_freq = 6e+5;
dens = 5e+18; %[m^-3]
freq = 8e+6;
w = 2*pi*freq;
L = 0.8; %[m]
n = 2;
k_parallel = n*pi/L;
N_para = k_parallel*c/w;
% electron cyclotron frequency
w_ce = e.*B/m_e; %[Hz]
% ion cyclotron frequency
w_ci = e.*B/m_i; %[Hz]
% plasma electron frequency
w_pe = sqrt((dens*e^2)/(eps_0*m_e));
% plasma ion frequency
w_pi = sqrt((dens*e^2)/(eps_0*m_i));
w_LH = sqrt(((w_ce.*w_ci).^2 + (w_pi.*w_ce).^2)./(w_pe^2 + w_ce.^2));
%%%%%
%%% check for dividing by zero %%%
e_denom = (w_ce.^2 - ((1 + 1i*e_coll_freq/w)^2)*w^2);
divZeroElements = find(abs(e_denom)<0.01);
e_denom(divZeroElements) = NaN;
i_denom = (w_ci.^2 - ((1 + 1i*ion_coll_freq/w)^2)*w^2);
divZeroElements2 = find(abs(i_denom)<0.01);
eps1 = (1 + ((1 + 1i*e_coll_freq/w)*w_pe^2)./e_denom + ((1 + 1i*ion_coll_freq/w)*w_pi^2)./i_denom);
eps2 = (-((w_pe^2).*(w_ce./w))./e_denom + ((w_pi^2)*(w_ci./w))./i_denom);
eps3 = (1 - w_pe^2/(((1 + 1i*e_coll_freq/w)^2)*w^2) - w_pi^2/(((1 + 1i*ion_coll_freq/w)^2)*w^2));
Alpha = eps1 - N_para^2 - (eps2.^2)./eps1;
Beta = eps3.*(1-((N_para^2)./eps1));
Delta = N_para.*eps2./eps1;
Gamma = (N_para.*eps2.*eps3)./eps1;
%coefficients of 4th order polynomial
C1 = 1;
C2 = 0;
C3 = (-1).*(Alpha+Beta);
C4 = 0;
C5 = (Alpha.*Beta) - (Gamma.*Delta);
ROOTS1(:,1) = sqrt((-C3 + sqrt(C3.^2 - 4.*C5))./2);
ROOTS2(:,1) = sqrt((-C3 - sqrt(C3.^2 - 4.*C5))./2);
ROOTS3(:,1) = -sqrt((-C3 + sqrt(C3.^2 - 4.*C5))./2);
ROOTS4(:,1) = -sqrt((-C3 - sqrt(C3.^2 - 4.*C5))./2);
% normalized x_axis
x_axis = w_LH./w;
figure()
plot(x_axis,imag(ROOTS1),'--k')
hold on
plot(x_axis,real(ROOTS1))
plot(x_axis,real(ROOTS2))
plot(x_axis,imag(ROOTS2),'--b')
ylim([-1.5e+4 1.5e+4])
xlabel('\omega_{LH}/\omega')
ylabel('p')
0 comentarios
Respuesta aceptada
Katey Stevenson
el 21 de Sept. de 2022
Editada: Katey Stevenson
el 21 de Sept. de 2022
0 comentarios
Más respuestas (0)
Ver también
Categorías
Más información sobre Bartlett en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!