Why are the output variables, that I expect to be real, complex?
1 view (last 30 days)
Show older comments
RITA PIA SESSA
on 16 Sep 2022
Commented: Walter Roberson
on 18 Sep 2022
Hello, I am writing a code to plot a mathematical model that depends on real number k and complex number omega=omega_r+i*omega_i.
In order to find the values of omega_r and omega_i that satisfy the condition: equation==0, I performed a for loop to change k with each iteration and find the corresponding needed values of omega_r and omega_i.
I am expecting to obtain two real-valued outputs, that are the real and imaginary coefficient of the complex number omega, but at the end I obtain that both omega_r and omega_i are actually complex numbers, too.
At the end, I would expect to plot a curve in the k - omega_i space as in the picture (the continuous line).
Can anybody help me understand what I did wrong? Thank you in advance
1 Comment
Torsten
on 16 Sep 2022
%Constant parameters%
a = 0.005; %sheet half thickness%
eps0 = 8.9*10^-12; %vacuum dielectric constant%
eps_2 = 80; %water relative dielectric constant%
epsg = 1.006; %air relative dielectric constant%
epsilon_G = eps0*epsg; %air dielectric constant%
epsilon_L = eps0*eps_2; %water dielectric constant%
gamma = 0.0728; %water-air surface tension%
mu = 8.94*10^-4; %water viscosity%
rho_2 = 1000; %water density%
rhog = 1.2250; %air density%
sigma = 500; %water conductivity%
%Plotted conditions%
We = 400; %Weber number%
Re = 1000; %Reynolds number%
rho = 0.001; %density ratio%
tau = 0.01; %hydrodynamic time/electrical relaxation time%
U_g = 0.3; %gas velocity ratio%
Eu = 0; %Euler number%
%Dispersion relation part%
syms omega_r omega_i
omega_r_arr = zeros(21,0);
omega_i_arr = zeros(21,0);
k = zeros(21,0);
k(1) = 0;
for j = 1:21
omega = omega_r+1i*omega_i; %complex number omega%
omega_ = k(j)*U_g-omega; %omega soprasegnato%
l_2 = k(j)^2+1i*k(j)*Re-1i*omega*Re; %l^2%
l = sqrt(l_2); %l%
DR = ((k(j)^3)/We)...
+((omega_^2)*(-1*rho))...
+((((k(j)^2+l_2)^2)*tanh(k(j)))-((4*k(j)^3)*l*tanh(l))/(Re^2))...
+((Eu*(k(j)^2)*((epsg^2*tanh(k(j))*(k(j)^2-l_2)^2)+(eps_2^2*tanh(k(j))*(k(j)^2-l_2)*(k(j)^2-tau*Re-l_2))+(epsg*eps_2*((k(j)^2)*tanh(k(j))*((k(j)^2)*(-2)+2*tau*Re+tanh(k(j))*tau*Re)-l*(2*tau*Re*k(j)*(1+tanh(k(j)))*tanh(l)-l*tanh(k(j))*(4*(k(j)^2)+tanh(k(j))*tau*Re-l_2*2))))))/((k(j)^2-l_2)*((k(j)^2-tau*Re-l_2)*eps_2+(k(j)^2-l_2)*tanh(k(j))*epsg))) == 0;
[wr(j), wi(j)] = vpasolve(DR, [omega_r, omega_i]);
omega_r_arr(j) = wr(j);
omega_i_arr(j) = wi(j);
k(j+1) = k(j)+0.01;
end
k=0:0.01:0.2;
plot(k,omega_i_arr)
Accepted Answer
Torsten
on 16 Sep 2022
%Constant parameters%
a = 0.005; %sheet half thickness%
eps0 = 8.9*10^-12; %vacuum dielectric constant%
eps_2 = 80; %water relative dielectric constant%
epsg = 1.006; %air relative dielectric constant%
epsilon_G = eps0*epsg; %air dielectric constant%
epsilon_L = eps0*eps_2; %water dielectric constant%
gamma = 0.0728; %water-air surface tension%
mu = 8.94*10^-4; %water viscosity%
rho_2 = 1000; %water density%
rhog = 1.2250; %air density%
sigma = 500; %water conductivity%
%Plotted conditions%
We = 400; %Weber number%
Re = 1000; %Reynolds number%
rho = 0.001; %density ratio%
tau = 0.01; %hydrodynamic time/electrical relaxation time%
U_g = 0.3; %gas velocity ratio%
Eu = 0; %Euler number%
%Dispersion relation part%
syms omega_r omega_i
omega_r_arr = zeros(21,0);
omega_i_arr = zeros(21,0);
k = zeros(21,0);
k(1) = 0;
for j = 1:21
omega = omega_r+1i*omega_i; %complex number omega%
omega_ = k(j)*U_g-omega; %omega soprasegnato%
l_2 = k(j)^2+1i*k(j)*Re-1i*omega*Re; %l^2%
l = sqrt(l_2); %l%
DR = ((k(j)^3)/We)...
+((omega_^2)*(-1*rho))...
+((((k(j)^2+l_2)^2)*tanh(k(j)))-((4*k(j)^3)*l*tanh(l))/(Re^2))...
+((Eu*(k(j)^2)*((epsg^2*tanh(k(j))*(k(j)^2-l_2)^2)+(eps_2^2*tanh(k(j))*(k(j)^2-l_2)*(k(j)^2-tau*Re-l_2))+(epsg*eps_2*((k(j)^2)*tanh(k(j))*((k(j)^2)*(-2)+2*tau*Re+tanh(k(j))*tau*Re)-l*(2*tau*Re*k(j)*(1+tanh(k(j)))*tanh(l)-l*tanh(k(j))*(4*(k(j)^2)+tanh(k(j))*tau*Re-l_2*2))))))/((k(j)^2-l_2)*((k(j)^2-tau*Re-l_2)*eps_2+(k(j)^2-l_2)*tanh(k(j))*epsg)));
[wr(j), wi(j)] = vpasolve([real(DR)==0,imag(DR)==0], [omega_r, omega_i]);
omega_r_arr(j) = wr(j);
omega_i_arr(j) = wi(j);
k(j+1) = k(j)+0.01;
end
k=0:0.01:0.2;
plot(k,omega_i_arr)
More Answers (1)
Walter Roberson
on 16 Sep 2022
syms omega_r omega_i real
might help. If not then pass in
[-inf inf; -inf inf]
after the list of variables.
2 Comments
Walter Roberson
on 18 Sep 2022
[wr_out, wi_out] = vpasolve(DR, [omega_r, omega_i], [-inf inf; -inf inf]);
if isempty(wr_out)
wr(j) = sym(nan);
wi(j) = sym(nan);
else
wr(j) = wr_out(1);
wi(j) = wi_out(1);
end
See Also
Categories
Find more on Lighting, Transparency, and Shading in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!