0
$\begingroup$

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

enter image description here

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ In fact, the initial numerical values for the solutions will depend on the values of K, findroot starts calculating the roots from K = 0, ie Wr = 0 and Wi = 0, by varying k from 0 to 5. $\endgroup$ – Betatron Nov 18 '14 at 11:11
  • 1
    $\begingroup$ If you use NSolve instead of FindRoot (changing code as necessary) you get the full set of solutions. THis can be useful if trying to track a particular algebraic function or trying to impose specific conditions e.g. 0<root<1. $\endgroup$ – Daniel Lichtblau Nov 18 '14 at 14:39
2
$\begingroup$

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.

enter image description here

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]

enter image description here

(To get lines instead of points :)

ListLinePlot[resRe // Transpose]

enter image description here

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]

enter image description here

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]

enter image description here

$\endgroup$
  • $\begingroup$ Please, How to select only the solutions for Wi<0 and Wr>>Wi ? $\endgroup$ – Betatron Nov 18 '14 at 22:09
  • $\begingroup$ I want to plot only the negatif Imaginary solutions Wi(K)<0 and positif Real solutions Wr(K)>0 Please, how do it? $\endgroup$ – Betatron Dec 11 '14 at 23:14
  • $\begingroup$ @betatron I remember another answer to your post which explained how to use Solve/NSolve which is much more adpated here. It seems it has been deleted ... ?? I will take a look at your problem. $\endgroup$ – SquareOne Dec 12 '14 at 13:37
  • $\begingroup$ @betatron, sorry for the delay, i have what you want and i'll post it in a few hours. Just to be sure : I understand that you want Wi(K)<0 and Wr(k)>0 at the same time ? $\endgroup$ – SquareOne Dec 13 '14 at 15:09
  • $\begingroup$ Yes @SquareOne, that's what I want. Thank you. $\endgroup$ – Betatron Dec 13 '14 at 15:29

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.