# Roots of a characteristic equation. Can Mathematica solve it?

I have been using Mathematica for a few years, but with basic things. Now I need to find the first 10 roots of the following equation:

Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) == 0


I have been trying to use the commands Solve, SolveAlways and FindRoot. But the only one that worked is the last one, but it needs a value near the root. Could you help me?

• What's wrong with FindRoot[Cos[y] - 2 E^(-y)/(1 + E^(-2 y)), {y, 1}]? Mar 9, 2016 at 21:00
• Note that 2 E^(-y)/(1 + E^(-2 y)) is the same as Sech[y]. Mar 9, 2016 at 21:59
• In particular, since the hyperbolic secant goes to zero very rapidly as $y \to \infty$, the roots (except perhaps for the first few) are going to be very very close to half-integer multiples of $\pi$. Mar 9, 2016 at 22:12
• Mar 10, 2016 at 1:19

NDSolve[ ] allows a fine grained search for events through the interval by using the clause WhenEvent[ ].
Here we use it to find the roots and to stop searching when ten roots have been found.

Sow[ ]and Reap[ ]are used to collect our findings.

dy = D[Cos[x] - 2 E^(-x)/(1 + E^(-2 x)), x];
sol = Reap@ NDSolveValue[{y'[x] == dy, y == 0, a == 0,
WhenEvent[y[x] == 0, {a[x] -> a[x] + 1, Sow[x]}],
WhenEvent[a[x] == 10, {tf = x, "StopIntegration"}]},
y, {x, 0, 33}, DiscreteVariables -> {a}
];
Plot[sol[][t], {t, 0, tf},
Epilog -> {PointSize[Large], Red, Point[Thread[{sol[[2, 1]], 0}]]}] • Not a solution I would have thought of, but I suppose it does the trick! Also note that sol[[2,1]] / Pi yields {1.50562, 2.49975, 3.50001, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5} (see my comment on the original question.) Mar 9, 2016 at 22:13
• @MichaelSeifert Yep, you're right of course and that was my first approach. I changed it because Runny posted his answer a few seconds before me and I thought this one is a "classic" approach Mar 9, 2016 at 22:24

Exploiting Michael's observation in the comments, we can use the somewhat regular root spacing along with the Delves-Lyness method to locate the roots precisely:

Table[Re[NIntegrate[# (Sech[#] Tanh[#] - Sin[#])/(Cos[#] - Sech[#]) &[
π (k + 1/2) + π Exp[I t]/4] Exp[I t], {t, -π, π},
Method -> {"Trapezoidal", "SymbolicProcessing" -> 0}]]/8,
{k, 10}]
{4.730040744862704, 7.853204624095836, 10.995607838001671, 14.137165491257463,
17.27875965739948, 20.420352245626063, 23.56194490204045, 26.703537555508188,
29.845130209103246, 32.98672286269282}


If more digits are needed, these can then be fed to Root[] for refinement, like so:

N[Root[{Cos[#] - Sech[#] &, #}] & /@ %, 25]
{4.730040744862704026024048, 7.853204624095837556477067,
10.99560783800167090666903, 14.13716549125746417710592,
17.27875965739948143809107, 20.42035224562606109093641,
23.56194490204045507539202, 26.70353755550818624841941,
29.84513020910325426700149, 32.98672286269281956154720}


Using findNextRoot from my answer to How to find the next root larger than a specified value, numerically?...

I'm not sure if the "first ten roots" starts includes 0 or just positive roots. The function findNextRoot uses NDSolve, like Dr. belisarius's answer, and the stiffness at y == 0 is an issue for the "Projection" method that findNextRoot uses by default as well as for FindRoot.

findNextRoot[Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) == 0, {y, 0}, 10,
"NDSolveOptions" -> {Method -> Automatic}]
(*
{{y -> 4.73004}, {y -> 7.8532}, {y -> 10.9956}, {y -> 14.1372},
{y -> 17.2788}, {y -> 20.4204}, {y -> 23.5619}, {y -> 26.7035},
{y -> 29.8451}, {y -> 32.9867}}
*)


Including zero:

findNextRoot[Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) == 0, {y, -1/10}, 10,
"NDSolveOptions" -> {Method -> Automatic}]
(*
{{y -> 0.0000245925}, {y -> 4.73004}, {y -> 7.8532},
{y -> 10.9956}, {y -> 14.1372}, {y -> 17.2788}, {y -> 20.4204},
{y -> 23.5619}, {y -> 26.7035}, {y -> 29.8451}}
*)


Should one feel that the root at y == 0 is not accurate enough, observe that Cos is 1, so at machine precision a residual of \$MachineEpsilon is pretty good:

Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) /. %
(*
{1.11022*10^-16, 1.41603*10^-9, 4.49942*10^-10, -1.46555*10^-7, -9.86477*10^-8,
-4.52107*10^-8, 5.93708*10^-7, 2.04556*10^-7, -2.59637*10^-10, 4.35094*10^-9}
*)


The residuals at the other roots are explained by the default AccuracyGoal of FindRoot being 8. Compare:

findNextRoot[Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) == 0,
{y, -1/10}, 10,
"NDSolveOptions" -> {Method -> Automatic},
"PostProcess" -> (FindRoot[##, AccuracyGoal -> 12] &)]
Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) /. %
(*
{{y -> 0.0000589511}, {y -> 4.73004}, {y -> 7.8532}, {y -> 10.9956},
{y -> 14.1372}, {y -> 17.2788}, {y -> 20.4204}, {y -> 23.5619},
{y -> 26.7035}, {y -> 29.8451}}

{0., -2.67147*10^-16, -4.14382*10^-16, 3.90326*10^-16, -2.89806*10^-16,
-1.71698*10^-15, 1.76219*10^-15, 1.43247*10^-16, -1.50236*10^-15, -1.17313*10^-15}
*)


Basically the question is about using Solve as the most straightforward tool for finding exact solutions however all the other existing answers demonstrate numerical techniques and find approximate solutions. It was observed in the comments that 2 E^(-y)/(1 + E^(-2 y)) == Sech[y] == 1/Cosh[y]. As Cosh[y] is an exponentially increasing function it is obvious that there are infinitely many solutions of the given equations very close to those of this trivial equation Cos[y] == 0. Solve will not give any solutions unless it is suplemented by an appropriate codition restricting the space of solutions to a finite number of them. There were many asnwers discussing this issue and it is easy to find them on this site. The exact solutions of transcendental equations in Mathematica are given in terms of the Root objects since the version 7, see e.g. this post by Roger Germundsson Mathematica 7, Johannes Kepler, and Transcendental Roots.

Since we would like to find the first 10 solutions and having said the remarks above we conclude that the condition supplementing the equation should be roughly 0 <= y <= 10 Pi but to be ensured one can put 0 <= y <= 11 Pi.

Now we can use Solve as well as Reduce, e.g.

sols = DeleteDuplicates[ y /. Solve[ Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) == 0
&& 0 <= y <= 35, y]];


We have used DeleteDuplicates because Solve produces four copies of the same solution for y == 0. In fact we have obtained all the desired solutions:

Length @ sols

 11


All the above is true if we asume we need only real roots. However if we are to find all complex solutions we need to do it this way:

solComp = ReIm[y]/.{ ToRules @ Reduce[ Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) == 0
&& 0 <= Abs[y] <= 12, y]};

ContourPlot[
{ Re[Cos[x + I y] - 2 E^(-x - I y)/(1 + E^(-2 (x + I y)))],
Im[Cos[x + I y] - 2 E^(-x - I y)/(1 + E^(-2 (x + I y)))]},
{x, -11, 11}, {y, -11, 11},
PlotPoints -> 75, MaxRecursion -> 5,
Epilog -> {Table[{Dashed, Gray, Circle[{0, 0}, i]}, {i, sols[[2;;4]]}],
Red, PointSize[0.015], Point[solComp]},
PlotLegends -> Placed["Expressions", Bottom]] This plot sheds some light why Solve finds four copies of trivial solution z == x + I y == 0. Moreover we can find that on the both axes roots are in equal distances from the origin exchanging x and y because the equation is equivalent to Cos[x + I y] == 1/Cosh[x + I y] and {Cos[I y], Cosh[I y]} yields {Cosh[y], Cos[y]} i.e. they exchange their roles. However on the level of the given function 1/Cosh[x + I y] becomes 1/Cos[y] for x == 0 and this is why there occurs a singularity when Cos[y] == 0.

Reduce[Cosh[x + I y] == 0 && -(1/4) < x < 1/4 && 5/4 < y < 7/4, {x, y}]

   x == 0 && y == Pi/2


This is a the first singular point of Sech[ x + I y].

In case of more sophisticated equations one should plot given functions to estimate the range 0 <= y <= yMax. Let's plot the first 10 real and positive roots:

Plot[ Cos[y] - 2 E^(-y)/(1 + E^(-2 y)), {y, 0, 35},
Epilog -> {Red, PointSize[0.02], Point[Tuples[{sols, {0}}]]}] Approximate values of y one finds with

N @ sols

{ 0., 4.73004, 7.8532, 10.9956, 14.1372, 17.2788, 20.4204,
23.5619, 26.7035, 29.8451, 32.9867}


and as I pointed out earlier differences between consecutive roots tend to Pi:

Differences @ %

 { 4.73004, 3.12316, 3.1424, 3.14156, 3.14159, 3.14159, 3.14159,
3.14159, 3.14159, 3.14159}


quite rapidly, what give us an estimation for the condition 0 <= y <= yMax to find much more exact solutions.

Here we demonstrate the structure of contours {Re[f[z]==0, Im[f[z]==0} near the first singular point on the imaginary axis.

cp = ContourPlot[
{ Re[Cos[x + I y] - 2 E^(-x - I y)/(1 + E^(-2 (x + I y)))] == 0,
Im[Cos[x + I y] - 2 E^(-x - I y)/(1 + E^(-2 (x + I y)))] == 0},
{x, -2.2, 2.2}, {y, -2.2, 2.2}, ContourStyle -> {Orange, Green},
Epilog -> {Red, PointSize[0.015], Point[solComp]}, PlotPoints -> 75,
MaxRecursion -> 5];

GraphicsRow[
Table[ Show[
ContourPlot[{ f[Cos[x + I y] - 2 E^(-x - I y)/(1 + E^(-2 (x + I y)))]},
{x, -2.2, 2.2}, {y, -2.2, 2.2}, PlotPoints -> 100, MaxRecursion -> 5,
ColorFunction -> ColorData["DeepSeaColors"]], cp], {f, {Re, Im}}]] Althought the green and orange lines intersect, there is no root since a singular point does not belong to the domain of the underlying function.

How about NSolve?

NSolve[Cos[y] - 2 E^(-y)/(1 + E^(-2 y)) == 0 && 0 < y < 11 Pi, y]
(*
{{y -> 4.73004}, {y -> 7.8532}, {y -> 10.9956}, {y -> 14.1372},
{y -> 17.2788}, {y -> 20.4204}, {y -> 23.5619}, {y -> 26.7035},
{y -> 29.8451}, {y -> 32.9867}}
*)


The following should give you about 10 roots of the equation:

FindRoot[Cos[y] - 2 E^(-y)/(1 + E^(-2 y)), {y, #}] & /@ Range[0, 30, 3]

(* {{y -> 0.}, {y -> -23.5619}, {y -> 4.73004}, {y -> 7.8532}, {y ->
10.9956}, {y -> 14.1372}, {y -> 17.2788}, {y -> 20.4204}, {y ->
23.5619}, {y -> 26.7035}, {y -> 29.8451}} *)