Numerically solving an algebraic equation

Question

I have asked a similar question earlier, but I am still stuck at the same problem. I would like to solve an algebraic equation of the form

Tan[d*Sqrt[-k^2 + w^2*(1 + (1 - w^2)^(-1))]] == 
Sqrt[k^2 - w^2]/Sqrt[-k^2 + w^2*(1 + (1 - w^2)^(-1))]

and plot the result w as a function of k. I am particularly interested in two limiting cases: small d and large d. For small d, the Tangent function can be approximated by its argument to lowest order and the resulting equation:

d*Sqrt[-k^2 + w^2*(1 + (1 - w^2)^(-1))] == 
Sqrt[k^2 - w^2]/Sqrt[-k^2 + w^2*(1 + (1 - w^2)^(-1))]

can be solved analytically. For larger d, this approximation does not work, and I need something else.

What I have tried so far:

1) Taking the derivative on both sides with respect to one of the variables while keeping the other as a function of this variable, and then use NDSolve[] in order to solve the differential equation.

2) Trying to substitute both variables to get an expression that can be separated.

3) Try out the answers from my previous question applied to this modified equation.

4) One comment on the previous question suggested to use methods from Solve symbolically a transcendental trigonometric equation and plot its solutions, which I tried to implement.

None of these have worked so far. I can see that for d = 0, there are exactly two solutions (two branches), so I would expect that these two solutions will remain but change with bigger d. Unfortunately, every attempt to solve this equation numerically either ends up with no solution or something that is clearly wrong (discontinuous curve, values in the $10^{30}$, etc.). For instance, the solution approach described by the first answer in my previous question would give something like

Can someone help me find a solution for the two cases d = 0.1 and d = 100, in the range $0 \leq k \leq 10$?

Answer 1

You can do series expansion around d=0 to see what' going on.

feq = Tan[d*Sqrt[-k^2 + w^2*(1 + (1 - w^2)^(-1))]] - 
        Sqrt[k^2 - w^2]/Sqrt[-k^2 + w^2*(1 + (1 - w^2)^(-1))];

ce0[k_, w_, d_] = 
   ComplexExpand[feq, TargetFunctions -> {Re, Im}] // 
      Simplify[#, 0 < w < 1 && 0 < k < 10] &;

ser1 = Series[feq, {d, 0, 1}] // Normal

(*   -(Sqrt[k^2 - w^2]/Sqrt[-k^2 + w^2 (1 + 1/(1 - w^2))]) + 
      d Sqrt[-k^2 + w^2 (1 + 1/(1 - w^2))]   *)

ser3 = Series[feq, {d, 0, 3}] // Normal;

wsol1[k_, d_] = Solve[0 < d && 0 < k < 10 && ser1 == 0, w, Reals]

(*   {{w -> ConditionalExpression[-\[Sqrt]Root[-k^2 + 
     d^2 k^4 + (1 + 2 k^2 - 4 d^2 k^2 - 2 d^2 k^4) #1 + (-2 + 
        4 d^2 - k^2 + 6 d^2 k^2 + d^2 k^4) #1^2 + (1 - 4 d^2 - 
        2 d^2 k^2) #1^3 + d^2 #1^4 &, 2], 
0 < k < 10 && d > 0]}, 
{w -> 
ConditionalExpression[\[Sqrt]Root[-k^2 + 
    d^2 k^4 + (1 + 2 k^2 - 4 d^2 k^2 - 2 d^2 k^4) #1 + (-2 + 
       4 d^2 - k^2 + 6 d^2 k^2 + d^2 k^4) #1^2 + (1 - 4 d^2 - 
       2 d^2 k^2) #1^3 + d^2 #1^4 &, 2], 0 < k < 10 && d > 0]}}   *)

Manipulate[
   Plot[Evaluate[(w /. wsol1[k, d])], {k, 0, 10}, 
     PlotRange -> All], {{d, 10^-3}, 0, 100}]

Solving to higher series order than 3 doesn't change very much.

wsol3[k_, d_] = Solve[0 < d && 0 < k < 10 && ser3 == 0, w, Reals]

Manipulate[
  Plot[Evaluate[{(w /. wsol1[k, d]) - (w /. wsol3[k, d])}], {k, 0, 10},
    PlotRange -> All], {{d, 10^-3}, 0, 100}]

Make a ContourPlot to see, the greater d, the more solutions appear, all starting at k==w and go as asymptotes to the series solution. (The dots in the lower right corner are artefacts due to WorkingPrecision).

ContourPlot[
  Evaluate[{k == w, ce0[k, w, 1/2] == 0, w == (w /. wsol1[k, 1/2]), 
  w == (w /. wsol3[k, 1/2])}], {k, .8, 2}, {w, .8, 1}, 
   PlotPoints -> 200, GridLines -> Automatic, 
   ContourStyle -> {Red, Blue, Green, Magenta}]

For given d and k, you get infinite many solutions near w->1.

NSolve[0 < w < 99999/100000 && ce0[15/10, w, 1/2] == 0, w, Reals]

(*   {{w -> 0.881366}, {w -> 0.989096}, {w -> 0.996958}, {w -> 
0.998618}, 
{w -> 0.999216}, {w -> 0.999497}, {w -> 0.99965}, {w -> 
0.999742}, {w -> 0.999803}, {w -> 0.999844}, {w -> 0.999874}, {w ->
0.999895}, {w -> 0.999912}, {w -> 0.999925}, {w -> 
0.999935}, {w -> 0.999944}, {w -> 0.999951}, {w -> 0.999956}, {w ->
0.999961}, {w -> 0.999965}, {w -> 0.999968}, {w -> 
0.999971}, {w -> 0.999974}, {w -> 0.999976}, {w -> 0.999978}, {w ->
0.99998}, {w -> 0.999981}, {w -> 0.999983}, {w -> 0.999984}, {w ->
0.999985}, {w -> 0.999986}, {w -> 0.999987}, {w -> 
0.999988}, {w -> 0.999989}}   *)

Edit

Let me note, that the curve beginning at about w=.95 and every second curve above is arificial due to branch cut. (Have no time to correct it now.)

Answer 2

More an exploration than an answer. Take your equation, multiply left and right side by the denominator of the right side. Call out each side. Left side...

left = Sqrt[-k^2 + w^2 (1 + 1/(1 - w^2))] Tan[d Sqrt[-k^2 + w^2 (1 + 1/(1 - w^2))]];

Right side

right = Sqrt[k^2 - w^2];

Both are complex numbers for values of d and range of k you are considering. Look at the norm of left-right for value d = 0.1.

norm = (Norm@(left - right)) /. d -> .1;

Manipulate[Plot[(norm /. k -> kk), {w, 0, 20}, 
          PlotRange -> {Automatic, {-10, 20}}], {kk, 0, 10}]

Eyeballing it, there is a singularity at w=1 always, and a minimum at w=k which never reaches zero, a necessary (but not sufficient) condition for your equation to hold.