Find extreme values in a polynomial expression

Question

I have been hunting down and could not find an answer to my question. As I've now been struggling for quite a while and believe it is a simple trick that's preventing me from getting over this problem, I'm posting my question.

I have the following expression:

k^2*c^2 == omega^2*(1 + Sum[4 Pi kappaR[[i]]/(1 -omega^2/omegaR[[i]]^2), {i, 1, 3}])

with

kappaR = List[0.07141914, 0.03246304, 0.05539915]; 
omegaR =  List[0.125285/ℏ, 10.6661/ℏ, 18.1252/ℏ];

and

 ℏ = 6.5821192815/10^16;(* this is hbar/e *);
 c = 3*10^8;(* Velocity of light in vacuum, in m.s^-1 *)

If I graphically solve this polynomial expression, I get the following:

To do this I simply did

Plot[ℏ omega/.Sort@Quiet@N@Solve[ k^2*c^2 == 
  omega^2*(1 + Sum[4 Pi kappaR[[i]]/(1 -omega^2/omegaR[[i]]^2), {i, 1, 3}]),omega], 
  {k, -.55 10^9, .55 10^9}]

and played a little bit around with the PlotRange.

Anyway, we see that this relation has eight branches solution to it (two are very close to the horizontal axis and cannot be seen easily on the picture, that's because of the value of the poles that are an order of magnitude different). I would like to find the extreme values in omega of the third branch from the top:

It has the above behavior and omega varies between the first and second poles of the relation (i.e. omegaR[[1]] < omega < omegaR[[2]])

I tried the following for the maximum and got the result below

Input:

 FindMaximum[{ k^2*c^2 == omega^2*(1 + Sum[4 Pi kappaR[[i]]/(1 -omega^2/omegaR[[i]]^2), {i, 1, 3}]), 
   omegaR[[1]] <= omega <= omegaR[[2]]}, omega]

warning/error message:

FindMaximum::nrnum: The function value -(90000000000000000 k^2==-3.11568*10^42) is not a real number at {omega} = {1.90341*10^14}. >>

output:

FindMaximum[{ k^2*c^2 == omega^2*(1 + Sum[4 Pi kappaR[[i]]/(1 -omega^2/omegaR[[i]]^2), {i, 1, 3}]), 
  omegaR[[1]] <= omega <= omegaR[[2]]}, omega]

From here I must say I don't quite know how to proceed and any help and suggestions would be very welcome ;)

Thanks in advance guys

Answer 1

I found it easier to create explicit polynomials.

\[HBar] = 6.5821192815/10^16;(*this is hbar/e*);
c = 3*10^8;(*Velocity of light in vacuum,in m.s^-1*)

kappaR = {kp1, kp2, kp3};
omegaR = {or1, or2, or3};

poly = Numerator[
   Together[
    k^2*c^2 - (omega^2*(1 + 
         Sum[4 Pi kappaR[[i]]/(1 - omega^2/omegaR[[i]]^2), {i, 1, 
           3}]))]];

npoly = poly /. 
   Thread[kappaR -> {0.07141914, 0.03246304, 0.05539915}] /. 
  Thread[omegaR -> {0.125285/\[HBar], 10.6661/\[HBar], 
     18.1252/\[HBar]}]

From here we solve for the eight branches as Root functions.

solns = omega /. Solve[npoly == 0, omega, Quartics -> False];

Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

Some plotting showed that the fourth of the solutions corresponded to the one you want. I will remark that it does not seem to have a maximizer. Instead it looks to have a horizontal asymptote it approaches from below. It can be approximated though.

FindMaximum[\[HBar] solns[[4]], {k, 10^6}]

During evaluation of In[385]:= FindMaximum::sdprec: Line search unable to find a sufficient increase in the function value with MachinePrecision digit precision. >>

(* Out[385]= {10.6661, {k -> 1.07066362485*10^15}} *)