1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving for complex phase speed eigenvalue

  1. Aug 9, 2011 #1
    Hey there, first timer poster here

    1. The problem statement, all variables and given/known data
    I'm working on a barotropic linear instability analysis and I've been having trouble getting an expression for the complex phase speed eigenvalue [itex]C = C_r + i*C_i[/itex] for the purpose of plotting a dispersion diagram ([itex]C_i[/itex] vs [itex]k[/itex] or [itex]C_r[/itex] vs [itex]k[/itex])


    2. Relevant equations

    The Normal mode solution:

    [itex]\psi(x,y,t) = \Psi(y)e^{ik(x - Ct)}[/itex]

    where k is the zonal wavenumber

    And the eigenvalue problem I'm trying to solve for C

    [itex]\frac{\partial^{2} \Psi}{\partial y^{2}} \Psi + (\frac{\beta - \frac{\partial^{2} U}{\partial y^2}}{U-C} -k^{2} )\Psi = 0[/itex]

    where U is a velocity profile thats given in gaussian form \e^{-y^2} and \beta is a given constant.

    3. The attempt at a solution

    The only things undetermined are C,k and y. The purpose is to find C for a given k and y. However I'm not too sure on the procedure in solving an eigenvalue expression like this so I tried to follow the Chebyshev method in Mathematica to solve for the eigenvalue C and I've ended up with a large output at the
    Eigenvalues[result] part of the code.
    that is too large to post full so here is a small part of it showing values and the general form:
    [itex]
    Root[......] \ , \ Root[-600. C^5 E^((3 y^2)/2000000000000) k^12 +
    150000. C^4 E^((3 y^2)/2500000000000) k^12 + ......] \ ,Root[....]
    [/itex]

    But I'm not sure whether I'm on the right track and if I am, how to go about plotting C vs k.

    Thanks

    Sorry if I'm not sounding clear
     
    Last edited: Aug 9, 2011
  2. jcsd
  3. Aug 17, 2011 #2
    Just thought I'd post the Mathematica code I'm working with:

    Code (Text):

    Uo = 100
    b = 6*10^(-11)
    R = 60268000
    r = R*Sin[76]
    U = Uo*Exp[-b*(y^2)/(2*Uo)]
    Ld = N[Sqrt[2*Uo/b]]
    \[CapitalOmega] = N[(2*\[Pi]/(39 + 60*39 + 10*60*60))]
    \[Beta] = 2*\[CapitalOmega]*Cos[(76) Degree]/R
    Clear[y]

    Uyy2 =  (b^{2} y^{2}/Uo)*Exp[-b*(y^2)/(2*Uo)] - b*Exp[-b*(y^2)/(2*Uo)]

    DE = \[Phi]''[
       y] + (- k^2 +  (\[Beta] - D[U, {y, 2}] )/(U - g)) \[Phi][y]


    t1 = DE /. {\[Phi][y] -> a[i]*ChebyshevT[i, y],
       D[\[Phi][y], {y, a1_}] -> D[a[i]*ChebyshevT[i, y], {y, a1}]}

    t2 = ExpandAll[t1]

    n = 12

    t3 = Sum[t2, {i, 0, n}]
    t4 = Simplify[t3]
    t5 = (t4 *(1/Sqrt[1 - y^2])*ChebyshevT[i, y] 2)/(\[Pi] c[i])
    t6 = Table[t5, {i, 0, n, 2}];
    t7 = (t6 /. c[0] -> 2) /. Table[c[i] -> 1, {i, 2, n}];
    t8 = Collect[Expand[Numerator[t7]], y];
    t9 = Transpose[Table[Coefficient[t8, y, i], {i, 0, 2 n, 2}]]
    t10 = Table[y^i/(\[Pi] Sqrt[1 - y^2]), {i, 0, 2 n, 2}]
    Timing[t11 = \!\(
    \*SubsuperscriptBox[\(\[Integral]\), \(-1\), \(1\)]\(t10 \
    \[DifferentialD]y\)\)]
    t12 = Simplify[Expand[t9.t11]]
    t13 = Transpose[Table[Coefficient[t12, a[i]], {i, 0, n, 2}]]
    b1 = Table[1, {i, 0, n, 2}]
    t14 = Append[Table[t13[[i]], {i, 1, n/2}], b1]
    s1 = Solve[t14[[n/2 + 1]].Table[a[i], {i, 0, n, 2}] == 0, a[n]]
    t15 = Table[t14[[i]].Table[a[i], {i, 0, n, 2}], {i, 1, n/2}]
    t16 = Expand[t15 /. s1]
    t17 = Transpose[
      Table[Table[Coefficient[t16[[1, i]], a[j]], {i, 1, n/2}], {j, 0,
        n - 2, 2}]]
    (*How would I go about geting the complex phase speed g from this?*)
    (*Like this? Am I approaching it wrong*)
    vals = Eigenvalues[-t17 /. g -> 0]

     
    I'm just not very sure on how to go about solving this complex eigenvalue problem numerically.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook