Finding all/most roots of a discontinuous function more consistently?

Question

I have the equation:

f[x_]:=α Tan[α*a] - Sqrt[c - α^2]

and ideally I want to find all of its positive zeros, given a and c, with variable alpha. The problem is that it's a discontinuous graph with large derivatives at the zeros. This is a plot of $\tan(\alpha a)$ vs $\sqrt{c-\alpha^2}/\alpha$, so that the intersection of the two graphs is a zero of the function:

Manipulate[
 Plot[{ Tan[α a], 
   Sqrt[c - α^2]/α}, {α, -0.1, 4}, PlotRange -> {-10, 10}]
 , {{a, 1.56}, -5, 5}, {{c, 10}, 0, 10}]

So, for that case there are five positive zeros (ignoring the single one at the end, I don't want to try to find that numerically).

To solve for them, I used FindRoot, and based off the properties of the tangent curve, I figure the best place to start would be somewhere well behaved in the tangent function, like when αa=[Pi](n+1/4) for integer n. So to get the roots I call FindRoot whenever this formula for alpha gives a value of Sqrt[c-alpha^2] that is real:

wellroots[a2_, c2_] :=
 Table[HoldForm[
     FindRoot[
      α Tan[α Rationalize[a] ] == 
       Sqrt[Rationalize[c] - α^2],
      {α, Pi *(n + 1/4)/(a)},
      WorkingPrecision -> 30
      ]
     ] /. {a -> a2, c -> c2, n -> n2},
   {n2, 0, Floor[Sqrt[c2] a2/Pi - 1/8]}] // ReleaseHold

(I needed the HoldForm and ReleaseHold commands to make sure FindRoot wasn't evaluated until the function was actually called. I used WorkingPrecision to try to smooth out the problems I was getting, and I used Rationalize so that it turned c into an exact value so I didn't get a "function precision less than working precision" error/warning from FindRoot. The Floor[...] term just puts a bound on n2.)

This works decently, and I can plot all approximated roots:

Manipulate[
 Show[
  Plot[{ Tan[α a], 
    Sqrt[c - α^2]/α}, {α, -0.1, 4}, 
   PlotRange -> {-1, 10}],
  Graphics[Join[{PointSize[0.02]},
    (Point[{α, Tan[α a]}]) /. wellroots[a, c]
    ]],
  Graphics[(Line[{{x, -5}, {x, 5}}])
    /. Table[{x -> Pi (n + 1/4)/(a)}, {n, 0, 
      Floor[Sqrt[c] a/Pi - 1/8]}]
   ]
  ]
 , {{a, 1.56}, -5, 5}, {{c, 10}, 0, 10}]

(the black lines are the initial approximations passed in to FindRoot, and the black circles are the approximated solutions.)

The problem is that as I drag the manipulate bars for a and c, sometimes the display of roots found flickers, showing that wellroots didn't find all possible roots. (and, the flickering isn't just a graphical thing, when I stop dragging the bar, some roots are still not shown)

So, I guess it's a bit much to ask to find all n roots with only n numerical attempts, but, how can I increase reliability so that I can expect to always (or, almost always) find all the roots?

And, only the positive roots matter, because the negative roots are just the negatives of all the positives.

(this question comes from "A First Course in Computational Physics" by DeVries, section "A little quantum mechanics problem", and is about time independent solutions to the schrodinger equation in a square well with finite side energies, V=c^2*hbar^2/(2m), where c is the value in the equation above. alpha is Sqrt[2 m E/hbar^2], where E is the energy in the schrodinger equation, same units as V. It's not homework/not for a course.)

Answer 1

Your equation:

$$x \tan (a x)-\sqrt{c-x^2}=0$$

implies

$$x^2 \tan ^2(a x)=c-x^2$$

So

$$x^2 \left(\tan ^2(a x)+1\right)=c$$

and

$$x^2 \sec ^2(a x)=c$$

Calling

$$\text{c1}=\sqrt{c}$$

we get

$$x^2 \sec ^2(a x)=\text{c1}^2$$

Taking square roots (remember both signs of c1should be considered):

$$\cos (a x)=\frac{x}{\text{c1}}$$

So all roots will have x <= c1

Now we could guess how many roots are there .... based on the period of the Cos[] and c1. Left as exercise :) (and remember to take in account both signs). Let's call that number n.
Now you can solve your problem by doing:

FindInstance[Cos[a x] - x/c1 == 0 && x > 0, {x}, Reals, n]

and the same for -c1, and the number of roots you expect for the negative part.

As an example, we take the values in your question:

Plot[{Cos[x 5], x/Sqrt@10, -x/Sqrt@10}, {x, 0, Sqrt@10}]

Column@N@FindInstance[Cos[5 x] - x/Sqrt@10 == 0 && x > 0, {x}, Reals, 5]

{x->0.295446}
{x->1.00732}
{x->1.47383}
{x->2.36841}
{x->2.63091}

Again, this was done for the positive c1 only

Answer 2

Mathematica doesn't yet have an approximate numerical solver for finding all/many roots as far as I know, but for 1D problems I once had very satisfactory results with Ted Ersek's RootSearch package: http://library.wolfram.com/infocenter/Demos/4482/.

It seems to work OK here:

Get["RootSearch`"]
Manipulate[
 Plot[{Tan[\[Alpha] a], 
   Sqrt[c - \[Alpha]^2]/\[Alpha]}, {\[Alpha], -0.1, 4}, 
  PlotRange -> {-10, 10}, 
  Epilog -> {PointSize[Large], Red, 
    Point[{#, Tan[# a]}] & /@ (\[Alpha] /. 
       Quiet@RootSearch[
         Tan[\[Alpha]*a] == 
          Sqrt[c - \[Alpha]^2]/\[Alpha], {\[Alpha], -0.1, 4}])}, 
  Exclusions -> Cos[\[Alpha] a] == 0, 
  ExclusionsStyle -> Dashed], {{a, 1.56}, -5, 5}, {{c, 10}, 0, 10}]

Answer 3

If you restrict the range of the solve variable, then in many cases NSolve will be able to compute all solutions within that range. For instance:

NSolve[Tan[α 5] == Sqrt[10 - α^2]/α && 0 < α < 10, α, Reals]

{{α -> 0.295446}, {α -> 0.885702}, {α -> 1.47383}, {α -> 2.05744}, {α -> 2.63091}, {α -> 3.15508}}

So, you can incorporate the intersection points by adding an epilog:

Manipulate[
    Plot[
        {Tan[α a], Sqrt[c-α^2]/α}, {α, -0.1, 4},
        Epilog->{
            PointSize[Large],
            Point[Values @ NSolve[z == Tan[α a] == Sqrt[c-α^2]/α && 0<α<c, {α,z}, Reals]]
        },
        PlotRange->{-10, 10}
    ],
    {{a, 5}, -5, 5},
    {{c, 10}, 0, 10}
]