2

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.

Amatya
  • 6,888
  • 3
  • 26
  • 35
  • 1
    Does this question somewhat help you? – xzczd Nov 08 '12 at 07: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. – xzczd Nov 10 '12 at 05: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. – Amatya Nov 10 '12 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. – Amatya Nov 10 '12 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. – xzczd Nov 10 '12 at 11:53

2 Answers2

3

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} *)
b.gates.you.know.what
  • 20,103
  • 2
  • 43
  • 84
  • 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? – Amatya Nov 08 '12 at 08:53
  • As long as neither of sol is 0 then yes, badFactor is finite (but possibly very large). – b.gates.you.know.what Nov 08 '12 at 09: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. :) – Amatya Nov 08 '12 at 09:08
2

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?

sebhofer
  • 2,741
  • 18
  • 25
  • I just now realized this is basically the same as Daniel suggests in the answer @xzczd linked to... – sebhofer Nov 08 '12 at 08: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 – Amatya Nov 08 '12 at 08:33
  • @sebhofer Try {f1, f2} /. fr though. – b.gates.you.know.what Nov 08 '12 at 08:43
  • @b.gatessucks whoa.. I got {-3.67838*10^-9, 208.209} when I tried {f1, f2} /. fr. Whats going on? – Amatya Nov 08 '12 at 08:56
  • Just what you said yourself, your function is hyper sensitive. A very good solution to a slightly changed f2 is not a good solution for the original function. – b.gates.you.know.what Nov 08 '12 at 09:05
  • @b.gatessucks I don't think you can deduce that from your observation. I guess that the discrepancy is caused by floating-point arithmetic. The point is that f1r and f2r are the same as f1 and f2 up to errors which are smaller than MachinePrecision. However, manipulation of the rationalized functions does not involve floating-point arithmetic and is thus exact. – sebhofer Nov 08 '12 at 09:19
  • 1
    囧…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… – xzczd Nov 08 '12 at 13:23
  • BTW, I found that the option MaxIterations isn't necessary here and we only need WorkingPrecision -> 9 for this case. – xzczd Nov 08 '12 at 13:30
  • @xzczd That's an interesting observation. I wonder why Rationalize and SetPrecision behave differently in this case. However I find that I need a WorkingPrecision of at least 20 to get decent accuracy for the result. A value of 9 is even less then MachinePrecision, so I don't understand how that should work out. (Btw for WP equal to 9, Mathematica even displays a warning about no significant digits!) – sebhofer Nov 08 '12 at 15:11
  • @sebhofer Er…maybe I forgot to Clear the variables or just mistook the result yesterday, I tried the code again and this time WorkingPrecision -> 17 is needed: though the output for f2r is framed by a red box (what's the exact meaning of it?), {f1, f2} is OK. And to avoid the red box, WorkingPrecision -> 19 is enough. – xzczd Nov 09 '12 at 05:29
  • And I found that the documentation has already given a explanation for SetPrecision[…, 0] and Rationalize[…, Infinity], not difficult to understand. A somewhat related question is this. At last, let me give the same sigh as @b.gatessucks : this function is so sensitive! – xzczd Nov 09 '12 at 08:44
  • @sebhofer Rationalize is not working out for me. Even though I am getting solutions, when I substitute the solutions into expressions whose values I have a good idea about I am getting non-sense answers. I can post the example in my question. – Amatya Nov 09 '12 at 11:10