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 ``Rationalizeand
Rationalize[..,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.
Rationalize[…, 10^-16]
butRationalize[…, 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$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$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 inContourPlot
. $\endgroup$WorkingPrecision
(for example 20) inFindRoot
is still needed, and the optionMaxIterations
isn't necessary. $\endgroup$