1
$\begingroup$

I'm solving this system of equations with Mathematica:

roots = Solve[
  k1*A*N - k2*C == 0 && k3*B*N - k4*D == 0 && N0 - C - D - N == 0 && 
   A0 - C - A == 0 && B0 - D - B == 0 && k1 >= 0 && k2 >= 0 && 
   k3 >= 0 && k4 >= 0 && A >= 0 && B >= 0 && N >= 0 && C >= 0 && 
   D >= 0 && A0 >= 0 && B0 >= 0 && N0 >= 0, {A, B, C, D, N}, Reals, 
  Method -> Reduce]

the solutions in roots contain ConditionalExpression and Root objects. Plugging in values to roots gives the correct answers:

roots /. {A0 -> 100, B0 -> 100 , N0 -> 50, k1 -> 0.1, k2 -> 0.1, 
  k3 -> 0.1, k4 -> 0.1}

yields: {{A -> 75.1652, B -> 75.1652, C -> 24.8348, D -> 24.8348, N -> 0.330403}}

however I want to get the full expression of the (real) solutions, without Root objects. But when I use ToRadicals it doesn't preserve the solutions. Plugging in:

ToRadicals[roots] /. {A0 -> 100, B0 -> 100 , N0 -> 50, k1 -> 0.1, 
  k2 -> 0.1, k3 -> 0.1, k4 -> 0.1}

gives imaginary solutions:

{{A -> -7.58604*10^13 - 8.17902*10^11 I, 
  B -> 7.58604*10^13 + 8.17902*10^11 I, 
  C -> 7.58604*10^13 + 8.17902*10^11 I, 
  D -> -7.58604*10^13 - 8.17902*10^11 I, N -> -1. - 1.42109*10^-14 I}}

which aren't the solutions I'm looking for. Is there a way to get the full expression of the correct solution? Assuming here that all the constraints are met (i.e., all variables are positive reals). Thanks.

Improve this question
$\endgroup$
6
  • $\begingroup$ Might be that the radical forms are numerically unstable for evaluation in the range of interest. $\endgroup$
    – Daniel Lichtblau
    Oct 19, 2020 at 22:55
  • $\begingroup$ @DanielLichtblau but then why does it work in the same ranges when we use Root objects? That's what I don't understand. The result of roots /. {A0 -> 100, B0 -> 100 , N0 -> 50, k1 -> 0.1, k2 -> 0.1, k3 -> 0.1, k4 -> 0.1} is perfectly correct. $\endgroup$
    – pythonuser
    Oct 19, 2020 at 23:40
  • 1
    $\begingroup$ Root objects do not have the same issues with numeric (in)stability of evaluation as radical reformulations.One will not, for example, encounter a near-vanishing denominator in numerically evaluating a Root object after parameter substitution. $\endgroup$
    – Daniel Lichtblau
    Oct 20, 2020 at 13:46
  • 1
    $\begingroup$ For your parameter range, is there reason to expect there might be only one positive result? And does the external software have a root-finding capability e.g. by Newton's method? $\endgroup$
    – Daniel Lichtblau
    Oct 20, 2020 at 13:53
  • 1
    $\begingroup$ Root objects, after substitution of parameters, are handled by validated methods for root isolation and refinement. I do not know what might be the analog in Python. Given that you know there will be only one positive solution, you might try finding the largest root of the derivative (a cubic) numerically, then giving a larger value as a starting point for Newton iterations. How to find that root of the derivative? Work recursively. Two more derivatives brings it to linear. Granted, this is a bit of code to put in place. But it is fairly low-level and should run fast. $\endgroup$
    – Daniel Lichtblau
    Oct 20, 2020 at 16:49

1 Answer 1

Reset to default
3
$\begingroup$

The problem is outlined in the documentation of ToRadicals with a minimal example similar to

Root[#^3 - a &, 1]
ToRadicals[%]
{%%, %} /. {a -> 1}

ToRadicals minimal example

The problem is that Root refers to the n-th root in its second argument. For roots involving parameters ToRadicals does not guarantee/preserve this ordering. That being said for a given set/ranges of parameters one can find the correct Radical. For the example above for a<0 ToRadicals[Root[#^3 - a &, 1]] yields the correct result while for a>0 ToRadicals[Root[#^3 - a &, 3]] is equal to Root[#^3 - a &, 3]. Depending on the parameter one has to find the correct branch/root.

Regarding the concrete system in at hand using roots/.Root[a_,n_]:>Root[a,1]//ToRadicals for the given parameter set results in a numerically consistent result.

Hack

As Daniel Lichtblau already speculated in the comments: the result using radicals seems to be numerically not as stable as the solution using Root. Furthermore there is no guarantee that this 'hack' (...Root[a_,n_]:>Root[a,1]...) works for all parameter sets. I tried a few other consistent sets (all parameters greater 0) and it worked for them as well (which is in no way a guarantee that the given expression in radicals is valid for all consistent parameter sets). The other roots result in negative values in A,B,C,D or N so it is pretty easy to choose the correct branch for a given parameter set.

I want to point out that I do not recommend using the presented 'hack' without proofing that it works for all parameter sets involving only positive quantities. Even if one would succeed in such a proof (which I think is unlikely due to the complicated nature of the involved expressions) I would not recommend using it because the resulting expression is extremely lengthy (its LeafCount is 16303 compared to 800 for the roots) and because it is numerically inferior to the expression using Root. The implementation of Root is obviously clever enough to work around the numerical problems and instabilities encountered when using ToRadicals.

Improve this answer
$\endgroup$
1
  • $\begingroup$ thank you. the issue is that I want to take the result of Mathematica and use it in another programming language, so having an expression that doesn't use Mathematica-specific function calls/objects like Root would have been crucial (ie only using sqrt(), multiplication, division, powers, etc.). It sounds like this is not possible due to numerical instability of expressing this solution with radicals. $\endgroup$
    – pythonuser
    Oct 20, 2020 at 1:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.