Why are the output variables, that I expect to be real, complex?

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

[Torsten]

%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)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

Answer 1

[Torsten]

%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)

[RITA PIA SESSA]

Thank you! I think I made some mistakes when writing the (very long) equation and that is why the shape is not quite the same :)

Answer 2

[Walter Roberson]

syms omega_r omega_i real

might help. If not then pass in

[-inf inf; -inf inf]

after the list of variables.

[RITA PIA SESSA]

I tried both your suggestions, but received the following error, always in the same code line, that is where the vpasolve is.

[Walter Roberson]

[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