Plotting $\omega(k_x, k_z)$ in $(\omega, k)$ plane with assumption of $k^2 = k_x^2 + k_z^2$

Question

I'm currently studying atmospheric-gravity waves and their dispersion relations.

For analysis I need to plot the dispersion function:

\begin{equation} \omega^4 - \omega^2\cdot C_s^2(k_x^2 + k_z^2 + \dfrac{1}{4H^2}) + (\gamma - 1)g^2\cdot k_x^2 = 0 \end{equation} where, $C_s^2$, $H$, $\gamma$ and $g$ are some constant values.

I used Solve to get $\omega$ as a function of $k_x^2$ and $k_z^2$ - $\omega = \omega(k_x^2, k_z^2)$.

Code for solving:

   Solve[ω^4 - ω^2*(Subscript[k, x]^2 + 
        Subscript[k, z]^2 + 1/(4*H^2)) + (γ - 1)*g^2*Subscript[k, x]^2 == 0, ω]

What I want:

Plot[f, {k, 0, 20}]

where f is any solution of Solve. It plot nothing, because Mathematica didn't know about $k^2 = k_x^2 + k_z^2$ relation.

Answer 1

I'm pretty unsure about this, because I don't know much about manifolds and other difficult words. So I hope I learn something here too.

I use your equation with all parameters set to simple values as Chris did

eqn = w^4 - w^2 ((kx^2 + kz^2) + 1/4) + 2 kx^2 == 0;

Your second equation k^2==kx^2+kz^2 defines a tube with different radii k. We could for instance look at both equations and choose k==3

contPlot = ContourPlot3D[{w^4 - w^2 ((kx^2 + kz^2) + 1/4) + 2 kx^2 == 0, 
   kx^2 + kz^2 == 9}, {kx, -5, 5}, {w, -5, 5}, {kz, -5, 5}, 
  ContourStyle -> {Directive[Opacity[0.3], Red], Automatic}]

We are interested in the intersection of both contours. What we can do here, is to parametrize k^2==kx^2+kz^2 differently. A tube along w with radius k can expressed as the parametric equation

$$f(w,\phi)=\{k \cos(\phi), w,k \sin(\phi)\}$$

therefore we define a transformation rule

rule = Thread[{kx, w, kz} :> {k*Cos[phi], w, k*Sin[phi]}]

(* {kx :> k Cos[phi], w :> w, kz :> k Sin[phi]} *)

We can now apply this transformation and solve your initial equation for w. With this we get solutions for w which only depend on phi and k. With those solutions, we have an explicit parametrization of the curve in 3d

sol = Solve[eqn /. rule, w];
paramPlot = 
 ParametricPlot3D[{kx, w, kz} /. rule /. sol /. k :> 3, {phi, 0, 
    2 Pi}, PlotStyle -> Red] /. Line[pts_] :> Tube[pts, 0.1]

And we can of course combine them to see whether it fits with our imagination

Show[{contPlot, paramPlot}]

What you want to have now is a plot where on the first axis is k and on the second w

ParametricPlot[Evaluate[{k, w} /. sol], {k, 0, 9}, {phi, -Pi, Pi}]

Answer 2

Let's give ad hoc values to your physical parameters

eqn = ω^4 - ω^2 ((kx^2 + kz^2) + 1/4) + 2 kx^2 

Let's now apply the extra condition on the modulus of $k$

eqn= eqn /. kz -> Sqrt[k^2 - kx^2]; 

then we can find zeros of the dispersion relation at fixed kx

 Table[ContourPlot[
 eqn == 0 /. kx -> i/4 // Release, {ω, 0, 2}, {k, 0, 2}, 
 ContourStyle -> ColorData[10][i],
 RegionFunction -> Function[{k, ω}, k > i/4], 
 FrameLabel -> {k, ω}], {i, 1, 4}] // Show   

which produces

where I have added a cut to impose k>kx as suggested by Rahul Narain.

Alternatively, following more closely his strategy (i.e. proceed at fixed angle between 'kx' and 'ky')

 eqn = ω^4 - ω^2 ((kx^2 + kz^2) + 1/4) + 2 kx^2 /. 
 kx -> k Cos[θ] /. kz -> k Sin[θ]
 Table[ContourPlot[eqn == 0 /.θ-> Pi i/8 // Release, 
 {ω, 0, 2}, {k, 0,2}, ContourStyle -> ColorData[10][i], 
 FrameLabel -> {k, ω}], {i, 1, 8}] // Show

Answer 3

Using the same choice as @chris (and assuming you only want $\omega>0$, you can get the solution (in this case 2 solutions actually) as :

eqn[cs2_, oofh2_, gg_, kx_, kz_] = \[Omega]^4 - \[Omega]^2 cs2 ((kx^2 + kz^2) + oofh2) + gg kx^2 ;    

sol[cs2_?NumericQ, oofh2_?NumericQ, gg_?NumericQ, kx_?NumericQ, kz_?NumericQ] := 
  Solve[{eqn[cs2, oofh2, gg, kx, kz] == 0, \[Omega] >= 0}, {\[Omega]}, Reals][[All, 1, 2]]

Then you can simply plot it as a function of the parameters :

ParametricPlot[Thread[{Sqrt[kx^2 + kz^2], sol[1., 1/4., 2, kx, kz]}], {kx, -2, 2}, {kz, -2, 2}, 
          AxesLabel -> {" k ", " \[Omega] "}] 

However you get possibly better information with a more canonical plot :

Plot3D[sol[1., 1/4., 2, kx, kz], {kx, 0, 5}, {kz, 0, 5}, AxesLabel -> {" kx ", " kz ", " \[Omega] "}]