Roots of a characteristic equation. Can Mathematica solve it?

Question

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?

Answer 1

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] == 0, a[0] == 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[[1]][t], {t, 0, tf}, 
     Epilog -> {PointSize[Large], Red, Point[Thread[{sol[[2, 1]], 0}]]}]

Answer 2

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}

Answer 3

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[0] 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}
*)

Answer 4

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.

Answer 5

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}}
*)

Answer 6

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}} *)