# Not getting expected result from numerical evaluation of solution given by Solve

Let $m$ be a large integer and $z\in (0,1)$. I want to find $0 < k < m$ such that

$\frac{\sqrt{5-\frac{4k}{m}}-1}{2} + \frac{m-\frac{k}{2}\sqrt{5-\frac{4k}{m}}}{m^2\sqrt{5-\frac{4k}{m}}}$ is closest to $z$. So I set up an equation by letting the preceding expression equal $z$ and solved it for $k$. That is,

Solve[(Sqrt[5 - 4 k/m] - 1)/2 + (m - k/2 Sqrt[5 - 4 k/m])/(Sqrt[5 - 4 k/m] m^2) == z, k]


and it gave some complicated symbolic solutions, as the equation can be easily reduced to a cubic equation. Then I plugged in the actual values $m=20000$ and $z=0.1$,

N[% /. {m -> 20000, z -> 0.1}]


which evaluates as

{{k -> 143373. + 0.0527344 I},
{k -> -3.2001*10^13 + 0. I},
{k -> -100572. - 0.0546875 I}}


which is clearly wrong. There should be a $k$ that is less than $m$, because the initial expression is at least $1/(2m)$ (plugging in $k = m$) and my choice of $z$ is greater than $1/(2m)$.

However, if I plug in the numbers directly into the equation,

Solve[(Sqrt[5 - 4 k/m] - 1)/2 + (m - k/2 Sqrt[5 - 4 k/m])/(Sqrt[5 - 4 k/m] m^2)
== z /. {m -> 20000, z -> 0.1}, k]


{{k -> 17800.5}, {k -> 25000.}}


So I am wondering whether this is a bug.

• Block[{z = 1/10, m = 20000}, Solve[eq == z, k, Reals]] or in longer form Block[{m = 20000, z = 1/10}, Solve[(Sqrt[5 - 4 k/m] - 1)/ 2 + (m - k/2 Sqrt[5 - 4 k/m])/(Sqrt[5 - 4 k/m] m^2) == z, k, Reals]] // N – ciao Sep 29 '15 at 23:35

You have not encountered a bug but a weakness in machine precision arithmetic. If you switch to Mathematica's more accurate arbitrary precision arithmetic, you can get away with an unconstrained (sloppy) problem statement.

sol =
Solve[
(Sqrt[5 - 4 k/m] - 1)/2 + (m - k/2 Sqrt[5 - 4 k/m])/(Sqrt[5 - 4 k/m] m^2) == z,
k];
N[sol /. {m -> 20000, z -> 1/10}, 25] // Chop

{{k -> 24999.99996527657202992317},
{k -> -3.200096001780046597611895*10^13},
{k -> 17800.46601084237600423031}}


But putting the problem to Solve with proper constraints (as others have observed) is a much better way to proceed because Solve, working under the constraints, give a solution that evaluates well even at machine precison.

sol =
Solve[
(Sqrt[5 - 4 k/m] - 1)/2 + (m - k/2 Sqrt[5 - 4 k/m]) /
(Sqrt[5 - 4 k/m] m^2) == z
&&
0 < k < m,
k, Reals]
N[sol /. {m -> 20000, z -> .1}]

{{k -> 17800.5}, {k -> Undefined}}


You need to build your assumptions into the problem. Reduce appears to work better than Solve for your problem. Let's look at both.

Solve

If one limits the solution to real numbers and wraps Solve with Assumptions which apply your inequalities a solution is found

sol = Assuming[
m > k > 0 && z > 0,
Solve[(Sqrt[5 - 4 k/m] - 1)/
2 + (m - k/2 Sqrt[5 - 4 k/m])/(Sqrt[5 - 4 k/m] m^2) == z, k,
Reals
]
]


Substituting {m, z} = {20000., 0.1} results in

sol /. {m -> 20000., z -> 0.1}

{{k -> Undefined}, {k -> 17800.5}, {k -> 25000.}}


The middle rule matches your expected result.

Reduce

It appears to go smoother using Reduce.

sol = Reduce[(Sqrt[5 - 4 k/m] - 1)/
2 + (m - k/2 Sqrt[5 - 4 k/m])/(Sqrt[5 - 4 k/m] m^2) == z &&
m > k > 0 && z > 0, k, Reals]


sol /. {m -> 20000., z -> 0.1}

k == 17800.5


Use Minimize. Since you want k such that the expression is "closest" you need to allow for the difference being positive or negative. I chose to square the difference to handle this.

diff[z_, m_] =
z - (Sqrt[5 - 4 k/m] - 1)/
2 + (m - k/2 Sqrt[5 - 4 k/m])/(Sqrt[5 - 4 k/m] m^2);

soln = With[{m = 20000, z = 0.1}, Minimize[{diff[z, m]^2, 0 < k < m}, k]]

(*  {6.83353*10^-19, {k -> 17799.5}}  *)


The difference is

diff[z, m] /. {z -> 0.1, m -> 20000, soln[[-1, 1]]}

(*  -8.26651*10^-10  *)


Compared with the original result

diff[z, m] /. {z -> 0.1, m -> 20000,
Solve[(Sqrt[5 - 4 k/m] - 1)/
2 + (m - k/2 Sqrt[5 - 4 k/m])/(Sqrt[5 - 4 k/m] m^2) ==
z /. {m -> 20000, z -> 0.1}, k][[1, 1]]}

(*  0.0000388349  *)


Minimize gives the "closest"