Message Boards Message Boards

Please Help : Error in plotting the solutions of FindRoot

Betatron Steve
Posted 10 years ago

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.

Attachments:
POSTED BY: Betatron Steve
1 Reply
Sort By:

One could do

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

In[76]:= (* The solutions W *)
Clear[W]
W[K_, \[Mu]_] := With[{preC = 100}, 
   FindRoot[SetPrecision[f[\[CapitalOmega], K], preC], {\[CapitalOmega], K/Sqrt[\[Mu]/(1 + K^2)]}, WorkingPrecision -> preC]
 ][[-1, -1]]

In[78]:= {Re[#], Im[#]} &[W[3.1, 100]]

During evaluation of In[78]:= FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than 100.` digits of working precision to meet these tolerances. >>

Out[78]=  {-1624.8869305883541028830136870262106382101189280581163075734345152591252489\80638534881282900723186087, 
34717.52852573805626771928431449226093876668088758334764286400786043682850613291079441186321008020593}

In[90]:= Off[FindRoot::lstol]
ListPlot[Table[{Re[#], Im[#]} &[W[x, 100]], {x, 0.1, 5, 0.001}], Frame -> True, PlotLabel -> "K\[Element]{0.1,5}  \[Mu]=100", 
 PlotRange -> All, PerformanceGoal -> "Quality", PlotMarkers -> {Automatic, 6}]
On[FindRoot::lstol]

to plot an Argand diagram of the roots found

enter image description here

One could also do a 

Off[FindRoot::lstol]
ParametricPlot[{Re[#], Im[#]} &[W[x, 100]], {x, 0.1, 5}, Frame -> True, PlotLabel -> "K\[Element]{0.1,5}  \[Mu]=100", 
 PlotRange -> All, PerformanceGoal -> "Quality"]
On[FindRoot::lstol]

to see some structure I'm unable to decipher. If one enters W[0,100] a divergency is produced, because K = 0 and the search of FindRoot starts at [CapitalOmega] = 0, so the denominator is really 0.
POSTED BY: Udo Krause
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract