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


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


 ℏ = 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: enter image description here

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: enter image description here

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


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


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

  • $\begingroup$ You will need to Maximize with respect to k, not omega. $\endgroup$ Oct 15, 2014 at 22:27

1 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[
    k^2*c^2 - (omega^2*(1 + 
         Sum[4 Pi kappaR[[i]]/(1 - omega^2/omegaR[[i]]^2), {i, 1, 

npoly = poly /. 
   Thread[kappaR -> {0.07141914, 0.03246304, 0.05539915}] /. 
  Thread[omegaR -> {0.125285/\[HBar], 10.6661/\[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}} *)
  • $\begingroup$ Hi @Daniel, thank you very much for this answer. It really looks like exactly what I needed :) $\endgroup$
    – MXJ
    Oct 16, 2014 at 9:28
  • $\begingroup$ Hi @DanielLichtblau, One comment: I find interesting and "unpredictable" the way in which Solve sorts its solutions. Usually to circumvent this I would use Sort[solns]... If I do so the branch of interest is #6 (the solutions are then obviously sorted in growing order of omega). But then executing FindMaximum[\[HBar] Sort[solns][[6]], {k, 10^6}] returns {11.6345, {k -> 1.1*10^7}}, which is obviously wrong! I wonder why that is and shall investigate this further. $\endgroup$
    – MXJ
    Oct 16, 2014 at 9:34
  • $\begingroup$ That 6th solution might be the branch above the one of interest? They do come kind of close as you go far in eaither direction. As for solution ordering, I think it will be by root branch number, with the two quartics (if I recall correctly) themselves selected in an unknown order. "Branch number" is actually a misnomer in that it really means "root order", but only gets resolved upon substituting actual values for the parameters. The branches can and do jump, as in the case of the one of interest (notice the cusp, where two algebraic function branches have crossed). $\endgroup$ Oct 16, 2014 at 13:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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