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

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

-

$$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

-
For the number of roots, a good estimate would be twice the number of whole periods until you get to where x/c1=1, plus one. (so, largest k for which $2k\pi<c1$). But, for example, if c1 is a bit larger than 1/(2 pi) and a=1, then you could have x be only a bit less than 2 pi and calculate 1 root when there are actually 3. But you can plug in this number +2, and, FindInstance doesn't seem to mind it if it can't find the requested number of solutions. (just wondering if you had something else in mind) –  NeuroFuzzy Feb 22 '13 at 10:29
Awesome! I'm not sure how I missed this... simply transforming the equation so that it didn't have tangent in it would have done. "FindInstance" is a bit of a magic function in this though. (oh, and in the previous comment, I was only considering c1 with positive sign) –  NeuroFuzzy Feb 22 '13 at 10:35
Although, FindInstance seems to be a lot slower than FindRoot, so I'll stick to the wellroots function in my post. Plus there are a lot of extraneous roots caused by squaring, in this solution. But, the transforming the equation to get rid of the nasty 1/cos(ax) is what I needed. Thanks! –  NeuroFuzzy Feb 22 '13 at 11:08
Isn't it cheating to use real math in responses? (Okay, I did in fact give this a vote.) –  Daniel Lichtblau Feb 22 '13 at 16:09
@DanielLichtblau I think they are valid only if you don't need a pencil to work them out :) –  belisarius Feb 22 '13 at 16:19

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

-