solve numerically a nonlinear (plasma wave dispersion) relation

Hi everyone,

I'm currently trying to solve numerically the following nonlinear (dispersion) relation (i.e. wave frequency f or pulsation omega, as a function of wave number k_B) with 8 fixed parameters:

that should have 3 solutions (red, light blue and orange curves), according to the following plot:

thus I wrote the following simple code :

Parameters

 c = 3e5; 
 we = 8e4; 
 wp = 1.3e3;
 wc = 3.1e2;
 Omegae = -4.4e2; 
 Omegap = 2.4e-1; 
 Omegac = 1.3e-2; 
 uc = 402;

 %%Set arrays
 k_array = 1e-5:1e-6:1e-2; % k values to solve for
 w_array = zeros(length(k_array),1); 3 solutions expected for each k
 i = 1; % loop counter

 %%Loop over k
 for k = k_array(1):k_array(end)

 clear w
 syms w

 eqn = c^2*k^2/w^2 + we^2/(w^2-Omegae^2)...
 +wp^2/(w^2-Omegap^2)+(wc^2/w^4)*(w^2+k^2*uc^2)...
 +(2*w^2/(wc^2*k^2*uc^2)*(Omegae*we^2/(w^2-Omegae^2)+Omegap*wp^2/(w^2-Omegap^2))^2/...
 (1+sqrt(1+4*w^6/(wc^4*k^4*uc^4)*(Omegae*we^2/(w^2-Omegae^2)+Omegap*wp^2/(w^2-Omegap^2))^2))) == 0;

 root = vpasolve(eqn,w);
 w_array(i,:) = double(root);
 i=i+1;
 end

but I get the error: Subscripted assignment dimension mismatch, because the solver does not find any solution (empty sym).

Where did I go wrong ?

Thanks for your help. HB

[Star Strider]

Before you attempt to evaluate ‘eqn’, explore:

c = 3e5;
we = 8e4;
wp = 1.3e3;
wc = 3.1e2;
Omegae = -4.4e2;
Omegap = 2.4e-1;
Omegac = 1.3e-2;
uc = 402;

syms k w
eqn = c^2*k^2/w^2 + we^2/(w^2-Omegae^2)...
     +wp^2/(w^2-Omegap^2)+(wc^2/w^4)*(w^2+k^2*uc^2)...
     +(2*w^2/(wc^2*k^2*uc^2)*(Omegae*we^2/(w^2-Omegae^2)+Omegap*wp^2/(w^2-Omegap^2))^2/...
     (1+sqrt(1+4*w^6/(wc^4*k^4*uc^4)*(Omegae*we^2/(w^2-Omegae^2)+Omegap*wp^2/(w^2-Omegap^2))^2)));
eqn = expand(eqn);
eqn = simplify(eqn, 'Steps', 50)
sw = solve(eqn,w, 'ReturnConditions',true)

w = sw.w
conds = sw.conditions

You will need to adjust your code to accommodate the roots.

[Hugo B]

Dear Star Strider,

thanks a lot for your quick answer!

However, when I compute the roots as you indicated, I see 3 roots in [0 Inf] (which is expected) but they are all identical (and thus are the conditions), so I'm not sure if it cannot solve or what I have to do to adapt the code (I cannot change the parameters as I'm trying to reproduce the results obtained in the figure above).

Or maybe I missed something?

Thanks in advance for your help.

[Torsten]

If you plot your function for the parameters given, you'll see that it has no real roots.

Best wishes

Torsten.

[Hugo B]

Dear Torsten,

Thanks for your answer. Indeed 2 solutions (omega_R and omega_L) should be complex (whose real parts correspond to orange and light blue curves) and one (omega_c, the one I'm interested in) is expected to be purely imaginary (i.e. its real part is 0 which corresponds to the red curve and its imaginary part is the black curve, i.e., gamma_c).

But when I try to solve the equation it gives no root at all (or maybe I should specify that I want complex and/or imaginary solutions and not only real roots?)

Thanks again for your help!

[David Goodmanson]

Hi Hugo, I can see the units on the plot, but what units are being used in the calculation? 3e5 would seem to imply sec and km.

[Hugo B]

Hi David,

thanks for your answer. In the calculation indeed units are in km and s (i.e., resulting k should be in km-1 and w in rad.s-1)