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:

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

graph of b tan(b*a) and sqrt(c-b^2)

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_] :=
      α 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:

  Plot[{ Tan[α a], 
    Sqrt[c - α^2]/α}, {α, -0.1, 4}, 
   PlotRange -> {-1, 10}],
    (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}]

Graph of zeros with starting approximations and calculated roots

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


3 Answers 3


Your equation:

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


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


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


$$x^2 \sec ^2(a x)=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}]

Mathematica graphics

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


Again, this was done for the positive c1 only

  • $\begingroup$ 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) $\endgroup$
    – user273
    Feb 22, 2013 at 10:29
  • $\begingroup$ 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) $\endgroup$
    – user273
    Feb 22, 2013 at 10:35
  • $\begingroup$ 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! $\endgroup$
    – user273
    Feb 22, 2013 at 11:08
  • 2
    $\begingroup$ Isn't it cheating to use real math in responses? (Okay, I did in fact give this a vote.) $\endgroup$ Feb 22, 2013 at 16:09
  • $\begingroup$ @DanielLichtblau I think they are valid only if you don't need a pencil to work them out :) $\endgroup$ Feb 22, 2013 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:

 Plot[{Tan[\[Alpha] a], 
   Sqrt[c - \[Alpha]^2]/\[Alpha]}, {\[Alpha], -0.1, 4}, 
  PlotRange -> {-10, 10}, 
  Epilog -> {PointSize[Large], Red, 
    Point[{#, Tan[# a]}] & /@ (\[Alpha] /. 
         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}]

enter image description here


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:

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

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.