# FindRoot errors

I'm trying to solve a system of two equations in two variables for some range of parameter values. I use ContourPlot to see where I should look for a solutions and then I use FindRoot to get the solutions. However, I am getting some "you need more than machine precision.." errors and sometimes when I rescale the function the error goes away. But I never thought to plug the solution back into the main expression to confirm that I get 0 until a few moments back and unfortunately I am not getting zeros.

My function is hyper sensitive to the values and the solutions I get from this FindRoot I want to use to compute other things later. So I am worried about picking up non-solutions. What are some good practices I can use in these conditions? How can I ensure that I do have a solution? Is ContourPlot reliable in terms of where to look?

Here is the exact problem I am solving:

My parameter is c0. I am looking to solve the system for "small" c0. The example below has c0=0.004 and I am looking to see how small I could push it between 0.01 > c0 >0.0001

The variables I am solving for are H0S and H0D.

 f1 = (-1.906273202485612 + 891.8237204611049 c0) H0D^1.7049946543060779 -    10.246354837464178 H0D^1.8549946543060778 - 0.3513303329336004 H0D^2.704994654306078 + (50.83395206628298 - 891.8237204611049 c0) H0S^1.7049946543060779 +   10.246354837464178 H0S^1.8549946543060778 - .0702660665867201 H0S^2.704994654306078;

f2 = (-1.093726797514388 + 511.6850514687195 c0)/H0D^2.9716613209727445 - 6.736078957579072/H0D^2.821661320972744 - 0.4820030003997329/H0D^1.9716613209727445 + (   29.166047933717014 - 511.6850514687195 c0)/H0S^2.9716613209727445 + 6.736078957579072/H0S^2.821661320972744 - 0.09640060007994662/H0S^1.9716613209727445;

c0 = 0.004
sub = 1088.96;
slb = 1088.94;
dub = 2.1774212*10^(-6);
dlb = 2.177421*10^(-6);

ContourPlot[{f1, f2}, {H0S, slb, sub}, {H0D, dlb, dub}]


This is what the ContourPlot looks like:

Then I try to find solutions:

{hs, hd, c0} = {H0S, H0D, c0} /. FindRoot[{f1 == 0, f2 == 0}, {H0S, (slb + sub)/2, slb,
sub}, {H0D, (dlb + dub)/2, dlb, dub}, MaxIterations -> Infinity]

{1088.95, 2.17742*10^-6, 0.004}

Scale1 = f1 /. {H0S -> 2000, H0D -> 0.000001};
Scale2 = f2 /. {H0S -> 2000, H0D -> 0.000001};

{hs, hd, c0} = {H0S, H0D, c0} /. FindRoot[{(f1)/Scale1 == 0, (f2)/Scale2 == 0}, {H0S, (slb + sub)/2,
slb, sub}, {H0D, (dlb + dub)/2, dlb, dub}, MaxIterations -> Infinity]

{1088.95, 2.17742*10^-6, 0.004}


When I do FindRoot without the scale, I get an error message about "..more than machine precision.." but when I do the FindRoot with the scale then I don't get any errors. But the solutions are the same.

However, when I do this:

  f1 /. {H0S -> hs, H0D -> hd}
f2 /. {H0S -> hs, H0D -> hd}
(f1/Scale1) /. {H0S -> hs, H0D -> hd}
(f2/Scale2) /. {H0S -> hs, H0D -> hd}


I get:

4.6912*10^-11
-15.7874
-1.80356*10^-18
-2.22429*10^-16


My f2=-15 at the proposed solution. So what's going on? Can I reliably use this solution going forward? How do I get it to give me the solution that makes f1 and f2 small..like zero... the solution ContourPlot seems to indicate that exists.

EDIT to include problems with RationalizeandRationalize[..,0]

fa=-1.90627 H0D^1.70499 + 52.3946 f H0D^1.70499 - 10.2464 H0D^1.85499 -0.35133 H0D^2.70499 + 11.1478 f H0S^1.70499 + 10.2464 H0S^1.85499 - 0.0702661 H0S^2.70499 + (0.523946 (9. + 100. f - 1. H0D) H0D^1.70499 H0S^1.70499 (877653. - 1. H0S^2.97166) (877653. - 0.000389055 H0S^4.67666))/((-877653. + H0S^2.97166) (877653. H0S^1.70499 - 1. H0S^4.67666 + H0D^4.67666 (1 - 0.000389055 H0S^1.70499) + H0D^1.70499 (-877653. + 0.000389055 H0S^4.67666))) + (0.523946 (9. + 100. f - 1. H0D) H0D^3.40999 (877653. -
1. H0S^2.97166) (-877653. + 0.000389055 H0S^4.67666))/((-877653. + H0S^2.97166) (877653. H0S^1.70499 - 1. H0S^4.67666 + H0D^4.67666 (1 - 0.000389055 H0S^1.70499) +
H0D^1.70499 (-877653. + 0.000389055 H0S^4.67666)))

fb=    (30.0615 f)/H0D^2.97166 - 1.09373/H0D^2.97166 + (-0.857219 -
0.280427 H0D^0.85)/H0D^2.82166 - 5.87886/H0D^2.82166 - \
0.201576/H0D^1.97166 + (6.39606 f)/H0S^2.97166 + (
0.857219 -
0.0560854 H0S^0.85)/H0S^2.82166 + 5.87886/H0S^2.82166 - \
0.0403152/H0S^1.97166 + (
0.300615 (9. + 100. f - 1. H0D) (877653. -
1. H0S^2.97166) (-877653. + 0.000389055 H0S^4.67666))/(
H0D^1.26667 (-877653. + H0S^2.97166) (877653. H0S^1.70499 -
1. H0S^4.67666 + H0D^4.67666 (1 - 0.000389055 H0S^1.70499) +
H0D^1.70499 (-877653. + 0.000389055 H0S^4.67666))) - (
0.300615 (9. + 100. f - 1. H0D) H0D^1.70499 (877653. -
1. H0S^2.97166) (-877653. + 0.000389055 H0S^4.67666))/(
H0S^2.97166 (-877653. + H0S^2.97166) (877653. H0S^1.70499 -
1. H0S^4.67666 + H0D^4.67666 (1 - 0.000389055 H0S^1.70499) +
H0D^1.70499 (-877653. + 0.000389055 H0S^4.67666)))

f = 0.3
ContourPlot[{fa == 0, fb == 0}, {H0S, 405, 408}, {H0D, 0.1, 3}]


According to ContourPlot the solution is near (406,1.50)

far = Rationalize[fa, 10^-16];
fbr = Rationalize[fa, 10^-16];
fa0 = Rationalize[fa, 0];
fb0 = Rationalize[fa, 0];
Clear[f]

In[21]:= {hsr, hdr} = {H0S, H0D} /.
FindRoot[{far == 0, fbr == 0} /. {f -> 0.3}, {H0S, 410, 360.5,
420}, {H0D, 1, 0.0001, 5}, MaxIterations -> Infinity]

Out[21]= {406.678, 0.316565}

In[434]:= {hs, hd} = {H0S, H0D} /.FindRoot[{fa == 0, fb == 0} /. {f -> 0.3}, {H0S, 400,    360.5, 420}, {H0D, 1, 0.0001, 5}, MaxIterations -> Infinity]

Out[434]= {406.43, 1.54421}

In[12]:= {hs0, hd0} = {H0S, H0D} /.
FindRoot[{fa0 == 0, fb0 == 0} /. {f -> 0.3}, {H0S, 400, 360.5,
420}, {H0D, 1, 0.0001, 5}, MaxIterations -> Infinity]

Out[12]= {406.175, 2.32183}


The solutions I get using rationalize are clearly incorrect if ContourPlot is to be believed. Now I will substitute these solutions in the expression below:

 exp=117.218 - 1754.39 f (0.047 + (0.00047 (9. + 100 f - H0D) (877653. -
H0S^2.97166) (-(877653./H0S^1.70499) + 0.000389055 H0S^2.97166))/(f (-877653. + H0S^2.97166) (877653./H0D^1.70499 - 0.000389055 H0D^2.97166 - 877653./H0S^1.70499 + H0D^2.97166/H0S^1.70499 + 0.000389055 H0S^2.97166 - H0S^2.97166/
H0D^1.70499))) + 877653. ((30.0615 f)/H0D^2.97166 - (0.336512 (0.833333 + 2.54737/H0D^0.85))/H0D^1.97166 - (0.122684 (2.97166 (3. + 16.9824 H0D^0.15 +          0.833333 H0D) H0D^1.97166 - (0.833333 + 2.54737/H0D^0.85) H0D^2.97166))/H0D^4.94332 + (
0.300615 (9. + 100 f - H0D) (877653. - H0S^2.97166) (-(877653./H0S^1.70499) +
0.000389055 H0S^2.97166))/(H0D^2.97166 (-877653. + H0S^2.97166) (877653./H0D^1.70499 - 0.000389055 H0D^2.97166 - 877653./H0S^1.70499 + H0D^2.97166/H0S^1.70499 + 0.000389055 H0S^2.97166 - H0S^2.97166/H0D^1.70499))) + 0.000389055 (52.3946 f H0D^1.70499 - (    0.213828 (2.97166 (3. + 16.9824 H0D^0.15 + 0.833333 H0D) H0D^1.97166 - (0.833333 + 2.54737/
H0D^0.85) H0D^2.97166))/H0D^0.266667 + (0.523946 (9. + 100 f - H0D) H0D^1.70499 (877653. - H0S^2.97166) (-(877653./H0S^1.70499) + 0.000389055 H0S^2.97166))/((-877653. + H0S^2.97166) (877653./H0D^1.70499 - 0.000389055 H0D^2.97166 - 877653./H0S^1.70499 +
H0D^2.97166/H0S^1.70499 + 0.000389055 H0S^2.97166 - H0S^2.97166/H0D^1.70499)))

In[22]:= exp /. {H0S -> hs, H0D -> hd, f -> 0.3}
exp /. {H0S -> hs0, H0D -> hd0, f -> 0.3}
exp /. {H0S -> hsr, H0D -> hdr, f -> 0.3}

Out[22]= 99.0462

Out[23]= -58748.5

Out[24]= 5.63457*10^7


99 is around where the answer should be whereas -58752 and 5*10^7 are just absurdly wrong. This is making me real hesitant about using Rationalize. At least in this instance I knew what the answer was gonna look like. For many other expressions I don't have a good idea and I could introduce errors and things could go wildly wrong as the errors compound.

Any help will be greatly appreciated.

• Does this question somewhat help you? Commented Nov 8, 2012 at 7:21
• If you had read the link and all the comments below carefully, you should have noticed that what you need isn't Rationalize[…, 10^-16] but Rationalize[…, 0] for your 2nd sample; SetPrecision[…, Infinity] will also work. BTW I really suggest you to clean up your sample a little (your 2nd sample is not even available for direct running! ) so it can attract more attention. Commented Nov 10, 2012 at 5:08
• @xzczd Thanks. I'll follow up on Rationalize[..,0]. Also, sorry about the example, I will correct it. I must've made a mistake after pasting it while I was pressing back space and delete to get the pasted stuff four spaces away from the margin. Commented Nov 10, 2012 at 10:44
• @xzczd I used Rationalize[..,0] in my second example and the solution I am getting is incorrect and is giving me absurd answers when I use it to compute other variables. It also contradicts what I see in ContourPlot. Commented Nov 10, 2012 at 11:13
• Oh, I forgot to mention that a higher WorkingPrecision (for example 20) in FindRoot is still needed, and the option MaxIterations isn't necessary. Commented Nov 10, 2012 at 11:53

Not a fully satisfactorily solution but more of a starting point. The problem is the function f2 because of the negative powers; if one stays away from zeros then you can try to factor them out.

(* Made the interval a bit wider *)
c0 = 0.004
sub = 1090;
slb = 1088;
dub = 2.178 10^(-6);
dlb = 2.177 10^(-6);

powerH0D = Cases[f2, Power[H0D, p_] -> p , 2];
powerH0S = Cases[f2, Power[H0S, p_] -> p , 2];

badFactor = Power[H0D, Min[powerH0D]] Power[H0S, Min[powerH0S]];

f2new = f2/badFactor // Expand ;

sol = FindRoot[{(f1 /. num_Real :> SetPrecision[num, \[Infinity]]) == 0,
(f2 /. num_Real :> SetPrecision[num, \[Infinity]]) == 0},
{H0S, (slb + sub)/2, slb, sub}, {H0D, (dlb + dub)/2, dlb, dub},
MaxIterations -> Infinity, WorkingPrecision -> 40];

{f1, f2new} /. sol
(* {4.6912*10^-11, 1.20205*10^-7} *)

{f1, f2} /. sol
(* {4.6912*10^-11, 0.212433} *)

• ok cool. so is the idea that f2(x)=good(x)*bad(x) then if bad(a)!=Infinity and good(a)=0 then a is a solution to f2 ? In my example, badFactor/.sol is around 10^16.. can we say that that's not (or not going to) Infinity? Commented Nov 8, 2012 at 8:53
• As long as neither of sol is 0 then yes, badFactor is finite (but possibly very large). Commented Nov 8, 2012 at 9:04
• ...phbbhtt.. my brains are fried. Ocfcourse.. I know what badFactor is.. literally ... I can't believe I asked the second part of the question. :) Commented Nov 8, 2012 at 9:08

Try Rationalize on your equations first

f1r = Rationalize[f1, 10^-16];
f2r = Rationalize[f2, 10^-16];


and then use a higher WorkingPrecision

fr = FindRoot[{f1r == 0, f2r == 0}, {H0S, (slb + sub)/2,slb,sub},
{H0D, (dlb + dub)/2, dlb, dub}, MaxIterations -> Infinity,
WorkingPrecision -> 32];


The results should be better:

f1r /. fr
(*0.*10^-25*)
f2r /. fr
(*0.*10^-15*)


Does that help?

• I just now realized this is basically the same as Daniel suggests in the answer @xzczd linked to... Commented Nov 8, 2012 at 8:03
• @sebhofer... that's pretty cool. I'm trying to read Daniel's answer but since it has details of the other question, it is not immediately clear to me why Rationalize reduces any numerical instability or other issues. If you can expand a bit on why Rationalize works then that would be awesome Thanks Commented Nov 8, 2012 at 8:33
• @sebhofer Try {f1, f2} /. fr though. Commented Nov 8, 2012 at 8:43
• @b.gatessucks whoa.. I got {-3.67838*10^-9, 208.209} when I tried {f1, f2} /. fr. Whats going on? Commented Nov 8, 2012 at 8:56
• 囧…after changing the Rationalize[…, 0] in sebhofer's code into SetPrecison[…, Infinity]`, I get the correct result as @b.gatessucks … I don't know the exact reason though… Commented Nov 8, 2012 at 13:23