Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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[\[Omega]^4 - \[Omega]^2*(Subscript[k, x]^2 + 
        Subscript[k, z]^2 + 1/(4*H^2)) + (\[Gamma] - 1)*g^2*Subscript[k, x]^2 == 0, \[Omega]]

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.

share|improve this question
Could you please share your Mathematica code and specify what are the problems you get ? – b.gatessucks 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.gatessucks 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
up vote 9 down vote accepted

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}]

Mathematica graphics

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]

Mathematica graphics

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

Show[{contPlot, paramPlot}]

Mathematica graphics

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}]

Mathematica graphics

share|improve this answer
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

 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

Mathematica graphics

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

Mathematica graphics

share|improve this answer

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] "}] 

enter image description here

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] "}]

enter image description here

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.