0
$\begingroup$

Numerical precision/accuracy question here.

I am currently working on solving nonlinear equations with Reduce. Oftentimes, Reduce returns values, or combinations of values, that are wild (e.g., 10^{-169}, 10^{334}). An example from my work is below:

ConstraintsPDKD3 = {0.0782082 f + 1/3 (0.0392699 - 1.02013 k1 - 1.78339 k2 - 2.10318 k3),
  0.0782082 f + 
  1/3 (-0.011781 - 0.653523 k1 - 1.14525 k2 - 1.36932 k3), 
 1/2 (0.54957 f - 1.2747 k1 - 2.20973 k2 - 2.52383 k3 - 
    Sqrt[(-0.54957 f + 1.2747 k1 + 2.20973 k2 + 2.52383 k3)^2 - 
     4 (-0.523076 f k1 + 0.304552 k1^2 - 0.89299 f k2 + 
        1.05458 k1 k2 + 0.909021 k2^2 - 0.984585 f k3 + 
        1.20109 k1 k3 + 2.05057 k2 k3 + 1.13055 k3^2)]), 
 1/2 (0.54957 f - 1.2747 k1 - 2.20973 k2 - 2.52383 k3 + 
    Sqrt[(-0.54957 f + 1.2747 k1 + 2.20973 k2 + 2.52383 k3)^2 - 
     4 (-0.523076 f k1 + 0.304552 k1^2 - 0.89299 f k2 + 
        1.05458 k1 k2 + 0.909021 k2^2 - 0.984585 f k3 + 
        1.20109 k1 k3 + 2.05057 k2 k3 + 1.13055 k3^2)])}

where f, k1, k2, and k3 are all parameters. I then used Reduce:

ParamSols3 = 
  Reduce[
    ConstraintsPDKD3[[1]] == 0 && ConstraintsPDKD3[[2]] == 0 && 
    ConstraintsPDKD3[[3]] < 0 && ConstraintsPDKD3[[4]] < 0, 
    {k1, k2, k3, f}]`

and got:

((k1 <= -20.339 && -3.00902*10^-219 (-3.20072*10^217 + 
          2.91008*10^218 k1) - 
       1.60253*10^-169 Sqrt[-8.267332762368689*10^337 - 
         7.702472465239018*10^336 k1 + 
         1.895553130068739*10^334 k1^2] < 
      k2 < -3.00902*10^-219 (-3.20072*10^217 + 2.91008*10^218 k1) + 
       1.60253*10^-169 Sqrt[-8.267332762368689*10^337 - 
         7.702472465239018*10^336 k1 + 
         1.895553130068739*10^334 k1^2] && 
     k3 == 6.17779*10^-63 (1.12604*10^61 - 8.08636*10^61 k1 - 
         1.40754*10^62 k2)) || (-20.339 < 
      k1 < -16.499 && ((-3.00902*10^-219 (-3.20072*10^217 + 
              2.91008*10^218 k1) - 
           1.60253*10^-169 Sqrt[-8.267332762368689*10^337 - 
             7.702472465239018*10^336 k1 + 
             1.895553130068739*10^334 k1^2] < 
          k2 <= -2.43563*10^-187 (-8.96998*10^185 + 
              3.42892*10^186 k1) - 
           4.18487*10^-153 Sqrt[-9.56392*10^304 - 1.04989*10^304 k1 - 
             2.85003*10^302 k1^2] && 
         k3 == 6.17779*10^-63 (1.12604*10^61 - 8.08636*10^61 k1 - 
             1.40754*10^62 k2)) || (-2.43563*10^-187 (-8.96998*10^185 \
+ 3.42892*10^186 k1) + 
           4.18487*10^-153 Sqrt[-9.56392*10^304 - 1.04989*10^304 k1 - 
             2.85003*10^302 k1^2] <= 
          k2 < -3.00902*10^-219 (-3.20072*10^217 + 
              2.91008*10^218 k1) + 
           1.60253*10^-169 Sqrt[-8.267332762368689*10^337 - 
             7.702472465239018*10^336 k1 + 
             1.895553130068739*10^334 k1^2] && 
         k3 == 6.17779*10^-63 (1.12604*10^61 - 8.08636*10^61 k1 - 
             1.40754*10^62 k2)))) || (-16.499 <= 
      k1 < -10.4639 && -3.00902*10^-219 (-3.20072*10^217 + 
          2.91008*10^218 k1) - 
       1.60253*10^-169 Sqrt[-8.267332762368689*10^337 - 
         7.702472465239018*10^336 k1 + 
         1.895553130068739*10^334 k1^2] < 
      k2 < -3.00902*10^-219 (-3.20072*10^217 + 2.91008*10^218 k1) + 
       1.60253*10^-169 Sqrt[-8.267332762368689*10^337 - 
         7.702472465239018*10^336 k1 + 
         1.895553130068739*10^334 k1^2] && 
     k3 == 6.17779*10^-63 (1.12604*10^61 - 8.08636*10^61 k1 - 
         1.40754*10^62 k2))) && 
 f == 1.04399*10^-36 (-1.6032*10^35 + 4.16471*10^36 k1 + 
     7.28071*10^36 k2 + 8.58629*10^36 k3)

I don't mind having conditional expressions, but the very large and very small exponentials make it a mess.

To ensure there are values that satisfy these equations/inequalities that are on a reasonable order of magnitude, I have used FindInstance and recovered many solutions such as

ParamSols3 = 
  FindInstance[
    ConstraintsPDKD3[[1]] == 0 && ConstraintsPDKD3[[2]] == 0 && 
    ConstraintsPDKD3[[3]] < 0 && ConstraintsPDKD3[[4]] < 0, 
    {k1, k2, k3, f}, 20]
{{k1 -> -20.1284, k2 -> 17.0906, k3 -> -4.73627, 
  f -> -0.234088}, {k1 -> -17.046, k2 -> 15.371, k3 -> -4.78079, 
  f -> -0.302194}, {k1 -> -17.4559, k2 -> 15.4754, k3 -> -4.66681, 
  f -> -0.269071}, {k1 -> -14.8976, k2 -> 12.7925, k3 -> -3.61189, 
  f -> -0.0824322}, {k1 -> -11.7108, k2 -> 10.62, k3 -> -3.3148, 
  f -> -0.0764295}, {k1 -> -19.4291, k2 -> 16.5602, k3 -> -4.62438, 
  f -> -0.22237}, {k1 -> -11.7319, k2 -> 10.3636, k3 -> -3.08135, 
  f -> -0.0240402}, {k1 -> -19.2433, k2 -> 16.4122, k3 -> -4.58848, 
  f -> -0.217884}, {k1 -> -19.2912, k2 -> 16.0814, k3 -> -4.27693, 
  f -> -0.147596}, {k1 -> -18.8736, k2 -> 17.2364, k3 -> -5.48986, 
  f -> -0.425605}, {k1 -> -1015., k2 -> 871., k3 -> -250.256, 
  f -> -36.1311}, {k1 -> -19.5651, k2 -> 15.9578, k3 -> -4.03263, 
  f -> -0.088018}, {k1 -> -17.6034, k2 -> 15.0429, k3 -> -4.21706, 
  f -> -0.166124}, {k1 -> -16.5728, k2 -> 14.0221, k3 -> -3.84428, 
  f -> -0.102562}, {k1 -> -441., k2 -> 384.25, k3 -> -113.75, 
  f -> -16.5665}, {k1 -> -18.9253, k2 -> 15.6425, k3 -> -4.07804, 
  f -> -0.11021}, {k1 -> -10.6265, k2 -> 9.49818, k3 -> -2.88101, 
  f -> -0.000293597}, {k1 -> -18.0594, k2 -> 16.1724, k3 -> -4.97143, 
  f -> -0.325519}, {k1 -> -19.1149, k2 -> 16.3089, k3 -> -4.56279, 
  f -> -0.214586}, {k1 -> -1092., k2 -> 939., k3 -> -270.92, 
  f -> -39.2793}}

To be clear, I am not looking for solutions within a specific region of order of magnitude, I just am looking for a way to have the results from Reduce presented in a manner where these small/large numbers have already been handled by Mathematica.

$\endgroup$
1
$\begingroup$

You could try something like:

ExpandAll[ParamSols3] /. a_Real x_^n_:>Sign[a] (Expand[x Abs[a]^(1/n)])^n

k1 < -10.5768 && 0.0968595 - 0.8754 k1 - Sqrt[-2.12282 - 0.198645 k1 + 0.000194755 k1^2] < k2 < 0.0968595 - 0.8754 k1 + Sqrt[-2.12282 - 0.198645 k1 + 0.000194755 k1^2] && k3 == 0.0695649 - 0.49956 k1 - 0.869566 k2 && f == -0.167373 + 4.34792 k1 + 7.60104 k2 + 8.96402 k3

Addendum

If you get General::munfl errors, then it makes sense to rationalize the input first, and then numericize after (as suggested by Bill). The numericized result can then be polished as suggested by my answer. To clarify the steps, I will rewrite your constraints as follows:

eqns = And[
    ConstraintsPDKD3[[1]]==0,
    ConstraintsPDKD3[[2]]==0,
    ConstraintsPDKD3[[3]]<0,
    ConstraintsPDKD3[[4]]<0
];

Then, putting everything together:

sols = Reduce[Rationalize[eqns, 0], {k1, k2, k3, f}] //
    N //
    ExpandAll //
    ReplaceAll[a_Real x_^n_ :> Sign[a] Expand[x Abs[a]^(1/n)]^n]

k1 < -10.5768 && 0.0968595 - 0.8754 k1 - Sqrt[-2.12282 - 0.198645 k1 + 0.000194755 k1^2] < k2 < 0.0968595 - 0.8754 k1 + Sqrt[-2.12282 - 0.198645 k1 + 0.000194755 k1^2] && k3 == 0.0695649 - 0.49956 k1 - 0.869566 k2 && f == -0.167373 + 4.34792 k1 + 7.60104 k2 + 8.96402 k3

As for differences between my answer and yours, I can only speculate that the ConsraintPDKD3 expression you put in your question is not the same expression that you used in your Reduce call.

$\endgroup$
  • $\begingroup$ Hey @Carl! Thanks so much for the help :) Just a quick question if that is alright? I have used the method you suggested above, but got a different answer and precision warning. The warning reads "1.60253^{-169^{2}} is too small to represent as a normalised machine number; precision may be lost." One of the conditional expressions I got was very close to yours, however the k1 and k2 expressions are different. (I would post the result but it puts me over the character limit). Are these possibly related or could it be something else? Thanks again! $\endgroup$ – Cameron F. May 14 '19 at 16:45
1
$\begingroup$
N[Reduce[
  Rationalize[{constraint1==0,constraint2==0,constraint3==0,constraint4==0}],
  {k1, k2, k3, f}]]

seems to solve this particular example. You might need to include a second ,0 argument for Rationalize to force rationalization of every decimal approximation in every case.

If you have other examples where this doesn't work and can show those then perhaps we can find a way to deal with those too.

$\endgroup$
  • $\begingroup$ Hey @Bill ! Thanks for the suggestion, unfortunately in my case Constraints[[3]] < 0 && Constraints[[4]] < 0, which gives the same result with or without the Rationalize. I will take a look into this function though! $\endgroup$ – Cameron F. May 14 '19 at 16:48
  • 1
    $\begingroup$ Your third and fourth constraints have lots of decimal approximations that are almost exactly equal. Solving ends up subtracting those two almost equal, but not necessarily exactly equal, things resulting in large or huge exponents. To really overcome this problem I suspect will require dealing with this problem. Look at Rationalize[third,0]- Rationalize[fourth,0]//Simplify and look at the size of the coefficients that Mathematica has to deal with. Maybe you can really analyze your problem and come up with much simpler constraints. Lots of decimal approximations are often a problem. $\endgroup$ – Bill May 14 '19 at 18:01

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.