Help Please : Error in plotting the solutions of FindRoot

Question

I want to find the roots of an equation using FindRoot as a function of a real parameter K and then plot the real and imaginary solutions vs. the parameter K. The function I'm trying to find the roots is:

         f[\[CapitalOmega]_,K_]:=1.`\[VeryThinSpace]+(0.00021600988594924262` K^10+(3.81005811091083`*^-6-5.617765774643865`*^-9 I) K^2 \[CapitalOmega]^8-(1.938898556339015`*^-8-4.736384338194926`*^-11 I) \[CapitalOmega]^10+K^4 \[CapitalOmega]^6 ((-0.00024705957314623457`+1.7313407860214411`*^-7 I)+1.3234889800848443`*^-23 \[CapitalOmega]^2)+K^6 \[CapitalOmega]^4 ((0.005208736011404687`\[VeryThinSpace]+8.502422744684232`*^-8 I)+4.235164736271502`*^-22 \[CapitalOmega]^2)+K^8 \[CapitalOmega]^2 ((0.0021316316776359315`\[VeryThinSpace]+1.0578397620930897`*^-8 I)+6.776263578034403`*^-21 \[CapitalOmega]^2))/(2.138711742071709`*^-6 K^12+0.000021082195553663905` K^10 \[CapitalOmega]^2+0.000051246570978881814` K^8 \[CapitalOmega]^4-3.481280809505708`*^-6 K^6 \[CapitalOmega]^6+8.665879028798112`*^-8 K^4 \[CapitalOmega]^8-9.45611694948034`*^-10 K^2 \[CapitalOmega]^10+3.831913991301718`*^-12 \[CapitalOmega]^12)

The real solutions Wr :

          Subscript[W, r][K_,\[Mu]_]:=Re[\[CapitalOmega]/.FindRoot[1.+(0.00021600988594924262 K^10+(3.81005811091083*^-6-5.617765774643865*^-9 I) K^2 \[CapitalOmega]^8-(1.938898556339015*^-8-4.736384338194926*^-11 I) \[CapitalOmega]^10+K^4 \[CapitalOmega]^6 ((-0.00024705957314623457+1.7313407860214411*^-7 I)+1.3234889800848443*^-23 \[CapitalOmega]^2)+K^6 \[CapitalOmega]^4 ((0.005208736011404687+8.502422744684232*^-8 I)+4.235164736271502*^-22 \[CapitalOmega]^2)+K^8 \[CapitalOmega]^2 ((0.0021316316776359315+1.0578397620930897*^-8 I)+6.776263578034403*^-21 \[CapitalOmega]^2))/(2.138711742071709*^-6 K^12+0.000021082195553663905 K^10 \[CapitalOmega]^2+0.000051246570978881814 K^8 \[CapitalOmega]^4-3.481280809505708*^-6 K^6 \[CapitalOmega]^6+8.665879028798112*^-8 K^4 \[CapitalOmega]^8-9.45611694948034*^-10 K^2 \[CapitalOmega]^10+3.831913991301718*^-12 \[CapitalOmega]^12),{\[CapitalOmega],K/Sqrt[\[Mu] (1+K^2)]}]];

And the imaginary solutions Wi :

          Subscript[W, i][K_,\[Mu]_]:=Im[\[CapitalOmega]/.FindRoot[1.`\[VeryThinSpace]+(0.00021600988594924262` K^10+(3.81005811091083`*^-6-5.617765774643865`*^-9 I) K^2 \[CapitalOmega]^8-(1.938898556339015`*^-8-4.736384338194926`*^-11 I) \[CapitalOmega]^10+K^4 \[CapitalOmega]^6 ((-0.00024705957314623457`+1.7313407860214411`*^-7 I)+1.3234889800848443`*^-23 \[CapitalOmega]^2)+K^6 \[CapitalOmega]^4 ((0.005208736011404687`\[VeryThinSpace]+8.502422744684232`*^-8 I)+4.235164736271502`*^-22 \[CapitalOmega]^2)+K^8 \[CapitalOmega]^2 ((0.0021316316776359315`\[VeryThinSpace]+1.0578397620930897`*^-8 I)+6.776263578034403`*^-21 \[CapitalOmega]^2))/(2.138711742071709`*^-6 K^12+0.000021082195553663905` K^10 \[CapitalOmega]^2+0.000051246570978881814` K^8 \[CapitalOmega]^4-3.481280809505708`*^-6 K^6 \[CapitalOmega]^6+8.665879028798112`*^-8 K^4 \[CapitalOmega]^8-9.45611694948034`*^-10 K^2 \[CapitalOmega]^10+3.831913991301718`*^-12 \[CapitalOmega]^12),{\[CapitalOmega],K/Sqrt[\[Mu](1+K^2)]}]];

When I plotted the real part Wr I got the following numerical roundoff error :

        Block[{\[Mu]=100},Plot[{Subscript[W, r][K,\[Mu]]},{K,0,5},Frame->True,PlotRange->{0,0.2},PlotRangePadding->0,WorkingPrecision->100]]

And When I plotted the imaginary part Wi I got also the numerical roundoff error :

       Block[{\[Mu]=100},Plot[{-Subscript[W, i][K,\[Mu]] },{K,0,5},Frame->True,PlotRange->{0,0.02},PlotRangePadding->0,WorkingPrecision->100]]

Please, Can you give me now some hint here? I confess that I' m not very savvy with Mathematica and don' t really know which methods are implemented.

Answer 1

According to my tests, your function (where only even powers of omega occur) have symmetric solutions. Depending on the initial starting value for the root search, FindRoot converges to a solution or the opposite. That's why you observe the oscillations in your plot.

Solution : To prevent FindRoot from choosing "randomly" one solution or the opposite, you just have to specify a better initial starting value in FindRoot, closer to one side (or the other) of the symmetric solutions.

Try to see what solutions gives you FindRoot for particular values of K and other initial values.

Also, try to give two initial starting values to FindRoot instead of one, this way it will prevent FindRoot to use symbolic derivatives (see the detailed docs).

Also, use the option WorkingPrecision inside FindRoot (actually i did that with your function because FindRoot could not converge for some values).

This is the plot I get if i choose negative solutions.

Edit

Getting the roots with Solve/NSolve

As indicated in the comments by @DanielLichtblau, Solve/NSolve can be much more useful here because you can get directly the full set of solutions, and you can also try to impose some constraints to get only a part of these solutions (which is actually what you ask for in your comments.)

Let's see how to find the roots of your function f[omega,k], for example when k=1 and given your conditions : Real part of the solution must be >0 and Imaginary part must be negative. You just input :

sol = NSolve[f[om, 1.] == 0  && Re[om] > 0 && Im[om] < 0, om,
WorkingPrecision -> 16]

which returns 3 solutions

{{om -> 71.4852755969679 - 0.08779418975395 I},
 {om -> 8.57845490027755 - 0.51910299658873 I},
 {om -> 0.0000580288402574 - 0.449329051149141 I}}

You can check that the conditions (&& Re[om] > 0 && Im[om] < 0) are satisfied here (try to remove these conditions to see the full set of solutions which in your case are just the symmetric ones).
Now let's just check if these are really roots :

f[#, 1.] & /@ (om /. sol )

returns

{-1.21465*10^-7 - 4.29844*10^-8 I, 1.45325*10^-10 - 1.77549*10^-10 I,
8.13978*10^-9 + 5.19932*10^-9 I}

not bad, but you can of course impose some higher working precision for better results.

Now, if you want to automatize and for example plot the real part of the solutions for all k's between 0 and 5, you can do this :

resRe = Table[Thread@{k, Re[om /. 
NSolve[f[om, k] == 0  && Re[om] > 0 && Im[om] < 0 , om ]]},
{k, 0.01, 5, 0.4}]

returns a list of results for each of the 3 family of solutions and for different k's :

{{{0.01, 71.1328}, {0.01, 0.0857852}, {0.01, 5.80234*10^-7}},
{{0.41, 71.1923}, {0.41, 3.51719}, {0.41, 0.00002379}},
{{0.81, 71.3643}, {0.81, 6.94857}, {0.81, 0.0000470019}},
...

which you can plot

ListPlot[resRe // Transpose]

(To get lines instead of points :)

ListLinePlot[resRe // Transpose]

In the same way, to get the plot for the imaginary part :

resIm = Table[Thread@{k, Im[om /. 
NSolve[f[om, k] == 0  && Re[om] > 0 && Im[om] < 0 , om ]]},
{k, 0.01, 5, 0.4}];

then

ListPlot[resIm // Transpose]

And if you want to plot the Real and Imaginary part at the same time :

res = Table[
om /. NSolve[f[om, k] == 0  && Re[om] > 0 && Im[om] < 0 , 
  om ], {k, 0.01, 5, 0.4}] /. x_Complex -> {Re[x], Im[x]} ;
ListPlot[res // Transpose]