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

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.

• Could you please share your Mathematica code and specify what are the problems you get ? – b.gates.you.know.what Oct 29 '12 at 9:51
• And how should I share code, when copying from Mathematica I get plenty of some strange additions like \FractioBox and \SuperscriptBox? – m0nhawk Oct 29 '12 at 10:04
• If you start with the equation I'll try to edit it so you can see how it's done. – b.gates.you.know.what Oct 29 '12 at 10:06
• so you want say contours of iso kx in the (omega,k) plane which are solution to F(omega,k,kx)=0 above? – chris Oct 29 '12 at 10:16
• I've updated the question, if it isn't clear enough now - I will update further. – m0nhawk Oct 29 '12 at 10:20

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}] • Its only a curve because I fixed arbitrarily one degree of freedom (I think anyway...) – chris Oct 29 '12 at 15:45
• and stop playing the trick of throwing lots of nice looking 3D curves: I invented it ;o) (just kidding nice plots and good answer!) – chris Oct 29 '12 at 15:47
• @chris Thanks. I see now that you set one dimension. I was just looking at your images and started trying it by myself. – halirutan Oct 29 '12 at 15:53

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[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[i],
FrameLabel -> {k, ω}], {i, 1, 8}] // Show 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] "}] 