2
$\begingroup$

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:

enter image description here

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

enter image description here

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.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Does this question somewhat help you? $\endgroup$ – xzczd Nov 8 '12 at 7:21
  • $\begingroup$ 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. $\endgroup$ – xzczd Nov 10 '12 at 5:08
  • $\begingroup$ @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. $\endgroup$ – Amatya Nov 10 '12 at 10:44
  • $\begingroup$ @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. $\endgroup$ – Amatya Nov 10 '12 at 11:13
  • $\begingroup$ Oh, I forgot to mention that a higher WorkingPrecision (for example 20) in FindRoot is still needed, and the option MaxIterations isn't necessary. $\endgroup$ – xzczd Nov 10 '12 at 11:53
3
$\begingroup$

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} *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$ – Amatya Nov 8 '12 at 8:53
  • $\begingroup$ As long as neither of sol is 0 then yes, badFactor is finite (but possibly very large). $\endgroup$ – b.gates.you.know.what Nov 8 '12 at 9:04
  • $\begingroup$ ...phbbhtt.. my brains are fried. Ocfcourse.. I know what badFactor is.. literally ... I can't believe I asked the second part of the question. :) $\endgroup$ – Amatya Nov 8 '12 at 9:08
2
$\begingroup$

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?

Improve this answer
$\endgroup$
12
  • $\begingroup$ I just now realized this is basically the same as Daniel suggests in the answer @xzczd linked to... $\endgroup$ – sebhofer Nov 8 '12 at 8:03
  • $\begingroup$ @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 $\endgroup$ – Amatya Nov 8 '12 at 8:33
  • $\begingroup$ @sebhofer Try {f1, f2} /. fr though. $\endgroup$ – b.gates.you.know.what Nov 8 '12 at 8:43
  • $\begingroup$ @b.gatessucks whoa.. I got {-3.67838*10^-9, 208.209} when I tried {f1, f2} /. fr. Whats going on? $\endgroup$ – Amatya Nov 8 '12 at 8:56
  • 1
    $\begingroup$ 囧…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… $\endgroup$ – xzczd Nov 8 '12 at 13:23

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.