0
$\begingroup$

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?

|improve this question
$\endgroup$
  • $\begingroup$ What's wrong with FindRoot[Cos[y] - 2 E^(-y)/(1 + E^(-2 y)), {y, 1}]? $\endgroup$ – Sjoerd C. de Vries Mar 9 '16 at 21:00
  • 1
    $\begingroup$ Note that 2 E^(-y)/(1 + E^(-2 y)) is the same as Sech[y]. $\endgroup$ – Michael Seifert Mar 9 '16 at 21:59
  • 1
    $\begingroup$ 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$. $\endgroup$ – Michael Seifert Mar 9 '16 at 22:12
  • $\begingroup$ Related: mathematica.stackexchange.com/questions/109012/… $\endgroup$ – Michael E2 Mar 10 '16 at 1:19
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Dr. belisarius Mar 10 '16 at 16:22
8
$\begingroup$

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

Mathematica graphics

|improve this answer
$\endgroup$
  • $\begingroup$ 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.) $\endgroup$ – Michael Seifert Mar 9 '16 at 22:13
  • $\begingroup$ @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 $\endgroup$ – Dr. belisarius Mar 9 '16 at 22:24
7
$\begingroup$

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}
|improve this answer
$\endgroup$
4
$\begingroup$

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}
*)
|improve this answer
$\endgroup$
4
$\begingroup$

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

enter image description here

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

enter image description here

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

enter image description here

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.

|improve this answer
$\endgroup$
3
$\begingroup$

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}}
*)
|improve this answer
$\endgroup$
2
$\begingroup$

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}} *)
|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.