0
$\begingroup$

I have this (rather long) transcendental equation, that has multiple singularities where it goes to infinity, then reappears at negative infinity. When I try and find the roots of the equation via NSolve (or Reduce), it gives me the positions of the singularities as well as the zero crossings, which is incorrect:

enter image description here

However, when I look for a single root with FindRoot, it will only give me the position of a zero crossing, even if I give FindRoot a singularity as the search point. My question is - what command can I use to get the correct results I am getting from FindRoot, except for all roots within a certain domain, rather than one at a time? And my second question is, why are NSolve and Reduce giving me singularities as solutions when they are not solutions?

Here is my code - transcendental equation:

resonanceCondition = (1.0215026378769585`*^-7 Cos[
    0.038659981280734605` f] + 
  0.012261290032792185` Sqrt[f] Cos[0.038659981280734605` f] + 
  367.9364783156771` f Cos[0.038659981280734605` f] - 
  1.7385441873343195`*^-12 f^(3/2) Cos[0.038659981280734605` f] + 
  1.9168458821735513`*^11 f^2 Cos[0.038659981280734605` f] - 
  0.000905985363774135` f^(5/2) Cos[0.038659981280734605` f] - 
  27.186785669459642` f^3 Cos[0.038659981280734605` f] + 
  3.5778253945592175`*^-17 f^(7/2) Cos[0.038659981280734605` f] - 
  7.294985015701372`*^9 f^4 Cos[0.038659981280734605` f] + 
  1.349642838920257` f^(9/2) Cos[0.038659981280734605` f] + 
  0.0005593326523029028` f^5 Cos[0.038659981280734605` f] + 
  7.617136885646142`*^-25 f^(11/2) Cos[0.038659981280734605` f] + 
  145698.31404595848` f^6 Cos[0.038659981280734605` f] + 
  3.968314731945664`*^-16 f^(13/2) Cos[0.038659981280734605` f] + 
  1.190810871788677`*^-11 f^7 Cos[0.038659981280734605` f] + 
  0.003101895368576676` f^8 Cos[
    0.038659981280734605` f] + (1.0215026378769585`*^-7 + 
     0.012261290032792185` Sqrt[f] + 367.9364783156771` f + 
     1.73854418733432`*^-12 f^(3/2) + 
     1.916845882173552`*^11 f^2 + 
     0.0008919683885642353` f^(5/2) + 26.766164635165886` f^3 - 
     3.577825394559219`*^-17 f^(7/2) + 
     7.075853251101013`*^9 f^4 - 1.349642838920257` f^(9/2) - 
     0.0005593326523029029` f^5 - 
     7.617136885646143`*^-25 f^(11/2) - 145698.31404595854` f^6 - 
     3.968314731945666`*^-16 f^(13/2) - 
     1.190810871788677`*^-11 f^7 - 
     0.0031018953685766773` f^8) Cosh[
    2.8222`*^-11 + 1.6937701111797905`*^-6 Sqrt[f]] - 
  3.9263268389666276`*^-8 f Sin[0.038659981280734605` f] - 
  0.00471284462232609` f^(3/2) Sin[0.038659981280734605` f] - 
  141.42292112410644` f^2 Sin[0.038659981280734605` f] + 
  9.035967129830372`*^-7 f^(5/2) Sin[0.038659981280734605` f] - 
  7.422693884237215`*^10 f^3 Sin[0.038659981280734605` f] + 
  470.7485518941833` f^(7/2) Sin[0.038659981280734605` f] + 
  0.04187499717395471` f^4 Sin[0.038659981280734605` f] - 
  1.8614725761418453`*^-11 f^(9/2) Sin[0.038659981280734605` f] + 
  2.1815699347445037`*^7 f^5 Sin[0.038659981280734605` f] - 
  0.009697750162420625` f^(11/2) Sin[0.038659981280734605` f] + 
  8.2863`*^-10 f^6 Sin[0.038659981280734605` f] + 
  0.4316929951106164` f^7 Sin[
    0.038659981280734605` f] + (2.866240292375598`*^-17 + 
     5.160603931704003`*^-12 Sqrt[f] + 
     3.0971854210038807`*^-7 f + 0.006196017403536449` f^(3/2) + 
     53.38370833435124` f^2 + 3.203873913991683`*^6 f^(5/2) + 
     0.020624423363802637` f^3 - 
     1.8346263052726576`*^-6 f^(7/2) + 
     1.0744746292968111`*^7 f^4 - 955.7888618948019` f^(9/2) - 
     4.248775821083848`*^-7 f^5 - 
     3.6303916368586526`*^-14 f^(11/2) - 
     221.34930665701606` f^6 - 
     0.00001891332246165412` f^(13/2)) Sinh[
    2.8222`*^-11 + 
     1.6937701111797905`*^-6 Sqrt[
      f]])/(f (5.329080691237708`*^-19 + 
    6.396596693986087`*^-14 Sqrt[f] + 1.9194891031013206`*^-9 f + 
    1.` f^2) (-1.2730400062331583`*^9 Cos[
      0.038659981280734605` f] + 
    3.007633699724958`*^-6 Sqrt[f] Cos[0.038659981280734605` f] + 
    0.0902529952868547` f Cos[0.038659981280734605` f] + 
    4.7019279839115955`*^7 f^2 Cos[0.038659981280734605` f] + 
    6.396596693986088`*^-14 f^(5/2)
      Cos[0.038659981280734605` f] + 
    1.9194891031013215`*^-9 f^3 Cos[0.038659981280734605` f] + 
    1.` f^4 Cos[
      0.038659981280734605` f] + (-1.2730400062331583`*^9 - 
       3.007633699724958`*^-6 Sqrt[f] - 0.0902529952868547` f - 
       4.7019279839115955`*^7 f^2 - 
       6.396596693986088`*^-14 f^(5/2) - 
       1.9194891031013215`*^-9 f^3 - 
       1.0000000000000002` f^4) Cosh[
      2.8222`*^-11 + 1.6937701111797905`*^-6 Sqrt[f]] + 
    4.893155395016759`*^8 f Sin[0.038659981280734605` f] - 
    3.1263949972942227` f^(3/2)
      Sin[0.038659981280734605` f] + (-0.3572030481737306` - 
       21437.87990287647` Sqrt[f] - 71359.37236924547` f^2) Sinh[
      2.8222`*^-11 + 1.6937701111797905`*^-6 Sqrt[f]])) == 0;

Finding and plotting the roots with NSolve:

sol = NSolve[(resonanceCondition ) && 20 < f < 500]
p1 = Plot[((resonanceCondition) // First) /. f -> freq, {freq, 20, 
    500}];
pts = Table[{f /. sol[[j]], (resonanceCondition // First) /. 
     sol[[j]]}, {j, 1, Length[sol]}]; 
p2 = ListPlot[pts, PlotStyle -> Red];
Show[p1, p2]

Comparing results with FindRoot:

FindRoot[(resonanceCondition // First) , {f, 165}]
FindRoot[(resonanceCondition // First) , {f, 164}]
$\endgroup$
  • 4
    $\begingroup$ Consider Plot[((resonanceCondition) // First) /. f -> freq, {freq, 164., 164.4}]. Note the range of your plot and this one. $\endgroup$ – Michael E2 Apr 18 at 3:55
  • 1
    $\begingroup$ @MichaelE2 this means NSolve was correct, right? But OP said where it goes to infinity, then reappears at negative infinity. so I was wondering what is happening here. It looks OP then was wrong in this statement. $\endgroup$ – Nasser Apr 18 at 4:03
  • 2
    $\begingroup$ @Nasser Yes, I believe NSolve is correct. The OP's plot shows the graph going to ±7, not to infinity actually. I think the OP misinterpreted the plot. $\endgroup$ – Michael E2 Apr 18 at 4:06
  • 1
    $\begingroup$ @Jeremiah, since you're finding the roots of a ratio of two functions, have you considered just applying NSolve[]/FindRoot[] to the numerator only? $\endgroup$ – J. M.'s discontentment Apr 18 at 4:34
  • 2
    $\begingroup$ Nevertheless, it would surely be less of a computational burden to just consider the numerator, unless you have good reason to believe that the numerator and denominator have common factors. $\endgroup$ – J. M.'s discontentment Apr 18 at 4:36
2
$\begingroup$

This is just a long comment. Answer is given in comment above by @MichaelE2. The answer by NSolve is correct.

But if you want, you could always filter out solutions found, which are not as close to zero as you wanted as follows

sol = NSolve[resonanceCondition && 20 < f < 500]

Mathematica graphics

sel = If[Abs[ First[resonanceCondition]/. #] >10*$MachineEpsilon, False,True] & /@ sol;
sol = Pick[sol, sel]

Mathematica graphics

p1 = Plot[((resonanceCondition) // First) /. f -> freq, {freq, 20, 500}];
pts = Table[{f /. sol[[j]], (resonanceCondition // First) /. 
     sol[[j]]}, {j, 1, Length[sol]}];
p2 = ListPlot[pts, PlotStyle -> Red];
Show[p1, p2]

Mathematica graphics

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Just to add: it is possible to use the MeshFunctions capability of Plot[] to find initial approximations for roots, which can then be passed to FindRoot[] for further polishing, cf. this thread. $\endgroup$ – J. M.'s discontentment Apr 18 at 4:25
  • $\begingroup$ @Nasser I was just using sol[[1;;;;2]] to filter results. Turns out that approach was VERY incorrect. $\endgroup$ – Jeremiah Rose Apr 18 at 4:38
  • $\begingroup$ @JeremiahRose I would suggest filtering based on sensitivity: obj = resonanceCondition /. Equal -> Subtract; Abs[obj] < Abs[D[obj, f] * (f*$MachineEpsilon)] /. sol (accept a solution if True) where $\Delta f = $ f*$MachineEpsilon allows for an error in the last bit of the computed solution for f. Make $\Delta f$ larger if you want more tolerance. $\endgroup$ – Michael E2 Apr 18 at 15:35

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.