# Why can Solve solve this system of expressions but not a similar system?

In a meteorology textbook (Hogan & Mason, 2012) the following system is specified: n isn't specified above, but it's just a+b+c+d.

I tried to use Mathematica to reproduce the results from this textbook, in which a b c d are replaced with H F s. The following code works out PC in a few seconds:

eqn = {p == (a + d)/(a + b + c + d), s == (a + c)/(a + b + c + d),
h == (a)/(a + c), f == b/(b + d)}

Solve[eqn, {p}, {a, b, c, d}]


Then I tried to work out DSS:

eqn = {p == (ad - bc)/(Sqrt[(a + b) (c + d) (a + c) (b + d)]),
s == (a + c)/(a + b + c + d), h == (a)/(a + c), f == b/(b + d)}

Solve[eqn, {p}, {a, b, c, d}]


Mathematica ran for hours and couldn't work out the solution.

Googling I found a forum thread that suggested to improve speed by putting assumptions on the variables. I can assume that a b c d are non-negative integers and that n is a positive integer and that p h f s are non-negative reals. I can further assume that h f s are >= 0 and <= 1, and that p is >= -1 and <= 1.

The documentation for Solve doesn't indicate any way to put assumptions on the variables. I tried the simple solution

Solve[eqn, {p}, {a, b, c, d}, Reals]


but that also just kept running for ages.

Then I tried adapting @MarcoB's solution from this thread, I tried:

Assuming[
{{p, f, h, s} \[Element] Reals, {a, b, c, d} \[Element] Integers,
p >= -1, p <= 1, f >= 0, f <= 1, h >= 0, h <= 1, s >= 0, s <= 1},

Simplify[
Solve[eqn, {p}, {a, b, c, d}]
]
]


This also ran for hours without reaching a solution. Can anyone explain why Mathematica is not reaching a solution? I get that the DSS system is more complicated than the PC system, but it surprised me that this difference was sufficiently large that it precluded Mathematica from reaching an answer.

Hogan, R. J., & Mason, I. B. (2012). Deterministic forecasts of binary events. Forecast Verification: A Practitioner's Guide in Atmospheric Science, Second Edition, 31-59. Chicago.

• ad is a single variable; a d is the same as a*d; which do you want? Adding the Sqrt doubles the underlying degree of the system, not to mention all the extra products, which makes it considerably harder to solve. Oct 14, 2016 at 2:15
• @MichaelE2 Sorry, that's a typo in my question. But the issue of extreme slowness also occurs when the code is correctly rendered as (a d - b c) Oct 14, 2016 at 6:37

Is this what you want?

Clear[a, b, c, d, e, f];
eqn = {p^2 == (a d - b c)^2/(Sqrt[(a + b) (c + d) (a + c) (b + d)])^2,
s == (a + c)/(a + b + c + d), h == (a)/(a + c), f == b/(b + d)};

Solve[eqn, {p}, {a, b, c, d}]
(*
{{p -> -(((f - h) Sqrt[-1 + s] Sqrt[s])/
Sqrt[-f + f^2 + f s - 2 f^2 s - h s + 2 f h s + f^2 s^2 - 2 f h s^2 + h^2 s^2])},
{p -> ((f - h) Sqrt[-1 + s] Sqrt[s])/
Sqrt[-f + f^2 + f s - 2 f^2 s - h s + 2 f h s + f^2 s^2 - 2 f h s^2 + h^2 s^2]}}
*)

• +1. Why doesn't MMA solve p==... but does p^2==...? I mean, is squaring the equation a well justified method or a trick? Oct 14, 2016 at 8:00
• @corey979 It "rationalizes" the equation (gets rid of Sqrt[]) and is taught in algebra/pre-calculus texts in the US. it's usually easier for symbolic solvers to deal with polynomial systems, which this is easily converted to. It possibly introduces an extraneous solution, since both branches ±Sqrt[..] will satisfy the equation. That also makes it easier, since Mathematica does not have to analyze the branches (with so many variables); I suspect that is what takes so long. I think the checking also happens if you specify the domain Reals. Oct 14, 2016 at 10:23

Even after fixing your error in specifying the equations (your equation for p involved the variables ad and bc instead of the products a d and b c), Solve[] still took its sweet time. We can help Solve[] out a bit by using GroebnerBasis[] as a preprocessor, which we can do since we have an algebraic system of equations (no transcendental functions involved). More precisely, we can remove the intermediate variables a, b, c, d like so:

neq = Simplify[First @
GroebnerBasis[{p == (a d - b c)/(Sqrt[(a + b) (c + d) (a + c) (b + d)]),
s == (a + c)/(a + b + c + d), h == (a)/(a + c),
f == b/(b + d)}, {p, h, s, f}, {a, b, c, d}]]
f^2 (p^2 (-1 + s) - s) (-1 + s) + h s (h - p^2 - h s + h p^2 s) -
f (-1 + s) (-2 h s + p^2 (-1 + 2 h s))


Solve[] can now be used to solve for p:

Solve[neq == 0, p] // FullSimplify
{{p -> ((-f + h) Sqrt[-1 + s] Sqrt[s])/Sqrt[(f (-1 + s) - h s) (1 + f (-1 + s) - h s)]},
{p -> ((f - h) Sqrt[-1 + s] Sqrt[s])/Sqrt[(f (-1 + s) - h s) (1 + f (-1 + s) - h s)]}}


where two solutions were returned. Since (you can check this yourself!) h (1 - f) - (1 - h) f == h - f, the first solution is the desired one.

Unfortunately, the expression within the square root in the original source seems iffy, since (f (-1 + s) - h s) (1 + f (-1 + s) - h s)/(s (s - 1)) doesn't quite simplify to (h/(1 - s) + f/s) ((1 - h)/(1 - s) + (1 - f)/s). Checking which of these formulae is sensible to your application is something I'll leave for you to determine.